/* @(#)spematch.c	17.1.1.1 (ES0-DMD) 01/25/02 17:56:14 */
/*===========================================================================
  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               */
/*                     printf --> SCTPUT 16-SEP-1999  SW         */
/* ------------------------------------------------------------- */

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

/* All Global Definitions for the accumulation array */

#define  ipos2D(col,row,npix) row*npix[0]+col
#define  MLEN 255

char     name[60];
char     Msg[MLEN];
float    *pntr, maxval;
int      imno, npix;
double   start, step, stend;

double          *pntr_buf;
int             nbrow;

int    tmp_pos[2000];
double tmp_iso[2000];
double tmp_res[2000];

/* Starting the main code */

main()
{

   int             linwav, catwav, linid, lindif, liny, linerr, liniso;
   char            line[60], lincat[60];
   double          lambda_lin, lambda_cat, error, error_min;
   int             row, rowcat, null, pass;
   int             order, ordref, bufsize, imno;
   float           thres[3], fraction;
   int             status, kunit, actvals, tidlin, tidcat;
   int             nbcol, nsort, allcol, allrow;
   int             nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat;
   int             rowlow, rowupp, rowref, index, number;
   double          resref, reslow, resupp, residual, finres;
   double          resref_abs, reslow_abs, resupp_abs, res_abs;
   double          lambda_ref, lambda_low, lambda_upp;

   SCSPRO("echmatch");


   status = SCKRDR("INPUTR", 1, 1, &actvals, thres, &kunit, &null);
   status = SCKGETC("IN_A", 1, 60, &actvals, line);
   status = SCKGETC("IN_B", 1, 60, &actvals, lincat);

   /* Read all keywords for the accumulator */
   status = SCKGETC("OUT_A",1, 60, &actvals, name);
   status = SCKRDI("INPUTI", 1, 1, &actvals, &npix, &kunit, &null);
   status = SCKRDD("INPUTD", 1, 1, &actvals, &start,&kunit, &null);
   status = SCKRDD("INPUTD", 2, 1, &actvals, &step, &kunit, &null);
   status = SCKRDR("INPUTR", 2, 1, &actvals, &fraction, &kunit, &null);

   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.;
          break;
	}
       if (lambda_low == lambda_upp) {
          printf("Error: Column :WAVE of the table %s contains duplicated wavelength : %f\n",
                  lincat, lambda_low);
          SCKWRD("OUTPUTD",&lambda_low,1,1,&kunit);
          thres[0] = -2.;
          break;
	}
       lambda_low = lambda_upp;
     }
   status = SCKWRR("INPUTR", &thres[0], 1, 1, &kunit);
   SCSEPI();
   exit(0);
 }

   /* 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, ":ERROR", &linerr);
   TCCSER(tidlin, ":DISTANCE", &liniso);

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


   bufsize = 11*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 = 2L;
   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);

      /* Read the error on wavelength if the column exists */

      if (linerr > 0) 
        TCERDD(tidlin, row, linerr,   &error, &null);
      else
        error = -1.;

      /* 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];
          if (error < 0.) error = reslow;
          finres = pntr_buf[ipos(2,row,nbrow)];
 
          pntr_buf[ipos(7,row,nbrow)] = lambda_lin;
          pntr_buf[ipos(8,row,nbrow)] = reslow;
          pntr_buf[ipos(9,row,nbrow)] = error;

          pntr_buf[ipos(10,row,nbrow)] = lambda_ref;
          pntr_buf[ipos(11,row,nbrow)] = 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 row = */


   /* Determines the standard error on a fraction of the most isolated lines           TCEWRD(tidlin, row, liniso, &reslow); */

   number = (int) ((float) nbrow / 10); 
   if (number < 10 ) number = 10;
   stdev(number,&error_min);
   sprintf(Msg,"Std Dev on %d lines: %g wav. units \n",number,error_min);
   SCTPUT (Msg);

   open_acc();

   for (row = 1; row <= nbrow; row++) {
     lambda_lin = pntr_buf[ipos(7,row,nbrow)];

     if (lambda_lin != 0.) {
       finres = pntr_buf[ipos(2,row,nbrow)];
       reslow = pntr_buf[ipos(8,row,nbrow)];

       error  = pntr_buf[ipos(9,row,nbrow)];

       accumulate(log(finres/error),log(reslow/error));
     }}

   findmax(fraction); 
   close_acc();

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

   lambda_lin = pntr_buf[ipos(7,row,nbrow)];

   if (lambda_lin != 0.) {

    finres = pntr_buf[ipos(2,row,nbrow)];
    reslow = pntr_buf[ipos(8,row,nbrow)];
    error  = pntr_buf[ipos(9,row,nbrow)]*exp(maxval);
    /* error = 2.5*error_min; */

    lambda_ref = pntr_buf[ipos(10,row,nbrow)];
    resref     = pntr_buf[ipos(11,row,nbrow)];

   /* printf("Comparing line %f with catalog %f.\n ",lambda_ref,lambda_lin);
   printf("Residual is %f and maximum residual is %f\n",finres,reslow); */

   TCEWRD(tidlin, row, liniso, &reslow);
   TCEWRD(tidlin, row, lindif, &resref);
   if (finres < error && error <= reslow) { 
          TCEWRD(tidlin, row, linid,  &lambda_ref);
          /* printf("**** Identified line: %f  %f\n",lambda_ref,resref); */
 	}

  }} /* End of if lambda_lin and for (row = */

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



open_acc()

{

char  ident[81], cunit[33];
int   i;

strcpy(ident,"Error proportionality coefficient accumulator");
strcpy(cunit,"Frequency                       ");

SCIPUT(name,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,1,&npix,
       &start,&step,ident,cunit, (char **)&pntr, &imno);

stend = start + (npix-1)*step;

for (i=0; i<npix; i++) pntr[i] = 0.;

}


accumulate(min,max)

double min,max;

{

int    pixmin, pixmax, i;
double vmin, vmax;

if (max < min) return;

if (min < start) vmin = start;
else             vmin = min;

if (max > stend)   vmax = stend;
else             vmax = max;

pixmin = (vmin-start)/step;
pixmax = (vmax-start)/step;

for (i=pixmin; i<=pixmax; i++) pntr[i] += 1;

}


close_acc()
{SCFCLO(imno);}


findmax(fraction)

float fraction;

{

int    pixmax, i, pixopt;
double vmax, vopt;

vmax   = pntr[0];
pixmax = 0;

for (i=1; i<npix; i++) {
    if (pntr[i] > vmax) {
        vmax = pntr[i];
        pixmax = i;
      }
  }

maxval = start + (pixmax-1)*step;
/* printf("Found optimum at %f . Value: %f\n",maxval,pntr[pixmax]); */

vopt = fraction*pntr[pixmax];
pixopt = pixmax;

/* printf("Fraction: %f , Ref level: %f\n",fraction,vopt); */

while (vopt < pntr[pixopt] && pixopt < npix)  {
      pixopt++;
      /* printf("Current pos. %d Level: %f\n",pixopt,pntr[pixopt]); */
    }

maxval = start + (pixopt-1)*step;
/* printf("Selection level at %f . Value: %f\n",maxval,pntr[pixopt]); */

}

stdev(number, error)

int    number;
double *error;

{

int row, posval, nb, flag, i;
double maxval, val, std, res, lambda;

for (nb=1; nb<=number; nb++) {

maxval = -1.; /* Sort positive values */
posval = 0;

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

    flag = 0;

    lambda = pntr_buf[ipos(7,row,nbrow)];
    if (lambda < 1.) flag = 1;

    for(i=1; i<nb; i++) {
        if (row == tmp_pos[i])  {
           flag = 1;
	 }}

    if (flag == 0) {
      val = pntr_buf[ipos(8,row,nbrow)];
      res = pntr_buf[ipos(2,row,nbrow)];
      if (val > res && val > maxval){
        maxval = val; posval = row;
      }}
  }

tmp_pos[nb] = posval;
tmp_iso[nb] = maxval;
tmp_res[nb] = pntr_buf[ipos(2,posval,nbrow)];

/* printf("Max nb %d Pos %d Isol. %f  Res. %f \n",
    nb,posval,maxval,tmp_res[nb]); */

}

/* Now sorting tmp array by decreasing residual value */ 

{
int flag, tmpp;
double tmpr, tmpi;

flag = 1;

while (flag != 0) {
  flag =0;
  for(row=2; row<=number; row++) {
      if (tmp_res[row-1] < tmp_res[row]) {
        flag += 1;
        tmpr = tmp_res[row-1], tmpi = tmp_iso[row-1], tmpp = tmp_pos[row-1];
        tmp_res[row-1] = tmp_res[row];
        tmp_iso[row-1] = tmp_iso[row];
        tmp_pos[row-1] = tmp_pos[row];
        tmp_res[row] = tmpr, tmp_iso[row] = tmpi, tmp_pos[row] = tmpp;
      }
    }
}}

/*
std = 0.;
for (row=1; row<=number; row++) {
  res = tmp_res[row]; 
  std += res*res;  
}
*error = sqrt(std/(nb-1));
*/

*error = tmp_res[number/2];

/* Writing the residuals in a small auxiliary ascii file */

{
FILE   *fp;
fp = fopen("dat.dat","w");
for (row=1; row<=number; row++) 
    fprintf(fp,"%f %f\n",tmp_res[row],tmp_iso[row]);
fclose(fp);
}

}