/* nag_2d_cheb_fit_lines (e02cac) Example Program.
 *
 * Copyright 2001 Numerical Algorithms Group.
 *
 * Mark 7, 2001.
 */

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

int main(void)
{
  /* Scalars */
  double ymax;
  Integer exit_status, i, j, k, l, mi, mj, n, r, t,
          na, one;
  NagError fail;

  /* Arrays */
  double *a = 0, *f = 0, *ff = 0, *w = 0,
         *x = 0, *xmax = 0, *xmin = 0, *y = 0;
  Integer *m = 0;

  INIT_FAIL(fail);
  exit_status = 0;
  Vprintf("e02cac Example Program Results\n");

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

  /* Input the number of lines Y = Y(I) on which data is given, */
  /* and the required degree of fit in the X and Y directions */
  while (scanf("%ld%ld%ld%*[^\n] ", &n, &k, &l) != EOF)
    {
      Vprintf("\n");
      if (n > 0)
        {
          /* Allocate arrays m, y, xmin and xmax */
          if ( !(m = NAG_ALLOC(n, Integer)) ||
               !(y = NAG_ALLOC(n, double)) ||
               !(xmin = NAG_ALLOC(n, double)) ||
               !(xmax = NAG_ALLOC(n, double)))
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

          mj = 0;
          /* Input Y(I), the number of data points on Y = Y(I) and the */
          /* range of X-values on this line, for I = 1,2,...N */
          for (i = 0; i < n; ++i)
            {
              Vscanf("%lf%ld%lf%lf%*[^\n] ", &y[i], &mi, &xmin[i], &xmax[i]);
              m[i] = mi;
              mj += mi;
            }

          /* Allocate arrays x, f, ff, w and a */
          na = (k + 1) * (l + 1);
          if ( !(x = NAG_ALLOC(mj, double)) ||
               !(f = NAG_ALLOC(mj, double)) ||
               !(ff = NAG_ALLOC(mj, double)) ||
               !(w = NAG_ALLOC(mj, double)) ||
               !(a = NAG_ALLOC(na, double)) )
            {
              Vprintf("Allocation failure\n");
              exit_status = -1;
              goto END;
            }

          /* Input the X-values and function values, F, together with */
          /* their weights, W. */
          for (i = 0; i < mj; ++i)
            Vscanf("%lf%lf%lf", &x[i], &f[i], &w[i]);
          Vscanf("%*[^\n] ");

          /* Evaluate the coefficients, A, of the fit to this set of data */
          one = 1;
          e02cac(m, n, k, l, x, y, f, w, a, xmin, xmax, y,  one, y, one, &fail);
          if (fail.code != NE_NOERROR)
            {
              Vprintf("Error from e02cac.\n%s\n", fail.message);
              exit_status = 1;
              goto END;
            }

          Vprintf("     Data Y     Data X     Data F   Fitted F   Residual\n");
          Vprintf("\n");

          mi = 0;
          for (r = 1; r <= n; ++r)
            {
              t = mi + 1;
              mi += m[r-1];
              ymax = y[n-1];
              if (n == 1)
                ymax += 1.0;

              /* Evaluate the fitted polynomial at each of the data points */
              /* on the line Y = Y(R) */
             e02cbc(t, mi, k, l, x, xmin[r-1], xmax[r-1],
                      y[r-1], y[0], ymax, ff, a, &fail);
             if (fail.code != NE_NOERROR)
               {
                 Vprintf("Error from e02cbc.\n%s\n", fail.message);
                 exit_status = 1;
                 goto END;
               }
             /* Output the data and fitted values on the line Y = Y(R) */
              for (i = t-1; i < mi; ++i)
                {
                  Vprintf("%11.4f%11.4f%11.4f%11.4f",
                          y[r-1], x[i], f[i], ff[i]);
                  Vprintf("%11.2e\n", ff[i] - f[i]);
                }
              Vprintf("\n");
            }

          /* Output the Chebyshev coefficients of the fit */
          Vprintf("Chebyshev coefficients of the fit\n");
          Vprintf("\n");

          for (j = 1; j <= k + 1; ++j)
            {
              for (i = (j - 1) * (l + 1); i < j * (l + 1); ++i)
                Vprintf("%11.4f ", a[i]);
              Vprintf("\n");
            }
        }
    }

 END:
  if (a) NAG_FREE(a);
  if (f) NAG_FREE(f);
  if (ff) NAG_FREE(ff);
  if (w) NAG_FREE(w);
  if (x) NAG_FREE(x);
  if (xmax) NAG_FREE(xmax);
  if (xmin) NAG_FREE(xmin);
  if (y) NAG_FREE(y);
  if (m) NAG_FREE(m);

  return exit_status;
}