/* @(#)necmatch.c	17.1.1.1 (ES0-DMD) 01/25/02 17:51:30 */
/*===========================================================================
  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) 1991 European Southern Observatory           */
/* .IDENT       echmatch.c                                       */
/* .AUTHOR      Pascal Ballester,  ESO - Garching                */
/* .KEYWORDS    Spectroscopy, Echelle,                           */
/* .PURPOSE     associates line catalog and computed wavelengths */
/* .VERSION     1.0    Creation    15-JAN-1991  PB               */
/* ------------------------------------------------------------- */

#include <tbldef.h>
#include <midas_def.h>
#define  ipos(col,row,nrow)   (col-1)*nrow + row - 1

main()
{

   int             linwav, catwav, linid, lindif, liny, null, kunit;
   char            line[60], lincat[60];
   double          ident, lambda_lin, lambda_cat;
   int             row, last, match_flag, rowcat, written = 0;
   int             order, ordref, error, bufsize, imno;
   float           minmax[2], thres[3];
   int             status, npar, iav, knull, actvals, tidlin, tidcat;
   int             nbcol, nbrow, nsort, allcol, allrow;
   int             nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat;
   int             rowlow, rowupp, rowref, index, nbsel;
   double          resref, reslow, resupp, residual;
   double          resref_abs, reslow_abs, resupp_abs, res_abs;
   double          lambda_ref, lambda_low, lambda_upp, lambda_inf, lambda_sup;
   double          *pntr_buf;

   SCSPRO("echmatch");

   status = SCKRDR("INPUTR", 1, 1, &actvals, thres, &kunit, &null);

   status = SCKGETC("IN_A", 1, 60, &actvals, line);
   if (status)
      SCTPUT("Error while reading IN_A");

   status = SCKGETC("IN_B", 1, 60, &actvals, lincat);
   if (status)
      SCTPUT("Error while reading IN_B");

   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 column :WAVE in line catalog                    */

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


   /* Verification loop */
   /* If the threshold parameter thres<0, the line catalog is explored to
      make sure that wavelengths are sorted by increasing value and no 
      wavelength is duplicated */

   if (thres[0] < 0) {
   thres[0] = 1.;
   TCERDD(tidcat, 1, catwav, &lambda_low, &null);
   for (row=2; row<=nbrowcat; row++) {
       TCERDD(tidcat, row, catwav, &lambda_upp, &null);
       if (lambda_low > lambda_upp) {
          printf("Warning: Column :WAVE of the table %s is not sorted by increasing wavelength\n",lincat);
          thres[0] = -1.;
	}
       if (lambda_low == lambda_upp) {
          printf("Error: Column :WAVE of the table %s contains duplicated wavelength : %f\n",lincat, lambda_low);
          thres[0] = -2.;
	}
       lambda_low = lambda_upp;
     }
   status = SCKWRR("INPUTR", &thres[0], 1, 1, &kunit);
   SCSEPI();
   exit();
 }

   /* Locates columns :IDENT, :WAVEC, :RESIDUAL in table LINE */

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

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

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

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


   bufsize = 6*nbrow;

   SCFCRE("buffer.bdf",D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,bufsize,&imno);
   SCFMAP(imno,F_X_MODE,1, bufsize, &actvals, (char **)&pntr_buf);


   /* Matching loop */

   SCTPUT ("Identifying lines...\n");

   rowcat = 2;

   for (row = 1; row <= nbrow; row++) {

      /* Reads current wavelength to be identified and relative order number */
      TCERDD(tidlin, row, linwav, &lambda_lin, &null);
      TCERDI(tidlin, row, liny,   &ordref, &null);
      /* printf ("Reading line.tbl row = %d line = %f\n",row,lambda_lin); */

      pntr_buf[ipos(6,row,nbrow)] = 0.;

      if (lambda_lin != 0.) {

          pntr_buf[ipos(6,row,nbrow)] = 1.;
          reslow_abs = 99999.9;
          resupp_abs = 99999.9;

          /* Finds the residual to the nearest lower line in the same order */
          if (row > 1) {
            rowlow = row - 1;
            TCERDI(tidlin, rowlow, liny,  &order, &null);
            if (order == ordref) {
	      TCERDD(tidlin, rowlow, linwav, &lambda_low, &null);
              reslow = lambda_lin - lambda_low;
              reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
	    }
	  }

          /* Finds the residual to the nearest upper line in the same order */
          if (row < nbrow) {
            rowupp = row + 1;
            TCERDI(tidlin, rowupp, liny,  &order, &null);
            if (order == ordref) {
	      TCERDD(tidlin, rowupp, linwav, &lambda_upp, &null);
              resupp = lambda_lin - lambda_upp;
              resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
	    }
	  }

          /* Compares both residuals and selects the minimum */
          if (row == 1) reslow_abs = resupp_abs;
          if (row == nbrow) resupp_abs = reslow_abs;
          /* printf("Table line: low,high %f %f\n",reslow_abs, resupp_abs); */

          if (reslow_abs < resupp_abs)
             pntr_buf[ipos(1,row,nbrow)] = reslow_abs;
          else
             pntr_buf[ipos(1,row,nbrow)] = resupp_abs;

          
          /* Reads the line catalog at the first entry point and initializes variables */
          /* printf ("Reading catalog thar.tbl rowcat = %d\n",rowcat); */
          if (rowcat > 1) {
          rowlow = rowcat - 1;
	  TCERDD(tidcat, rowlow, catwav, &lambda_low, &null);
          rowupp = rowcat;
	  TCERDD(tidcat, rowupp, catwav, &lambda_upp, &null);
          reslow = lambda_lin - lambda_low;
          reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
          resupp = lambda_lin - lambda_upp;
          resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
	}

          /* The line catalog is assumed to be ordered. Depending on which direction
          is found the minimum residual, the catalog will be searched further in this
          direction */

          if (reslow_abs < resupp_abs) {
              index = -1;
              resref = reslow;
              resref_abs = reslow_abs;
              lambda_ref = lambda_low;
              rowcat = rowlow;
	    }
          else {
              index = 1;
              resref = resupp;
              resref_abs = resupp_abs;
              lambda_ref = lambda_upp;
              rowcat = rowupp;
	    }

          res_abs = resref_abs;
          residual = resref;
          lambda_cat = lambda_ref;

          /* Explores the catalog in the direction indicated by index and finds the
          position of the minimum residual */
          while (res_abs <= resref_abs) {
             resref = residual;
             resref_abs = res_abs;
             lambda_ref = lambda_cat;
             rowcat += index;
             if (rowcat<=1 || rowcat >= (nbrowcat-1)) {
                 printf ("Error: reached limit of line catalog\n");
                 printf ("Wavelength: %f Position: %d\n",lambda_lin,rowcat);
                 break;
	       }
             else {
	     TCERDD(tidcat, rowcat, catwav, &lambda_cat, &null);
  /* printf ("Matching loop at rowcat = %d cat = %f\n",rowcat,lambda_cat); */
             residual = lambda_lin - lambda_cat;
             res_abs = (residual < 0) ?  -residual : residual; /* abs. value */
	   }
	   }

             rowref = rowcat - index;
             /* printf("Ref: Read wavelength %f at position %d\n",lambda_ref,rowref); */
             /* rowref is now the position of the minimum residual */
             pntr_buf[ipos(2,row,nbrow)] = resref_abs;
             pntr_buf[ipos(3,row,nbrow)] = lambda_ref;
             pntr_buf[ipos(4,row,nbrow)] = lambda_lin;

          if (rowref > 1) {
            rowlow = rowref - 1;
            TCERDD(tidcat, rowlow, catwav,  &lambda_low, &null);
             /* printf("Low: Read wavelength %f at position %d\n",lambda_low,rowlow); */
             reslow = lambda_ref - lambda_low;
             reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
	    }

          if (rowref < nbrowcat && rowref > 1) {
            rowupp = rowref + 1;
            TCERDD(tidcat, rowupp, catwav,  &lambda_upp, &null);
            /* printf("Upp: Read wavelength %f at position %d\n",lambda_upp,rowupp); */
            resupp = lambda_ref - lambda_upp;
            resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
	    }

          if (rowref == 1) reslow_abs = resupp_abs;
          if (rowref == nbrowcat) resupp_abs = reslow_abs;
          /* printf("Catalog: low,high %f %f\n",reslow_abs, resupp_abs); */

          if (reslow_abs < resupp_abs)
             pntr_buf[ipos(5,row,nbrow)] = reslow_abs;
          else
             pntr_buf[ipos(5,row,nbrow)] = resupp_abs;


          if (pntr_buf[ipos(1,row,nbrow)] < pntr_buf[ipos(5,row,nbrow)])
             reslow = pntr_buf[ipos(1,row,nbrow)];
          else
             reslow = pntr_buf[ipos(5,row,nbrow)];
          
          reslow = reslow * thres[0];

          /* printf("Comparing line %f with catalog %f.\nResidual is %f and maximum residual is %f\n",lambda_ref,lambda_lin,pntr_buf[ipos(2,row,nbrow)],reslow); */

          if (pntr_buf[ipos(2,row,nbrow)] < reslow) {
	          TCEWRD(tidlin, row, linid,  &lambda_ref);
	          TCEWRD(tidlin, row, lindif, &resref);
      /* printf("**** Identified line: %f  %f\n",lambda_ref,resref); */
		}

      }				/* End of if lambda_lin != 0.  */

/*    printf("Row %d Order %d - %f  %f  %f  %f  %f\n",row,ordref, */
/*    pntr_buf[ipos(1,row,nbrow)],pntr_buf[ipos(2,row,nbrow)],  */
/*    pntr_buf[ipos(3,row,nbrow)],pntr_buf[ipos(4,row,nbrow)],  */
/*    pntr_buf[ipos(5,row,nbrow)]);  */

   }				/* End of for */


   TCTCLO(tidlin);
   TCTCLO(tidcat);
   SCSEPI();
}