/* nag_2d_spline_fit_scat(e02ddc) Example Program
 *
 * Copyright 1991 Numerical Algorithms Group.
 *
 * Mark 2, 1991.
 *
 * Mark 6 revised, 2000.
 */

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

#define MMAX 30
#define NMAX 7
int main(void)
{

  double f[MMAX], weights[MMAX], x[MMAX], y[MMAX], fg[NMAX*NMAX],
    px[NMAX], py[NMAX];
  Integer i, j, m, rank, nx, ny, npx, npy, nxest, nyest;
  double  s, delta, fp, xhi, yhi, xlo, ylo;
  Nag_Start start;
  Nag_2dSpline spline;
  double warmstartinf;

  nxest = 14;
  nyest = 14;
  Vprintf("e02ddc Example Program Results\n");
  Vscanf("%*[^\n]"); /* Skip heading in data file */
  /* Input the number of data-points m. */
  Vscanf("%ld",&m);
  if (m>0 && m<=MMAX)
    {
      /* Input the data-points and the weights. */
      for (i=0; i<m;  i++)
        Vscanf("%lf%lf%lf%lf",&x[i],&y[i],&f[i],&weights[i]);
      start = Nag_Cold;
      if (scanf("%lf",&s) != EOF)
        {
          /* Determine the spline approximation. */

          e02ddc(start, m, x, y, f, weights, s, nxest, nyest, &fp,
                 &rank, &warmstartinf, &spline, NAGERR_DEFAULT);
          nx = spline.nx;
          ny = spline.ny;
          Vprintf("\nCalling with smoothing factor s = %11.4e, nx = %1ld,"
                  " ny = %1ld\n",s,nx,ny);

          Vprintf("rank deficiency = %1ld\n\n",(nx-4)*(ny-4)-rank);

          /* Print the knot sets, lamda and mu. */
          Vprintf("Distinct knots in x direction located  at\n");
          for (j=3; j<spline.nx-3; j++)
            Vprintf("%12.4f%s",spline.lamda[j],((j-3)%5==4 || j==spline.nx-4)
                    ? "\n" : " ");
          Vprintf("\nDistinct knots in y direction located  at\n");
          for (j=3; j<spline.ny-3; j++)
            Vprintf("%12.4f%s",spline.mu[j],((j-3)%5==4 || j==spline.ny-4)
                    ? "\n" : " ");
          Vprintf("\nThe B-spline coefficients:\n\n");
          for (i=0; i<ny-4; i++)
            {
              for (j=0; j<nx-4; j++)
                Vprintf("%9.2f",spline.c[i+j*(ny-4)]);
              Vprintf("\n");
            }

          Vprintf("\n Sum of squared residuals fp = %11.4e\n",fp);
          if (nx==8 && ny==8)
            Vprintf("The spline is the least-squares bi-cubic polynomial\n");

          /* Evaluate the spline on a rectangular grid at npx*npy points
           * over the domain (xlo to xhi) x (ylo to yhi).
           */

          Vscanf("%ld%lf%lf",&npx,&xlo,&xhi);
          Vscanf("%ld%lf%lf",&npy,&ylo,&yhi);

          if (npx<=NMAX && npy<=NMAX)
            {
              delta = (xhi-xlo)/(npx-1);
              for (i=0; i<npx;  i++)
                px[i] = MIN(xhi,xlo+i*delta);
              for (i=0; i<npy;  i++)
                py[i] = MIN(yhi,ylo+i*delta);

              e02dfc(npx, npy, px, py, fg, &spline, NAGERR_DEFAULT);

              Vprintf("\nValues of computed spline:\n\n");
              Vprintf("          x");
              for (i=0; i<npx; i++)
                Vprintf("%8.2f ",px[i]);
              Vprintf("\n     y\n");
              for (i=npy-1; i>=0; i--)
                {
                  Vprintf("%8.2f   ",py[i]);
                  for (j=0; j<npx; j++)
                    Vprintf("%8.2f ", fg[npy*j+i]);
                  Vprintf("\n");
                }
              /* Free memory used by spline */
              NAG_FREE(spline.lamda);
              NAG_FREE(spline.mu);
              NAG_FREE(spline.c);
            }
          else
            {
              Vfprintf(stderr,"npx or npy is out of range: npx = %4ld,"
                       "npy = %4ld\n",npx,npy);
              /* Free memory used by spline */
              NAG_FREE(spline.lamda);
              NAG_FREE(spline.mu);
              NAG_FREE(spline.c);
              return EXIT_FAILURE;
            }
        }
      return EXIT_SUCCESS;
    }
  else
    {
      Vfprintf(stderr,"m is out of range: m = %4ld,\n",m);
      return EXIT_FAILURE;
    }
}