/* nag_1d_cheb_interp_fit(e02afc) Example Program
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 */

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

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;
}