/* @(#)extract_spec.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:58 */
/* @(#)extract_spec.c	14.1.1.1 (ESO-IPG) 09/16/99 09:47:00 */
/* 
   Otmar Stahl, Anton Malina

   extract_spec.c  

   extract spectra with 
   standard or optimum extraction

*/

/* system includes */

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

/* general Midas includes */

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

/* FEROS specific includes */

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

int extract_spec
#ifdef __STDC__
(
 float imageo1[], int fit_deg,
 char mode, float v0, float gain, float thres, int niter,
 int nact, int nlines, int npts, int *npoints,
 int straight1, int xpos1, int dy1, char coeffile[],
 BOOLEAN get_spatial_profile, BOOLEAN write_imgs,
 char fit_image[], char mask_image[]
 )
#else
     (
      imageo1, fit_deg,
      mode, v0, gain, thres, niter,
      nact, nlines, npts, npoints, straight1, xpos1, dy1, coeffile,
      get_spatial_profile, write_imgs, fit_image, mask_image
      )
     int  nact, fit_deg, niter, nlines, npts, straight1, xpos1, dy1;
     float imageo1[], v0, gain, thres;
     char mode;
     int *npoints;
     char coeffile[], fit_image[], mask_image[];
     BOOLEAN get_spatial_profile, write_imgs;
#endif
{
  double x;
  double starto[2], stepo[2];
  double **p_coeff1;
  float *fxx, *ypos1;

  float mask_pix, profile_pix;
  float x_sum, raw_sum;
  float lhcuts[4];
  float *x1, *y1, *z1, *vari;
  float *raw, *msk_row;
  int i, j, jj, k, offset, rawoffset, xposoffset, actvals;
  int where;
  int imnomask, imnofit, tid1;
  int kunit, knul, naxis;
  int npixo[2], npixout[2], null[10], icol[20]; 
  char line[80];
  char cunit[64], ident[72];
  BOOLEAN do_optext;

  int masked_counter[25];

  /* array for counting how often a pixel in a specific position along the
     cross-order-profile has been masked */
  
  int total_masked;

  /* total number of masked pixels (at any position along the COP) */


  npixout[0] = npts; npixout[1] = nlines * nact;
  naxis = 2;

  lhcuts[0] = lhcuts[2] = 0.0;
  lhcuts[1] = lhcuts[3] = 10000.0; 

  msk_row = vector(0, nlines);
  raw = vector(0, npts * nlines);
  fxx = vector(0, npts);
  ypos1 = vector(0, npts);
  x1 = vector(0, nlines);                 /* alloc mem for var x1 */
  y1 = vector(0, nlines);                  /* alloc mem for var y1 */
  z1 = vector(0, data_per_pixel * nlines);       /* alloc mem for var z1 */
  vari = vector(0, nlines);
  p_coeff1 = dmatrix(1, data_per_pixel, 1, fit_deg * nlines * nact);


  strcpy(cunit, "ADU             x-axis [mm]     y-axis [mm]");

  /* create images that contain the fitted profile and the masked pixels */
  /* and write descriptors */

  SCDRDD(straight1, "START", 1, 2, &actvals, starto, &kunit, &knul);
  SCDRDD(straight1, "STEP", 1, 2, &actvals, stepo, &kunit, &knul);
  
  strcpy(ident, "fitted FEROS frame ");

  SCFCRE(fit_image, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE, nlines * npts * nact,
	 &imnofit);
  SCDWRI(imnofit, "NAXIS", &naxis, 1, 1, &kunit);
  SCDWRI(imnofit, "NPIX", npixout, 1, 2, &kunit);
  SCDWRD(imnofit, "START", starto, 1, 2, &kunit);
  SCDWRD(imnofit, "STEP", stepo, 1, 2, &kunit);
  SCDWRR(imnofit, "LHCUTS", lhcuts, 1, 4, &kunit);
  SCDWRC(imnofit, "IDENT", 1, ident, 1, 72, &kunit);
  SCDWRC(imnofit, "CUNIT", 1, cunit, 1, 64, &kunit);
  strcpy(ident, "masked pixels in FEROS spectrum ");

  SCFCRE(mask_image, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE, nlines * npts * nact,
	 &imnomask);
  SCDWRI(imnomask, "NAXIS", &naxis, 1, 1, &kunit);
  SCDWRI(imnomask, "NPIX", npixout, 1, 2, &kunit);
  SCDWRD(imnomask, "START", starto, 1, 2, &kunit);
  SCDWRD(imnomask, "STEP", stepo, 1, 2, &kunit);
  SCDWRR(imnomask, "LHCUTS", lhcuts, 1, 4, &kunit);
  SCDWRC(imnomask, "IDENT", 1, ident, 1, 13, &kunit);
  SCDWRC(imnomask, "CUNIT", 1, cunit, 1, 64, &kunit);


  /* populate arrays */

  total_masked = 0;
  for(k = 0; k < nlines; k++)
    masked_counter[k] = 0;

  if(toupper(mode == 'S')) /* standard extraction */
    {
      sprintf(line, "standard extraction without cosmic detection");
      SCTPUT(line);
      for (j = 0; j < nact; j++)           /* orders */
	{
	  rawoffset = npts * nlines * j + 1;

	  SCFGET(straight1, rawoffset, npts * nlines, &actvals, (char *) raw);

	  sprintf(line, "order %i of %i", j + 1, nact);
	  SCTPUT(line);

	  /* reads frame into mem at *raw */

	  offset = npts * (nact - j - 1);
	  xposoffset = npts * j + 1 ;
	      
	  SCFGET(xpos1, xposoffset, npts, &actvals, (char *) fxx); 

	  for (i = 0; i < npts; i++) 
	    imageo1[offset + i] = 0.0;

	  for (i = 0; i < npoints[j+1]; i++)         /* rows */
	    {
	      for (jj = 0; jj < nlines; jj++)
		{
		  imageo1[offset + ((int) fxx[i])-1] += raw[jj * npts + i];
		}
	    }
	}
    }
  else  /* mode == 'O' || mode == 'M' */
    {
      do_optext = (toupper(mode) == 'O');
      if (do_optext)
	sprintf(line, "optimum extraction with cosmic detection");
      else
	sprintf(line, "standard extraction with cosmic detection");
      SCTPUT(line);

      /* determine spatial profile */
      
      if (get_spatial_profile) 	                  
	{
	  /* initialize table */
	  TCTINI (coeffile, F_TRANS, F_O_MODE, 3*fit_deg, nact, &tid1); 

	  /* initialize columns for coefficients of the fit polynomial */

	  for (i = 1; i <= fit_deg; i++)         
	    {
	      /* 0.0 +- 0.05 Pixels off center */
	      sprintf (line, "FIT%04i", i);
	      TCCINI (tid1, D_R8_FORMAT, 1, "E14.7", "None", line, &icol[i-1]);
	    }
	  /* 0.3 +- 0.05 Pixels off center */
	  for (i = fit_deg + 1; i <= 2 * fit_deg; i++)
	    {
	      sprintf (line, "FIT2_%04i", i - fit_deg);
	      TCCINI (tid1, D_R8_FORMAT, 1, "E14.7", "None", line, &icol[i-1]);
	    }
	  /* -0.3 +- 0.05 Pixels off center */
	  for (i = 2 * fit_deg + 1; i <= 3 * fit_deg; i++) 
	    {
	      sprintf (line, "FIT3_%04i", i - 2 * fit_deg);
	      TCCINI (tid1, D_R8_FORMAT, 1, "E14.7", "None", line, &icol[i-1]);
	    }


	  extract_profiles(straight1, xpos1, dy1,
			   nact, nlines, npts, npoints,
			   v0, gain, thres,
			   tid1, icol, fit_deg);
	  TCTCLO(tid1);             /* close table */
	}

      /* here for reading cross-dispersion data */
      
      TCTOPN (coeffile, F_I_MODE, &tid1);         /* open table file */

      for (i = 1; i <= data_per_pixel * fit_deg; i++)
	icol[i - 1] = i;                           /* populate array */

      for (i = 0; i < nact; i++)                   /* orders */  
	{
	  offset = i * fit_deg * nlines + 1;
	  for (j = 1; j <= nlines; j++)
	    {
	      TCRRDD(tid1, i * nlines + j, fit_deg, icol, 
		     &p_coeff1[1][offset + (j - 1) * fit_deg], null);
	      TCRRDD(tid1, i * nlines + j, fit_deg, &icol[fit_deg], 
		     &p_coeff1[2][offset + (j - 1) * fit_deg], null);
	      TCRRDD(tid1, i * nlines + j, fit_deg, &icol[2 * fit_deg],
		     &p_coeff1[3][offset + (j - 1) * fit_deg], null);

	    } /* read table row as double into arrays p_coeff1 */
	}

      TCTCLO(tid1);             /* close table */
      npixo[0] = nlines;
      npixo[1] = nact * (data_per_pixel * (nlines + 1));


      for(j = 0; j < nact; j++)        /* orders */
	{
          sprintf(line, "order %i of %i", j + 1, nact);
	  SCTPUT(line);

	  rawoffset = npts * nlines * j + 1;
	  SCFGET(straight1, rawoffset, npts * nlines, &actvals, (char *) raw);

	  /* reads frame into mem at *raw */

	  for (i = 0; i < npts; i++)
	    imageo1[npts * (nact - j - 1) + i] = 0.0;

	  offset = npts * j + 1 ;

	  SCFGET(xpos1, offset, npts, &actvals, (char *) fxx);

	      /* read ypos1 from file */

	  SCFGET(dy1, offset, npts, &actvals, (char *) ypos1); 

	  for (i = 0; i < npoints[j+1]; i++)                   /* rows */
	    {
	      /* deviation from pixel center */

	      x = (double) fxx[i];
	      offset = j * fit_deg * nlines + 1;
	      for (jj = 0; jj < nlines; jj++)
		{
		  for(k = 0; k < data_per_pixel; k++)
		    {
		      z1[jj * data_per_pixel + k] = 
			(float) eval_dpoly(x, 
					   &p_coeff1[data_per_pixel-k]
					   [offset + jj * fit_deg - 1], 
					   fit_deg);		    
		    }
		  y1[jj] = raw[jj * npts + i];
		}

	      /* interpolate spatial profile, z1 in x1 out */

	      interpolate(z1, data_per_pixel * nlines, ypos1[i], x1, nlines);  

	      /* normalize interpolated profile */

	      x_sum = 0.0;
	      for(jj = 0; jj < nlines; jj++)
		{
		  x_sum += x1[jj];
		}
	      for(jj = 0; jj < nlines; jj++)
		x1[jj] /= x_sum;
	      
	      /* now do optimum extraction */
 
	      imageo1[npts * (nact-j-1) + ((int) fxx[i]) - 1 ] = 
		opt_ext(x1, y1, vari, nlines, niter,
			v0, gain, thres,
			0.2, do_optext, msk_row, 
			masked_counter, total_masked);
	      
	      if (write_imgs)
		{
		  x_sum = 0.0;
		  raw_sum = 0.0;
		  for (jj = 0; jj < nlines; jj++)
		    {
		      mask_pix = 2000.0 * msk_row[jj];
		      profile_pix = 
			x1[jj] * imageo1[npts * (nact-j-1) + 
					((int) fxx[i]) - 1];
		      where = npts * (nlines * j + jj) + i + 1;
		      SCFPUT(imnofit, where, 1, (char *) &profile_pix);
		      SCFPUT(imnomask, where, 1, (char *) &mask_pix);
		    }

		}
	    } /* end for i */
	}   /* end for j */
    }     /* end else */

  free_vector(raw, 0, npts * nlines);
  free_vector(msk_row, 0, nlines);
  free_vector(fxx, 0, npts);
  free_vector(ypos1, 0, npts);
  free_vector(x1, 0, nlines);                 
  free_vector(y1, 0, nlines);                 
  free_vector(z1, 0, data_per_pixel * nlines); 
  free_vector(vari, 0, nlines);
  free_dmatrix(p_coeff1, 1, data_per_pixel, 1, fit_deg * nlines * nact);

  SCFCLO(imnofit);
  SCFCLO(imnomask);

  if(toupper(mode != 'S')) /* standard extraction */
  {
    for(k = 0; k < nlines; k++)
      {
	sprintf(line, "pixel %i has been masked %i times", 
		k, masked_counter[k]);
	SCTPUT(line);
      }
  }
  return 0;
}