/* nag_2d_cheb_fit_lines (e02cac) Example Program.
 *
 * Copyright 2014 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;
  printf("nag_2d_cheb_fit_lines (e02cac) Example Program Results\n");

  /* Skip heading in data file */
  scanf("%*[^\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)
    {
      printf("\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)))
            {
              printf("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)
            {
              scanf("%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)))
            {
              printf("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)
            scanf("%lf%lf%lf", &x[i], &f[i], &w[i]);
          scanf("%*[^\n] ");

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

          printf(
                  "     Data Y     Data X     Data F   Fitted F   Residual\n");
          printf("\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) */
              /* nag_2d_cheb_eval (e02cbc).
               * Evaluation of fitted polynomial in two variables
               */
              nag_2d_cheb_eval(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)
                {
                  printf("Error from nag_2d_cheb_eval (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)
                {
                  printf("%11.4f%11.4f%11.4f%11.4f",
                          y[r-1], x[i], f[i], ff[i]);
                  printf("%11.2e\n", ff[i] - f[i]);
                }
              printf("\n");
            }

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

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

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

  return exit_status;
}