# include # include # include # include # include # include "chebyshev_series.h" /******************************************************************************/ double echebser0 ( double x, double coef[], int nc ) /******************************************************************************/ /* Purpose: ECHEBSER0 evaluates a Chebyshev series. Discussion: This function implements a modification and extension of Maess's algorithm. Table 6.5.1 on page 164 of the reference gives an example for treating the first derivative. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ECHEBSER0, the value of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; x2 = 2.0 * x; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; } value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ double echebser1 ( double x, double coef[], int nc, double *y1 ) /******************************************************************************/ /* Purpose: ECHEBSER1 evaluates a Chebyshev series and first derivative. Discussion: This function implements a modification and extension of Maess's algorithm. Table 6.5.1 on page 164 of the reference gives an example for treating the first derivative. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ECHEBSER1, the value of the Chebyshev series at X. Output, double *Y1, the value of the 1st derivative of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; c0 = coef[nc-1]; x2 = 2.0 * x; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + x2 * c1; } } *y1 = c0 - c2; value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ double echebser2 ( double x, double coef[], int nc, double *y1, double *y2 ) /******************************************************************************/ /* Purpose: ECHEBSER2 evaluates a Chebyshev series and two derivatives. Discussion: This function implements a modification and extension of Maess's algorithm. Table 6.5.1 on page 164 of the reference gives an example for treating the first derivative. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ECHEBSER2, the value of the Chebyshev series at X. Output, double *Y1, the value of the 1st derivative of the Chebyshev series at X. Output, double *Y2, the value of the 2nd derivative of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; double d0; double d1 = 0.0; double d2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; c0 = coef[nc-1]; d0 = coef[nc-1]; x2 = 2.0 * x; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + x2 * c1; } if ( 1 < i ) { d2 = d1; d1 = d0; d0 = c0 - d2 + x2 * d1; } } *y2 = ( d0 - d2 ) * 4.0; *y1 = c0 - c2; value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ double echebser3 ( double x, double coef[], int nc, double *y1, double *y2, double *y3 ) /******************************************************************************/ /* Purpose: ECHEBSER3 evaluates a Chebyshev series and three derivatives. Discussion: This function implements a modification and extension of Maess's algorithm. Table 6.5.1 on page 164 of the reference gives an example for treating the first derivative. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ECHEBSER3, the value of the Chebyshev series at X. Output, double *Y1, the value of the 1st derivative of the Chebyshev series at X. Output, double *Y2, the value of the 2nd derivative of the Chebyshev series at X. Output, double *Y3, the value of the 3rd derivative of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; double d0; double d1 = 0.0; double d2 = 0.0; double e0; double e1 = 0.0; double e2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; c0 = coef[nc-1]; d0 = coef[nc-1]; e0 = coef[nc-1]; x2 = 2.0 * x; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + x2 * c1; } if ( 1 < i ) { d2 = d1; d1 = d0; d0 = c0 - d2 + x2 * d1; } if ( 2 < i ) { e2 = e1; e1 = e0; e0 = d0 - e2 + x2 * e1; } } *y3 = ( e0 - e2 ) * 24.0; *y2 = ( d0 - d2 ) * 4.0; *y1 = c0 - c2; value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ double echebser4 ( double x, double coef[], int nc, double *y1, double *y2, double *y3, double *y4 ) /******************************************************************************/ /* Purpose: ECHEBSER4 evaluates a Chebyshev series and four derivatives. Discussion: This function implements a modification and extension of Maess's algorithm. Table 6.5.1 on page 164 of the reference gives an example for treating the first derivative. Licensing: This code is distributed under the MIT license. Modified: 29 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ECHEBSER3, the value of the Chebyshev series at X. Output, double *Y1, *Y2, *Y3, *Y4, the value of the first four derivatives of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; double d0; double d1 = 0.0; double d2 = 0.0; double e0; double e1 = 0.0; double e2 = 0.0; double f0; double f1 = 0.0; double f2 = 0.0; int i; double y0; b0 = coef[nc-1]; c0 = coef[nc-1]; d0 = coef[nc-1]; e0 = coef[nc-1]; f0 = coef[nc-1]; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + 2.0 * x * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + 2.0 * x * c1; } if ( 1 < i ) { d2 = d1; d1 = d0; d0 = c0 - d2 + 2.0 * x * d1; } if ( 2 < i ) { e2 = e1; e1 = e0; e0 = d0 - e2 + 2.0 * x * e1; } if ( 3 < i ) { f2 = f1; f1 = f0; f0 = e0 - f2 + 2.0 * x * f1; } } y0 = ( b0 - b2 ) / 2.0; *y1 = c0 - c2; *y2 = ( d0 - d2 ) * 2.0 * 2.0; *y3 = ( e0 - e2 ) * 6.0 * 4.0; *y4 = ( f0 - f2 ) * 24.0 * 8.0; return y0; } /******************************************************************************/ double evenchebser0 ( double x, double coef[], int nc ) /******************************************************************************/ /* Purpose: EVENCHEBSER0 evaluates an even Chebyshev series. Discussion: This function implements Clenshaw's modification of his algorithm for even series. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double EVENCHEBSER0, the value of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; x2 = 4.0 * x * x - 2.0; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; } value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ double evenchebser1 ( double x, double coef[], int nc, double *y1 ) /******************************************************************************/ /* Purpose: EVENCHEBSER1 evaluates an even Chebyshev series and first derivative. Discussion: This function implements a modification and extension of Maess's algorithm. Table 6.5.1 on page 164 of the reference gives an example for treating the first derivative. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double EVENCHEBSER1, the value of the Chebyshev series at X. Output, double *Y1, the value of the 1st derivative of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; c0 = coef[nc-1]; x2 = 4.0 * x * x - 2.0; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + x2 * c1; } } *y1 = ( c0 - c2 ) * 4.0 * x; value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ double evenchebser2 ( double x, double coef[], int nc, double *y1, double *y2 ) /******************************************************************************/ /* Purpose: EVENCHEBSER2 evaluates an even Chebyshev series and first two derivatives. Discussion: This function implements a modification and extension of Maess's algorithm. Table 6.5.1 on page 164 of the reference gives an example for treating the first derivative. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double EVENCHEBSER2, the value of the Chebyshev series at X. Output, double *Y1, the value of the 1st derivative of the Chebyshev series at X. Output, double *Y2, the value of the 2nd derivative of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; double d0; double d1 = 0.0; double d2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; c0 = coef[nc-1]; d0 = coef[nc-1]; x2 = 4.0 * x * x - 2.0; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + x2 * c1; } if ( 1 < i ) { d2 = d1; d1 = d0; d0 = c0 - d2 + x2 * d1; } } *y2 = ( d0 - d2 ) * 64.0 * x * x + ( c0 - c2 ) * 4.0; *y1 = ( c0 - c2 ) * 4.0 * x; value = 0.5 * ( b0 - b2 ); return value; } /******************************************************************************/ double oddchebser0 ( double x, double coef[], int nc ) /******************************************************************************/ /* Purpose: ODDCHEBSER0 evaluates an odd Chebyshev series. Discussion: This function implements Clenshaw's modification of his algorithm for odd series. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ODDCHEBSER0, the value of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; int i; double value; double x2; b0 = coef[nc-1]; x2 = 4.0 * x * x - 2.0; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; b0 = coef[i] - b2 + x2 * b1; } value = ( b0 - b1 ) * x; return value; } /******************************************************************************/ double oddchebser1 ( double x, double coef[], int nc, double *y1 ) /******************************************************************************/ /* Purpose: ODDCHEBSER1 evaluates an odd Chebyshev series and the first derivative. Discussion: This function implements a modification and extension of Clenshaw's algorithm. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ODDCHEBSER1, the value of the Chebyshev series at X. Output, double *Y1, the value of the 1st derivative of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; double coefi; int i; double value; double x2; coefi = 2.0 * coef[nc-1]; b0 = coefi; c0 = coefi; x2 = 4.0 * x * x - 2.0; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; coefi = 2.0 * coef[i] - coefi; b0 = coefi - b2 + x2 * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + x2 * c1; } } *y1 = ( c0 - c2 ) * 4.0 * x * x + ( b0 - b2 ) * 0.5; value = ( b0 - b2 ) * 0.5 * x; return value; } /******************************************************************************/ double oddchebser2 ( double x, double coef[], int nc, double *y1, double *y2 ) /******************************************************************************/ /* Purpose: ODDCHEBSER2 evaluates an odd Chebyshev series and first two derivatives. Discussion: This function implements a modification and extension of Clenshaw's algorithm. Licensing: This code is distributed under the MIT license. Modified: 21 April 2014 Author: Manfred Zimmer Reference: Charles Clenshaw, Mathematical Tables, Volume 5, Chebyshev series for mathematical functions, London, 1962. Gerhard Maess, Vorlesungen ueber Numerische Mathematik II, Analysis, Berlin, Akademie_Verlag, 1984-1988, ISBN: 978-3764318840, LC: QA297.M325.   Parameters: Input, double X, the evaluation point. -1 <= X <= +1. Input, double COEF[NC], the Chebyshev series. Input, int NC, the number of terms in the Chebyshev series. 0 < NC. Output, double ODDCHEBSER2, the value of the Chebyshev series at X. Output, double *Y1, the value of the 1st derivative of the Chebyshev series at X. Output, double *Y2, the value of the 2nd derivative of the Chebyshev series at X. */ { assert ( 0 < nc && fabs ( x ) <= 1.0 ); double b0; double b1 = 0.0; double b2 = 0.0; double c0; double c1 = 0.0; double c2 = 0.0; double d0; double d1 = 0.0; double d2 = 0.0; double coefi; int i; double value; double x2; coefi = 2.0 * coef[nc-1]; b0 = coefi; c0 = coefi; d0 = coefi; x2 = 4.0 * x * x - 2.0; for ( i = nc - 2; 0 <= i; i-- ) { b2 = b1; b1 = b0; coefi = 2.0 * coef[i] - coefi; b0 = coefi - b2 + x2 * b1; if ( 0 < i ) { c2 = c1; c1 = c0; c0 = b0 - c2 + x2 * c1; } if ( 1 < i ) { d2 = d1; d1 = d0; d0 = c0 - d2 + x2 * d1; } } value = ( b0 - b2 ) * 0.5 * x; x2 = x * x; *y1 = ( c0 - c2 ) * 4.0 * x2 + ( b0 - b2 ) * 0.5; *y2 = ( ( d0 - d2 ) * 64.0 * x2 + ( c0 - c2 ) * 12.0 ) * x; return value; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }