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

#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);

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

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

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

      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));
            }
        }
      /* 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)
        {
          Vprintf("Error from nag_1d_cheb_interp_fit (e02afc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      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)
        {
          /* 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)
            {
              Vprintf("Error from nag_1d_cheb_eval (e02aec).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

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