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

   Otmar Stahl

   fit_back.c  

   interpolate background 

   1) compute median between orders
   2) smoothing splines along lines, compute background values at grid-points
   3) bi-cubic splines over the grid-points
   4) evaluate splines at all points

*/

/* general Midas includes */

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

/* FEROS specific includes */

#include <glsp.h>
#include <proto_nrutil.h>
#include <proto_mutil.h>
#include <echelle.h>

int fit_back
#ifdef __STDC__
(
 float xval[], float yval[], float **ya, float **y2a, 
 int npix[], int imno, int m, int n, 
 int nact, float **centers, int xysize[], int bgdist, int fibmode
 )
#else
     (
      xval, yval, ya, y2a, npix, imno, m, n, nact, 
      centers, xysize, bgdist, fibmode
      )
     int npix[], imno, nact, m, n, xysize[], bgdist, fibmode; 
     float **ya, **y2a, xval[], yval[], **centers;
#endif
{
  int i, ii, j, jj, k, iarr1, iy, actvals, center;
  int splwidth[2], splwin[2];
  float *arr1, *xpos, *ypos, *imagex1;
  double xxd, *xn, *fn, *w, *a, *b, *c, *d, ausg[3];
  double xposs, weight;
  int status;

  weight = 1.0e-06; /* smoothing parameter, this is a magic number ? */  

  splwin[0] = xysize[0];     
  splwin[1] = xysize[1];
  splwidth[0] = splwin[0] * 2 + 1;
  splwidth[1] = splwin[0] * 2 + 1;

  xn = dvector(0, 2*nact + 2);                  /* alloc mem for vars */
  fn = dvector(0, 2*nact + 2);
  w = dvector(0, 2*nact + 2);
  a = dvector(0, 2*nact + 2);
  b = dvector(0, 2*nact + 2);
  c = dvector(0, 2*nact + 2);
  d = dvector(0, 2*nact + 2);

  imagex1 = vector(0, npix[0] * splwidth[1]);

  arr1 = vector (0, splwidth[0] * splwidth[1] + 1);
  xpos = vector (1, 2*nact + 1);
  ypos = vector (1, 2*nact + 1);

  /* init various arrays */

  for(k = 1; k <= n; k++ )                         /* for each row */
    {
      j = (int) yval[k];

      /* read one line  */

      SCFGET(imno, (j - splwin[1]) * npix[0] + 1, npix[0] * splwidth[1], 
	     &actvals,  (char *) imagex1); /* data are at *imagex1 */

      /* loop over the orders */

      ii = 0;
      for (i = 1; i <= nact; i++)           /* for all orders  */
	{

	  /* between two orders if sufficient space */

	  if ( (i > 1) && (centers[i][j] - centers[i-1][j]) > bgdist)
	    {
	      xposs = (centers[i][j] + centers[i-1][j]) / 2;
	      center = (int) xposs;
	      
	      if (center > splwin[0] && center < npix[0] - splwin[0])
		{
		  ii++;
		  iarr1 = 0;
		  for(jj = -splwin[0]; jj <= splwin[0]; jj++)
		    {
		      for (iy = 0; iy < splwidth[1]; iy++)
			arr1[++iarr1] = imagex1[center + jj + iy * npix[0]];
		      
		    }
		  xpos[ii] = xposs;
		  ypos[ii] = select_pos((int) iarr1/2 , iarr1 , arr1);
		}
	    }

	  /* between the fibers for fiber mode 2 */
	  
	  if (fibmode > 1)
	    {
	      xposs = (centers[i][j]);
	      
	      center = (int) xposs;
	      
	      /* 
		 populate array for median filtering,
		 this needs more work, e.g.
		 median over a number of lines
	      */ 
	      
	      if (center > splwin[0] && center < npix[0] - splwin[0])
		{
		  ii++;
		  iarr1 = 0;
		  for(jj = -splwin[0]; jj <= splwin[0]; jj++)
		    {
		      for (iy = 0; iy < splwidth[1]; iy++)
			arr1[++iarr1] = imagex1[center + jj + iy * npix[0]];
		      
		    }
		  xpos[ii] = xposs;
		  ypos[ii] = select_pos((int) iarr1/2 , iarr1 , arr1);
		}
	    }
	}
      
      for(i = 0; i < ii; i++)
	{
	  xn[i] = (double) xpos[i + 1];
	  fn[i] = (double) ypos[i + 1];
	  w[i] = weight;
	}

      /* smoothing spline in X-direction and compute new grid:
	 see Engeln-Muellges & Reuther, 1990; BI Wissenschaftsverlag Mannheim;
	 Formelsammlung zur numerischen Mathematik mit C-Programmen, page 236f
      */
    
      status = glspnp(ii - 1, xn, fn, w, 2, 0.0, 0.0, a, b, c, d);

      for(i = 1; i <= m; i++)
	{
	  xxd =  spval(ii - 1, (double) xval[i], a, b, c, d, xn, ausg);
	  ya[i][k] = (float) xxd;

	  /* smoothed function values at x[i = 1 to m] are in ya[i][k] for any row 
	     k (= 1 .. n) */
	}
    }

  /* to compute second derivatives, see Recipes, pg. 128 - we need
     an rectangular grid.*/

  splie2(xval, yval, ya, m, n, y2a); 
  
  free_dvector(xn, 0, 2*nact + 2);
  free_dvector(fn, 0, 2*nact + 2);
  free_dvector(w, 0, 2*nact + 2);
  free_dvector(a, 0, 2*nact + 2);
  free_dvector(b, 0, 2*nact + 2);
  free_dvector(c, 0, 2*nact + 2);
  free_dvector(d, 0, 2*nact + 2);

  free_vector(imagex1, 0, npix[0] * splwidth[1]);
  free_vector(arr1, 0, (splwidth[0] * splwidth[1] + 1));
  free_vector(xpos, 1,2*nact + 1);
  free_vector(ypos, 1, 2*nact + 1);

  return 0;
}