/* This file contains:
	doBlev
	FitToOverscan
*/

# include <stdio.h>
# include <math.h>		/* sqrt */
# include <string.h>

# include <hstio.h>

# include "../acs.h"
# include "../acsinfo.h"
# include "../acserr.h"
# include "../acsdq.h"		/* for CALIBDEFECT */

static void FitToOverscan (SingleGroup *, int, int, int *, double, short);

/* This routine subtracts the bias determined from the overscan region
   on each line.  Only a portion of the trailing overscan region will be
   used to find the bias level.

   If an output text file for bias levels (third command-line argument)
   was specified, this routine also creates this file, writes the bias
   level for each image line, and then closes the file.

   The CRPIX and LTV keyword values in the output extension headers will be
   updated.

   NOTE:  This task now performs the subtractions in-place on x, resulting
   in output which still has the overscan regions.  

   Warren Hack, 1998 June 2:
   	Revised for ACS data... 
*/

int doBlev (ACSInfo *acs, SingleGroup *x, int chip, float *meanblev, 
		int *overscan, int *driftcorr) {

/* arguments:
ACSInfo *acs    i: calibration switches, etc
SingleGroup *x    io: image to be calibrated
int chip        i: chip number
float meanblev   o: mean value of bias levels that were subtracted
int *overscan     io: true if bias determined from an overscan region in x
int *driftcorr    o: true means correction for drift along lines was applied
*/

	extern int status;

	double biaslevel;	/* value to subtract from a line of image */
	double averagedrift;	/* drift averaged over all (output) columns */
	double sumbias;	    /* sum of bias levels (for getting mean) */
	int binx, biny;		/* bin factors */
	int trimx1, trimx2;	/* width to trim off ends of each line */
	int trimy1, trimy2;	/* amount to trim off ends of each column */
	int biassect[2];	/* section to use for finding bias level */
	int i, j;
	double di, dj;
	int endx, endy;		/* boundaries of actual image */
	int begx, begy;		/* boundaries of actual image */
	int sizex, sizey;
	int amp;			/* Counter for amps used */
	int numamps;		/* Number of AMPS used to readout chip */
	char *ccdamp; /* Amps which are used for this chip */
	int dodrift = YES;		/* use virtual overscan to get blevdrift? */
	int avgbias = NO;		/* average together trailing and leading bias */
	char *amploc; 
	int bias_loc;
	int bias_ampx, bias_ampy;
	int bias_orderx[4] = {0,1,0,1};
	int bias_ordery[4] = {0,0,1,1};
	double deval, biaslevela, biaslevelb;
	int npix;
	int too_few;
    int too_fewa, too_fewb;
		
	/* Function definitions */
	int BlevDrift (SingleGroup *, int *, int *, int, int *, int *, short);
	double DriftEval (double);
	double DriftMean (double);
	double BlevEval (double);
	int AvgBiasSect (char *, int, int *, int *);
	int FindBlev (SingleGroup *, int, int *, short, double *, int *);
	void parseWFCamps (char *, int, char *);

	/* Allocate space for ccdamp... */
	ccdamp = (char *) calloc (NAMPS+1, sizeof(char));

	/* initial values */
	*driftcorr = 0;
	*meanblev = 0.;
	
	binx = acs->bin[0];
	biny = acs->bin[1];

	if ((binx != 1 && binx != 2 && binx != 4) ||
	    (biny != 1 && biny != 2 && biny != 4)) {
	    trlerror ("(doBlev) bin size must be 1, 2, or 4.");
        free (ccdamp);
	    return (status = 1001);
	}

	/* Get biassect, the location of the region to use for determining
	   the bias level, and vx & vy, the location in the virtual
	   overscan to use for determining the drift with column number.
	   Also get the widths of the trim regions, the amount to remove
	   from each axis of the input when copying to output.
	  ** This routine is new to ACS...  It reads in overscan and trim
	  **	region specifications from an OSCNTAB, rather than having
	  **	it hardwired into the code as it was in CALSTIS. 
	*/
	if (*overscan != YES) {

		/* If no overscan region, subtract BIAS level from CCDTAB */
	    trlwarn ("Overscan region is too small to do BLEVCORR; ");
	    trlmessage ("           bias from CCDTAB will be subtracted.");
	    biaslevel = acs->ccdbias;
	    for (j = 0;  j < x->sci.data.ny;  j++) {
		  for (i = 0;  i < x->sci.data.nx;  i++)
		    Pix (x->sci.data,i,j) = Pix (x->sci.data,i,j) - biaslevel;
	    }
	    *meanblev = biaslevel;
        free (ccdamp);
	    return (status);
	}

	/* Determine bias levels from BIASSECT columns specified in OSCNTAB */
	
	/* Copy out overscan info for ease of reference in this function*/
	trimx1 = acs->trimx[0];
	trimx2 = acs->trimx[1];
	trimy1 = acs->trimy[0];
	trimy2 = acs->trimy[1];

	/* 
	biassect[0] = acs->biassecta[0];
	biassect[1] = acs->biassecta[1];
	*/
	/* Establish which amps from ccdamp string are appropriate 
		for this chip */
	ccdamp[0] = '\0';
	
	if (acs->detector == WFC_CCD_DETECTOR) {
		/* Set up the 2 AMPS used per line */ 
		parseWFCamps (acs->ccdamp, chip, ccdamp);

	} else {
		/* HRC observation; use full CCDAMP string */
		strcpy (ccdamp, acs->ccdamp);
	}

	/* How many amps are used for this chip */	
	numamps = strlen (ccdamp);
	
	/* Are we going to calculate drift in bias from virtual overscan? */
	/* If the end points of vx and vy are zero, no section was specified */
	if (acs->vx[1] == 0 && acs->vy[1] == 0 ) dodrift = NO;	

	/* For each amp used, determine bias level and subtract from 
		appropriate region of chip */
	for (amp = 0; amp < numamps; amp++) {
		bias_loc = 0;
		/* Determine whether leading and trailing sections
			are to be averaged for this amp... */
		avgbias = AvgBiasSect (ccdamp, amp, acs->biassecta, acs->biassectb);
		
		/* determine which section goes with the amp */
		/* bias_amp = 0 for BIASSECTA, = 1 for BIASSECTB */
		amploc = strchr (AMPSORDER, ccdamp[amp]);
		bias_loc = *amploc - *ccdamp;
		bias_ampx = bias_orderx[bias_loc];
		bias_ampy = bias_ordery[bias_loc];
		
		/* Requirement: at least 1 section should be specified!	 */
		/* If both bias sections are specified,... */
		if (acs->biassecta[1] != 0 && acs->biassectb[1] != 0) {
			/* select section nearest the amp based on bias_amp */
			biassect[0] = (bias_ampx == 0) ? acs->biassecta[0] : acs->biassectb[0];
			biassect[1] = (bias_ampx == 0) ? acs->biassecta[1] : acs->biassectb[1];
		} else {
			/* otherwise, select the non-zero bias section */
			biassect[0] = (acs->biassecta[1] == 0) ? acs->biassectb[0] : acs->biassecta[0];
			biassect[1] = (acs->biassecta[1] == 0) ? acs->biassectb[1] : acs->biassecta[1];
		}
	
		/* Compute range of pixels affected by each amp */	
		/* Original designations...
		** endx = x->sci.data.nx - trimx2;
		** endy = x->sci.data.ny - trimy2;
		** sizex = x->sci.data.nx - (trimx1 + trimx2);
		** sizey = x->sci.data.ny - (trimy1 + trimy2);
		*/
		begx = trimx1 + acs->ampx * bias_ampx;
		endx = (bias_ampx == 0 && acs->ampx != 0) ? acs->ampx + trimx1 : x->sci.data.nx - trimx2;
		begy = trimy1 + acs->ampy* bias_ampy;
		endy = (bias_ampy == 0 && acs->ampy != 0) ? acs->ampy +trimy1 : x->sci.data.ny - trimy2; 
		sizex = endx - begx;
		sizey = endy - begy;		
		
		/* Make sure that endx and endy do not extend beyond the bounds of the
			image... WJH 8 Sept 2000
		*/
		if (endx > x->sci.data.nx) endx = x->sci.data.nx;
		if (endy > x->sci.data.ny) endy = x->sci.data.ny;
		
		/* At this point, we decide how we are to determine the bias level...
			average two sections across the line for 1 amp, or use only
			one section... 
		*/

		if (dodrift == YES) {
			/* Fit a line to the virtual overscan region as a function of
			   column number.
			*/
			if (BlevDrift (x, acs->vx, acs->vy, trimx1, biassect, driftcorr,
			acs->sdqflags)) {
                free (ccdamp);
	    		return (status);
            }

			/* Evaluate the fit for each line, and subtract from the data. */
			averagedrift = DriftMean ((double)sizex);
			
		} else {
			averagedrift = 0;
		}
		dj = (double)begy;
		sumbias = 0.;
	
		if (avgbias == NO) {
		
			/* For each image line, determine the bias level from the overscan
			   in x, and fit to these values as a function of line number in
			   the output image.
			*/
			FitToOverscan (x, sizey, trimy1, biassect, acs->ccdbias,
			acs->sdqflags);

			for (j = begy;  j < endy;  j++) {
				dj = dj + 1.0;
	    		biaslevel = BlevEval (dj);
	    		/* bias for current line plus average drift (constant) */
	    		sumbias += (biaslevel + averagedrift);

				di = (double)begx;
	    		for (i = begx;  i < endx;  i++, di++){
					di = di + 1.0;
				if (dodrift == YES) {
					deval = DriftEval(di);
				} else {
					deval = 0;
				}
				Pix (x->sci.data,i,j) =
					Pix (x->sci.data,i,j) - biaslevel - deval;
				}
			}
		} else {
			/* Average two sections... */

			/* For each image line, determine the bias level from the overscan
			   in x.
			*/		
			too_few = 0;
            too_fewa = NO;
            too_fewb = NO;
			for (j = begy;  j < endy;  j++) {
			    if (FindBlev (x, j, acs->biassecta, acs->sdqflags, &biaslevela, &npix)) {
			   		too_fewa = YES;
			   		biaslevela = acs->ccdbias;
				}
				if (FindBlev (x, j, acs->biassectb, acs->sdqflags, &biaslevelb, &npix) ) {
					too_fewb = YES;
					biaslevelb = acs->ccdbias;
				}
                
                if ( (too_fewa == YES) || (too_fewb == YES) ) {
                    too_few++;
                    too_fewa = NO;
                    too_fewb = NO;
                }
                
				biaslevel = (biaslevela + biaslevelb)/2.;	
								 
				dj = dj + 1.0;
	    		/* bias for current line plus average drift (constant) */
	    		sumbias += (biaslevel + averagedrift);

				di = (double)begx;
	    		for (i = begx;  i < endx;  i++, di++){
					di = di + 1.0;

					if (dodrift == YES) {
						deval = DriftEval(di);
					} else {
						deval = 0;
					}
					Pix (x->sci.data,i,j) =
						Pix (x->sci.data,i,j) - biaslevel - deval;
				}
			}
			if (too_few > 0) {
	    		sprintf (MsgText, "(blevcorr) %d image line", too_few);
	    		if (too_few == 1)
					strcat (MsgText," has");
	    		else
					strcat (MsgText, "s have");
				strcat (MsgText," too few usable overscan pixels.");
				trlwarn (MsgText);
			}
		}	
	} /* End loop over AMPs */

	*overscan = YES;
	
	/* This is the mean value of all the bias levels subtracted. */
	*meanblev = sumbias / (double)sizey;
	
	/* free ccdamp space allocated here... */
    free (ccdamp);

	return (status);
}

/* This function determines the bias level from the overscan in the input
   image x, and it fits a straight line to these values as a function of
   line number in the output image.
*/

static void FitToOverscan (SingleGroup *x, int ny, int trimy1,
	int *biassect, double ccdbias, short sdqflags) {

/* arguments:
SingleGroup *x    i: input image
int ny            i: size of second axis after overscan has been removed
int trimy1        i: offset between input and output line numbers
int biassect[2]   i: section to use for bias level determination
double ccdbias    i: bias level to subtract if we can't get it from overscan
*/

	extern int status;

	double biaslevel;	/* bias level in one line of input image */
	int too_few = 0;	/* number of lines with too few good pixels */
	int npix;		/* number of pixels used to compute bias */
	int j;
	void BlevInit (int);
	void BlevAccum (int, double);
	int BlevFit (void);
	void BlevSet (double);
	int FindBlev (SingleGroup *, int, int *, short, double *, int *);

	BlevInit (ny / 2);			/* initialize for fitting */

	/* For each line, determine the bias level from the overscan in x.
	   Note that we loop over the number of pixels in the output image.
	   The argument j+trimy1 to FindBlev is the line number in the
	   input image x, with trimy1 the offset to the illuminated region.
	*/
	for (j = 0;  j < ny;  j++) {
	    if (FindBlev (x, j+trimy1, biassect, sdqflags, &biaslevel, &npix)) {
            too_few++;
            status = 0;			/* not fatal */
	    } else {
            /* Note:  j is line number in output image. */
            BlevAccum (j, biaslevel);	/* increment sums */
	    }
	}
	if (too_few > 0) {
	    sprintf (MsgText, "(blevcorr) %d image line", too_few);
	    if (too_few == 1)
			strcat (MsgText, " has");
	    else
			strcat (MsgText, "s have");
	    strcat (MsgText, " too few usable overscan pixels.");
		trlwarn (MsgText);
	}

	/* Fit a curve to the bias levels found. */
	if (BlevFit()) {
	    trlwarn ("No bias level data, or singular fit; ");
	    trlmessage ("            bias from CCDTAB will be subtracted.");
	    BlevSet (ccdbias);		/* assign the default value */
	}
}

/* This function determines whether one amp is used to read out an
	entire line or not.  If only one amp is used and both bias sections
	are specified, then avg the two regions together.   Otherwise, compute
	the bias from each section and apply it to the data separately.
*/
int AvgBiasSect (char *ccdamp, int amp, int *biassecta, int *biassectb) {

	int avgbias;
	char ampab[]="AB";
	char ampcd[]="CD";
	char amppair[2];
	int sum, i	;
	
	if (strchr(ampab,ccdamp[amp]) != NULL)
		strcpy(amppair, ampab);
	else
		strcpy(amppair, ampcd);

	sum = 0;
	for (i = 0; i < 2; i++) {
		if (strchr(ccdamp, amppair[i]) != NULL)
			sum++;
	}
	
	if (sum > 1) { 
		avgbias =  NO;
	 } else {
		/* 1 AMP only: if both sections are specified, avg them... */
		avgbias = (biassecta[1] !=0 && biassectb[1] != 0) ? YES: NO;
	}
	
	return(avgbias);
}