/* @(#)specorrel.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) 1992 European Southern Observatory           */
/* .IDENT       specorrel.c                                      */
/* .AUTHOR      Pascal Ballester,  ESO - Garching                */
/* .KEYWORDS    Spectroscopy,                                    */
/* .PURPOSE     Finds position of the correlation peak between   */
/*              the columns of two tables                        */
/* .VERSION     1.0    Creation    12-AUG-1992  PB               */
/* ------------------------------------------------------------- */

#include <tbldef.h>
#include <midas_def.h>
#define  MAXSIZE  10000

main()

{

    char   table1[61], table2[61], outima[61], xcol[21];
    char   refcol[21], wgtcol[21],text[81];
    int    tid1, tid2, xcol1, xcol2, sel1, sel2;
    int    refcol1, refcol2, wgtcol1, wgtcol2;
    int    pnb1, pnb2;
    double px1[MAXSIZE], px2[MAXSIZE], pref1[MAXSIZE];
    double pref2[MAXSIZE], pwgt1[MAXSIZE], pwgt2[MAXSIZE];
    int    reflog=0, wgtlog=0, vallog=0;
    int    cnt1, cnt2, pos;
    int    nbcol1, nbrow1, nsort1, allcol1, allrow1;
    int    nbcol2, nbrow2, nsort2, allcol2, allrow2;
    int    actvals, kunit, null;
    float  inputr[4], center, range;
    int    npix, imno;
    double start, step, shift;
    double dx, tol;
    char   ident[72], cunit[16], cheval[3];
    float  *pntr, maxcor, xpos;
    int    maxpos, pass;
    double pref, valref[2];

   SCSPRO("specorrel");
    
   SCKGETC("IN_A",   1,  60, &actvals, table1);
   SCKGETC("IN_B",   1,  60, &actvals, table2);
   SCKGETC("OUT_A",  1,  20, &actvals, outima);
   SCKGETC("INPUTC", 1,  20, &actvals, xcol);
   SCKGETC("INPUTC", 21, 20, &actvals, refcol);
   SCKGETC("INPUTC", 41, 20, &actvals, wgtcol);
   SCKGETC("INPUTC", 61,  1, &actvals, cheval);
   SCKRDR ("INPUTR", 1,  4,  &actvals, inputr, &kunit, &null);

   center = inputr[0];
   tol    = inputr[1];
   range  = inputr[2];
   step   = inputr[3];

   if (cheval[0] == '+')  {
     SCKRDD("INPUTD", 1, 2, &actvals, valref, &kunit, &null);
     vallog = 1;
   }

   TCTOPN(table1, F_I_MODE, &tid1);
   TCTOPN(table2, F_I_MODE, &tid2);
   TCIGET(tid1, &nbcol1, &nbrow1, &nsort1, &allcol1, &allrow1);
   TCIGET(tid2, &nbcol2, &nbrow2, &nsort2, &allcol2, &allrow2);
   TCCSER(tid1, xcol, &xcol1);
   TCCSER(tid2, xcol, &xcol2);

   /* reflog = 1 indicates that a reference column is used to 
      compare the items */
   if (refcol[0] != '+') {
       reflog = 1;
       TCCSER(tid1, refcol, &refcol1);
       TCCSER(tid2, refcol, &refcol2);
     }

   /* wgtlog = 1 indicates that a weight column is used for
      the incrementation */
   if (wgtcol[0] != '+') {
       wgtlog = 1;
       TCCSER(tid1, wgtcol, &wgtcol1);
       TCCSER(tid2, wgtcol, &wgtcol2);
     }


   /* Creates output image */

   strcpy(ident,"Cross-correlation image");
   strcpy(cunit,"Unitless");
   start  = center - range;
   npix   = 2.*range/step + 1;

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

    /* Initialisations */

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

    pnb1 = 0;
    for (cnt1 = 1; cnt1<=nbrow1; cnt1++) {
        TCSGET(tid1, cnt1, &sel1);
        if (sel1)  {
            if (reflog == 1) TCERDD(tid1, cnt1, refcol1, &pref1[pnb1], &null);
            if ((vallog == 1 && pref1[pnb1] >= valref[0] && 
                      pref1[pnb1] <= valref[1]) || (vallog == 0)) {
             TCERDD(tid1, cnt1, xcol1, &px1[pnb1], &null);
             if (wgtlog == 1) TCERDD(tid1, cnt1, wgtcol1, &pwgt1[pnb1], &null);
             pnb1++;
	   }
	  }
      }

    pnb2 = 0;
    for (cnt2 = 1; cnt2<=nbrow2; cnt2++) {
        TCSGET(tid2, cnt2, &sel2);
        if (sel2)  {
            if (reflog == 1) TCERDD(tid2, cnt2, refcol2, &pref2[pnb2], &null);
            if ((vallog == 1 && pref2[pnb2] >= valref[0] &&
                      pref2[pnb2] <= valref[1]) || (vallog == 0)) {
             TCERDD(tid2, cnt2, xcol2, &px2[pnb2], &null);
             if (wgtlog == 1) TCERDD(tid2, cnt2, wgtcol2, &pwgt2[pnb2], &null);
             pnb2++;
	   }
	  }
      }

    if (reflog != 1) {
       for (cnt1 = 0; cnt1<pnb1; cnt1++) pref1[cnt1] = 1.;
       for (cnt2 = 0; cnt2<pnb2; cnt2++) pref2[cnt2] = 1.;
     }

    if (wgtlog != 1) {
       for (cnt1 = 0; cnt1<pnb1; cnt1++) pwgt1[cnt1] = 1.;
       for (cnt2 = 0; cnt2<pnb2; cnt2++) pwgt2[cnt2] = 1.;
     }

    /* Cross-correlation loop */

    pref = pref1[0];

    for (cnt1 = 0; cnt1<pnb1; cnt1++) {

        pass = 0;

        for (cnt2 = 0; cnt2<pnb2; cnt2++) {

  	  if (pref1[cnt1] == pref2[cnt2]) {

            if (pass == 0) pass = 1;
            if (reflog == 1 && pref1[cnt1] != pref) {
                   pref = pref1[cnt1];
                   sprintf(text,"Reference value : %f\n",pref1[cnt1]);
                   SCTPUT(text);
		 }

            for (pos = 0; pos<npix; pos++) {

              shift = pos*step + start;
              dx = px1[cnt1] - px2[cnt2] + shift;
              if (dx < 0.) dx = (-1.)*dx;
              if (dx <= tol) pntr[pos] += pwgt1[cnt1]*pwgt2[cnt2];

	      }
	  }
          else
            if (pass == 1) break;
      }}

   /* Now find the maximum in the cross-correlation image */

   maxcor = pntr[0];
   maxpos = 0;
   for (pos = 1; pos<npix; pos++) {
       if (pntr[pos] > maxcor) {
          maxcor = pntr[pos];
          maxpos = pos;
	}
     }

   xpos = maxpos*step + start;

   /* Writes results and exit */

   SCKWRR("OUTPUTR",&xpos,1,1,&null);
   SCKWRR("OUTPUTR",&maxcor,2,1,&null);
   SCKWRI("OUTPUTI",&maxpos,1,1,&null);

   TCTCLO(tid1);
   TCTCLO(tid2);
   SCSEPI();

}