/* @(#)lnmatch.c	10.1.1.2 (ES0-DMD) 12/18/95 18:27:23 */
/*===========================================================================
  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
===========================================================================*/

/* @(#)lnmatch.c	10.1.1.2  (ESO)  12/18/95  18:27:23 */
/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++       */
/* .COPYRIGHT   (C) 1993 European Southern Observatory           */
/* .IDENT       echmatch.c                                       */
/* .KEYWORDS    Spectroscopy, Echelle,                           */
/* .AUTHORS     Pascal Ballester (ESO/Garching)                  */
/* .PURPOSE     associates line catalog and computed wavelengths */
/* .VERSION     1.0  Package Creation  17-MAR-1993               */
/* -------------------------------------------------------       */
#include <moscalib.h>
#include <proto_mos.h>

int match(verif, linid, linwav, liny, lindif, nbrow, 
        lcat, nbrowcat, alpha, stdres, dnull, reject)

double    lcat[], linid[], linwav[], liny[], lindif[], alpha;
double    dnull, *stdres;
int       reject[], verif, nbrow, nbrowcat;

{

   double          lambda_lin, lambda_cat;
   int             row, rowcat;
   int             order, ordref, bufsize;
   int             rowlow, rowupp, rowref, index;
   double          resref, reslow, resupp, residual;
   double          resref_abs, reslow_abs, resupp_abs, res_abs;
   double          lambda_ref, lambda_low, lambda_upp, lambda_next;
   int             nmatch=0, ldir, ldir_ref;
   int             not_bottom, not_top, debug=0, nbump=0;
   double          minlin, mincat;
   char            text[100];

   /* Verification loop */
   /* If the verification parameter verif<0, the line catalog is explored to
      make sure that wavelengths are sorted by increasing value and no 
      wavelength is duplicated. The value returned is >0 if no violation 
      was found */

   if (verif < 0) {
   verif = 1;
   lambda_low = lcat[1];
   lambda_next = lcat[2];

   ldir_ref = (lambda_next > lambda_low) ? 1 : -1;

   for (row=2; row<=nbrowcat; row++) {
       lambda_upp = lcat[row];
       ldir = (lambda_upp > lambda_low) ? 1 : -1;
       if (ldir != ldir_ref) {
          sprintf(text,"Warning: Column :WAVE of the line catalog is not sorted by increasing wavelength");
          SCTPUT(text);
          verif = -1;
	}
       if (lambda_low == lambda_upp) {
          sprintf(text,"Error: Column :WAVE of the line catalog contains duplicated wavelength : %f",lambda_low);
          SCTPUT(text);
          verif = -2;
	}
       lambda_low = lambda_upp;
     }
   return(verif);
 }

   /* Matching loop */

   rowcat  = 1;
   *stdres = 0.;

   /* Some assumptions are made concerning the input data :
        - Index of arrays are in the range 1 to N
        - All N values in the input arrays are valid
        - Line catalog is sorted by increasing wavelength and contain no duplicated lines 
   */

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

      reject[row] = NOT_PROCESSED; 
      linid[row]  = dnull;
      lindif[row] = dnull;

      /* Reads current wavelength to be identified and relative order number */
      lambda_lin = linwav[row];
      ordref = liny[row];
       if (debug == 1) {
           sprintf(text,"Reading line.tbl row = %d line = %f\n",row,lambda_lin);
          SCTPUT(text); 
	 }

          if (not_bottom = (row > 1 && liny[row-1] == ordref)) { 
              rowlow = row-1, lambda_low = linwav[rowlow], reslow = lambda_lin - lambda_low;
              reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
	    }

          if (not_top = (row < nbrow && liny[row+1] == ordref)) {
              rowupp = row+1, lambda_upp = linwav[rowupp], resupp = lambda_lin - lambda_upp;
              resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
	    }

          if (not_bottom && not_top)  minlin = (reslow_abs<resupp_abs) ? reslow_abs : resupp_abs;
          if (! not_bottom) minlin = resupp_abs;
          if (! not_top)    minlin = reslow_abs;

          /* Reads the line catalog at the first entry point and initializes variables */
          if (debug == 1) {
           sprintf(text,"Reading catalog thar.tbl rowcat = %d\n",rowcat);
           SCTPUT(text); 
	 }

          rowupp = rowcat,     lambda_upp = lcat[rowupp], resupp = lambda_lin - lambda_upp;
          resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */

          if (rowcat > 1) {
          rowlow = rowcat - 1, lambda_low = lcat[rowlow], reslow = lambda_lin - lambda_low;
          reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
          }
          else reslow_abs = resupp_abs + 1.;


          /* The line catalog is assumed to be sorted by increasing wavelength. 
             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) {

                ++nbump;
                 /* reject[row] = REACHED_LIMIT; */
                 break;
	       }
             else {
             lambda_cat = lcat[rowcat];
               if (debug == 1) {
           sprintf(text,"Matching loop at rowcat = %d lambda_cat = %f, lambda_lin = %f\n",
                      rowcat,lambda_cat,lambda_lin);
          SCTPUT(text); 
	 }
             residual = lambda_lin - lambda_cat;
             res_abs = (residual < 0) ?  -residual : residual; /* abs. value */
	   }
	   }

             rowref = rowcat - index;
               if (debug == 1) {
           sprintf(text,"Minimum found at %d, lambda_cat = %f, lambda_lin = %f, Residual = %f\n",
                         rowref, lambda_ref, lambda_lin, resref_abs);
          SCTPUT(text); 
	 }

          if (not_bottom = (rowref > 1)) { 
              rowlow = rowref-1, lambda_low = lcat[rowlow], reslow = lambda_ref - lambda_low;
              reslow_abs = (reslow < 0) ? -reslow : reslow; /* abs. value */
	    }

          if (not_top = (rowref < nbrowcat)) {
              rowupp = rowref+1, lambda_upp = lcat[rowupp], resupp = lambda_ref - lambda_upp;
              resupp_abs = (resupp < 0) ? -resupp : resupp; /* abs. value */
	    }

          if (not_bottom && not_top)  
             mincat = (reslow_abs<resupp_abs) ? reslow_abs : resupp_abs;
          if (! not_bottom) mincat = resupp_abs;
          if (! not_top)    mincat = reslow_abs;

          /* Find minimum of all minima */

          reslow = (minlin < mincat) ? minlin : mincat; 
          reslow *= alpha;

             if (debug == 1) {
           sprintf(text,"Comparing line %f with catalog %f.\nResidual is %f and maximum residual is %f\n",
             lambda_ref,lambda_lin,resref_abs,reslow);
          SCTPUT(text); 
	 }

      if (reject[row] == NOT_PROCESSED) {

          if (resref_abs < reslow ) {
                  linid[row] = lambda_ref, lindif[row] = resref;
                  reject[row] = SUCCESSFULL_MATCH, *stdres += resref_abs, nmatch++;
                   if (debug == 1) {
           sprintf(text,"**** Identified line: cat = %f  line = %f  Residual = %f\n",
                       lambda_ref,lambda_lin,resref);
          SCTPUT(text); 
	 }
		}
          else   reject[row] = MATCH_FAILED;

	} /* End of if (reject[row] == NOT_PROCESSED) */

     if (debug == 1) {
           sprintf(text,"Row %d Order %d - %f  %f  %f  %f  %f\n",row,ordref,minlin,resref_abs, lambda_ref, lambda_lin,mincat);
          SCTPUT(text); 
	 }


   }				/* End of for */

if (nmatch > 0)  *stdres /= (double) nmatch;

/*
if (nbump > 0) {
    if (nbump > 1) 
       sprintf(text,"Error: reached limit of line catalog (%d times)",nbump);
     else
       sprintf(text,"Error: reached limit of line catalog (%d time)",nbump);
  
     SCTPUT(text);
                }
*/

return(nmatch);

}