/* @(#)spffit.c	17.1.1.1 (ES0-DMD) 01/25/02 17:54:01 */
/*===========================================================================
  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
===========================================================================*/

/* @(#)spffit.c	17.1.1.1  (ESO)  01/25/02  17:54:01 */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        spffit.c
.MODULE       main program -- spffit.exe
.LANGUAGE     C
.AUTHOR       Cristian Levin - ESO La Silla
.PURPOSE      This programs makes a polynomial fitting over
              three posible spaces:
                 - ratio / wave
                 - colour / wave
                 - colour / frequency
.KEYWORDS     response curve, flux calibration
.VERSION 1.0  1-Apr-1991   Implementation
.ENVIRONMENT  UNIX
------------------------------------------------------------*/

#include <stdio.h>
#include <midas_def.h>

#define WAVE_RATIO	0
#define WAVE_COLOUR	1
#define FREQ_COLOUR	2

float *fvector();
double *dvector();
double eval_dpoly();
void dpoly();

main()
{
    int  unit;		/* useless */
    int nulval, actval;	/* useless */
    int  sortcol, aw, ar, ncols, nrows;
    float *af, xout, *yout;
    double *ad, *xin, *yin;
    int i, npix;
    int c_xin, c_yin;
    int idf, idt;
    float start, step;
    float fstart, fstep, flast;
    int param[2], fit_deg, fit_type;
    char image[80], table[80];

    SCSPRO( "SPFFIT" );

    SCKGETC( "IN_A",  1, 80, &actval, table );
    SCKGETC( "OUT_A", 1, 80, &actval, image );
    SCFOPN( image, D_R4_FORMAT, 0, F_IMA_TYPE, &idf );
    TCTOPN( table, F_I_MODE, &idt );

    SCKRDI( "INPUTI", 1, 2, &actval, param, &unit, &nulval );
    fit_type = param[0];
    fit_deg  = param[1]+1;
    SCDRDI( idf, "NPIX", 1, 1, &actval, &npix, &unit, &nulval );
    SCDRDR( idf, "START", 1, 1, &actval, &start, &unit, &nulval );
    SCDRDR( idf, "STEP", 1, 1, &actval, &step, &unit, &nulval );

    TCIGET( idt, &ncols, &nrows, &sortcol, &aw, &ar );

    if ( fit_type == WAVE_RATIO )
        TCCSER( idt, ":RATIO", &c_yin );
    else 
        TCCSER( idt, ":COLOUR", &c_yin );
    if ( fit_type == FREQ_COLOUR )
    	TCCSER( idt, ":FREQ", &c_xin );
    else
    	TCCSER( idt, ":WAVE", &c_xin );

    xin  = dvector(0, nrows - 1);
    yin  = dvector(0, nrows - 1);
    for ( i = 0; i < nrows; i++ ) {
	TCERDD( idt, i+1, c_xin, &xin[i], &nulval );
	TCERDD( idt, i+1, c_yin, &yin[i], &nulval );
    }
    ad = dvector( 1, fit_deg );
    af = fvector( 1, fit_deg );
    lfit( xin-1, yin-1, nrows, ad, fit_deg, dpoly );	/* fitting */
    for ( i = 1; i <= fit_deg; i++ )
	af[i] = ad[i];

    yout  = fvector(0, npix - 1);

    fstart = 1.0 / start;
    flast  = 1.0 / (start + step * (npix-1));
    fstep  = (flast - fstart) / (npix-1);
    for ( i = 0; i < npix; i++ ) {
	if ( fit_type == FREQ_COLOUR )
	    xout = fstart + fstep * i;
	else
	    xout = start + step * i;
	yout[i] = eval_dpoly( xout, ad, fit_deg );
    }

    if ( fit_type == FREQ_COLOUR ) {
    	SCDWRR( idf, "START", &fstart, 1, 1, &unit );
    	SCDWRR( idf, "STEP", &fstep, 1, 1, &unit );
    }
    SCFPUT( idf, 1, npix, (char *)yout );

	/* save coefficients */
    SCDWRR( idf, "COEFF", af+1, 1, fit_deg, &unit );

    free_fvector(xout, 0, npix - 1);
    free_fvector(yout, 0, npix - 1);

    SCSEPI();
}
 
void dpoly( x, p, np )  /* (p[i] <-- x**(i-1)), i = 1,np) */
double x, p[];
int np;
{
    int j;

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