/* @(#)dispersion.c	17.1.1.1 (ES0-DMD) 01/25/02 17:55:50 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* .COPYRIGHT   (C) 1993 European Southern Observatory     */
/* .IDENT       mosdisp.c                                  */
/* .AUTHORS     Pascal Ballester (ESO/Garching)            */
/*              Cristian Levin   (ESO/La Silla)            */
/* .KEYWORDS    Spectroscopy, Long-Slit, MOS               */
/* .PURPOSE                                                */
/* .VERSION     1.0  Package Modification 21-SEP-1993      */
/* ------------------------------------------------------- */

/* Purpose:
   This module is related to the dispersion relation and its representation.
   The upper modules (lncalib.c, lnrebin.c) do not know which representation is
   used (Chebyshev polynomials, or any other), nor how it is stored (which
   descriptor, meaning of coefficients). This design enables to replace 
   this module lndisp.c by any other based on a different representation.
 */

/* Note:

   Because the module is based on Numerical Recipes routines, which arrays
   start at index 1, the general rule is that arrays start at index 1. However
   the module can accomodate any starting index, which is specified at 
   initialisation time (function initdisp).
 */

/* Algorithms:

   The module provides ten top-level routines related to the evaluation of the
   dispersion relation. Only those six routines can be used in applications 
   related to the dispersion relation (wavelength calibration, rebinning, ...).
   The lower level routines must be considered as private.

   --- void mos_savedisp (save) ---

   Stores dispersion coefficients in array <save>.

   --- double fit_disp (ndata, degree, x, l) ---

   This routine evaluates the dispersion relation on the basis of ndata 
   couples of d values (x,l). Array indexes start at 0. The degree is 
   indicative. If ndata is too small, a smaller degree will be fitted 
   to the couples. Minimum value for ndata is 2. The procedure returns 
   a rough estimate of the average pixel size.

   --- double mos_fit_disp (ndata, degree, x, l, chi) ---

   This routine evaluates the dispersion relation on the basis of ndata 
   couples of d values (x,l). Array indexes start at 0. The degree is 
   indicative. If ndata is too small, a smaller degree will be fitted 
   to the couples. Minimum value for ndata is 2. The procedure returns 
   a rough estimate of the average pixel size and the error of the 
   dispersion relation.

   --- void eval_disp (x, l, ndata) ---

   Applies the last previously estimated (or loaded) dispersion relation 
   to ndata values of x and estimates the corresponding wavelength l. 
   Arrays start at index 0.

   --- void initdisp (name, mode, start) ---

   Open (mode = OLD) or Create (mode = NEW) a storage table (named <name>).

   --- void mos_initdisp (name, mode, start) ---

   Open (mode = OLD) or Create (mode = NEW) a storage table (named <name>).
   In addition to initdisp also the column :SLIT is created.

   --- void finishdisp() ---

   Closes the storage table

   --- void writedisp(line, pix_y, world_y) ---

   Store in the storage table the dispersion relation corresponding to the row
   number <line> and world coordinate <y>.

   --- void mos_writedisp(line, slit, pix_y, world_y, nyrow) ---

   Store in the storage table the dispersion relation corresponding to the 
   slitlet <slit>, row number <line>, and world coordinate <y>.

   --- double readdisp(line) ---

   Retrieves from the storage table the dispersion relation corresponding 
   to the row number <line>.  The procedure returns a rough estimate of 
   the average pixel size.

   --- void setdisp(degree, coefs) ---

   Sets the coefficients of the dispersion relation. coefs is an array
   which element 0 is the constant term and other elements represent
   the coefficients of higher degree

 */

#include <math.h>
#include <midas_def.h>
#include <tbldef.h>
#include <lndisp.h>
int FIT_SUCCESS = 0; /* Sucess indicator for the fit */

#include <proto_spec.h>
#include <numerical.h>

#ifndef PI  
#define PI 3.141592653589793
#endif

#define TOL 1.0e-5
#define MAXNCOE 100

#define  FUNCTION        fpoly   
/* #define  mosfunc         fcheb   */

/* void     (*function) (), function = fpoly;   */
/* Another way to define a variable pointing to a function */

double pixbin;			/* Pixel bin in wavelength units */
double coef[MAXNCOE];		/* Array of coefficients  */
int fdeg, ncoef, refdeg, maxcoef;	/* Number of coefficients, Degree */

int tide, colslit, colline, coly, colrms, colcoef[MAXNCOE];
int linpix, linrms;
int start_index, nbline;

void fpoly(), fcheb();

/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
mos_savedisp (double save[])
#else
void 
mos_savedisp (save)
     double save[];
#endif

{

  int i;

  /* printf("Dispersion Relation. Degree: %d. Refdeg: %d. MaxCoef:%d\n",
     fdeg,refdeg,maxcoef);
     printf("Coefficients: ");  */

  for (i = 0; i <= ncoef-1; i++)
    {
      save[i] = coef[i+1];
      /*     printf("save[%d] =  %f ",i, save[i]);  */
    }

}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
printdisp (void)
#else
void 
printdisp ()
#endif

{

  int i;

  printf ("Dispersion Relation. Degree: %d. Refdeg: %d. MaxCoef:%d\n",
	  fdeg, refdeg, maxcoef);
  printf ("Coefficients: ");
  for (i = 1; i <= ncoef; i++)
    printf (" %f ", coef[i]);
  printf ("\n");

}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
setrefdeg (int deg)
#else
void 
setrefdeg (deg)
     int deg;
#endif

{

  refdeg = deg;
  maxcoef = refdeg + 1;

}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
set_zero (int deg)
#else
void 
set_zero (deg)
     int deg;
#endif


{

  int i;

  fdeg = deg;
  refdeg = deg;
  ncoef = deg + 1;
  maxcoef = deg + 1;

  for (i = 1; i <= ncoef; i++)
    coef[i] = 0;
}
/* -------------------------------------------------------------------------- */
#ifdef __STDC__
void 
setdisp (int deg, double coefs[])
#else
void 
setdisp (deg, coefs)
     int deg;
     double coefs[];
#endif

{
  int i;

  fdeg = deg;
  refdeg = deg;
  ncoef = deg + 1;
  maxcoef = deg + 1;

  for (i = 1; i <= ncoef; i++)
    coef[i] = coefs[i - 1];

  FIT_SUCCESS = 1;

}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
double 
fit_disp (int *ndata, int *deg, double x[], double l[])
#else
double 
fit_disp (ndata, deg, x, l)

     int *ndata;
     int *deg;
     double x[], l[];
#endif

{

  double *sig, **covar, *chisq, *a, *xsh, *lsh;
  int *lista;
  int i;

  int maxdeg;
  /* printf ("NData : %d   ---   Deg value : %d\n", *ndata, *deg); */

  refdeg = *deg;
  maxdeg = *ndata - 1;
  if (refdeg > maxdeg)
    fdeg = maxdeg;
  else
    fdeg = refdeg;
  ncoef = fdeg + 1;
  maxcoef = refdeg + 1;

  FIT_SUCCESS = 0;

  for (i=0; i<MAXNCOE; i++) coef[i] = 0.;

  if (*ndata < 2)
    {
      printf ("Not enough lines (minimum is 2). \nNo dispersion relation computed\n");
      FIT_SUCCESS = -2;
      return (0.);
    }

  if (fdeg < 1)
    {
      printf ("Degree : %d. No dispersion relation fitted\n", *deg);
      FIT_SUCCESS = -1;
      return (0.);
    }

  covar = dmatrix (1, *ndata, 1, *ndata);
  chisq = dvector (0, *ndata);
  sig = dvector (1, *ndata);

  xsh = x + start_index - 1;
  lsh = l + start_index - 1;

  lista = ivector (1, ncoef);

  for (i = 1; i <= ncoef; i++)
    lista[i] = i;
  for (i = 1; i <= *ndata; i++)
    sig[i] = 1.;

  lfit2 (x, l, sig, *ndata, coef, ncoef, lista, ncoef, covar, chisq, FUNCTION);

  free_dmatrix (covar, 1, *ndata, 1, *ndata);
  free_dvector (chisq, 0, *ndata);
  free_dvector (sig, 1, *ndata);
  free_ivector (lista, 1, ncoef);

  FIT_SUCCESS = 1;
  pixbin = coef[2];
  return (pixbin);

}


/* -------------------------------------------------------------------------- */

#ifdef __STDC__
double 
mos_fit_disp (int *ndata, int *deg, double x[],
	      double l[], double *chi)
#else
double 
mos_fit_disp (ndata, deg, x, l, chi)
     int *ndata, *deg;
     double x[], l[], *chi;
#endif

{

  double *sig, **u, **v, *w, *chisq, *a, *xsh, *lsh, **covar;
  int    *lista;
  int    actvals, maxdeg, i; 
  char   poltyp[8];

/* printf("NData : %d   ---   Deg value : %d\n",*ndata,*deg); */

  refdeg = *deg; 
  maxdeg = *ndata - 1;
  if (refdeg > maxdeg)
    fdeg = maxdeg;
  else
    fdeg = refdeg;
  ncoef = fdeg + 1;
  maxcoef = refdeg + 1;

  if (*ndata < 2)
    {
      printf ("Not enough lines (minimum is 2). \nNo dispersion relation computed\n");
      return (-2.);
    }

  else if (fdeg < 1)
    {
      printf ("Degree : %d. No dispersion relation fitted\n", *deg);
      return (-1.);
    }

  else
    {
      covar = (double **) dmatrix (1, *ndata, 1, *ndata);
      chisq = (double *) dvector (0, *ndata);
      sig = (double *) dvector (1, *ndata);

      xsh = x + start_index - 1;
      lsh = l + start_index - 1;

      lista = (int *) ivector (1, ncoef);
      for (i = 1; i <= ncoef; i++)
	lista[i] = i;
      for (i = 1; i <= *ndata; i++)
        sig[i] = 1.;

      SCKGETC("POLTYP", 1, 8, &actvals, poltyp);
      if (toupper(poltyp[0]) == 'L')
        lfit2(x,l,sig,*ndata,coef,ncoef,lista,ncoef,covar,chisq,fleg);
      else
        lfit2(x,l,sig,*ndata,coef,ncoef,lista,ncoef,covar,chisq,fcheb);

      *chi = *chisq;

      free_dmatrix (covar, 1, *ndata, 1, *ndata);
      free_dvector (chisq, 0, *ndata);
      free_dvector (sig, 1, *ndata);
      free_ivector (lista, 1, ncoef);

      pixbin = coef[2];
      return (pixbin);
    }
}
/* -------------------------------------------------------------------------- */
#ifdef __STDC__
void 
mos_eval_disp (double x[], double l[], int n)
#else
void 
mos_eval_disp (x, l, n)
     double x[], l[];
     int n;
#endif

{

  int actvals, i, icoef;
  double xp[MAXNCOE];
  char   poltyp[8];

  SCKGETC("POLTYP", 1, 8, &actvals, poltyp);
  for (i = start_index; i < (n + start_index); i++)
    {
      l[i] = 0.;
      if (toupper(poltyp[0]) == 'L')
        fleg (x[i], xp, ncoef);
      else 
        fcheb (x[i], xp, ncoef);
      for (icoef = 1; icoef <= ncoef; icoef++)
	l[i] += coef[icoef] * xp[icoef];
    }
}
/* -------------------------------------------------------------------------- */
#ifdef __STDC__
void 
eval_disp (double x[], double l[], int n)
#else
void 
eval_disp (x, l, n)
     double x[], l[];
     int n;
#endif

{

  int i, icoef;
  double xp[MAXNCOE];

  if (FIT_SUCCESS <= 0) {
         printf("No dispersion relation fitted. No evaluation.\n");
         return;
         }

  for (i = start_index; i < (n + start_index); i++)
    {
      l[i] = 0.;
      FUNCTION (x[i], xp, ncoef);
      for (icoef = 1; icoef <= ncoef; icoef++)
	l[i] += coef[icoef] * xp[icoef];
    }
}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
initdisp (char name[], char mode[], int start)
#else
void 
initdisp (name, mode, start)
     char name[], mode[];
     int start;
#endif

{

  int icoef;
  char colnam[30], numb[10];
  int actvals, null;
  int kunit;
  int ncol, nrow, nsort, allcol, allrow;


  start_index = start;

  if (toupper (mode[0]) == 'N')
    {
      if (TCTINI (name, F_TRANS, F_IO_MODE, 5L, 10L, &tide))
	SCTPUT ("**** Error while creating output table");
      nbline = 0;
    }
  else
    {
      if (TCTOPN (name, F_IO_MODE, &tide))
	SCTPUT ("**** Error while opening output table");
      SCDRDD (tide, "LNPIX", 1L, 1L, &actvals, &pixbin, &kunit, &null);
      SCDRDI (tide, "LNDEG", 1L, 1L, &actvals, &refdeg, &kunit, &null);
      SCDRDI (tide, "LNCOE", 1L, 1L, &actvals, &maxcoef, &kunit, &null);
      fdeg = refdeg, ncoef = maxcoef;
      TCIGET (tide, &ncol, &nrow, &nsort, &allcol, &allrow);
      nbline = nrow;
    }

  TCCSER (tide, ":ROW", &colline);
  if (colline == (-1))
    TCCINI (tide, D_I4_FORMAT, 1, "I6", "Row Number",
	    "ROW", &colline);

  TCCSER (tide, ":Y", &coly);
  if (coly == (-1))
    TCCINI (tide, D_R8_FORMAT, 1, "F8.2", "Y Value",
	    "Y", &coly);

  for (icoef = 1; icoef <= maxcoef; icoef++)
    {
      strcpy (colnam, ":COEF_");
      sprintf (numb, "%d", icoef);
      strcat (colnam, numb);
      TCCSER (tide, colnam, &colcoef[icoef]);
      if (colcoef[icoef] == (-1))
	TCCINI (tide, D_R8_FORMAT, 1, "F16.10", "Coefficients",
		colnam, &colcoef[icoef]);
    }


    TCCSER(tide, ":PIXEL", &linpix);
       if (linpix == (-1)) TCCINI(tide, D_R8_FORMAT, 1, "F10.3", 
                                 "Angstrom/pixel", "PIXEL", &linpix);
 
    TCCSER(tide, ":RMS", &linrms);
       if (linrms == (-1)) TCCINI(tide, D_R8_FORMAT, 1, "F10.3", "Angstrom",
                                "RMS", &linrms);

}

/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
mos_initdisp (char name[], char mode[], int start)
#else
void 
mos_initdisp (name, mode, start)
     char name[], mode[];
     int start;
#endif

{

  int icoef;
  char colnam[30], numb[10];
  int actvals, null;
  int kunit;
  int ncol, nrow, nsort, allcol, allrow;


  start_index = start;

  if (toupper (mode[0]) == 'N')
    {				/*  NEW table */
      if (TCTINI (name, F_TRANS, F_IO_MODE, 5L, 10L, &tide))
	SCTPUT ("**** Error while creating output table");
      nbline = 0;
    }
  else
    {
      if (TCTOPN (name, F_IO_MODE, &tide))
	SCTPUT ("**** Error while opening output table");
      SCDRDD (tide, "LNPIX", 1L, 1L, &actvals, &pixbin, &kunit, &null);
      SCDRDI (tide, "LNDEG", 1L, 1L, &actvals, &refdeg, &kunit, &null);
      SCDRDI (tide, "LNCOE", 1L, 1L, &actvals, &maxcoef, &kunit, &null);
      /* fdeg=refdeg, ncoef=maxcoef; */
      TCIGET (tide, &ncol, &nrow, &nsort, &allcol, &allrow);
      nbline = nrow;
    }

  TCCSER (tide, ":SLIT", &colslit);
  if (colslit == (-1))
    TCCINI (tide, D_I4_FORMAT, 1, "I6", "Slit Number",
	    "SLIT", &colslit);

  TCCSER (tide, ":ROW", &colline);
  if (colline == (-1))
    TCCINI (tide, D_I4_FORMAT, 1, "I6", "Row Number",
	    "ROW", &colline);

  TCCSER (tide, ":Y", &coly);
  if (coly == (-1))
    TCCINI (tide, D_R8_FORMAT, 1, "F8.2", "Y Value",
	    "Y", &coly);

  TCCSER (tide, ":RMS", &colrms);
  if (colrms == (-1))
    TCCINI (tide, D_R8_FORMAT, 1, "F8.4", "Angstrom",
	    "RMS", &colrms);

  for (icoef = 1; icoef <= maxcoef; icoef++)
    {
      strcpy (colnam, ":COEF_");
      sprintf (numb, "%d", icoef);
      strcat (colnam, numb);
      TCCSER (tide, colnam, &colcoef[icoef]);
      if (colcoef[icoef] == (-1))
	TCCINI (tide, D_R8_FORMAT, 1, "F16.10", "Coefficients",
		colnam, &colcoef[icoef]);
    }
}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
finishdisp (void)
#else
void 
finishdisp ()
#endif

{

  int kunit;

  SCDWRD (tide, "LNPIX", &pixbin, 1L, 1L, &kunit);
  SCDWRI (tide, "LNDEG", &refdeg, 1L, 1L, &kunit);
  SCDWRI (tide, "LNCOE", &maxcoef, 1L, 1L, &kunit);


  TCSINI(tide);
  TCTCLO (tide);

}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
double 
readdisp (int y)
#else
double 
readdisp (y)
     int y;
#endif

/* long int  y;     Pixel number of the row */
{

  int line, yval, linext, ydif, ynext, null, icoef;
  int mindif;

  mindif = -1;

  for (line = 1; line <= nbline; line++)
    {
      TCERDI (tide, line, colline, &yval, &null);
      if (!null)
	{
	  ydif = y - yval;
	  if (ydif < 0)
	    ydif = (-1) * ydif;
	  if (mindif < 0)
	    mindif = ydif;
	  if (ydif <= mindif)
	    mindif = ydif, linext = line;
	}
    }

  fdeg = refdeg, ncoef = maxcoef;

  for (icoef = 1; icoef <= ncoef; icoef++)
    TCERDD (tide, linext, colcoef[icoef], &coef[icoef], &null);

  FIT_SUCCESS = 1;

  return (pixbin);

}

/*----------------------------------------------------------------------------*/
#ifdef __STDC__
int 
mos_readdisp (int y, int slit)
#else
int 
mos_readdisp (y, slit)
     int y, slit;
#endif

/* long int  y;     Pixel number of the row */
{

  int line, yval, linext, ydif, ynext, null, icoef;
  int mindif, actslit;

  mindif = -1;

  for (line = 1; line <= nbline; line++)
    {
      TCERDI (tide, line, colline, &yval, &null);
      TCERDI (tide, line, colslit, &actslit, &null);
      if (!null && actslit == slit)
	{
	  ydif = y - yval;
	  if (ydif < 0)
	    ydif = (-1) * ydif;
	  if (mindif < 0)
	    mindif = ydif;
	  if (ydif <= mindif)
	    mindif = ydif, linext = line;
	}
    }

  if (mindif == -1)
    return (-1);

  fdeg = refdeg, ncoef = maxcoef;

  for (icoef = 1; icoef <= ncoef; icoef++)
    TCERDD (tide, linext, colcoef[icoef], &coef[icoef], &null);

  return (0);

}



/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
writedisp (int line, int ypix, double y, double pixel, double rms)
#else
void 
writedisp (line, ypix, y, pixel, rms)
     int line, ypix;
     double y, pixel, rms;
#endif

{

  int icoef;

  TCEWRI (tide, line, colline, &ypix);
  TCEWRD (tide, line, coly, &y);

  if (nbline < line)
    nbline = line;

  for (icoef = 1; icoef <= maxcoef; icoef++)
    TCEWRD (tide, line, colcoef[icoef], &coef[icoef]);

 TCEWRD(tide, line, linpix,  &pixel);
 TCEWRD(tide, line, linrms,  &rms);

}
/* -------------------------------------------------------------------------- */

#ifdef __STDC__
void 
mos_writedisp (int line, int slit, int ypix,
	       double y, int numrow, double rms)
#else
void 
mos_writedisp (line, slit, ypix, y, numrow, rms)
     int line, slit, ypix, numrow;
     double y, rms;
#endif

{

  int icoef;

  TCEWRI (tide, line, colslit, &slit);
  TCEWRI (tide, line, colline, &ypix);
  TCEWRD (tide, line, coly, &y);
  TCEWRD (tide, line, colrms, &rms);

/*      if (nbline<line) nbline=line;   */
  if (nbline < line)
    nbline = numrow;

  for (icoef = 1; icoef <= maxcoef; icoef++)
    TCEWRD (tide, line, colcoef[icoef], &coef[icoef]);
}

/**********************************************************************/
/*                   Private section of module lndisp.c               */
/* (In all programs from Numerical Recipes, arrays start at index 1)  */
/**********************************************************************/

static float sqrarg;
#define SQR(a) (sqrarg=(a),sqrarg*sqrarg) 

#ifdef __STDC__
void 
lfit2 (double x[], double y[], double sig[], int ndata, double a[],
       int ma, int lista[], int mfit, double **covar, double *chisq,
       void (*funcs) (double, double *, int))
#else
void 
lfit2 (x, y, sig, ndata, a, ma, lista, mfit, covar, chisq, funcs)
     int ndata, ma, lista[], mfit;
     double x[], y[], sig[], a[], **covar, *chisq;
     void (*funcs) ();		/* ANSI: void (*funcs)(double, double *,int); */
#endif

{
  int k, kk, j, ihit, i;
  double ym, wt, sum, sig2i, **beta, *afunc;

  beta = (double **) dmatrix (1, ma, 1, 1);
  afunc = (double *) dvector (1, ma);
  kk = mfit + 1;
  for (j = 1; j <= ma; j++)
    {
      ihit = 0;
      for (k = 1; k <= mfit; k++)
	if (lista[k] == j)
	  ihit++;
      if (ihit == 0)
	lista[kk++] = j;
      else if (ihit > 1)
	nrerror ("Bad LISTA permutation in LFIT-1");
    }
  if (kk != (ma + 1))
    nrerror ("Bad LISTA permutation in LFIT-2");
  for (j = 1; j <= mfit; j++)
    {
      for (k = 1; k <= mfit; k++)
	covar[j][k] = 0.0;
      beta[j][1] = 0.0;
    }
  for (i = 1; i <= ndata; i++)
    {
      (*funcs) (x[i], afunc, ma);
      ym = y[i];
      if (mfit < ma)
	for (j = (mfit + 1); j <= ma; j++)
	  ym -= a[lista[j]] * afunc[lista[j]];
      sig2i = 1.0 / SQR (sig[i]);
      for (j = 1; j <= mfit; j++)
	{
	  wt = afunc[lista[j]] * sig2i;
	  for (k = 1; k <= j; k++)
	    covar[j][k] += wt * afunc[lista[k]];
	  beta[j][1] += ym * wt;
	}
    }
  if (mfit > 1)
    for (j = 2; j <= mfit; j++)
      for (k = 1; k <= j - 1; k++)
	covar[k][j] = covar[j][k];
  gaussj (covar, mfit, beta, 1);
  for (j = 1; j <= mfit; j++)
    a[lista[j]] = beta[j][1];
  *chisq = 0.0;
  for (i = 1; i <= ndata; i++)
    {
      (*funcs) (x[i], afunc, ma);
      for (sum = 0.0, j = 1; j <= ma; j++)
	sum += a[j] * afunc[j];
      *chisq += SQR ((y[i] - sum) / sig[i]);
    }
  covsrt (covar, ma, lista, mfit);
  free_dvector (afunc, 1, ma);
  free_dmatrix (beta, 1, ma, 1, 1);
}

#undef SQR 

#ifdef __STDC__
void 
fpoly (double x, double p[], int np)
#else
void 
fpoly (x, p, np)
     double x, p[];
     int np;
#endif

{
  int j;

  p[1] = 1.0;
  for (j = 2; j <= np; j++)
    p[j] = p[j - 1] * x;
}


#ifdef __STDC__
void 
fleg (double x, double pl[], int nl)
#else
void 
fleg (x, pl, nl)
     double x, pl[];
     int nl;
#endif

{
  int j;
  double twox, f2, f1, d;

  pl[1] = 1.0;
  pl[2] = x;
  if (nl > 2)
    {
      twox = 2.0 * x;
      f2 = x;
      d = 1.0;
      for (j = 3; j <= nl; j++)
	{
	  f1 = d;
	  d += 1.0;
	  f2 += twox;
	  pl[j] = (f2 * pl[j - 1] - f1 * pl[j - 2]) / d;
	}
    }
}

/*----------------------------------------------------------------------------*/
#ifdef __STDC__
void 
fcheb (double x, double pl[], int nl)
#else
void 
fcheb (x, pl, nl)
     double x, pl[];
     int nl;
#endif
{
  int j;
  double twox, f2, f1, d;

  pl[1] = 1.0;
  pl[2] = x;
  if (nl > 2)
    {
      twox = 2.0 * x;
      for (j = 3; j <= nl; j++)
	  pl[j] = twox * pl[j - 1] -  pl[j - 2];
    }
}