/* @(#)spintegr.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
===========================================================================*/

/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        spintegr.c
.MODULE       main program -- spintegr.exe
.LANGUAGE     C
.AUTHOR       Cristian Levin - ESO La Silla
.PURPOSE      This program makes the flux integration.
.KEYWORDS     flux integration, flux calibration.
.VERSION 1.0  1-Apr-1991   Implementation
.ENVIRONMENT  UNIX
------------------------------------------------------------*/

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

#define EPS		0.1

int	Line = 1;	/* suppose 1-dim spectrum */
int	Npix[2];
int 	Ncols, Nrows;

char 	Image[81], Table[81], Outtab[81]; 
int 	IdIma, IdTab;

double	Start, Step;
float 	*Flux, *Sumf, *Wavet, *Intens, *Wavei;
float	*Binw;

float 	Rnull;

struct colids {
	int wavet, flux, binw;
} col;

void init_midas();
void init_purify();
void end_midas();
void read_values();
void calculate_flux();
void update_table();
float linint();

float *fvector();
int   *ivector();

main()
{
    init_purify();
    init_midas();

    read_values();
    calculate_flux();
    update_table();

    end_midas();
}

void calculate_flux()
{
    int i, j;
    int binmid, jmid = 0;
    int jlo, jhi;
    float histep, lostep, hiint;

    for ( i = 0; i < Nrows; i++ ) {
	binmid = Binw[i] / 2.0;

        for ( j = jmid; j < Npix[0]-1; j++ )
	    if ( Wavet[i] < Wavei[j+1] )
		break;

	if ( Wavet[i] + binmid > Wavei[Npix[0]-1] ||
	     Wavet[i] - binmid < Wavei[0] )
	    continue;

        jmid = j+1;

        for ( j = jmid+1; j < Npix[0]; j++ )
	    if ( Wavet[i] + binmid <= Wavei[j] ) {
		jhi = j-1;
		break;
	    }

        for ( j = jmid; j >= 0; j-- )
	    if ( Wavet[i] - binmid >= Wavei[j] ) {
		jlo = j+1;
		break;
	    }
 
        histep  = linint( Wavet[i] + binmid, Wavei[jhi], 0.0,
			 Wavei[jhi+1], Step );
        lostep  = linint( Wavet[i] - binmid, Wavei[jlo-1], Step,
			 Wavei[jlo], 0.0 );
	hiint   = linint( Wavet[i] + binmid, Wavei[jhi], Intens[jhi],
			 Wavei[jhi+1], Intens[jhi+1] );

        Sumf[i] = (Intens[jlo] * lostep +  hiint * histep) / Binw[i];
	for ( j = jlo+1; j <= jhi; j++ )
            Sumf[i] = Sumf[i] + Intens[j] * Step / Binw[i];

    }
}

float linint( x, x1, y1, x2, y2 )	/* Linear interpolation */
float x, x1, y1, x2, y2;
{
    float y;
    
    y = y1 + (y2 - y1)*(x - x1) / (x2 - x1);

    return(y);
}

void read_values()
{
    int i, j;
    int actval; 	/* actual values returned */
    int unit;		/* useless */
    int nulval;		/* useless */
    int  sortcol, aw, ar;
    char msg[80];
    float flux0, fluxn, wave_end;
    float first_wave, last_wave; 

    	/* read parameters */
    SCKGETC( "IN_A",  1, 80, &actval, Image );
    SCKGETC( "IN_B",  1, 80, &actval, Table );
    SCKGETC( "OUT_A", 1, 80, &actval, Outtab );

    	/* read info of table */
    if ( TCTOPN( Table, F_I_MODE, &IdTab ) != 0 ) {
	sprintf( msg, "Table %s invalid. Stop.", Table );
	SCTPUT( msg );
	end_midas();
    }
    TCIGET( IdTab, &Ncols, &Nrows, &sortcol, &aw, &ar );

    TCCSER( IdTab, ":WAVE", &col.wavet );
    TCCSER( IdTab, ":FLUX_W", &col.flux );
    TCCSER( IdTab, ":BIN_W",  &col.binw );
    if ( col.wavet == -1 || col.flux == -1 || col.binw == -1 ) {
	SCTPUT( "**Error**  Wrong table columns. Correct names are:" );
	SCTPUT( "         :WAVE    (wavelength)" );
	SCTPUT( "         :FLUX_W  (flux in wavelength units)" );
	SCTPUT( "         :BIN_W   (Bin width)" );
	end_midas();
    }
    Flux  = fvector(0, Nrows + 2);
    Wavet = fvector(0, Nrows + 2);
    Binw  = fvector(0, Nrows + 2);
    Sumf  = fvector(0, Nrows + 2);
    for ( i = 0; i < Nrows; i++ ) {
	TCERDR( IdTab, i+1, col.flux,  &Flux[i], &nulval );
	TCERDR( IdTab, i+1, col.wavet, &Wavet[i], &nulval );
	TCERDR( IdTab, i+1, col.binw,  &Binw[i], &nulval );
	Sumf[i] = Rnull;
        /* if (i == 0 || Binw[i] < minbin)  minbin = Binw[i]; */
    }

    	/* read info of image */
    if ( SCFOPN( Image, D_R4_FORMAT, 0, F_IMA_TYPE, &IdIma ) != 0 ) {
	sprintf( msg, "Image %s invalid. Stop.", Image );
	SCTPUT( msg );
	end_midas();
    }
    SCDRDI( IdIma, "NPIX", 1, 2, &actval, Npix, &unit, &nulval );
    SCDRDD( IdIma, "START", 1, 1, &actval, &Start, &unit, &nulval );
    SCDRDD( IdIma, "STEP", 1, 1, &actval, &Step, &unit, &nulval );

    Intens = fvector(0, Npix[0] - 1);
    Wavei  = fvector(0, Npix[0] - 1);

    SCFGET( IdIma, (Line-1)*Npix[0]+1, Npix[0], &nulval, (char *)Intens );

    for ( i = 0; i < Npix[0]; i++ )
	Wavei[i] = Start + Step * i;


	/* interpolation of the table */

    if ( Start + Binw[0] / 2.0 > Wavet[0] ) 
	for ( i = 1; i < Nrows; i++ ) {
	    /* EPS is used to ensure the value will be inside the interval */
    	    first_wave = Start + Binw[i] / 2.0 + EPS;
	    if ( first_wave < Wavet[i] ) {
		flux0 = linint( first_wave, Wavet[i-1], Flux[i-1],
				Wavet[i], Flux[i] ); 
		for ( j = Nrows; j > i; j-- ) {
		    Wavet[j] = Wavet[j-1];
		    Flux[j] = Flux[j-1];
		    Binw[j] = Binw[j-1];
		}
		Wavet[i] = first_wave;
		Flux[i] = flux0;
		break;
	    }
	}

    wave_end = Wavei[Npix[0]-1];
    if ( wave_end - Binw[Nrows-1] / 2.0 < Wavet[Nrows-1] ) 
	for ( i = Nrows-2; i > 0; i-- ) {
    	    last_wave = wave_end - Binw[i] / 2.0 - EPS;
	    if ( last_wave > Wavet[i] ) {
		fluxn = linint( last_wave, Wavet[i], Flux[i],
				Wavet[i+1], Flux[i+1] ); 
		for ( j = Nrows; j > i+1; j-- ) {
		    Wavet[j] = Wavet[j-1];
		    Flux[j] = Flux[j-1];
		    Binw[j] = Binw[j-1];
		}
		Wavet[i+1] = last_wave;
		Flux[i+1] = fluxn;
		break;
	    }
	}
}

void update_table()
{
    int    unit;		/* useless */
    int i, idcol[6], id=0, j = 1;
    float divf, wfirst, wlast, kk;
    double dnum;
    int nwaves;
    char str[80];

    if ( TCTOPN(Outtab, F_IO_MODE, &id ) != 0 ) {
	sprintf( str, "Table %s couldn't be opened. Stop.", Outtab );
	SCTPUT( str );
	end_midas();
    }
    TCCINI( id, D_R4_FORMAT, 1, "F10.1", "Angstrom", "WAVE", &idcol[0] );
    TCCINI( id, D_R4_FORMAT, 1, "F13.5", " ", "FLUX",   &idcol[1] );
    TCCINI( id, D_R4_FORMAT, 1, "F13.5", " ", "SUM",    &idcol[2] );
    TCCINI( id, D_R4_FORMAT, 1, "F13.5", " ", "RATIO",  &idcol[3] );
    TCCINI( id, D_R4_FORMAT, 1, "F13.5", " ", "COLOUR", &idcol[4] );
    TCCINI( id, D_R4_FORMAT, 1, "F13.5", " ", "FREQ",   &idcol[5] );
	
    for ( i = 0; i < Nrows; i++ )
	if ( Sumf[i] != Rnull ) {
	    divf   = Sumf[i] / Flux[i];
	    TCEWRR( id, j, idcol[0], &Wavet[i] );
	    TCEWRR( id, j, idcol[1], &Flux[i] );
	    TCEWRR( id, j, idcol[2], &Sumf[i] );
	    TCEWRR( id, j, idcol[3], &divf );
	    if ( j == 1 ) 
		wfirst = Wavet[i];
	    wlast = Wavet[i];
	    j++;
	}

    /* save first wave and no. of waves for doing b-spline fitting */

    /* This way eventually cuts the borders, so we will force the size of the image
    SCDWRR( id, "WSTART", &wfirst, 1, 1, &unit ); 
    nwaves = (wlast - wfirst + 1) / Step;
    SCDWRI( id, "NWAVES", &nwaves, 1, 1, &unit );
    */

    SCDWRD( id, "WSTART", &Start, 1, 1, &unit );
    SCDWRD( id, "WSTEP", &Step, 1, 1, &unit );
    SCDWRI( id, "NWAVES", &Npix[0], 1, 1, &unit );

    TCTCLO( id );
}

void init_midas()
{
    int inull;
    double dnull;

    SCSPRO( "FLUX" );
    TCMNUL( &inull, &Rnull, &dnull ); /* NULL values */
}

void end_midas()
{
    free_fvector(Flux,   0, Nrows + 2);
    free_fvector(Wavet,  0, Nrows + 2);
    free_fvector(Binw,   0, Nrows + 2);
    free_fvector(Sumf,   0, Nrows + 2);
    free_fvector(Intens, 0, Npix[0] - 1);
    free_fvector(Wavei,  0, Npix[0] - 1);

    SCSEPI();
}
void init_purify()
{
    strcpy(Outtab,"");
    strcpy(Image,"");
    strcpy(Table,"");
}