/* nag_1d_cheb_interp_fit (e02afc) Example Program.
*
* Copyright 2017 Numerical Algorithms Group.
*
* Mark 26.2, 2017.
*
*/

#include <nag.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nage02.h>
#include <nagx01.h>
#include <math.h>

int main(void)
{
Integer exit_status = 0;
double *an = 0, d1, *f = 0, fit, pi, piby2n, *xcap = 0;

Integer i, j, n;
Integer r;
NagError fail;

INIT_FAIL(fail);

printf("nag_1d_cheb_interp_fit (e02afc) Example Program Results \n");

/* Skip heading in data file */
scanf("%*[^\n]");

/* nag_pi (x01aac).
* pi
*/
pi = nag_pi;
while ((scanf("%" NAG_IFMT "", &n)) != EOF)
{
if (n > 0) {
if (!(an = NAG_ALLOC(n + 1, double)) ||
!(f = NAG_ALLOC(n + 1, double)) ||
!(xcap = NAG_ALLOC(n + 1, double)))
{
printf("Allocation failure\n");
exit_status = -1;
goto END;
}
}
else {
printf("Invalid n.\n");
exit_status = 1;
return exit_status;
}

piby2n = pi * 0.5 / (double) n;
for (r = 0; r < n + 1; ++r)
scanf("%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));
}
}
/* nag_1d_cheb_interp_fit (e02afc).
* Computes the coefficients of a Chebyshev series
* polynomial for interpolated data
*/
nag_1d_cheb_interp_fit(n + 1, f, an, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_1d_cheb_interp_fit (e02afc).\n%s\n",
fail.message);
exit_status = 1;
goto END;
}

printf("\n          Chebyshev \n");
printf("   J   coefficient A(J) \n");
for (j = 0; j < n + 1; ++j)
printf(" %3" NAG_IFMT "%14.7f\n", j + 1, an[j]);
printf("\n   R    Abscissa   Ordinate      Fit \n");
for (r = 0; r < n + 1; ++r) {
/* nag_1d_cheb_eval (e02aec).
* Evaluates the coefficients of a Chebyshev series
* polynomial
*/
nag_1d_cheb_eval(n + 1, an, xcap[r], &fit, &fail);
if (fail.code != NE_NOERROR) {
printf("Error from nag_1d_cheb_eval (e02aec).\n%s\n", fail.message);
exit_status = 1;
goto END;
}

printf(" %3" NAG_IFMT "%11.4f%11.4f%11.4f\n", r + 1, xcap[r], f[r],
fit);
}
END:
NAG_FREE(an);
NAG_FREE(f);
NAG_FREE(xcap);
}
return exit_status;
}