/* @(#)lncalib.c	17.1.1.1 (ES0-DMD) 01/25/02 17:53:59 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  This program is free software; you can redistribute it and/or 
  modify it under the terms of the GNU General Public License as 
  published by the Free Software Foundation; either version 2 of 
  the License, or (at your option) any later version.
 
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
 
  You should have received a copy of the GNU General Public 
  License along with this program; if not, write to the Free 
  Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
/* .COPYRIGHT   (C) 1993 European Southern Observatory     */
/* .IDENT       lncalib.c                                  */
/* .AUTHORS     Pascal Ballester (ESO/Garching)            */
/*              Cristian Levin   (ESO/La Silla)            */
/* .KEYWORDS    Spectroscopy, Long-Slit                    */
/* .PURPOSE     Wavelength Calibration for 1D, Long and MOS*/
/* .VERSION     1.0  Package Creation  17-MAR-1993         */
/*                                                         */
/*
   .INPUT/OUTPUT:
       IN_A         :  line table   (def. line.tbl)
       IN_B         :  line catalog
       OUT_A        :  output coefficients table (def. coerbr.tbl)
       OUT_B        :  optional input guess table name
       INPUTC       :  mode (IDENT or GUESS)
       INPUTD/D/1/6 :  start(1),start(2),step(1),step(2),startw,pixel
       INPUTR/R/1/4 :  alpha,maxdev,tol,shift
       INPUTI/I/1/5 :  degx,degl,miniter,maxiter,ystart
       INPUTI(10)   :  Debug level
   -------------------------------------------------------
*/

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

#include <lncalib.h>
#include <lndisp.h>
#include <stdio.h>
#include <math.h>

int       inull;   /* Representation of NULL values */
float     rnull;
double    dnull;


main()

{

   char      line[60], lincat[60], outtab[60], mode[10]; 
             /* Table names */
   char      guess[60], text[120];
   int       tidlin, tidcat, kunit;    /* Table id numbers */
   int       nbcol, nbrow, nsort, allcol, allrow; /* Descriptors of line */
   int       nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat; 
             /* Descriptors of lincat */
   int       catwav, catint,  linid,  linwav,  linwavc;
   int       linx,  liny,  lindif;         /* Column numbers   */
   double    *colcat, *colid, *colwav, *colwavc, *colx;
   double    *coly, *coldif, *yval;          /* Column arrays    */
   int       linrej, *reject, *select, *posline, *yrow, nyrow, nbsel;
   double    *xid, *lid, *xline, *xshift, *yline;   /* Additional arrays */
   double    *linear, *coldel, *coldelc;
   int       lindel, lindelc, linpix, linrms, deg1=1L;
   int       ipar[10], degx, degl, minter, maxter, ystart, wrang[2];  
             /* Integer parameters */
   float     rpar[10], alpha, maxdev, tol, shift, imin;   /* Real parameters */
   int       status, iav, knull, actvals, null; 
            /* Standard Interfaces variables */

   int       ystart_row, rowstart, rowend, direction;
   double    ystart_world, delta, mindelta, startw, stepw, ord2;

   int       nid, degree, disp;
   int       iter, prevmatch, end_of_iter, ungraceful_exit;
   int       match(), nmatch;
   double    stdres, stdpix, pixel, tolwav;
   int       npos, ncat, verif=0;
   int       y_index, row, nblines=0;
   int       current_y_pix;
   double    current_y;
   double    dpar[10], start[2], step[2];
   int       linb;

   int      i,maxpos;
   double   resabs, maxres;

   SCSPRO("lncalib");

/* Read name of table line in keyword IN_A, lincat in IN_B 
   and mode in INPUTC(1:9). Mode can be IDENT or GUESS */

   status = SCKGETC("IN_A",  1, 60, &actvals, line);
   status = SCKGETC("IN_B",  1, 60, &actvals, lincat);
   status = SCKGETC("OUT_A", 1, 60, &actvals, outtab);
   status = SCKGETC("OUT_B", 1, 60, &actvals, guess);
   status = SCKGETC("INPUTC",1, 9,  &actvals, mode);

   status = SCKRDR("INPUTR", 1,  5, &actvals, rpar, &kunit, &null);
   status = SCKRDI("INPUTI", 1,  7, &actvals, ipar, &kunit, &null);
   status = SCKRDI("INPUTI", 10, 1, &actvals, &disp,&kunit, &null);  /* Display level */
   status = SCKRDD("INPUTD", 1,  7, &actvals, dpar, &kunit, &null);

   degx = ipar[0],  degl = ipar[1],   minter = ipar[2], maxter = ipar[3], ystart=ipar[4],   wrang[0] = ipar[5], wrang[1] = ipar[6];
   alpha = rpar[0], maxdev = rpar[1], tol = rpar[2],    shift = rpar[3],
   imin = rpar[4];
   start[0] = dpar[0], start[1] = dpar[1], step[0] = dpar[2], step[1] = dpar[3];
   startw = dpar[4], stepw = dpar[5], ord2= dpar[6];

/* Open table line and search for columns :WAVE, :WAVEC, :RESIDUA, 
    :REJECT. Column :IDENT is searched for only in mode IDENT. All
    table are mapped to double precision arrays */

   if (TCTOPN(line, F_IO_MODE, &tidlin))
      SCTPUT("**** Error while opening table line.tbl");

   TCIGET(tidlin, &nbcol, &nbrow, &nsort, &allcol, &allrow);

   if (TCTOPN(lincat, F_I_MODE, &tidcat))
      SCTPUT("**** Error while opening line catalog");

   TCIGET(tidcat, &nbcolcat, &nbrowcat, &nsortcat, &allcolcat, &allrowcat);

   /* Finds columns :WAVE and :INTENS in line catalog                    */

   TCCSER(tidcat, ":WAVE", &catwav);
      if (catwav == (-1)) SCTPUT("**** Column :WAVE not found");

   TCCSER(tidcat, ":INTENSITY", &catint);
      if (catint == (-1)) TCCSER(tidcat, ":STRENGTH", &catint);

   /* Locates columns :IDENT, :WAVEC, :X, :Y, :RESIDUAL in table line */

   TCCSER(tidlin, ":X", &linx);
      if (linx == (-1)) SCTPUT("**** Column :X not found");

   TCCSER(tidlin, ":Y", &liny);
      if (liny == (-1)) SCTPUT("**** Column :Y not found");

   if (toupper(mode[0]) == 'I')  {
        TCCSER(tidlin, ":IDENT", &linid);
        if (linid == (-1)) SCTPUT("**** Column :IDENT not found");
      }

   TCCSER(tidlin, ":WAVE", &linwav);
      if (linwav == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
	                       "WAVE", &linwav);

   TCCSER(tidlin, ":WAVEC", &linwavc);
      if (linwavc == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
	                       "WAVEC", &linwavc);

   TCCSER(tidlin, ":DELTA", &lindel);
      if (lindel == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
	                       "DELTA", &lindel);

   TCCSER(tidlin, ":DELTAC", &lindelc);
      if (lindelc == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
	                       "DELTAC", &lindelc);

   TCCSER(tidlin, ":RESIDUAL", &lindif);
      if (lindif == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom",
	                       "RESIDUAL", &lindif);

   TCCSER(tidlin, ":REJECT", &linrej);
      if (linrej == (-1)) TCCINI(tidlin, D_I4_FORMAT, 1, "I6", "Rejection Code",
	                       "REJECT", &linrej);

/* Read other parameters: 
             INPUTI(1)   Degree in X
             INPUTI(2)   Degree in Y
             INPUTI(3)   Minimal number of iterations
             INPUTI(4)   Maximal number of iterations

             INPUTR(1)   Parameter Alpha for matching
             INPUTR(2)   Maximal standard deviation of residuals (pixels)
             INPUTR(3)   Tolerance (pixels if >0, wavel. units if <0)
             INPUTR(4)   Shift in world coordinates (mode GUESS only)
*/

   if (disp >= 100) {
     sprintf(text,"Line table   : %s",line), SCTPUT(text);
     sprintf(text,"Line catalog : %s",lincat), SCTPUT(text);
     sprintf(text,"Mode         : %s\n",mode), SCTPUT(text);

     sprintf(text,"Nb of iterations (min,max) : %d , %d",minter,maxter), SCTPUT(text);
     sprintf(text,"Degree                     : %d",degx), SCTPUT(text);
     sprintf(text,"Tolerance (pixels)         : %f",tol), SCTPUT(text);
     sprintf(text,"Rejection parameter Alpha  : %f",alpha), SCTPUT(text);
     sprintf(text,"Maximum deviation (pixels) : %f",maxdev), SCTPUT(text);
     sprintf(text,"Shift (pixels)             : %f\n",shift), SCTPUT(text);
   }

/* Initialisations */

   /* Get representation of NULL values and allocate memory for auxiliary arrays */

   TCMNUL(&inull, &rnull, &dnull);

   /* Load line table and line catalog in memory. Table columns are stored in  double 
      although they can be of any type (I, R*4 or R*8) */
   
   colcat = (double *) osmmget((nbrowcat+1)*sizeof(double));
   ncat = read_col_cat(tidcat, nbrowcat, catwav, catint, wrang, imin, colcat);

   select = (int      *) osmmget((nbrow+1)*sizeof(int     ));
   nbsel = read_select(tidlin, nbrow, select);
 
   sprintf(text,"Number of lines in table line (total, selected) : %d , %d", nbrow, nbsel), SCTPUT(text);

   colx = (double *) osmmget((nbrow+1)*sizeof(double));
   read_col(tidlin, nbrow, linx,   colx);

   coly = (double *) osmmget((nbrow+1)*sizeof(double));
   read_col(tidlin, nbrow, liny,   coly);

   if (toupper(mode[0]) == 'I')  {
      colid = (double *) osmmget((nbrow+1)*sizeof(double));
      read_col(tidlin, nbrow, linid,  colid);
   }
 
   coldif = (double *) osmmget((nbrow+1)*sizeof(double));
   read_col(tidlin, nbrow, lindif, coldif);

   xid     = (double *)   osmmget((nbrow+1)*sizeof(double));
   lid     = (double *)   osmmget((nbrow+1)*sizeof(double));
   xline   = (double *)   osmmget((nbrow+1)*sizeof(double));
   xshift  = (double *)   osmmget((nbrow+1)*sizeof(double));
   yline   = (double *)   osmmget((nbrow+1)*sizeof(double));
   posline = (int      *) osmmget((nbrow+1)*sizeof(double));
   colwav  = (double *)   osmmget((nbrow+1)*sizeof(double));
   colwavc = (double *)   osmmget((nbrow+1)*sizeof(double));
   coldel  = (double *)   osmmget((nbrow+1)*sizeof(double));
   coldelc = (double *)   osmmget((nbrow+1)*sizeof(double));
   linear  = (double *)   osmmget((nbrow+1)*sizeof(double));
   reject  = (int      *) osmmget((nbrow+1)*sizeof(int     ));
   yval    = (double *)   osmmget((nbrow+1)*sizeof(double));
   yrow    = (int      *) osmmget((nbrow+1)*sizeof(int     ));

/* Write null values in all output columns */

    for (row=0; row<=nbrow; row++) {
         posline[row] = row;
         colwav[row]  = dnull;
         colwavc[row] = dnull;
         coldel[row]  = dnull;
         coldelc[row] = dnull;
         coldif[row]  = dnull;
         reject[row]  = inull;
       }

    write_dcol(tidlin, nbrow, posline, linwav, colwav);
    write_dcol(tidlin, nbrow, posline, linwavc,colwavc);
    write_dcol(tidlin, nbrow, posline, lindel, coldel);
    write_dcol(tidlin, nbrow, posline, lindelc,coldelc);
    write_dcol(tidlin, nbrow, posline, lindif, coldif);
    write_icol(tidlin, nbrow, posline, linrej, reject);

     /*  Check line catalog  */

     verif = -1;
     nmatch =  match(verif, colwav, colwavc, yline, coldif, nblines, colcat, 
                     ncat, alpha, &stdres, dnull, reject);
     verif = 0;

     if (nmatch < 0) {

         sprintf(text,"Error: Line catalog is not properly organized\n Apply command: SORT/TABLE {lincat} :WAVE\n And check that no wavelength is duplicated");
         SCETER(9,text);
       }

/* Go through the line table to identify the starting row number 
   for the different y values */

   y_index = 1, row=0;
   while (y_index <= nbsel) {
      yval[++row] = coly[y_index];
      yrow[row]   = y_index;
      while (y_index <= nbsel && coly[y_index] == yval[row]) ++y_index;
   }
   nyrow = row;
   yrow[nyrow+1] = nbsel + 1;

/* Finds the position ystart_row corresponding to the position of 
   ystart */

   ystart_world = (ystart-1)*step[1] + start[1];
   for (row=1; row<=nyrow; row++) {
        delta = yval[row] - ystart_world;
        if (delta < 0) 
            delta = delta*(-1.);
        if (row == 1)  mindelta = delta;
        if (delta <= mindelta)
            mindelta = delta, ystart_row = row;
      }

/* Loop on y values */

   for (direction=1; direction>=-1; direction-=2) {

/* Estimate a first dispersion relation in mode ident or read it
   in mode guess */

   switch (toupper(mode[0])) {

   case ('I') : {
        read_ident(colx, colid, nbsel, xid, lid, &nid);
        if (nid < 2) {
            sprintf (text,"Not enough identifications... Exiting."), SCTPUT(text);
            goto ciao;
	  }
        pixel = fit_disp (&nid, &degx, xid , lid);
        pixel *= step[0];
        if (pixel < 0.) pixel *= (-1.);
        break;
      }

    case ('G') :  {
        initdisp(guess,"OLD",1);
        pixel = readdisp(ystart);
        finishdisp();
        pixel *= step[0];
        if (pixel < 0.) pixel *= (-1.);
        break;
      }

   case ('L') : {
        double coefs[3];
        coefs[0] = startw;
        coefs[1] = stepw;
        coefs[2] = ord2;
        setdisp(2,coefs);
        pixel = stepw*step[0];
        if (pixel < 0.) pixel *= (-1.);
        break;
      }

   default: {
        sprintf(text,"Error in module lncalib.c:\nUnknown calibration method %s",mode), SCTPUT(text);
        SCETER(9, "Exiting...");
        break;
      }
   }


   if (direction == 1) {
        setrefdeg(degx);
        initdisp(outtab,"NEW",1);
        rowstart = ystart_row;
        rowend   = nyrow + 1;
      }

   if (direction == -1) {
        initdisp(outtab,"OLD",1);
        rowstart = ystart_row - 1;
        rowend   = 0;
      }

   for (row=rowstart; row != rowend; row+=direction) {

/* Read a set of values at constant y */

   current_y = yval[row];
   current_y_pix = 1 + (current_y - start[1])/step[1];

   nblines = 0;
   for (y_index=yrow[row]; y_index<yrow[row+1]; y_index++) {
         xline[++nblines]  =  colx[y_index];
         xshift[nblines]   = xline[nblines];
         if (toupper(mode[0]) == 'G' && row == ystart_row)  
                xshift[nblines] -= shift;
         yline[nblines]    =  coly[y_index];
          posline[nblines]  =  select[y_index];
       }

/* Iteration Loop */

   iter=0, prevmatch=0;
   /* if (row>1) disp=60; */

   end_of_iter     = FALSE;
   ungraceful_exit = FALSE;

   if (disp >= 100)  sprintf(text,"Calibrating for Y = %d",current_y_pix), SCTPUT(text);

   if (disp >= 100) {
      sprintf(text,"----------------------------------------------------------------------"), SCTPUT(text);
      sprintf(text,"--     Number of   ---      Standard Error      ---  Average Pixel ---"), SCTPUT(text);
      sprintf(text,"-- Ident. -- Lines --- Wav. Units --   Pixels   ---     Size       ---"), SCTPUT(text);
      sprintf(text,"----------------------------------------------------------------------"), SCTPUT(text);
    }

   while (!end_of_iter && !ungraceful_exit) {

     /* Compute wavelength for all lines */
     /* printdisp(); */
     if (iter == 0)  eval_disp(xshift, colwavc, nblines);
     else            eval_disp(xline, colwavc, nblines);

     /* match lines against the catalog */
     nmatch =  match(verif, colwav, colwavc, yline, coldif, nblines, colcat, 
                     ncat, alpha, &stdres, dnull, reject);

     stdpix = stdres/pixel;

     /* Test end of iteration */
     iter++;
     end_of_iter = ((iter >= maxter) || (iter > minter &&  nmatch == prevmatch));

     prevmatch = nmatch;

     /* Display intermediate results */
     if (disp >= 100) {
       sprintf(text,"--   %d   --  %d   ---  %f  --  %f  ---    %f    ---",
                    nmatch, nblines, stdres, stdpix,pixel), SCTPUT(text);
     }

     /* Estimate the new dispersion relation on the matched lines */

    if (!(ungraceful_exit = (stdpix > maxdev))) {
          read_ident(xline, colwav, nblines, xid, lid, &nid);
          pixel = fit_disp (&nid, &degx, xid, lid);
          pixel *= step[0];
         if (pixel < 0.) pixel *= (-1.);
	}
   }

/* Final Fit and Ciao... */

   if (ungraceful_exit) 
     sprintf(text,"Sorry, wrong identifications..."), SCTPUT(text);

   else {

    tolwav = (double) ((tol>0.) ?  tol*pixel : (0.-tol));
    maxres = tolwav;
    
    while (maxres >= tolwav) {

    maxres = 0.;

    for (i=1; i<=nblines; i++) {
         if (colwav[i] != dnull) {
             resabs = (coldif[i] < 0.) ? -coldif[i] : coldif[i];
             if (resabs > maxres) maxpos = i, maxres = resabs;
	   }
       }

    if (disp >= 150) 
       printf ("Max residual, tolerance, position : %f %f %d\n",
       maxres, tolwav,maxpos);

    if (maxres > tolwav) {
        colwav[maxpos] = dnull, reject[maxpos] = RESIDUAL_GT_TOL;
        read_ident (xline, colwav, nblines, xid, lid, &nid);
        fit_disp   (&nid, &degx, xid, lid);
        eval_disp  (xline, colwavc, nblines);
        stdres = 0. , actvals=0;
        for (i=1; i<=nblines; i++) {
             if (colwav[i] != dnull) {
                 coldif[i] = colwav[i] - colwavc[i];
                 actvals++, stdres += coldif[i]*coldif[i];
	       }
	   }
        stdres = sqrt(stdres/(double) actvals);
      }
  } /* End of while */

    if (disp >= 50)   
    sprintf(text,"Y = %3d -- %2d lines out of %2d -- RMS = %f wav. units",
    current_y_pix, nid, nblines, stdres), SCTPUT(text);
/** current_y --> yline
    sprintf(text,"Y = %3d -- %2d lines out of %2d -- RMS = %f wav. units",
    (int)current_y, nid, nblines, stdres), SCTPUT(text);
**/

    read_ident (xline, colwav, nblines, xid, lid, &nid);

    /* Now fitting a linear relation to get delta and deltac */
    pixel = fit_disp(&nid,&deg1,xid,lid);
    if (pixel < 0.) pixel *= (-1.);
    eval_disp(xline,linear,nblines);
    for(i=1; i<= nblines; i++) {
        if (colwav[i] != dnull) coldel[i]  = linear[i] - colwav[i];
        else coldel[i] = dnull;
        coldelc[i] = linear[i] - colwavc[i];
      }

    if (nid >= 2)    fit_disp   (&nid, &degx, xid, lid);

    linb = (int) ((current_y - start[1])/step[1] + 1.5);
    writedisp  (row, linb, current_y, pixel, stdres);

    write_dcol(tidlin, nblines, posline, linwav, colwav);
    write_dcol(tidlin, nblines, posline, linwavc,colwavc);
    write_dcol(tidlin, nblines, posline, lindel, coldel);
    write_dcol(tidlin, nblines, posline, lindelc,coldelc);
    write_dcol(tidlin, nblines, posline, lindif, coldif);

    write_icol(tidlin, nblines, posline, linrej, reject);

  }  /* End of else ... */

 } /* End of loop on y_index */

  finishdisp();

 } /* End of loop on direction */

 ciao:

  TCTCLO(tidcat);
  TCTCLO(tidlin);

  osmmfree((char *)colcat);
  if (toupper(mode[0]) == 'I') osmmfree((char *)colid);
  osmmfree((char *)colwav);
  osmmfree((char *)colx);
  osmmfree((char *)coly);
  osmmfree((char *)coldif);
  osmmfree((char *)yrow);
  osmmfree((char *)yval);

  osmmfree((char *)xid);
  osmmfree((char *)lid);
  osmmfree((char *)reject);

  SCSEPI();
}


write_dcol(tid, nb, select, colnb, col)

int      tid, nb, colnb, select[];
double   col[];

{

int      i;

for (i=1; i<=nb; i++)
    TCEWRD(tid, select[i], colnb, &col[i]);

}

write_icol(tid, nb, select, colnb, col)

int      tid, nb, colnb, select[], col[];

{

int      i;

for (i=1; i<=nb; i++)
    TCEWRI(tid, select[i], colnb, &col[i]);

}

fit_select (colx, colid, coldif, nb, tolwav, reject, xid, lid, nid)

double   tolwav, colid[], colx[], coldif[], xid[], lid[];
int      nb, *nid, reject[];

{

int       i;
double    resabs;

*nid = 0;

for (i=1; i<=nb; i++) {
     if (colid[i] != dnull) {
           resabs = (coldif[i] < 0) ? -coldif[i] : coldif[i];
           if (resabs < tolwav) {
               xid[++(*nid)] = colx[i];
               lid[*nid] = colid[i];
               }
           else reject[i] = RESIDUAL_GT_TOL;
	 }
   }

}

read_col(tid, nb, colnb, col)

  int      tid, colnb, nb;
  double   col[];

{
int      i, k=0, null, sel;

for (i=1; i<=nb; i++) {
    TCSGET(tid, i, &sel);
    if (sel) {
       TCERDD(tid, i, colnb, &col[++k], &null);
       if (null) col[k] = dnull;
     }
  }
return(k);
}

read_col_cat(tid, nb, colwav, colint, wrang, imin, col)

  int      wrang[];
  int      tid, colwav, colint, nb;
  double   col[];
  float imin;

{
int      i, k=0, null, null_intens=0, sel;
double wave;
float intens=0.;
char  text[60];

for (i=1; i<=nb; i++) {
    TCSGET(tid, i, &sel);
    if (sel) {
       TCERDD(tid, i, colwav, &wave, &null);
       if ( colint != -1 )
           TCERDR(tid, i, colint, &intens, &null_intens);
       if ( wave >= wrang[0] && wave <= wrang[1] && 
            (null_intens || intens >= imin || colint == -1) )
	   col[++k] = wave;	
     }
  }
sprintf(text,"Number of lines in catalog (total, selected)    : %d , %d",nb,k);
SCTPUT(text);

return(k);
}

read_select(tid, nb, col)

int      tid, nb, col[];

{
int      i, k=0, null, sel;

for (i=1; i<=nb; i++) {
    TCSGET(tid, i, &sel);
    if (sel)  col[++k] = i;
  }
return(k);
}

read_ident(colx, colid, nbsel, xid, lid, nid)

double   *colx, *colid, *xid, *lid;
int      nbsel, *nid;

{

int      i;

*nid = 0;

for (i=1; i<=nbsel; i++) {
    if (colid[i] != dnull) {
        xid[++(*nid)]     = colx[i];
        lid[*nid] = colid[i];
      }
}
}