/* @(#)mosdisp.c	8.3 (ESO-IPG) 11/23/94 10:06:03 */
/* @(#)mosdisp.c	5.1.1.5 (ESO-IPG) 4/20/93 18:31:39 */
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* .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 (moscalib.c, mosrebin.c) do not know which representation
   is used (Legendre polynomials, or any other), nor how it is stored (which
   descriptor, meaning of coefficients). This design enables to replace 
   this module mosdisp.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>

#define PI 3.141592653589793
#define MAXNCOE 100

/* #define  FUNCTION        fpoly */
#define  FUNCTION        fleg
void fleg(double, double [], int);
void fleg_2D(double, double, double [], int);
void lfit2(double [], double [], double [], int, double [],
           int, int [], int, double **, double *,
           void (*funcs)(double, double *, int));
void     fpoly(double, double [], int); 
int      *ivector(int, int);
double   **dmatrix(int, int, int, int), *dvector( int, int);
void     free_dmatrix(double **, int, int, int, int);
void     free_dvector(double *, int, int), free_ivector(int *, int, int);

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

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

int	 tide, colslit, colline, coly, colrms, colcoef[100];
int	 start_index, nbline;

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

void mos_savedisp(double save[])

{

  int i;

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

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

void printdisp( void )

{

  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");

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

void setrefdeg(int deg)

{

refdeg = deg;
maxcoef = refdeg+1;

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

void setdisp(int deg, double coefs[])


{

int i;

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

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

double fit_disp(int *ndata, int *deg, double x[], double l[])

{

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;


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

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

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);

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

}
*/

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

double mos_fit_disp(int *ndata, int *deg, double x[], 
                    double l[], double *chi)

{

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;


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

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

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);
*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];
    }
}
/* --------------------------------------------------------------------------*/

void eval_disp(double x[], double l[], int n)

{

int   i, icoef;
double     xp[100];

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];
     }
}
/* -------------------------------------------------------------------------- */

void initdisp(char name[], char mode[], int start)

{

int      icoef;
char     colnam[30], numb[10];
int actvals, null;
char     *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);
  TCIGET(tide, &ncol, &nrow, &nsort, &allcol, &allrow);
  fdeg=refdeg; 
  ncoef=maxcoef; 
  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]);
   }
}

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

void mos_initdisp(char name[], char mode[], int  start)

{

int      icoef;
char     colnam[30], numb[10];
int actvals, null;
char     *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]);
   }
}
/* -------------------------------------------------------------------------- */

void finishdisp( void )

{

char 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);

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

double readdisp(int y)

/* 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);

return(pixbin);

}
/* -------------------------------------------------------------------------- */
/*
void writedisp(int line, int ypix, double y)

{

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]);
}
*/
/* -------------------------------------------------------------------------- */

void mos_writedisp(int line, int slit, int ypix, 
                   double y, int numrow, double rms) 

{

  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 mosdisp.c               */
/* (In all programs from Numerical Recipes, arrays start at index 1)  */
/**********************************************************************/

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

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))
{
	int k,kk,j,ihit,i;
	double ym,wt,sum,sig2i,**beta,*afunc;
	void gaussj(double **, int, double **, int ); 
     	void covsrt(double **, int, int [], int );
        void nrerror(char []);
        void free_dmatrix(double **, int, int, int, int);
        void free_dvector(double *, int, int);
	double **dmatrix(int, int, int, int),*dvector(int, int);

	beta=dmatrix(1,ma,1,1);
	afunc=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

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

void lfit2(double x1[], x2[], double y[], double sig[], int ndata, double a[],
           int ma, int lista[], int mfit, double **covar, double *chisq,
           void (*funcs)(double, double, double *, int))
{
	int k,kk,j,ihit,i;
	double ym,wt,sum,sig2i,**beta,*afunc;
	void gaussj(double **, int, double **, int ); 
     	void covsrt(double **, int, int [], int );
        void nrerror(char []);
        void free_dmatrix(double **, int, int, int, int);
        void free_dvector(double *, int, int);
	double **dmatrix(int, int, int, int),*dvector(int, int);

	beta=dmatrix(1,ma,1,1);
	afunc=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)(x1[i],x2[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)(x1[i],x2[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

void fpoly(double x, double p[], int np)

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

void fleg( double x, double pl[], int nl)

{
	int j;
	float 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;
		}
	}
}

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

double mos_fit_disp_2D(int *ndata, int *deg, double x[], double y[],
		       double l[], double *chi)

{

  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;
  

  if (*ndata < 2) {
    printf("Not enough lines (minimum is 2). \nNo dispersion relation computed\n");
    return(-2.);
  }
  
  if (fdeg < 1)  {
    printf("Degree : %d. No dispersion relation fitted\n",*deg);
    return(-1.);
  }
  
  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_2D(x,y,l,sig,*ndata,coef,ncoef,lista,ncoef,covar,chisq,fleg_2D);
  *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_2D (double x[], double y[], double l[], int n)
#else
void 
mos_eval_disp_2D (x, y, l, n)
     double x[], y[], l[];
     int n;
#endif

{

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

  SCKGETC("POLTYP", 1, 8, &actvals, poltyp);
  for (i = start_index; i < (n + start_index); i++)
    {
      l[i] = 0.;
      if (toupper(poltyp[0]) == 'L')
        fleg_2D (x[i], y[i], xp, ncoef);
      else 
	{
	  sprintf(text,"Chebichev is not yet implemented");
	  SCTPUT(text);
	}
	/*        fcheb (x[i], xp, ncoef);*/
      for (icoef = 1; icoef <= ncoef; icoef++)
	l[i] += coef[icoef] * xp[icoef];
    }
}

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

void mos_writedisp_2D(int line, double ystart, int slit, int ypix, 
		      double y, int numrow, double rms) 
{
  int icoef,iexp;
  double coef2D[maxcoef-5];
  
  /* make 1D coeficients from 2D coeficients -- last 5 depend on (y-ystart)*/
  for (icoef = 1; icoef <= maxcoef - 5; icoef++)
    coef2D = &coef[icoef];
  for (iexp = 1; iexp < 5; iexp++)
    coef2D[1] = coef2D[1] + &coef[maxcoef-5+iexp]*(y-ystart);
  coef2D[2] = coef2D[2] + &coef[maxcoef]*(y-ystart);
  
  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-5; icoef++)
    TCEWRD(tide, line, colcoef[icoef], coef2D[icoef]);
}

void fleg_2D( double x1, double x2, double pl[], int nl)
  
{
  /* Legendre polynom in x1 - polynom in y and one mixed polynom x1*x2*/
  int j;
  float twox,f2,f1,d;
  
  pl[1]=1.0;
  pl[2]=x1;
  if (nl > 2) {
    twox=2.0*x1;
    f2=x1;
    d=1.0;
    for (j=3;j<=nl-5;j++) {
      f1=d;
      d += 1.0;
      f2 += twox;
      pl[j]=(f2*pl[j-1]-f1*pl[j-2])/d;
    }
    pl[nl-5]=1.0;
    for (j=nl-4;j<nl;j++)
      pl[j] = x2*pl[j-1];
    pl[nl] = x1*x2;
  }
  
}