/* nag_1d_cheb_fit(e02adc) Example Program
 *
 * Copyright 1998 Numerical Algorithms Group.
 *
 * Mark 5, 1998.
 *
 */

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

int main(void)
{
#define MMAX 200
#define KP1MAX 50
#define NROWS KP1MAX

  double d1;
  double xarg;
  double a[NROWS][KP1MAX];
  double s[KP1MAX], w[MMAX], x[MMAX], y[MMAX], xcapr;
  double  x1, ak[KP1MAX], xm, fit;

  Integer i, j, k, m, r;
  Integer iwght;
  Integer tdim;

  tdim = KP1MAX;
  Vprintf("e02adc Example Program Results \n");

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

  while ((scanf("%ld",&m)) != EOF)
    {
      if (m > 0 && m <= MMAX)
        {
          Vscanf("%ld",&k);
          Vscanf("%ld",&iwght);
          if (k + 1 <= KP1MAX)
            {
              for (r = 0; r < m; ++r)
                {
                  if (iwght != 1)
                    {
                      Vscanf("%lf",&x[r]);
                      Vscanf("%lf",&y[r]);
                      Vscanf("%lf",&w[r]);
                    }
                  else
                    {
                      Vscanf("%lf",&x[r]);
                      Vscanf("%lf",&y[r]);
                      w[r] = 1.0;
                    }
                }
              e02adc(m, k+1, tdim, x, y, w, &a[0][0], s, NAGERR_DEFAULT);
              for (i = 0; i <= k; ++i)
                {
                  Vprintf("\n");
                  Vprintf(" %s%4ld%s%12.2e\n","Degree",i,"   R.M.S. residual =",s[i]);
                  Vprintf("\n   J  Chebyshev coeff A(J) \n");
                  for (j = 0; j < i+1; ++j)
                    Vprintf(" %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];
              Vprintf("\n %s%4ld\n","Polynomial approximation and residuals for degree",k);
              Vprintf("\n   R   Abscissa     Weight   Ordinate  Polynomial  Residual \n");
              for (r = 1; r <= m; ++r)
                {
                  xcapr = (x[r-1] - x1 - (xm - x[r-1])) / (xm - x1);
                  e02aec(k+1, ak, xcapr, &fit, NAGERR_DEFAULT);
                  d1 = fit - y[r-1];
                  Vprintf(" %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);
                      e02aec(k+1, ak, xcapr, &fit, NAGERR_DEFAULT);
                      Vprintf("    %11.4f                      %11.4f\n",xarg,fit);
                    }
                }
            }
        }
      else
        {
          Vprintf( "Incorrect input value of n.\n");
          return EXIT_FAILURE;
        }

    }
  return EXIT_SUCCESS;
}