# include "data.h"

int fit (int argc, char **argv) {
  
  double **c, **b, *s, X, Y, dY;
  int i, j, nterm, mterm, order;
  int N, Weight, Quiet;
  Vector *xvec, *yvec, *dyvec;
  float *x, *y, *dy;
  char name[64];

  Quiet = FALSE;
  if (N = get_argument (argc, argv, "-q")) {
    Quiet = TRUE;
    remove_argument (N, &argc, argv);
  }
  if (N = get_argument (argc, argv, "-quiet")) {
    Quiet = TRUE;
    remove_argument (N, &argc, argv);
  }

  Weight = FALSE;
  if (N = get_argument (argc, argv, "-dy")) {
    remove_argument (N, &argc, argv);
    if ((dyvec = SelectVector (argv[N], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
    remove_argument (N, &argc, argv);
    Weight = TRUE;
  }

  if (argc != 4) {
    fprintf (stderr, "USAGE: fit x y order [-dy wt]\n");
    return (FALSE);
  }

  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);    

  order = atof (argv[3]);
  nterm = order + 1;
  mterm = 2*order + 1;

  if (xvec[0].Nelements != yvec[0].Nelements) {
    fprintf (stderr, "vectors must have same length\n");
    return (FALSE);
  }
  if (Weight && xvec[0].Nelements != dyvec[0].Nelements) {
    fprintf (stderr, "vectors must have same length\n");
    return (FALSE);
  }
  
  x = xvec[0].elements;
  y = yvec[0].elements;
  if (Weight) dy = dyvec[0].elements;

  ALLOCATE (s, double, 2*mterm);
  bzero (s, 2*mterm*sizeof(double));
  ALLOCATE (b, double *, nterm);
  ALLOCATE (c, double *, nterm);
  for (i = 0; i < nterm; i++) {
    ALLOCATE (c[i], double, nterm);
    bzero (c[i], nterm*sizeof(double));
    ALLOCATE (b[i], double, 1);
    bzero (b[i], sizeof(double));
  }

  for (i = 0; i < xvec[0].Nelements; i++, x++, y++) {
    if (!(finite(*x) && finite(*y))) continue;
    dY = 1.0;
    if (Weight) { 
      dY = 1.0 / SQ(*dy);
      dy ++;
    }
    X = 1*dY;
    Y = *y*dY;
    for (j = 0; j < nterm; j++) {
      s[j] += X;
      b[j][0] += Y;
      X = X * (*x);
      Y = Y * (*x);
    }
    for (j = nterm; j < mterm; j++) {
      s[j] += X;
      X = X * (*x);
    }
  }

  for (i = 0; i < nterm; i++) {
    for (j = 0; j < nterm; j++) {
      c[i][j] = s[i + j];
    }
  }
  gaussj (c, nterm, b, 1);

  /* print & save basic fit parameters */
  if (!Quiet) fprintf (stderr, "y = ");
  for (i = 0; i < nterm; i++) {
    sprintf (name, "C%d", i);
    set_variable (name, b[i][0]);
    if (!Quiet) fprintf (stderr, "%f x^%d ", b[i][0], i);
  }
  sprintf (name, "Cn");
  set_variable (name, (double) order);
  
  /* print & save basic fit parameters */
  if (!Quiet) fprintf (stderr, "\n");
  if (!Quiet) fprintf (stderr, "    ");
  for (i = 0; i < nterm; i++) {
    sprintf (name, "dC%d", i);
    set_variable (name, sqrt(c[i][i]));
    if (!Quiet) fprintf (stderr, "%f     ", sqrt(c[i][i]));
  }
  if (!Quiet) fprintf (stderr, "\n");

  return (TRUE);

}




