/* @(#)mosrebin.c	8.8 (ESO-IPG) 12/13/94 09:47:01 */
/* @(#)mosrebin.c	6.1.1.4  (ESO)  11/26/93  10:48:49 */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        mosrebin.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 <proto_spec.h>

#include <proto_nrutil.h>
#include <proto_splint.h>

#include <mutil.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[61], FramIn[61], FramOut[61]; 
int 	IdIn, IdOut;

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

int 	NpixOut;
int	Method;

#ifdef __STDC__
void read_parameters(void);
void update_frame(void);
void rebin( double *win, double *wout, float *yin, float *yout, 
            int nin, int nout );
void closest_index( double, double *, double *, int, 
                    int, double *, int * );
void read_coefs(void);
#else
void read_parameters();
void update_frame();
void rebin();
void closest_index();
void read_coefs();
#endif

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

/*-----------------for table MOS-------------------------------------------*/
    int icol[3],tid,iav, slit[100], upper[100], lower[100], null[3];
    int col, ncol, nslits, acol, arow, nsort, select, ii, status;
    float moscol[3];
    char table[60];

    SCSPRO("sprebin");
    read_parameters();

/*-----------------read table MOS------------------------------------------*/

  SCKGETC ("IN_C", 1, 60, &iav, table);

  TCTOPN (table, F_I_MODE, &tid);
  TCIGET (tid, &col, &ncol, &nsort, &acol, &arow);
 
  TCCSER (tid, ":slit", &icol[0]);
  TCCSER (tid, ":ystart", &icol[1]);
  TCCSER (tid, ":yend", &icol[2]);

  nslits = 0;
  for (i = 1; i <= ncol; i++)
    {
      TCSGET(tid, i, &select);
      if (select) {
        TCRRDR (tid, i, 3, icol, moscol, null);
        slit[nslits] = (int) moscol[0];
        lower[nslits] =  (int) ((moscol[1]-Starti[1])/Stepi[1]) + 1;
        upper[nslits] =  (int) ((moscol[2]-Starti[1])/Stepi[1]) + 1;
        nslits++;
      }
    }

  TCTCLO (tid);

/*-------------------end table MOS-----------------------------------*/

    read_coefs();  /* Now read only descriptor REBPAR */
    mos_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[0] + Stepi[0] * i);
    for ( i = 0; i < NpixOut; i++ )
	wout[i] = (Starto + Stepo * i);

    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++ ) {                */

      for ( ii = 0; ii < nslits; ii++) {
	for (i = lower[ii]; i <= upper[ii]; i++) {

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

	  status = mos_readdisp(i,slit[ii]);
	  
	  if (status == 0) {

	  mos_eval_disp(x, win, Npix[0]);

	  rebin( win, wout, yin, yout, Npix[0], 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, NpixOut - 1);
    free_fvector(yin, 0, Npix[0] - 1);
    free_fvector(yout, 0, NpixOut - 1);

    SCSEPI();
}

#ifdef __STDC__
void rebin( double *win, double *wout, float *yin, float *yout, 
            int nin, int nout )
#else
void rebin( win, wout, yin, yout, nin, nout )
double     *win, *wout;
float      *yin, *yout;
int        nin, nout;
#endif

{
    int i, j, k, idx, ix = 0;
    int j1, j2;
    double x1, x2, rx = -0.5;
    double a, b, c, y;
    double d, dd, ddd;
/*    double idx_real, dx; */
    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 */
}
#ifdef __STDC__
void closest_index( double rx, double *win, double *wout, int nin, 
                   int nout, double *x, int *ix )
#else
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 */
#endif

{
    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);
}
#ifdef __STDC__
void update_frame(void)
#else
void update_frame()
#endif
{
    int unit;		/* useless */
    float cuts[4];

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

#ifdef __STDC__
void read_parameters(void)
#else
void read_parameters()
#endif
{
    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, 2, &actval, &Starti[0], &unit, &nulval );
    SCDRDD( IdIn, "STEP", 1, 2, &actval, &Stepi[0], &unit, &nulval );
    SCDRDR( IdIn, "LHCUTS", 1, 2, &actval, Cutsi, &unit, &nulval );
}
#ifdef __STDC__
void read_coefs(void)
#else
void read_coefs()
#endif

{
    int id;
    int actval, nulval;	/* useless */
    float 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", 1L, 3L, &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 );

}