/* @(#)lnmatch.c	17.1.1.1 (ESO-IPG) 01/25/02 17:51:59 */
/* 

   Pascal Ballester (ESO/Garching)

   lnmatch.c

   line identification

*/

/* system includes */

#include <stdio.h>

/* general Midas includes */

#include <midas_def.h>

/* FEROS specific includes */

#include <lnmatch.h>

int match
#ifdef __STDC__
(
 int verif, double linid[], double linwav[], double liny[], double lindif[],
 int nbrow, double lcat[], int nbrowcat, double alpha, double *stdres,
 double dnull, int reject[]
 )
#else
     (
      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;
#endif
     
{
  double          lambda_lin, lambda_cat;
  int             row, rowcat;
  int             ordref;
  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 duplicate 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 = %7.1f, lambda_lin = %7.1f, Residual = %4.1f\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;
  
  return(nmatch);

}