/* nag_1d_cheb_interp_fit (e02afc) Example Program.
 *
 * Copyright 2014 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);

  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("%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)))
            {
              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(" %3ld%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(" %3ld%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;
}