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

/* @(#)sprebin.c	17.1.1.1  (ESO)  01/25/02  17:56:15 */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        sprebin.c
.MODULE       main program -- sprebin.exe
.LANGUAGE     C
.AUTHOR       Cristian Levin - ESO La Silla
.PURPOSE      linear, quadratic an spline rebinning.
.KEYWORDS     rebinning.
.VERSION 1.0  1-Apr-1991   Implementation
.ENVIRONMENT  UNIX
------------------------------------------------------------*/

#include <stdio.h>
#include <math.h>

#include <midas_def.h>

#include <gl_defs.h>
#include <lndisp.h>

#define MAXINTENS	99999.0
#define MININTENS	-99999.0

#define LINEAR		0
#define QUADRATIC	1
#define SPLINE		2

/* Keywords */
int	Fitd;
int	Npix[2];

char 	Coerbr[81], FramIn[81], FramOut[81]; 
int 	IdIn, IdOut;

double	Starti, Stepi,
	Starto, Endo, Stepo;
float	Cutsi[2];
float	Imax = MININTENS,
	Imin = MAXINTENS;

int 	NpixOut;
int	Method;

void read_parameters();
void   update_frame();
double splint();
void rebin();
void closest_index();
void read_coefs();
void flip_double(), flip_float();

double  *dvector(), **dmatrix();
float   *fvector();

main()
{
    int  nulval; /* garbage */
    double *x,    *win, *wout;
    float  *yin,  *yout;
    int done = 1, slice;
    int i, j, fsize,f_in, f_out;
    char msg[80];
    double    pixel;

    SCSPRO("sprebin");
    read_parameters();

    read_coefs();  /* Now read only descriptor REBPAR */
    initdisp(Coerbr, "OLD",0);

    x    = dvector(0, Npix[0] - 1);
    win  = dvector(0, Npix[0] - 1);
    wout = dvector(0, NpixOut - 1);
    yin  = fvector(0, Npix[0] - 1);
    yout = fvector(0, NpixOut - 1);

    for ( i = 0; i < Npix[0]; i++ )
	x[i]  = (Starti + Stepi * i);
    for ( i = 0; i < NpixOut; i++ )
	wout[i] = (Starto + Stepo * i);


    f_out  = (Stepo < 0.) ? -1 : 1;
    if (f_out == -1) flip_double(wout, NpixOut);

    fsize = NpixOut * Npix[1];
    SCFCRE( FramOut, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE, fsize, &IdOut );

    slice = Npix[1]/5;

    for ( i = 1; i <= Npix[1]; i++ ) {
    	SCFGET( IdIn, (i-1)*Npix[0]+1, Npix[0], &nulval, (char *)yin );

        pixel = readdisp(i);
        eval_disp(x, win, Npix[0]);

        f_in  = (win[0] > win[Npix[0]-1]) ? -1 : 1;
        if (f_in == -1) {flip_double(win,Npix[0]); flip_float(yin,Npix[0]);}

	rebin( win, wout, yin, yout, Npix[0], NpixOut);

        if (f_out == -1) flip_float(yout, NpixOut);

    	SCFPUT( IdOut, (i-1)*NpixOut+1, NpixOut, (char *)yout );

        /* printf("Now processing line %d\n",i); */
        if ( i == slice * done && done != 5 ) {
	    sprintf( msg, "%3d%% done...", done * 20 );
	    SCTPUT(msg);
	    done++;
        }
    }
    sprintf( msg, "100%% done..." );
    SCTPUT(msg);

    finishdisp();

    update_frame();

    free_dvector(x, 0, Npix[0] - 1);
    free_dvector(win, 0, Npix[0] - 1);
    free_dvector(wout, 0, Npix[0] - 1);
    free_fvector(yin, 0, Npix[0] - 1);
    free_fvector(yout, 0, Npix[0] - 1);

    SCSEPI();
}

void rebin( win, wout, yin, yout, nin, nout)
double     *win, *wout;
float      *yin, *yout;
int        nin, nout;
{
    int i, j, k, idx, ix = 0;
    int j1, j2;
    double x1, x2, dx;
    double a, b, c, y;
    double d, dd, ddd;
    double idx_real, rx = -0.5;
    double yval[3];

    int istart = 2;

    if ( Method == SPLINE ) {	/* No flux conservation ! */
    	for ( i = 0; i < nout; i++ ) {
    	    yout[i] = splint( wout[i], win-1, yin-1, nin, &istart );
	    if ( yout[i] < Imin ) Imin = yout[i];
	    if ( yout[i] > Imax ) Imax = yout[i];
	}
	return;
    }

    closest_index( rx, win, wout, nin, nout, &x1, &ix );
    j1 = NINT(x1);

    for ( i = 0; i < nout; i++ ) {
	rx += 1.0;
    	closest_index( rx, win, wout, nin, nout, &x2, &ix );
    	j2 = NINT(x2);

	dx = x2 - x1;
	d = 0.0;
	if ( Method == QUADRATIC ) {
	    if ( i == 0 ) {
		for ( k = 0; k < 3; k++ ) {
		    idx = j1 + k - 1;
		    if ( idx < 0 || idx > nin - 1 )
			yval[k] = 0.0;
		    else
			yval[k] = yin[idx];
		}
		a = 0.5*(yval[0] + yval[2]);
		b = 0.5*(a - yval[0]);
		c = (13.0*yval[1] - a) / 12.0;
		a = (a - yval[1]) / 3.0;
		y = x1 - j1;
		dd = ((a*y + b)*y + c)*y - 0.25*b + 0.125*a + 0.5*c;
	    }

	    for ( k = 0; k < 3; k++ ) {
	        idx = j2 + k - 1;
	        if ( idx < 0 || idx > nin - 1 )
		    yval[k] = 0.0;
	        else
		    yval[k] = yin[idx];
	    }
	    a = 0.5*(yval[0] + yval[2]);
	    b = 0.5*(a - yval[0]);
	    c = 1.083333333333*yval[1] - 0.083333333333*a;
	    a = 0.333333333333*(a - yval[1]);
	    y = x2 - j2;
	    d -= dd;
	    dd = ((a*y + b)*y + c)*y - 0.25*b;
	    ddd = 0.125*a + 0.5*c;
	    d += dd - ddd;
	    dd += ddd;
	}
	else if ( Method == LINEAR ) {
	    if ( i == 0 ) {
		if ( j1 < 0 || j1 > nin - 1 )
		    dd = 0.0;
	 	else
		    dd = (j1 - x1 - 0.5)*yin[j1];
	    }

	    if ( j2 < 0 || j2 > nin - 1 )
	        ddd = 0.0;
	    else
	        ddd = yin[j2];

	    d += dd;
	    dd = (j2 - x2 - 0.5)*ddd;
	    d = d - dd - ddd;
	}

	for ( j = j1; j <= j2; j++ )
	    if ( j >= 0 && j < nin )
		d += yin[j];

        yout[i] = d;
        x1 = x2;
	j1 = j2;

	if ( yout[i] < Imin ) Imin = yout[i];
	if ( yout[i] > Imax ) Imax = yout[i];

    } /* end of loop */
}

void closest_index( rx, win, wout, nin, nout, x, ix )
double rx;
double *win, *wout;
int nin, nout;
double *x;
int *ix; 		/* index of initial guess */
{
    double wave, waux;
    int irx, ix2;

    irx = rx; /* integer part of rx */

    if ( irx < 0 )
	irx = 0;
    else if ( irx >= nout - 1 )
	irx = nout - 2;

    wave = wout[irx] + (rx - irx) * Stepo;

    waux = win[*ix];
    if ( win[0] <= win[nin-1] ) {  /* increasing wavelength */
        if ( waux > wave ) { /* initial guess is too high */
	    while ( *ix > 0 && waux > wave ) {
	        (*ix)--;
	        waux = win[*ix];
	    }
	    ix2 = *ix + 1;
        }
        else {		/* initial guess is too low */
	    while ( *ix < nin - 1 && waux < wave ) {
	        (*ix)++;
	        waux = win[*ix];
	    }
	    ix2 = *ix - 1;
        }
    }
    else { /* decreasing wavelength */
        if ( waux > wave ) { /* initial guess is too high */
	    while ( *ix < nin - 1 && waux > wave ) {
	        (*ix)++;
	        waux = win[*ix];
	    }
	    ix2 = *ix - 1;
        }
        else {		/* initial guess is too low */
	    while ( *ix > 0 && waux < wave ) {
	        (*ix)--;
	        waux = win[*ix];
	    }
	    ix2 = *ix + 1;
        }
    }
	    
    /* interpolate to find X */
    *x = *ix + (wave - waux) / (win[ix2]-waux) * (ix2 - *ix);
}

void update_frame()
{
    int   unit;		/* useless */
    float cuts[4];
    int npix_output;

    cuts[0] = Cutsi[0];
    cuts[1] = Cutsi[1];
    cuts[2] = Imin;
    cuts[3] = Imax;

    SCDCOP( IdIn, IdOut, ALL_DESC, NULL ); 
    SCDWRD( IdOut, "START", &Starto, 1, 1, &unit );
    SCDWRD( IdOut, "STEP", &Stepo, 1, 1, &unit );
    SCDWRI( IdOut, "NPIX", &NpixOut, 1, 1, &unit );
    SCDWRR( IdOut, "LHCUTS", cuts, 1, 4, &unit );
}

void read_parameters()
{
    int  unit;		/* useless */
    char msg[80];
    int actval,nulval;	/* useless */
    int naxis;

    SCKGETC( "IN_A", 1, 60, &actval, FramIn );
    SCKGETC( "IN_B", 1, 60, &actval, Coerbr );
    SCKGETC( "OUT_A", 1, 60, &actval, FramOut );
    SCKRDI( "INPUTI", 1, 1, &actval, &Method, &unit, &nulval );
    SCKRDI( "FITD", 1, 1, &actval, &Fitd, &unit, &nulval );

    if ( SCFOPN( FramIn, D_R4_FORMAT, 0, F_IMA_TYPE, &IdIn ) != 0 ) {
	sprintf( msg, "Frame %s invalid...", FramIn );
	SCTPUT( msg );
	return;
    }

    SCDRDI( IdIn, "NAXIS", 1, 1, &actval, &naxis, &unit, &nulval );
    SCDRDI( IdIn, "NPIX", 1, 2, &actval, Npix, &unit, &nulval );
    if ( naxis == 1 )
	Npix[1] = 1;
    SCDRDD( IdIn, "START", 1, 1, &actval, &Starti, &unit, &nulval );
    SCDRDD( IdIn, "STEP", 1, 1, &actval, &Stepi, &unit, &nulval );
    SCDRDR( IdIn, "LHCUTS", 1, 2, &actval, Cutsi, &unit, &nulval );
}

void read_coefs()
{
    int i, j, id;
    int  actval, nulval;	/* useless */
    float value, rebpar[3];
    int   unit;

    if ( ! file_exists(Coerbr, ".tbl") )
	SCETER( 9 , "Coefficients table couldn't be opened. Stop.\n");

    TCTOPN(Coerbr, F_IO_MODE, &id);

    SCDRDR( id, "REBPAR", 1, 3, &actval, rebpar, &unit, &nulval );
    Starto = (double) rebpar[0];
    Endo   = (double) rebpar[1];
    Stepo  = (double) rebpar[2];
    NpixOut = (int) ((Endo - Starto) / Stepo + 0.5);

    TCTCLO( id );

}


void flip_double (w,n)

double *w;
int    n;

  {

    int middle, i;
    double tmp;

    middle = (int) (n-1)/2;

    for (i=0; i<=middle; i++) {
         tmp=w[i];
         w[i] = w[n-i-1];
         w[n-i-1] = tmp;
       }

  }

void flip_float (w,n)

float *w;
int    n;

  {

    int middle, i;
    float tmp;

    middle = (int) (n-1)/2;

    for (i=0; i<=middle; i++) {
         tmp=w[i];
         w[i] = w[n-i-1];
         w[n-i-1] = tmp;
       }

  }