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

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

int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;
  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);

  /* Check for command-line IO options */
  fpin = nag_example_file_io(argc, argv, "-data", NULL);
  fpout = nag_example_file_io(argc, argv, "-results", NULL);

  fprintf(fpout, "nag_1d_cheb_interp_fit (e02afc) Example Program Results \n");

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

  /* nag_pi (x01aac).
   * pi
   */
  pi = nag_pi;
  while ((fscanf(fpin, "%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)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }
        }
      else
        {
          fprintf(fpout, "Invalid n.\n");
          exit_status = 1;
          return exit_status;
        }

      piby2n = pi * 0.5 / (double) n;
      for (r = 0; r < n+1; ++r)
        fscanf(fpin, "%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)
        {
          fprintf(fpout, "Error from nag_1d_cheb_interp_fit (e02afc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      fprintf(fpout, "\n          Chebyshev \n");
      fprintf(fpout, "   J   coefficient A(J) \n");
      for (j = 0; j < n+1; ++j)
        fprintf(fpout, " %3ld%14.7f\n", j+1, an[j]);
      fprintf(fpout, "\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)
            {
              fprintf(fpout, "Error from nag_1d_cheb_eval (e02aec).\n%s\n",
                      fail.message);
              exit_status = 1;
              goto END;
            }

          fprintf(fpout, " %3ld%11.4f%11.4f%11.4f\n", r+1, xcap[r], f[r], fit);
        }
 END:
      if (fpin != stdin) fclose(fpin);
      if (fpout != stdout) fclose(fpout);
      if (an) NAG_FREE(an);
      if (f) NAG_FREE(f);
      if (xcap) NAG_FREE(xcap);
    }
  return exit_status;
}