/* @(#)echcalibrate.c	17.1.1.1 (ESO-IPG) 01/25/02 17:52:20 */
/* 

   Thomas Szeifert, Anton Malina, Otmar Stahl

   echcalibrate.c 

   compute globale fit for echelle dispersion

   Global formula completely revised:
   
   Now a 2-D polynomial (in x and m) is used to fit lambda * m
   Works much better and gives residuals < 0.003 A rms

*/

/* 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 <proto_nrutil.h>
#include <proto_fitnol.h>
#include <proto_mutil.h>
#include <lnmatch.h>
#include <2Dstuff.h>
#include <misc.h>
#include <echelle.h>

int main()
{

  /* general */

  int i, j, actvals, kun, knul, reset, poltyp;
  char line[80];

  /* io-related */

  int tid1, tid2, tid3, icol1[12], icol2[3];
  int ncol, nrow, nsort, acol, arow;
  int mcol, mrow, msort;

  int inull;
  float fnull;
  double dnull;
  
  int *rej;
  char do_wlc[2];
  char in_table1[60], in_table2[60], in_table3[60];
  char lbl[7];
  double *x1, *yabs, *yy1, *w1, *wc1, *r1, *w2, *f2, xrot, xcos, xsin;
  float fparam[8], xbuff, ybuff, ybuff1;
  double wbuff1;

  int nparams, buffdeg, norder, isel, iall;  
  int nmatch, prevnmatch = -1;
  double *xtmp, *itmp, *ytmp;
  double *dcoef, *buffcoef;
  int verif = 1;
  double alpha, stdres = 0, tol;
  double lambda_central, lambda, beta0, beta;
  double angle, xoffset, thres, g, f, theta_blaze;
  float peak;

  SCSPRO("echcalibrate");
  TCMNUL(&inull, &fnull, &dnull);
  SCKGETC("IN_A",   1, 60, &actvals, in_table1);
  SCKGETC("IN_C",   1, 60, &actvals, in_table2);
  SCKGETC("IN_D",   1, 60, &actvals, in_table3);
  SCKGETC("INPUTC", 1, 1, &actvals, do_wlc);
  SCKRDR("INPUTR",  1, 8,  &actvals, fparam, &kun, &knul);

  alpha = (double) fparam[0];
  tol = (double) fparam[1];
  thres = (double) fparam[2];
  angle = (double) fparam[3];          /* rotation angle (degrees) */
  xoffset = (double) fparam[4];        /* xoffset (mm) */
  g = (double) fparam[5];              /* grating constant (lines per mm) */  
  f = (double) fparam[6];              /* focal length of camera (mm) */
  theta_blaze = (double) fparam[7];    /* blaze angle (degrees) */
  nparams = 25;                        /* #parameters for OS formula */
  dcoef = dvector (1, nparams);

  SCTPUT ("identify spectral lines");
  SCTPUT ("-----------------------");
  sprintf (line, "observed line table:     %s", in_table1);
  SCTPUT (line);
  sprintf (line, "laboratory line catalog: %s", in_table2);
  SCTPUT (line);
  sprintf (line, "calibration method:     %c", toupper(do_wlc[0]));
  SCTPUT (line);
  SCTPUT ("------------");


  /* READ THE REFERENCE CATALOG 
     of laboratory wavelengths and peak intensities */

  TCTOPN(in_table2, F_I_MODE, &tid2);
  TCIGET(tid2, &ncol, &nrow, &nsort, &acol, &arow);
  TCCSER(tid2, "WAVE", &icol2[0]);
  TCCSER(tid2, "PEAK", &icol2[1]);

  if (icol2[0] == -1)
    {
      sprintf(line, "Error: missing input column in table %s", in_table2);
      SCTPUT(line);
      SCETER(9, "Exiting...");
    }

  w2 = dvector(1, nrow);           /* laboratory wavelength */
  f2 = dvector(1, nrow);           /* laboratory flux       */
  
  for (j = 1; j <= nrow; j++)
    {
      TCERDD(tid2, j, icol2[0], &w2[j], &knul);
    }
  TCTCLO(tid2);

  /* READ THE COMPLETE INPUT TABLE of observed line centers */

  TCTOPN (in_table1, F_IO_MODE, &tid1);
  TCIGET(tid1, &mcol, &mrow, &msort, &acol, &arow);
  TCCSER(tid1, "XABS", &icol1[0]);
  TCCSER(tid1, "MABS", &icol1[1]);
  TCCSER(tid1, "IDENT", &icol1[2]);
  TCCSER(tid1, "WAVEP", &icol1[3]);
  TCCSER(tid1, "WAVE", &icol1[4]);
  TCCSER(tid1, "WAVEC", &icol1[5]);
  TCCSER(tid1, "RESIDUAL", &icol1[6]);
  TCCSER(tid1, "RESIDP", &icol1[7]);
  TCCSER(tid1, "WAVEG", &icol1[8]);
  TCCSER(tid1, "PEAK", &icol1[9]);
  TCCSER(tid1, "YABS", &icol1[10]);
  TCCSER(tid1, "XROT", &icol1[11]);

  if (icol1[0] == -1 || icol1[1] == -1 || icol1[9] == -1)
  {
    sprintf(line, "Error: missing input column in table %s", in_table1);
    SCTPUT(line);
    SCETER(9, "Exiting...");
  }

  if (icol1[3] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "Angstroem", ":WAVEP", &icol1[3]);
  if (icol1[7] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "Angstroem", ":RESIDP", &icol1[7]);
  if (icol1[8] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "Angstroem", ":WAVEG", &icol1[8]);
  if (icol1[4] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "Angstroem", ":WAVE", &icol1[4]);
  if (icol1[5] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "Angstroem", ":WAVEC", &icol1[5]);
  if (icol1[6] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "Angstroem", ":RESIDUAL", &icol1[6]);
  if (icol1[10] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "mm", ":YABS", &icol1[10]);
  if (icol1[11] == -1)
    TCCINI (tid1,D_R8_FORMAT, 1, "F8.2", "mm", ":XROT", &icol1[11]);

  /* initialize coefficients */

  for(j = 1; j <= nparams; j++)
    dcoef[j] = 0.;

  switch(toupper(do_wlc[0]))
    {
    case ('G'): /* grating relation */

      dcoef[1] = fabs( 2. / g * sin(deg_to_rad(theta_blaze)) * 1.e7);
      break;
      
    case ('C'): /* read wlc-coeffs from IN_D */
      
      if (TCTOPN(in_table3, F_I_MODE, &tid3))
	SCTPUT ("Error opening wavelength coefficient table");
      for(j = 1; j <= nparams; j++)
	{
	  TCERDD(tid3, 1, j, &dcoef[j], &knul);
	}
      TCTCLO(tid3);
      break;
      
    default:
      dcoef[1] = fabs( 2. / g * sin(deg_to_rad(theta_blaze)) * 1.e7);

    }

  x1 = dvector(1, mrow);       /* X-AXES */
  yy1 = dvector(1, mrow);      /* Y-AXES = ECHELLE-ORDER */
  yabs = dvector(1, mrow);     /* Y-AXES on chip */
  w1 = dvector(1, mrow);       /* WAVE-LENGTH from column WAVE  (lab) */
  wc1 = dvector(1, mrow);      /* WAVE-LENGTH from column WAVEC (fit) */
  r1 = dvector(1, mrow);       /* residual */
  rej = ivector(1, mrow);
  xtmp = dvector(1, mrow);     /* selected X-AXES */
  ytmp = dvector(1, mrow);     /* selected echelle order */
  itmp = dvector(1, mrow);     /* selected WAVE-LENGTH from column IDENT */

  xsin = sin(deg_to_rad(angle)); 
  xcos = cos(deg_to_rad(angle));

  for (i = 1; i <= mrow; i++)
    {
      TCERDR(tid1, i, icol1[0], &xbuff, &knul);
      TCERDR(tid1, i, icol1[1], &ybuff, &knul);
      TCERDR(tid1, i, icol1[9], &peak, &knul);
      TCERDR(tid1, i, icol1[10], &ybuff1, &knul);
      TCERDD(tid1, i, icol1[4], &wbuff1, &knul);

      if (knul == 0)
	{
	  w1[i] = wc1[i] = (double) wbuff1;
	}
      else
	{
	  w1[i] = wc1[i] = dnull;
	}

      x1[i] = (double) xbuff;
      yy1[i] = (double) ybuff;
      yabs[i] = (double) ybuff1;

      r1[i] = dnull;
      
      if (peak > thres) /* only strong lines */
	{

	  switch(toupper(do_wlc[0])) 
	    {
	    case ('G'): 	      /* Grating relation */
	      
	      lambda_central = 2.0 * sin(deg_to_rad(theta_blaze)) *
		1.E+7 / (yy1[i] * g);
	      beta0 = deg_to_rad(theta_blaze); 
	      
	      /*   rotate coordinate system !! */ 
	      
	      xrot = (x1[i] + xoffset) * xcos - yabs[i]  * xsin; 
	      TCEWRD(tid1, i, icol1[11], &xrot);
	      beta = beta0 + xrot / f;
	      lambda =  (sin(deg_to_rad(theta_blaze)) + sin(beta)) * 
		1.E+7 / (yy1[i] * g); 
	      w1[i] = wc1[i] = closest_to(w2, lambda, 1.0 * tol, dnull); 
	      break;
	      
	    case ('C'): /* guess Coefficients from table */
	      
	      lambda = eval_OS (x1[i],yy1[i],dcoef,nparams);
	      w1[i] = wc1[i] = closest_to(w2, lambda, 3.0 * tol, dnull);  
	      break;

	    default:
	      ;
	      /* use identified lines */
	    }
	}
      else
	{
	  w1[i] = wc1[i] = dnull;
	}
      TCEWRD(tid1, i, icol1[8], &w1[i]); 
    }

  iall = mrow;

  norder = (int) (fabs(yy1[1] - yy1[mrow]) + 1);

  /* identify lines; first try */

  reset = 1;
  isel = 0;

  /* fit the wavelengths of the observed lines */

  fit_select(itmp, xtmp, ytmp, &isel, w1, wc1, x1, yy1, 
	     iall, nparams, dcoef, dnull, 0);

  /* search for lines in the selected reference catalog */

  isel = 0;
  stdres = 0.0;

  for (i = 1; i <= mrow; i++)
    {
      r1[i] = w1[i] - wc1[i];
      TCEWRD(tid1, i, icol1[3], &w1[i]);
      TCEWRD(tid1, i, icol1[7], &r1[i]);
      if (w1[i] != dnull)
	{
	  stdres += fabs (r1[i]);
	  isel++;
	}
    }
  stdres /= isel;
  
  reset = 0;
  isel = 0;
  nmatch = 0;

  do
    {

      /* fit the wavelengths of the observed lines */

      fit_select(itmp, xtmp, ytmp, &isel, w1, wc1, x1, yy1, iall, 
		 nparams, dcoef, dnull, 0);

      /* search for lines in the complete reference catalog */

      prevnmatch = nmatch;
      nmatch = match(verif, w1, wc1, yy1, r1, iall, w2, nrow, alpha, &stdres,
                     dnull, rej);
      for (i = 1; i <= iall; i++)
	{
	  if (fabs(r1[i]) >  tol)
	    w1[i] = r1[i] = dnull;
	}
    } while (prevnmatch != nmatch);

  free_dvector(w2, 1, nrow);
  free_dvector(f2, 1, nrow);

  /* fit the selected identified lines a last time */

  fit_select(itmp, xtmp, ytmp, &isel, w1, wc1, x1, yy1, iall, 
	     nparams, dcoef, dnull, 1);

  /* finally write out the data to the table tid1 */

  isel = 0;
  stdres = 0.;
  for (i = 1; i <= iall; i++)
    { 
      r1[i] = w1[i] - wc1[i];
      if (fabs(r1[i]) >  tol)
	w1[i] = r1[i] = dnull; 
      TCEWRD(tid1, i, icol1[4], &w1[i]);
      TCEWRD(tid1, i, icol1[5], &wc1[i]);
      TCEWRD(tid1, i, icol1[6], &r1[i]);
      if (w1[i] != dnull)
	{
	  stdres += fabs (w1[i] - wc1[i]);
	  isel++;
	}
    }
  stdres /= isel;

  /* write the dispersion coefficients to table */

  if (TCTINI(in_table3, F_TRANS, F_O_MODE, nparams, mrow, &tid3))
    SCTPUT ("Error creating table");
  for(j = 1; j <= nparams; j++)
    {
      sprintf(lbl, "coeff%02i", j);
      TCCINI(tid3, D_R8_FORMAT, 1, "E14.7", "", lbl, &j);
    }
  for(j = 1; j <= nparams; j++)
    {
      TCEWRD(tid3, 1, j, &dcoef[j]);
    }
  TCTCLO(tid3);

  /* write the dispersion coefficients to descriptors of table tid1 */

  buffdeg = 5;
  buffcoef = dvector(1, buffdeg * norder);
  
  global_to_poly 
    (dcoef, buffcoef, (int) (yy1[mrow]+0.5), 
     (int) (yy1[1]+0.5), buffdeg);

  SCDWRD(tid1, "DCOEF", &buffcoef[1], 1, norder * buffdeg, &kun); 
  SCDWRI(tid1, "NPARAMS", &nparams, 1, 1, &kun); 
  SCDWRI(tid1, "NORDER", &norder, 1, 1, &kun);
  SCDWRI(tid1, "FITDEG", &buffdeg, 1, 1, &kun);
  /* AK */
  poltyp = 0;
  SCDWRI(tid1, "POLTYP", &poltyp, 1, 1, &kun);

  /* deallocate memory */

  free_dvector(x1, 1, mrow);
  free_dvector(yy1, 1, mrow);
  free_dvector(w1, 1, mrow);
  free_dvector(wc1, 1, mrow);
  free_dvector(r1, 1, mrow);
  free_dvector(xtmp, 1, mrow);
  free_dvector(ytmp, 1, mrow);
  free_dvector(itmp, 1, mrow);
  free_ivector(rej, 1, mrow);
  free_dvector(dcoef, 1, nparams);
  free_dvector(buffcoef, 1, buffdeg * norder);

  TCTCLO(tid1);  
  SCSEPI();
  return 0;
} 

/* end of main */


/*---------------------------------------------------------------------------*/
void fit_select
#ifdef __STDC__
(
 double itmp[], double xtmp[], double ytmp[], int *isel,
 double w1[], double wc1[], double x1[], double yy1[], int iall,
 int nparams, double dcoef[], double dnull,
  int printout
 )
#else
     (
      itmp, xtmp, ytmp, isel, w1, wc1, x1, yy1,iall, nparams, dcoef,
      dnull, printout
      )
     int *isel, iall, nparams,  printout;
     double itmp[], xtmp[], ytmp[], w1[], wc1[], x1[], yy1[];
     double dcoef[], dnull;
#endif

{
  int i;
  int sel = 0;

  for (i = 1; i <= iall; i++) 
    {
      if (w1[i] != dnull)
	{
	  ++sel;
	  xtmp[sel] = x1[i];
	  ytmp[sel] = yy1[i];
	  itmp[sel] = w1[i];
	}
    }

  /* FIT the dispersion coefficients and compute 
     the wavelengths for all lines */

  *isel = sel;
  if (sel >= nparams + 1)
    {
      fit_wavelength(itmp, xtmp, ytmp, sel, wc1, x1, yy1, iall, 
			nparams, dcoef, printout);
    }
}

double closest_to
#ifdef __STDC__
(
 double *reference, double what, double howclose, double dnull
 )
#else
     (
      reference, what, howclose, dnull
      )
     double *reference, what, howclose, dnull;
#endif
{
  int i, j;
  double  x1, x2;

  x1 = what - howclose;
  x2 = what + howclose;

  i = 1;
  
  while (reference[i] < x1)
    i++;
  
  j = i-1;
  while (reference[j] < x2)
    j++;

  if ((j-i) == 1) /* identifications should be unique */
    {
      return reference[i];
    }
  else
    {
      return dnull;
    }	
}