/* nag_fit_dim2_cheb_lines (e02cac) Example Program.
 *
 * Copyright 2023 Numerical Algorithms Group.
 *
 * Mark 29.1, 2023.
 */

#include <nag.h>
#include <stdio.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_fit_dim2_cheb_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("%" NAG_IFMT "%" NAG_IFMT "%" NAG_IFMT "%*[^\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%" NAG_IFMT "%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_fit_dim2_cheb_lines (e02cac).
       * Least squares surface fit by polynomials, data on lines
       */
      nag_fit_dim2_cheb_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_fit_dim2_cheb_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_fit_dim2_cheb_eval (e02cbc).
         * Evaluation of fitted polynomial in two variables
         */
        nag_fit_dim2_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_fit_dim2_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;
}