/* @(#)mosrebin1d.c	8.6 (ESO-IPG) 11/30/94 15:56:47 */
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
.IDENT        mosrebin1d.c
.MODULE       main program -- mosrebin1D.exe
.LANGUAGE     C
.AUTHOR       Cristian Levin - ESO La Silla
              Otmar Stahl, LSW, rebin "image" of extracted spectra !!!
.PURPOSE      linear, quadratic an spline rebinning.
.KEYWORDS     rebinning.
.VERSION 1.0  1-Apr-1991   Implementation
.ENVIRONMENT  UNIX
------------------------------------------------------------*/

#include <tbldef.h>
#include <midas_def.h>
#include <gl_defs.h>
#include <proto_spec.h>

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

#define MAXINTENS   99999.0
#define MININTENS  -99999.0 

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

char    Coerbr[61];
double  Starti, Stepi;
double  Starto, Endo, Stepo;
int     Method, NpixOut, Npix[2], kun, knul;

float   Cutsi[2];
float   Imax = MININTENS, Imin = MAXINTENS;
/* float   **Coefs; */
       

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

main ()

{

  char inframe[60], rootima[40], outima[60], output[61];
  int dunit;
  char number[6];

  int status, actvals, null;	/* Std Interfaces variables */
  int naxis, ypos[100], slit[100], i, id;
  int  dnull, imno, imno2;

  double *x, *win, *wout;
  float *yin, *yout, rebpar[3];

  SCSPRO ("mosrebin1d");

  status = SCKGETC ("IN_A", 1, 60, &actvals, inframe);
  if (status)
    SCTPUT ("Error while reading IN_A");

  status = SCKGETC ("IN_B", 1, 60, &actvals, Coerbr);
  if (status)
    SCTPUT ("Error while reading IN_B");

  status = SCKGETC ("OUT_A", 1, 40, &actvals, rootima);
  if (status)
    SCTPUT ("Error while reading OUT_A");

  status = SCKRDI ("INPUTI", 1, 1, &actvals, &Method, &kun, &knul);

  /* open "line-by-line" frame and read descriptors */

    if ( SCFOPN( inframe, D_R4_FORMAT, 0, F_IMA_TYPE, &imno ) != 0 ) {
        sprintf( output, "Frame %s invalid...", inframe );
        SCTPUT( output );
        return;
    }

    SCDRDI( imno, "NAXIS", 1, 1, &actvals, &naxis, &dunit, &dnull );
    SCDRDI( imno, "NPIX", 1, 2, &actvals, Npix, &dunit, &dnull );
    if ( naxis == 1 )
        Npix[1] = 1;
    SCDRDD( imno, "START", 1, 1, &actvals, &Starti, &dunit, &dnull );
    SCDRDD( imno, "STEP", 1, 1, &actvals, &Stepi, &dunit, &dnull );
    SCDRDR( imno, "LHCUTS", 1, 2, &actvals, Cutsi, &dunit, &dnull );

  /* ypos is the original y-position of the spectra (before extraction) */

  SCDRDI (imno, "YPOS", 1, 100, &actvals, ypos, &dunit, &dnull);
  SCDRDI (imno, "SLIT", 1, 100, &actvals, slit, &dunit, &dnull);

  if (naxis == 1)
    {
      Npix[1] = 1;
    }

/*    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, &actvals, rebpar, &dunit, &dnull );
    Starto = (double) rebpar[0];
    Endo   = (double) rebpar[1];
    Stepo  = (double) rebpar[2];
    NpixOut = (int) ((Endo - Starto) / Stepo + 0.5);

    TCTCLO( id );

    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]  = (double) (Starti + Stepi * i);
    for ( i = 0; i < NpixOut; i++ )
        wout[i] = (double) (Starto + Stepo * i);

     /* loop over all rows, and images  */

  for (i = 0; i < Npix[1]/2; i++)
    {

      strcpy (outima, rootima);
      sprintf (number, "%04i", i + 1);
      strcat (outima, number);
      sprintf (output,"%s", outima);
      SCTPUT(output);

    SCFGET(imno, i*Npix[0]+1, Npix[0], &null, (char *)yin);

    status = mos_readdisp(ypos[i],slit[i]);

      if (status == 0) {
	mos_eval_disp(x, win, Npix[0]);      
	
	SCIPUT(outima, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE, 
	       1, &NpixOut, &Starto, &Stepo, 
	       "artificial Image", "Angstrom", (char **)&yout, &imno2);

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

  finishdisp();

  /* Close everything and good bye */

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

  SCFCLO (imno), 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;
    /* double dx; */
    double a, b, c, y;
    double d, dd, ddd;
    double 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 */
}

#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,*win,*wout,*x;
  int    nin,nout,*ix;
#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 ( 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;
    }
	    
    /* interpolate to find X */
    *x = *ix + (wave - waux) / (win[ix2]-waux) * (ix2 - *ix);
}