/* nag_1d_cheb_interp_fit(e02afc) Example Program * * Copyright 1998 Numerical Algorithms Group. * * Mark 5, 1998. * */ #include #include #include #include #include #include int main(void) { #define NMAX 199 #define NP1MAX NMAX+1 double d1; double xcap[NP1MAX]; double f[NP1MAX]; double piby2n, pi, an[NP1MAX], fit; Integer i, j, n; Integer r; Vprintf("e02afc Example Program Results \n"); /* Skip heading in data file */ Vscanf("%*[^\n]"); pi = X01AAC; while((scanf("%ld",&n)) != EOF) { if (n > 0 && n <= NMAX) { piby2n = pi * 0.5 / (double) n; for (r = 0; r < n+1; ++r) Vscanf("%lf",&f[r]); for (r = 0; r < n+1; ++r) { i = r; /* The following method of evaluating xcap = cos(pi*i/n) * ensures that the computed value has a small relative error * and, moreover, is bounded in modulus by unity for all * i = 0, 1, ..., n. (It is assumed that the sine routine * produces a result with a small relative error for values * of the argument between -PI/4 and PI/4). */ if (2*i <= n) { d1 = sin(piby2n * i); xcap[i] = 1.0 - d1 * d1 * 2.0; } else if (2*i > n * 3) { d1 = sin(piby2n * (n - i)); xcap[i] = d1 * d1 * 2.0 - 1.0; } else { xcap[i] = sin(piby2n * (n - 2*i)); } } e02afc(n+1, f, an, NAGERR_DEFAULT); Vprintf("\n Chebyshev \n"); Vprintf(" J coefficient A(J) \n"); for (j = 0; j < n+1; ++j) Vprintf(" %3ld%14.7f\n",j+1,an[j]); Vprintf("\n R Abscissa Ordinate Fit \n"); for (r = 0; r < n+1; ++r) { e02aec(n+1, an, xcap[r], &fit, NAGERR_DEFAULT); Vprintf(" %3ld%11.4f%11.4f%11.4f\n",r+1,xcap[r],f[r],fit); } } else { Vprintf( "Incorrect input value of n.\n"); return EXIT_FAILURE; } } return EXIT_SUCCESS; }