/* @(#)rebin.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:59 */
/*

  Thomas Szeifert, Anton Malina

  rebin.c              

  rebin to constant wavelength steps

*/

/* system includes */

#include <math.h>

/* FEROS specific includes */

#include <proto_dsplint.h>
#include <rebin.h>

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

{
    int i, j, k, idx, ix = 0;
    int j1, j2;
    double x1, x2;
    double a, b, c, y;
    double d, dd, ddd;
    double rx = -0.5;
    double yval[3];

    int istart = 2;
    int Imin = MAXINTENS;
    int Imax = MININTENS;

    if ( Method == SPLINE ) {	/* No flux conservation ! */
    	for ( i = 0; i < nout; i++ ) {
    	    yout[i] = dsplint( 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, step, scaling);
    j1 = NINT(x1);

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

	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 
#ifdef __STDC__
(
 double rx, double *win, double *wout, int nin, 
 int nout, double *x, int *ix , double step, char scaling
 )
#else
(
 rx, win, wout, nin, nout, x, ix, step, scaling
 )
     double rx, *win, *wout, *x, step;
     int    nin, nout, *ix;
     char scaling;
#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;

    switch (scaling)
      {
      case 'I':
	{
	  wave = wout[irx] + (rx - irx) * step;
	  break;
	}
      case 'O':
	{
	  wave = wout[irx] * exp((rx - irx) * step);
	  break;
	}
      }
    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);
}