/* nag_1d_cheb_fit (e02adc) 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>

int main(int argc, char *argv[])
{
  FILE     *fpin, *fpout;
#define A(I, J) a[(I) *tda + J]

  Integer  exit_status = 0, i, iwght, j, k, m, r, tda;
  NagError fail;
  double   *a = 0, *ak = 0, d1, fit, *s = 0, *w = 0, *x = 0, x1, xarg, xcapr,
            xm, *y = 0;

  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_fit (e02adc) Example Program Results \n");

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

  while ((fscanf(fpin, "%ld", &m)) != EOF)
    {
      if (m >= 2)
        {
          if (
            !(x = NAG_ALLOC(m, double)) ||
            !(y = NAG_ALLOC(m, double)) ||
            !(w = NAG_ALLOC(m, double)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }
        }
      else
        {
          fprintf(fpout, "Invalid m.\n");
          exit_status = 1;
          return exit_status;
        }

      fscanf(fpin, "%ld", &k);
      if (k >= 1)
        {
          if (!(a = NAG_ALLOC((k+1)*(k+1), double)) ||
              !(s = NAG_ALLOC(k+1, double)) ||
              !(ak = NAG_ALLOC(k+1, double)))
            {
              fprintf(fpout, "Allocation failure\n");
              exit_status = -1;
              goto END;
            }
          tda = k+1;
        }
      else
        {
          fprintf(fpout, "Invalid k.\n");
          exit_status = 1;
          return exit_status;
        }

      fscanf(fpin, "%ld", &iwght);
      for (r = 0; r < m; ++r)
        {
          if (iwght != 1)
            {
              fscanf(fpin, "%lf", &x[r]);
              fscanf(fpin, "%lf", &y[r]);
              fscanf(fpin, "%lf", &w[r]);
            }
          else
            {
              fscanf(fpin, "%lf", &x[r]);
              fscanf(fpin, "%lf", &y[r]);
              w[r] = 1.0;
            }
        }
      /* nag_1d_cheb_fit (e02adc).
       * Computes the coefficients of a Chebyshev series
       * polynomial for arbitrary data
       */
      nag_1d_cheb_fit(m, k+1, tda, x, y, w, a, s, &fail);
      if (fail.code != NE_NOERROR)
        {
          fprintf(fpout, "Error from nag_1d_cheb_fit (e02adc).\n%s\n",
                  fail.message);
          exit_status = 1;
          goto END;
        }

      for (i = 0; i <= k; ++i)
        {
          fprintf(fpout, "\n");
          fprintf(fpout, " %s%4ld%s%12.2e\n", "Degree", i,
                  "   R.M.S. residual =", s[i]);
          fprintf(fpout, "\n   J  Chebyshev coeff A(J) \n");
          for (j = 0; j < i+1; ++j)
            fprintf(fpout, " %3ld%15.4f\n", j+1, A(i, j));
        }
      for (j = 0; j < k+1; ++j)
        ak[j] = A(k, j);
      x1 = x[0];
      xm = x[m-1];
      fprintf(fpout, "\n %s%4ld\n\n",
              "Polynomial approximation and residuals for degree", k);
      fprintf(fpout,
              "   R   Abscissa     Weight   Ordinate  Polynomial  Residual \n");
      for (r = 1; r <= m; ++r)
        {
          xcapr = (x[r-1] - x1 - (xm - x[r-1])) / (xm - x1);
          /* nag_1d_cheb_eval (e02aec).
           * Evaluates the coefficients of a Chebyshev series
           * polynomial
           */
          nag_1d_cheb_eval(k+1, ak, xcapr, &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;
            }

          d1 = fit - y[r-1];
          fprintf(fpout, " %3ld%11.4f%11.4f%11.4f%11.4f%11.2e\n", r, x[r-1],
                  w[r-1], y[r-1], fit, d1);
          if (r < m)
            {
              xarg = (x[r-1]  + x[r]) * 0.5;
              xcapr = (xarg - x1 - (xm - xarg)) / (xm - x1);
              /* nag_1d_cheb_eval (e02aec), see above. */
              nag_1d_cheb_eval(k+1, ak, xcapr, &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, "    %11.4f                      %11.4f\n", xarg,
                      fit);
            }
        }
 END:
      if (fpin != stdin) fclose(fpin);
      if (fpout != stdout) fclose(fpout);
      if (a) NAG_FREE(a);
      if (x) NAG_FREE(x);
      if (y) NAG_FREE(y);
      if (w) NAG_FREE(w);
      if (s) NAG_FREE(s);
      if (ak) NAG_FREE(ak);
    }
  return exit_status;
}