/* @(#)rectify.c	17.1.1.1 (ESO-IPG) 01/25/02 17:52:00 */
/*  
    Otmar Stahl

    rectify.c 

    extract rectified profiles and fractional pixel offsets 
    there is no rebinning!

*/

/* system includes */

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

/* general Midas includes */

#include <midas_def.h>

/* FEROS specific includes */

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

int rectify
#ifdef __STDC__
(
 float **centers,  int imno, int npix[], int nact, double start0, 
 double step0, int width,  int off1, int off2,  int mode, char* out_image
)
#else
     (
      centers, imno, npix, nact, start0, step0, width, off1, off2, 
      mode, out_image
      )
     char *out_image;
     double start0, step0;
     float **centers;
     int imno, npix[], nact, width, off1, off2, mode;
#endif

{
  char c_name1[61], c_name2[61], d_name1[61], d_name2[61];
  char x_name1[61], x_name2[61];

  int actvals, i, jj, k, offset1, offset2, yposs, center, win;
  int *number1, *number2, npixo[2], npixl[2], dun;
  int imnoc1, imnoc2, imnod1, imnod2, imnox1, imnox2; 

  size_t name_len;

  float xpos, ypos, *line, ypos1;
  float *imagec1, *imagec2, *imaged1, *imaged2, *imagex1, *imagex2;

  double start[2], step[2];

  line = vector(0, npix[0]);

  number1 = ivector(1, nact);
  number2 = ivector(1, nact);

  start[0] = start0;
  step[0] = step0;
  start[1] = step[1] = 1.0;

  for (i = 1; i <= nact; i++)
    {
      number1[i] = number2[i] = 0;    /* initialize  arrays */
    }
  win = width / 2;

  npixo[0] = npix[1];
  npixo[1] = width * nact;

  npixl[0] = npixo[0];
  npixl[1] = nact;

  /* construct output file names */

  name_len = strcspn(out_image,".");
  out_image[name_len] = '\0';
  
  if (mode == 1)  /* one fiber */
    {
      strncpy (c_name1, out_image, name_len+1);
      strcat (c_name1, "_C");
      
      strncpy (d_name1, out_image, name_len+1);
      strcat (d_name1, "_D");

      strncpy (x_name1, out_image, name_len+1);
      strcat (x_name1, "_X");
    }
  else /* mode != 1, two fibers */
    {
      strncpy (c_name1, out_image, name_len+1);
      strcat (c_name1, "_C1");

      strncpy (d_name1, out_image, name_len+1);
      strcat (d_name1, "_D1");
      
      strncpy (x_name1, out_image, name_len+1);
      strcat (x_name1, "_X1");

      strncpy (c_name2, out_image, name_len+1);
      strcat (c_name2, "_C2");

      strncpy (d_name2, out_image, name_len+1);
      strcat (d_name2, "_D2");
      
      strncpy (x_name2, out_image, name_len+1);
      strcat (x_name2, "_X2");
    }
  

  /* create output frames */

  SCIPUT (c_name1, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	  2, npixo, start, step, "Cross", "no units", 
	  (char **) &imagec1, &imnoc1);
  SCIPUT (d_name1, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	  2, npixl, start, step, "Xdiff", "no units", 
	  (char **) &imaged1, &imnod1);
  SCIPUT (x_name1, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	  2, npixl, start, step, "Xpos", "no units", 
	  (char **) &imagex1, &imnox1);

  if (mode != 1)
    {
      SCIPUT (c_name2, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	      2, npixo, start, step, "Cross", "no units", 
	      (char **) &imagec2, &imnoc2);
      SCIPUT (d_name2, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	      2, npixl, start, step, "Xdiff", "no units", 
	      (char **) &imaged2, &imnod2);
      SCIPUT (x_name2, D_R4_FORMAT, F_O_MODE, F_IMA_TYPE,
	      2, npixl, start, step, "Xpos", "no units", 
	      (char **) &imagex2, &imnox2);
    }

  for (k = 0; k < npix[1]; k++)               /* for each row */
    {
      /* read one line  */

      SCFGET(imno, k * npix[0] + 1, npix[0], &actvals, (char *)line);

      for (i = 1 ; i <= nact; i++)       /* for each order */    
	{
	  offset1 = (i-1) * width * npixo[0];    /* row offset */
	  offset2 = (i-1) * npixl[0];

	  ypos = centers[i][k + 1];               /* y-pos in pixels */
	  xpos = k + 1;                           /* x-pos in pixels */

	  yposs = (int) (ypos + 0.5);       /* y-pos of next  pixel center */
	  ypos1 = ypos - yposs;             /* deviation from pixel center */
      
	  center =  yposs + off1;
	  if (center > win && center < npix[0] - win)
	    {
	      for (jj = -win ; jj <= win; jj++)
		{
		  imagec1[offset1 + number1[i] + (jj + win) * npixo[0]]
		    = line[center + jj];
		}
	      imaged1[offset2+number1[i]] = ypos1;
	      imagex1[offset2+number1[i]] = xpos;
	      number1[i]++;
	    }

	  /* two spectra */

	  if (mode != 1)
	    {
	      center =  yposs + off2;
	      if (center > win && center < npix[0] - win)
		{
		  for (jj = -win ; jj <= win; jj++)
		    {
		      imagec2[offset1 + number2[i] + (jj + win) * npixo[0]]
			= line[center + jj];
		    }
		  imaged2[offset2+number2[i]] = ypos1;
		  imagex2[offset2+number2[i]] = xpos;
		  number2[i]++;
		}
	    }
	}
    }
  
  SCDWRI(imnoc1, "NORDER", &nact, 1, 1, &dun);
  SCDWRI(imnoc1, "NLINES", &width, 1, 1, &dun);
  SCDWRI(imnoc1, "NPOINTS", &number1[1], 1, nact, &dun);

  SCDWRI(imnod1, "NPOINTS", &number1[1], 1, nact, &dun);
  SCDWRI(imnox1, "NPOINTS", &number1[1], 1, nact, &dun);
  
  SCFCLO(imnoc1);
  SCFCLO(imnod1);
  SCFCLO(imnox1);
  
  if (mode != 1)
    {
      SCDWRI(imnoc2, "NORDER", &nact, 1, 1, &dun);
      SCDWRI(imnoc2, "NLINES", &width, 1, 1, &dun);
      SCDWRI(imnoc2, "NPOINTS", &number2[1], 1, nact, &dun);
      
      SCDWRI(imnod2, "NPOINTS", &number2[1], 1, nact, &dun);
      SCDWRI(imnox2, "NPOINTS", &number2[1], 1, nact, &dun);

      SCFCLO(imnoc2);
      SCFCLO(imnod2);
      SCFCLO(imnox2);
    }
  
/* free memory */

  free_ivector(number1, 1, nact);          
  free_ivector(number2, 1, nact);
  free_vector(line, 0, npix[0]); 

  return 0;
}