Basic infrastructure changes needed to support saturation map. (#563)

* Basic infrastructure changes needed to support saturation map.

* Implementation of the use of a full image to flag full-well saturated pixels.
The new calibration reference file is associated with the SATUFILE FITS keyword.
Full-welll saturated pixels are no longer flagged in the DoDQI step.  The new
function, doFullWellSat(), which controls the flagging is executed after the
BLEVCORR and BIASCORR steps have been done.

* Added parentheses to decision equations to make the operator
precedence explicit.

* Replace some hardcoded values in dofwsat.c with defined variables, and
update the date of changes to correspond to a time close to the actual
release date.

* Forget to update the date in wf3version.h.

* Very minor changes to a comment and added a space for pound include.

pkg/wfc3/calwf3/Dates

@@ -1,3 +1,10 @@
+17-Apr-2023 - MDD - Version 3.7.0
+    -  Implementation to use an image to detect and flag full-well saturation versus a simple scalar.
+       The new routine, dofwsat.c, is similar to what has been done for calacs. The WFC3 implementation
+       is more complicated in that there are serial virtual overscan columns, as well as binned images,
+       to accommodate.  The detection/flagging occurs after blev and bias correction while the output is
+       still in counts.
+
 27-May-2021 - MDD - Version 3.6.2
     - Bug fix to address calwf3.e crashing (Abort trap: 6) when taking an existing *_ima.fits (IR) file
       (produced with FLATCORR=PERFORM and CRCORR=OMIT) with CRCORR set to PERFORM. This reentrant

pkg/wfc3/calwf3/History

@@ -1,3 +1,10 @@
+### 17-Apr-2023 - MDD - Version 3.7.0
+    -  Implementation to use an image to detect and flag full-well saturation versus a simple scalar.
+       The new routine, dofwsat.c, is similar to what has been done for calacs. The WFC3 implementation
+       is more complicated in that there are serial virtual overscan columns, as well as binned images,
+       to accommodate.  The detection/flagging occurs after blev and bias correction while the output is
+       still in counts.
+
 ### 27-May-2021 - MDD - Version 3.6.2
     - Bug fix to address calwf3.e crashing (Abort trap: 6) when taking an existing *_ima.fits (IR) file
       (produced with FLATCORR=PERFORM and CRCORR=OMIT) with CRCORR set to PERFORM. This reentrant

pkg/wfc3/calwf3/Updates

@@ -1,3 +1,10 @@
+Updates for Version 3.7.0 17-Apr-2023 - MDD
+    -  Implementation to use an image to detect and flag full-well saturation versus a simple scalar.
+       The new routine, dofwsat.c, is similar to what has been done for calacs. The WFC3 implementation
+       is more complicated in that there are serial virtual overscan columns, as well as binned images,
+       to accommodate.  The detection/flagging occurs after blev and bias correction while the output is
+       still in counts.
+
 Updates for Version 3.6.2 27-May-2021 - MDD
     - Bug fix to address calwf3.e crashing (Abort trap: 6) when taking an existing *_ima.fits (IR) file
       (produced with FLATCORR=PERFORM and CRCORR=OMIT) with CRCORR set to PERFORM. This reentrant

pkg/wfc3/calwf3/include/wf3dq.h

@@ -17,7 +17,7 @@
 # define UNSTABLE         32    /* IR unstable pixel */
 # define WARMPIX          64    /* warm pixel */
 # define BADBIAS         128    /* bad bias value */
-# define SATPIXEL        256    /* full-well or a-to-d saturated pixel */
+# define SATPIXEL        256    /* full-well saturated pixel */
 # define BADFLAT         512    /* bad flatfield value */
 # define TRAP           1024    /* UVIS charge trap, SINK pixel */
 # define SPIKE          1024    /* CR spike detected during cridcalc IR */ 

pkg/wfc3/calwf3/include/wf3info.h

@@ -86,6 +86,10 @@
     M. De La Pena 2020 March:
     Read the PCTERNOI value from the raw image header for possible 
     use in the CTE reduction.
+
+    M. De La Pena 2022 February:
+    Added the satmap variable which rendered the "saturate" variable
+    obsolete.  Removed "saturate" as a cleanup operation.
  
 */
 
@@ -153,7 +157,6 @@ typedef struct {
     float mean_gain;    /* mean actual gain of all amps */
     int ampx;           /* first column affected by amps on 2/4amp readout*/
     int ampy;           /* first row affected by amps on 2/4amp readout*/
-    float saturate;     /* CCD saturation level */
     int trimx[4];       /* Width of overscan to trim off ends of each line */
     int trimy[2];       /* Amount of overscan to trim off ends of each col */
     int vx[4];
@@ -227,6 +230,7 @@ typedef struct {
     RefImage nlin;	    /* non-linearity coefficients image */
     RefImage zoff;      /* super zero  */
     RefTab pctetab;     /*uvis CTE parameter table*/
+    RefImage satmap;    /* full-well saturation image */
 
 } WF3Info;
 

pkg/wfc3/calwf3/include/wf3version.h

@@ -4,7 +4,7 @@
 /* This string is written to the output primary header as CAL_VER. */
 
 
-# define WF3_CAL_VER "3.6.2 (May-27-2021)"
-# define WF3_CAL_VER_NUM "3.6.2"
+# define WF3_CAL_VER "3.7.0 (Apr-17-2023)"
+# define WF3_CAL_VER_NUM "3.7.0"
 
 #endif /* INCL_WF3VERSION_H */

pkg/wfc3/calwf3/lib/dodqi.c

@@ -121,6 +121,12 @@ static void FirstLast (double *, double *, int *, int *, int *, int *,
    H.Bushouse, 2007 Dec 21:
     Updated to use new rewrite of FirstLast routine provided by P. Hodge
     from calstis.
+
+   M. De La Pena, 2022 February
+    Removed flagging of full-well saturated pixels based upon a science
+    pixel value being greater than a defined scalar value.  Use of a
+    new full-well saturation image supersedes the functionality previously
+    done in this routine.
 */
 
 int doDQI (WF3Info *wf3, SingleGroup *x) {
@@ -154,7 +160,6 @@ SingleGroup *x    io: image to be calibrated; DQ array written to in-place
 	int i, j, i0, j0;	/* indexes for scratch array ydq */
 	int m, n;		/* indexes for data quality array in x */
 	short sum_dq;		/* for binning data quality array */
-	float sat;
 
 	int row;		/* loop index for row number */
 	int dimx, dimy;
@@ -178,25 +183,18 @@ SingleGroup *x    io: image to be calibrated; DQ array written to in-place
 	    return (status);
 
 	/* For the CCD, check for and flag saturation. */
-	sat = wf3->saturate;
 	if (wf3->detector != IR_DETECTOR) {
             dimx = x->sci.data.nx;
             dimy = x->sci.data.ny;
 
 	    for (j = 0;  j < dimy;  j++) {
-		 for (i = 0;  i < dimx;  i++) {
-		      /* Flag a-to-d saturated pixels with 2048 dq bit */
-		      if (Pix (x->sci.data, i, j) > ATOD_SATURATE) {
-			  sum_dq = DQPix (x->dq.data, i, j) | ATODSAT;
-			  DQSetPix (x->dq.data, i, j, sum_dq); /* atod sat */
-		      }
-		      /* Flag full-well or atod saturated pixels with 256 bit */
-		      if (Pix (x->sci.data, i, j) > sat || 
-			  Pix (x->sci.data, i, j) > ATOD_SATURATE) {
-			  sum_dq = DQPix (x->dq.data, i, j) | SATPIXEL;
-			  DQSetPix (x->dq.data, i, j, sum_dq);	/* saturated */
-		      }
-		 }
+			for (i = 0;  i < dimx;  i++) {
+				/* Flag a-to-d saturated pixels with 2048 dq bit */
+				if (Pix (x->sci.data, i, j) > ATOD_SATURATE) {
+					sum_dq = DQPix (x->dq.data, i, j) | ATODSAT;
+					DQSetPix (x->dq.data, i, j, sum_dq); /* atod sat */
+				}
+			}
 	    }
 	}
 

pkg/wfc3/calwf3/lib/getccdtab.c

@@ -23,7 +23,6 @@ typedef struct {
 	IRAFPointer cp_readnoise[NAMPS];
 	IRAFPointer cp_ampx;
 	IRAFPointer cp_ampy;
-	IRAFPointer cp_saturate;
 	IRAFPointer cp_pedigree;
 	IRAFPointer cp_descrip;
 	int nrows;			/* number of rows in table */
@@ -40,7 +39,6 @@ typedef struct {
 	float readnoise[NAMPS];
 	int ampx;
 	int ampy;
-	float saturate;
 } TblRow;
 
 static int OpenCCDTab (char *, TblInfo *);
@@ -112,6 +110,10 @@ static int CloseCCDTab (TblInfo *);
     21 Oct 2009: H.Bushouse:
 	Added computation of mean_gain to GetCCDTab. mean_gain is now used
 	in flatcorr steps for doing gain conversion.
+
+    16 Feb 2022: M. De La Pena:
+    The "saturate" variable became obsolete once the full-well saturation
+    map was implemented. Removed "saturate" as a cleanup operation.
 */
 
 int GetCCDTab (WF3Info *wf3, int dimx, int dimy) {
@@ -188,15 +190,14 @@ int     dimy      i: number of lines in exposure
 		}
 		wf3->mean_gain /= NAMPS;
 
-                /* For WFC3/UVIS exposures, correct ampx to match the actual
-		   size of the image for more seamless processing of subarrays.
-		   Leave ampx as is for WFC3/IR exposures. */
+        /* For WFC3/UVIS exposures, correct ampx to match the actual
+           size of the image for more seamless processing of subarrays.
+           Leave ampx as is for WFC3/IR exposures. */
 		if (wf3->detector == CCD_DETECTOR)
 		    wf3->ampx = (tabrow.ampx > dimx) ? dimx : tabrow.ampx;
 		else
 		    wf3->ampx = tabrow.ampx;
 		wf3->ampy = tabrow.ampy;		
-		wf3->saturate = tabrow.saturate;
 
 		break;
 	    }
@@ -276,7 +277,6 @@ static int OpenCCDTab (char *tname, TblInfo *tabinfo) {
 	c_tbcfnd1 (tabinfo->tp, "READNSED", &tabinfo->cp_readnoise[3]);
 	c_tbcfnd1 (tabinfo->tp, "AMPX", &tabinfo->cp_ampx);
 	c_tbcfnd1 (tabinfo->tp, "AMPY", &tabinfo->cp_ampy);
-	c_tbcfnd1 (tabinfo->tp, "SATURATE", &tabinfo->cp_saturate);
 
 	/* Initialize counters here... */
 	missing = 0;
@@ -308,7 +308,6 @@ static int OpenCCDTab (char *tname, TblInfo *tabinfo) {
 	if (tabinfo->cp_readnoise[3] == 0 ) { missing++; nocol[i] = YES;} i++;
 	if (tabinfo->cp_ampx == 0 ) { missing++; nocol[i] = YES;} i++;
 	if (tabinfo->cp_ampy == 0 ) { missing++; nocol[i] = YES;} i++;
-	if (tabinfo->cp_saturate == 0) { missing++; nocol[i] = YES;} i++;
 
 	if (PrintMissingCols (missing, NUMCOLS, nocol, colnames, "CCDTAB", tabinfo->tp) )
 		return(status);
@@ -437,10 +436,6 @@ static int ReadCCDTab (TblInfo *tabinfo, int row, TblRow *tabrow) {
 	if (c_iraferr())
 	    return (status = TABLE_ERROR);
 
-	c_tbegtr (tabinfo->tp, tabinfo->cp_saturate, row, &tabrow->saturate);
-	if (c_iraferr())
-	    return (status = TABLE_ERROR);
-
 	return (status);
 }
 

pkg/wfc3/calwf3/lib/wf3info.c

@@ -71,6 +71,11 @@
   M. De La Pena 2020 March:
     Added reading of PCTERNOI keyword from the raw primary header for 
     possible use in the CTE reduction
+
+  M. De La Pena 2022 February:
+    Added variable, satmap - the reference image for full-well
+    saturation. Use of this image rendered wf3->saturate variable
+    obsolete.  Removed "wf3->saturate" as part of the cleanup operation.
 */
 
 void WF3Init (WF3Info *wf3) {
@@ -128,7 +133,6 @@ void WF3Init (WF3Info *wf3) {
 	}
 	wf3->ampx = 0;
 	wf3->ampy = 0;
-	wf3->saturate = 0.;
 	wf3->trimx[0] = 0;
 	wf3->trimx[1] = 0;
 	wf3->trimx[2] = 0;
@@ -185,6 +189,7 @@ void WF3Init (WF3Info *wf3) {
 	InitRefTab (&(wf3->oscn));
 	InitRefTab (&(wf3->atod));
     InitRefTab (&(wf3->pctetab));
+    InitRefImg (&(wf3->satmap));
 
 	/* Initialize reference images and tables for WF32D */
 	InitRefImg (&(wf3->dark));

pkg/wfc3/calwf3/lib/wscript

@@ -30,7 +30,7 @@ def build(bld):
             ../wf3ccd/doatod.c   ../wf3ccd/dobias.c ../wf3ccd/doblev.c
             ../wf3ccd/doccd.c   ../wf3ccd/doflash.c ../wf3ccd/findblev.c
             ../wf3ccd/findover.c ../wf3ccd/getflags.c
-            ../wf3ccd/getccdsw.c ../wf3ccd/sink.c
+            ../wf3ccd/getccdsw.c ../wf3ccd/sink.c ../wf3ccd/dofwsat.c
 
             ../wf3ir/blevcorr.c   ../wf3ir/cridcalc.c   ../wf3ir/darkcorr.c
             ../wf3ir/doir.c   ../wf3ir/dqicorr.c   ../wf3ir/flatcorr.c

pkg/wfc3/calwf3/wf3ccd/doccd.c

@@ -55,11 +55,15 @@
     Moved the subarray codes to a general function and changed the calculation
     of where the subarray is. Also moved the DQ bit to the singlegroup extension
 
+    M. De La Pena February 2022
+    Modified to apply the full-well saturation flags stored as an image
+    to the data in doFullWellSat() instead of in the doDQI step.
+
  */
 
 # include <string.h>
 # include <stdio.h>
-#include "hstcal.h"
+# include "hstcal.h"
 # include "hstio.h"
 # include "wf3.h"
 # include "wf3info.h"
@@ -235,7 +239,6 @@ int DoCCD (WF3Info *wf3, int extver) {
         if (dqiHistory (wf3, x.globalhdr))
             return (status);
 
-
     /* ANALOG TO DIGITAL CORRECTION. */
     AtoDMsg (wf3, extver);
     if (wf3->atodcorr == PERFORM) {
@@ -313,6 +316,27 @@ int DoCCD (WF3Info *wf3, int extver) {
         if (biasHistory (wf3, x.globalhdr))
             return (status);
 
+    /* Apply the saturation image.
+       Strictly speaking, the application of the full-well saturation image is
+       not a calibration step (i.e., there is no SATCORR), but the application
+       of a 2D image to flag pixels versus using a single scalar to flag
+       saturated pixels as previously done in DQICORR will be done in doFullWellSat()
+       after BLEVCORR and BIASCORR.  This correction should only be done if both
+       BLEVCORR and BIASCORR have been performed.  This flagging is only applicable
+       for the UVIS. */
+
+    if (wf3->biascorr == PERFORM && wf3->blevcorr == PERFORM) {
+        SatMsg (wf3, extver);
+        sprintf(MsgText, "\nFull-well saturation flagging being performed.");
+        trlmessage(MsgText);
+        if (doFullWellSat(wf3, &x)) {
+            return (status);
+        }
+    } else {
+        sprintf(MsgText, "\nNo Full-well saturation flagging being performed.\n");
+        trlwarn(MsgText);
+    }
+
     /*UPDATE THE SINK PIXELS IN THE DQ MASK OF BOTH SCIENCE IMAGE SETS
      IT'S DONE HERE WITH ONE CALL TO THE FILE BECAUSE THEY NEED TO BE
      PROCESSED IN THE RAZ FORMAT JAY USES, THOUGH ONLY ONE CHIP DONE HERE
@@ -340,8 +364,8 @@ int DoCCD (WF3Info *wf3, int extver) {
                 */
 
                 /* Make temporary full group
-                 use the sinkfile as the reference
-                sub2full will reset the pixel values
+                   use the sinkfile as the reference
+                   sub2full will reset the pixel values
                 */
                 initSingleGroup(&fullarray);
                 allocSingleGroup(&fullarray,RAZ_COLS/2,RAZ_ROWS, True);
@@ -551,6 +575,20 @@ static void BiasMsg (WF3Info *wf3, int extver) {
     }
 }
 
+static void SatMsg (WF3Info *wf3, int extver) {
+
+    int OmitStep (int);
+    void PrSwitch (char *, int);
+    void PrRefInfo (char *, char *, char *, char *, char *);
+
+    trlmessage ("");
+    if (extver == 1 && !OmitStep (wf3->biascorr)) {
+
+        PrRefInfo ("satufile", wf3->satmap.name, wf3->satmap.pedigree,
+                wf3->satmap.descrip, "");
+    }
+}
+
 static void FlashMsg (WF3Info *wf3, int extver) {
 
     int OmitStep (int);

pkg/wfc3/calwf3/wf3ccd/doccd.h

@@ -4,6 +4,7 @@
 /*USEFUL LIB FUNCTIONS*/
 static void AtoDMsg (WF3Info *, int);
 static void BiasMsg (WF3Info *, int);
+static void SatMsg (WF3Info *, int);
 static void FlashMsg (WF3Info *, int);
 static void BlevMsg (WF3Info *, int);
 static void dqiMsg  (WF3Info *, int);
@@ -13,6 +14,7 @@ int GetCorner (Hdr *, int , int *, int *);
 int doAtoD (WF3Info *, SingleGroup *);
 int atodHistory (WF3Info *, Hdr *);
 int doBias (WF3Info *, SingleGroup *);
+int doFullWellSat(WF3Info *, SingleGroup *);
 int biasHistory (WF3Info *, Hdr *);
 int doFlash (WF3Info *, SingleGroup *, float *);
 int flashHistory (WF3Info *, Hdr *);