/* @(#)necdef.c	14.2 (ESO-IPG) 05/23/00 18:07:19 */
/*===========================================================================
  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
===========================================================================*/

/* Program  : echdef.c                           */
/* Author   : P. Ballester  -  ESO Garching      */
/* Date     : 29.04.91                           */
/*                                               */
/* Purpose  :                                    */
/* Follow orders. Create and fill in table ORDER */
/*                                               */
/* Input    :                                    */
/*      IN_A      Name of Flat-Field frame       */
/*      IN_B      Name of input table            */
/*      OUT_A     Name of order table            */
/*      INPUTI(1) Step (pixels)                  */
/*      INPUTI(2) Lower scan limit               */
/*      INPUTI(3) Upper scan limit               */
/*      INPUTR(1) Hot threshold                  */
/*                                               */

#include <math.h>
#include <tbldef.h>
#include <midas_def.h>

#define  MAXORD    1000
#define  MID       500     /* MID = MAXORD/2  */
#define  NBROW     10000
#define  NBCOL     10

#define  ipos(col,row,siz)   row * siz + col  /* Pointer value */
#define  nint(f)             (int) (f + 0.49999)  /* Nearest integer, 0.4999 for HP-machines */

float    *pntra;
int      nrow, ncol, stkmin=MID, stkmax=MID;
int      tbrow=1, scan[2];
float    slope[MAXORD], intercept[MAXORD], fwhm[MAXORD], userthres[MAXORD];
int      ordsta[MAXORD],ordend[MAXORD], ordernum[MAXORD];
float    xstack[MAXORD], ystack[MAXORD];

main()
{
      char      frame[TEXT_LEN+1];
      char      inptab[TEXT_LEN+1], ordtab[TEXT_LEN+1];
      char      ident[TEXT_LEN+1], cunit[16*3 + 1], text[TEXT_LEN+1];

      int  order, nb_order, stepix, stkdif, stklimit, counter, ordlen;

      int  imnoa, naxis, npix[2], null, status, npar, kunit;
      int  iav, actvals, tid;
      int  ordcol, xcol, ycol, xcent;

      extern   int    nrow, ncol;
      extern   float  *pntra;
      extern   float  slope[], intercept[], fwhm[];
      extern   int ordsta[],ordend[];
      extern   int scan[];

      float     hot_thres, thres[2];
      double    ord_thres, mini_thres;
      double    EstimThresh();

      double    start[2], step[2];

      SCSPRO("echdef");

      SCKGETC ("IN_A", 1, 60,  &actvals, frame);
      SCKGETC ("IN_B", 1, 60,  &actvals, inptab);
      SCKGETC ("OUT_A", 1, 60, &actvals, ordtab);

      SCKRDI ("INPUTI", 1, 1, &actvals, &stepix,     &kunit, &null);
      SCKRDI ("INPUTI", 2, 2, &actvals, &scan[0],    &kunit, &null);  
      SCKRDR ("INPUTR", 1, 1, &actvals, &thres[0],   &kunit, &null);

      hot_thres = thres[0];

      if (scan[1]<scan[0]) {
             nrow = scan[1];
             scan[1] = scan[0];
             scan[0] = nrow;
	   }

      /* Fixed on  23/05/2000 to prevent going out of the frames boundaries */
      scan[0] = scan[0] + 2;
      scan[1] = scan[1] - 2;

      strcpy(ident, " ");
      strcpy(cunit, " ");
      SCIGET (frame, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 2,
              &naxis, npix, start, step, ident, cunit, (char **)&pntra,
              &imnoa);

      nrow  =  npix[1];
      ncol  =  npix[0];

      TCTINI (ordtab, F_TRANS, F_O_MODE, NBCOL, NBROW, &tid);

      TCCINI (tid, D_R4_FORMAT, 1, "I6",   "  ", "ORDER", &ordcol);
      TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "  ", "X",     &xcol);
      TCCINI (tid, D_R4_FORMAT, 1, "F7.1", "  ", "Y",     &ycol);

      nb_order = ReadInput (inptab);

      /* printf("Number of orders : %d\n",nb_order); */

      for (order=1; order<=nb_order; order++) {

           /* Defines the geometrical central position of the order */
           xcent     = WhereStart (slope[order], intercept[order], 
                       &ordsta[order], &ordend[order], &ordlen);


           /* Defines the absolute threshold of the order */
           if (userthres[order]>0.001)  {
               ord_thres  = userthres[order];
               mini_thres = userthres[order];
	     }
           else
               ord_thres = EstimThresh (order, nb_order, hot_thres, 
                       xcent, &mini_thres);


           /* printf ("xcent, ordlen, thres, mini: %d %d %f %f\n",xcent,
               ordlen,ord_thres,mini_thres); */

           /* The threshold must be low enough to find at least */
           /* three quarters of the positions aint  each order  */

           stklimit = 9*ordlen/stepix/10;
           stkdif   = stklimit - 1;       /* Initial value */
           counter  = 0;                  /* Number of loops (maxi 10) */

           while (stkdif < stklimit && counter < 3) {

             /* Follows the order to estimate order center and */
             /* background center                              */
             Follow     (order, slope[order], intercept[order], ord_thres,
                       hot_thres, stepix, xcent);

             stkdif = stkmax - stkmin + 1;
             counter++;

             if (stkdif < stklimit)      ord_thres = 0.75 * ord_thres;
             if (ord_thres < mini_thres) ord_thres = mini_thres;
	   }

           /* Display information line */
           sprintf (text,"Order: %3d  Threshold: %2f  Numb. of values: %d",
                    ordernum[order], ord_thres, stkdif);
           SCTPUT  (text);  


           /* Writes data in the output table */
           UpdateOut  (tid, ordernum[order], ordcol, xcol, ycol, stepix);

	 }

      SCDWRI(tid,"ORDSTA",&ordsta[1],1,nb_order,&kunit);
      SCDWRI(tid,"ORDEND",&ordend[1],1,nb_order,&kunit);

      TCTCLO(tid);

      SCSEPI();
    }


ReadInput (name)

char      name[];

{
     extern  float     slope[], intercept[], fwhm[], userthres[];

     int   tid, slope_col, inter_col, fwhm_col;
     int   order_col, thres_col;
     int   nbcol, nbrow, nsort, allcol, allrow;
     int   order, null, nb_order;

     TCTOPN (name, F_I_MODE, &tid);

     TCIGET (tid, &nbcol, &nbrow, &nsort, &allcol, &allrow);

     TCCSER (tid, ":SLOPE", &slope_col);
     TCCSER (tid, ":ORIG",  &inter_col);
     TCCSER (tid, ":FWHM",  &fwhm_col);
     TCCSER (tid, ":THRES", &thres_col);
     TCCSER (tid, ":ORDER", &order_col);

     nb_order = nbrow;

     for (order=1; order<=nb_order; order++) {
          TCERDR (tid, order, slope_col, &slope[order],     &null);
          TCERDR (tid, order, inter_col, &intercept[order], &null);
          TCERDR (tid, order, fwhm_col,  &fwhm[order],      &null);
          TCERDR (tid, order, thres_col, &userthres[order], &null);
          TCERDI (tid, order, order_col, &ordernum[order],  &null);
	}

     TCTCLO(tid);
     return(nb_order);
   } 


WhereStart (slope, interc, ordsta, ordend, ordlen)

    double  slope, interc;
    int     *ordsta, *ordend;
    int     *ordlen;

{

    int     xsta, ysta, xend, yend, xcent;
    extern  int   nbrow, nbcol;
    extern  int scan[2];

    ysta  = nint(interc);
    xsta  = (ysta > scan[0])    ?  1   : nint((scan[0]-interc)/slope);

    yend  = nint(slope*ncol+interc);
    xend  = (yend < scan[1])    ? ncol : nint((scan[1]-interc)/slope);

    xcent = (xsta + xend)/2;
    *ordlen = xend - xsta;
    *ordsta = xsta;
    *ordend = xend;
    return(xcent);

  }


double EstimThresh (order, nbord, hot_th, xcent, minith)

    int order, nbord, xcent;
    double   hot_th, *minith;

{

    extern   int    nrow, ncol;
    extern   float  *pntra;
    extern  float   slope[], intercept[], fwhm[];
    extern  int  scan[];

    int    row, col, pos, ylow, yupp, ycent, ymax, ymin;
    float       mini, maxi;
    double      threshold;

    ycent = nint (slope[order] * xcent + intercept[order] - 1);
    ymax  = scan[1] - 1;
    ymin  = scan[0] - 1;

    if (order < nbord) {
        yupp  = nint (slope[order+1] * xcent + intercept[order+1])-1;
        yupp  = (ycent + yupp)/2;
      }
    else
        yupp  = ymax;   /* order == nbord */

    if (yupp > ymax) yupp = ymax;

    if (order > 1) {
        ylow  = nint (slope[order-1] * xcent + intercept[order-1])-1;
        ylow  = (ylow + ycent)/2;
      }
    else
        ylow  = ymin;  /* order == 1 */

    if (ylow < ymin)  ylow = ymin;

/* Init */

    mini = maxi = pntra[ipos(xcent,ylow,ncol)];

/* Find the minimum and maximum value */

    for (row=ylow; row<yupp; row++) {
               pos = ipos(xcent,row,ncol);
               if (pntra[pos] > maxi) maxi   = pntra[pos];
               if (pntra[pos] < mini) mini = pntra[pos];
	    }


    threshold = (maxi - mini)*0.60 + mini;
    *minith   = (threshold - mini)*0.5 + mini;
    return(threshold);

}


Follow (order, slope, intercept, ord_thres, hot_thres, step, xprev)

    int  order, step, xprev;
    double    slope, intercept, ord_thres, hot_thres;

{

    int    cnt=0, xnext, xini, eof;
    int         direct=1, pos;
    double      ycent, yprev, ynext, yini;
    double      slope_ini, inter_ini;
    extern int  ncol;

    eof   = FindCenter (xprev, slope, intercept, ord_thres, hot_thres, 
                        direct, &ycent);

    Store(xprev, ycent, MID);

    xini = xprev;
    yini = ycent;
    slope_ini = slope;
    inter_ini = intercept;

    for (direct = -1; direct<=1; direct += 2) {

        xprev = xini;
        yprev = yini;
        pos   = MID;
        eof   = 0;      /* Init of EndOfFrame flag */
        slope     = slope_ini;
        intercept = inter_ini;

        while (!eof) {
           xnext = xprev + direct*step;
           eof = FindCenter (xnext, slope, intercept, ord_thres, hot_thres, 
                             direct, &ynext);
           if (!eof) {
              pos += direct;
              Store (xnext, ynext, pos);
              slope     = (ynext - yprev)/(xnext - xprev);
              intercept = (xnext*yprev - xprev*ynext)/(xnext - xprev);
              xprev     = xnext;
              yprev     = ynext;
	    }
         }
      }
   }


FindCenter (x, slop, inter, low_th, high_th, direct, ycent)

int x, direct;
double   slop, inter, low_th, high_th;
double   *ycent;

{

   float   ymid;
   int     pos, row, col, ylow, yupp, pass;
   double  Center();

   extern   int    nrow, ncol;
   extern   int  scan[2];
   extern   float  *pntra;

   ymid = slop*x+inter;

   if (1<=x && x<=ncol && scan[0]<=ymid && ymid<=scan[1]) {
       col = x - 1;
       row = nint(ymid) - 1;

       /* Find row of maximum signal connex to the center */
       pass=1;
       while(pass != 0) {
       pos  = ipos(col, row, ncol);
       pass = 0;
       if (pntra[ipos(col, (row+1), ncol)] > pntra[pos])  {pass=1; row++;}
       if (pntra[ipos(col, (row-1), ncol)] > pntra[pos])  {pass=1; row--;}
     }

        if (pntra[pos] > (float) low_th) {

            *ycent = Center(col,row,low_th,high_th,&ylow,&yupp)+1.;

            if (ylow > scan[0] && yupp < scan[1])
                 return(0); /* Centering correctly done */
            else return(1);
	  }

        else return(1); /* No more signal above the threshold */
}

    else  return(1); /* Reached edge of frame */

 }


double Center  (col, row, low_th, high_th, ylow, yupp)

    double  low_th, high_th;
    int     row, col, *ylow, *yupp;

{

    int     pos;
    double  pixel, ycent;
    double  ysum, sum;

    extern   int    nrow, ncol;
    extern   float  *pntra;

    /* Moves to the lower limit of the order */
    while (pntra[ipos(col, row, ncol)] >= low_th  && row > 0 )  row--;

    /* Computes center of gravity between lower and upper limit */
    sum  = 0.;
    ysum = 0.;
    *ylow = row++;
    while ((pixel = pntra[ipos(col, row, ncol)]) >= low_th && row < nrow) {

          if (pixel < high_th) {
                ysum += row*(pixel-low_th);
                sum  += (pixel-low_th);
	      }
          *yupp = row++;
	}

     ycent = ysum/sum;

     *ylow += 1;   /* Convert from array index to pixel location */
     *yupp += 1;

     return(ycent);
  }



Store (x, y, pos)

     int    x,pos;
     double y;

  {
     extern int   stkmin, stkmax;
     extern float xstack[], ystack[];

     xstack[pos] = x;
     ystack[pos] = y;

     if (pos < stkmin) stkmin = pos;
     if (pos > stkmax) stkmax = pos;

   }

UpdateOut (tid, order, ordcol, xcol, ycol, stepix)

     int  tid, order, ordcol, xcol, ycol, stepix;

{
     extern int  tbrow;
     extern int   stkmin, stkmax;
     extern float xstack[], ystack[];

     int    stkpntr, xempty=1;

     while(xempty < xstack[stkmin]) {
       TCEWRI  (tid, tbrow,   ordcol, &order);
       TCEWRI  (tid, tbrow++, xcol,   &xempty);
       xempty += stepix;
     }

     for (stkpntr=stkmin; stkpntr<=stkmax; stkpntr++) {

          TCEWRI (tid, tbrow, ordcol, &order);
          TCEWRR (tid, tbrow, xcol, &xstack[stkpntr]);
          TCEWRR (tid, tbrow, ycol, &ystack[stkpntr]);
          tbrow++;
	}

     xempty = xstack[stkmax] + stepix;
     while (xempty < ncol) {
        TCEWRI  (tid, tbrow,   ordcol, &order);
        TCEWRI  (tid, tbrow++, xcol,   &xempty);
        xempty += stepix;
      }

     TCEWRI  (tid, tbrow,   ordcol, &order);
     TCEWRI  (tid, tbrow++, xcol,   &ncol);

     stkmin = stkmax = MID;

   }

EstimThres2 (slope, interc)

    double  slope, interc;

{

    int     xsta, ysta, xend, yend, xcent;
    extern  int   nrow, ncol;
    extern  int  scan[2];

    ysta  = nint(interc);
    xsta  = (ysta > 1)    ?  1   : nint((scan[0]-interc)/slope);

    yend  = nint(slope*ncol+interc);
    xend  = (yend < nrow) ? ncol : nint((scan[1]-interc)/slope);

    xcent = (xsta + xend)/2;
    return(xcent);

  }