# include <stdio.h>
# include <stdlib.h>
# include <math.h>
# include <float.h>

# include <xtables.h>
# include <hstio.h>

# include "../stis.h"
# include "../stisdq.h"
# include "calstis6.h"


/* 
   Find the region in the image that is contaminated by geocoronal Lya.
   This region must be avoided by the cross-correlation procedure. This
   routine computes the avoidance window limits, in physical (image)
   units.

   The Lya approximate position is found using the CRVAL1, CRPIX1 and
   CD1_1 values read from the SCI extension header. The slit width in the
   dispersion direction is computed (approximately) from the slit width
   in arcsec stored in the slit structure, corrected by the plate scale
   in the cross-dispersion direction and the ratio of the two LTM values
   (the plate scale in the dispersion direction is in A/pixel).

   The Lya position is then refined by cross-correlating a rectangular
   function with same width as the slit, with the image data collapsed
   in cross-dispersion direction. The collapsing uses the DQ flag value 
   applicable for the cross-correlation procedure. The brightest peak
   in the cross-correlation result is accepted as the final Lya position,
   except if it differs from the original position by more than 50 pixels.
   In that case the second brightest peak is used, and the procedure is
   repeated iteratively. If the procedure fails down to the 5th brightest
   peak, the original Lya position is used.

   The Lya avoidance window is set to be 20 pixels wider than the projected
   slit width.






   Revision history:
   ----------------
   17 Dec 98  -  Implemented (I.Busko)

*/

void FindLya (StisInfo6 *sts, SingleGroup *in, ApInfo *slit, 
              int *avoid1, int *avoid2) {

/* arguments:
StisInfo6 *sts         i: calibration switches and info
SingleGroup *in	       i: input image
ApInfo *slit           i: slit width taken from here
int avoid1, avoid2;    o: avoidance region limits (in physical pixels)
*/

	float *image, *crosscor, norm, pval;
	int center, width, pindex;
	int i, j, image_index, box_index;

	/* Initialize to default. */
	*avoid1 = 0;
	*avoid2 = 0;

	/* Compute Lya approximate center and projected slit width. */
	center = (int)((1216. - sts->crval[0]) / sts->cd[0] + sts->crpix[0]);
	width  = (int)(slit->width[0] / sts->cd[1] / 3600. / sts->ltm[1] * 
                                                             sts->ltm[0]);

	/* If Lya not present, return immediately. */
	if ((center + width/2) < 0 ||
	    (center - width/2) > in->sci.data.nx)
	    return;

	/* Build collapsed version of image data. */
	image = (float *) calloc (in->sci.data.nx, sizeof (float));
	for (j = 0; j < in->sci.data.ny; j++) {
	    for (i = 0; i < in->sci.data.nx; i++) {
                if (!(DQPix (in->dq.data, i, j) & sts->cc_sdqflags))
	            image[i] += Pix (in->sci.data, i, j);
	    }
	}

	/* Cross-correlate. */
	crosscor = (float *) calloc (in->sci.data.nx, sizeof (float));
	for (i = 0; i < in->sci.data.nx; i++) {
	    norm = 0.0;
	    for (box_index = 0; box_index < width; box_index++) {
	        image_index = box_index - width/2 + i;
                if (image_index > -1 && image_index < in->sci.data.nx) {
	            crosscor[i] += image[image_index];
	            norm        += 1.0;
	        }
	    }
	    if (norm > 0.0)
	        crosscor[i] /= norm;
	}

	free (image);

	/* Find peak. */
	for (j = 1; j <= 5; j++) {
	    pval = - FLT_MAX;
	    for (i = 0; i < in->sci.data.nx; i++) {
	        if (crosscor[i] > pval) {
	            pval = crosscor[i];
	            pindex = i;
	        }
	    }
	    if (abs (pindex - center) > 50) {
	        crosscor[pindex] = - FLT_MAX;
	        continue;
	    } else {

	        /* Found acceptable peak. */
	        *avoid1 = pindex - width/2 - 10;
	        *avoid2 = pindex + width/2 + 10;
	        free (crosscor);
	        return;
	    }
	}

	/* Didn't find acceptable peak; use approximate center then. */
	*avoid1 = center - width/2 - 10;
	*avoid2 = center + width/2 + 10;

	free (crosscor);
}