From c0235e90802b975a2b59a18f1a38dbd19d8d0c9c Mon Sep 17 00:00:00 2001 From: Pey Lian Lim <2090236+pllim@users.noreply.github.com> Date: Tue, 10 Mar 2026 18:05:33 -0400 Subject: [PATCH 1/5] Implement X-CTE based on code from Dr Jay Anderson. Bump calwf3 version and add history notes. Minor clean-ups for Y-CTE code to make it easier to read and have consistent behavior with the rest of the pipeline. --- pkg/wfc3/Dates | 3 + pkg/wfc3/History | 3 + pkg/wfc3/Updates | 3 + pkg/wfc3/include/cte.h | 1 + pkg/wfc3/include/wf3version.h | 4 +- pkg/wfc3/lib/CMakeLists.txt | 1 + pkg/wfc3/lib/wf3cte/cte_serial_corr.c | 3032 +++++++++++++++++++++++++ pkg/wfc3/lib/wf3cte/wf3cte.c | 468 ++-- 8 files changed, 3289 insertions(+), 226 deletions(-) create mode 100644 pkg/wfc3/lib/wf3cte/cte_serial_corr.c diff --git a/pkg/wfc3/Dates b/pkg/wfc3/Dates index de3acdf9e..1c0ba0b7b 100644 --- a/pkg/wfc3/Dates +++ b/pkg/wfc3/Dates @@ -1,3 +1,6 @@ +12-Mar-2026 - PLL - Version 3.8.0 + - New serial (X) CTE correction for UVIS fullframe data. + 07-Jan-2026 - PLL - Version 3.7.3 - Bug fix to remove per-second in BUNIT for CRREJ output for IR data. diff --git a/pkg/wfc3/History b/pkg/wfc3/History index 751e0b30b..d86ec3d85 100644 --- a/pkg/wfc3/History +++ b/pkg/wfc3/History @@ -1,3 +1,6 @@ +### 12-Mar-2026 - PLL - Version 3.8.0 + - New serial (X) CTE correction for UVIS fullframe data. + ### 07-Jan-2026 - PLL - Version 3.7.3 - Bug fix to remove per-second in BUNIT for CRREJ output for IR data. diff --git a/pkg/wfc3/Updates b/pkg/wfc3/Updates index 2ef336487..aa75525e8 100644 --- a/pkg/wfc3/Updates +++ b/pkg/wfc3/Updates @@ -1,3 +1,6 @@ +Updates for Version 3.8.0 12-Mar-2026 - PLL + - New serial (X) CTE correction for UVIS fullframe data. + Updates for Version 3.7.3 07-Jan-2026 - PLL - Bug fix to remove per-second in BUNIT for CRREJ output for IR data. diff --git a/pkg/wfc3/include/cte.h b/pkg/wfc3/include/cte.h index e30e34fbb..1a4e31dee 100644 --- a/pkg/wfc3/include/cte.h +++ b/pkg/wfc3/include/cte.h @@ -94,6 +94,7 @@ int cteHistory (WF3Info *, Hdr *); int free_array(float **ptr, int rows, int columns); int GetCTESwitch (WF3Info *, Hdr *); int initCTETrl (char *, char *); +int sub_xctecor(SingleGroup *, double); int makeRAZ(SingleGroup *, SingleGroup *, SingleGroup *); int undoRAZ(SingleGroup *, SingleGroup *, SingleGroup *); diff --git a/pkg/wfc3/include/wf3version.h b/pkg/wfc3/include/wf3version.h index 3a7c5e050..34d43e3ee 100644 --- a/pkg/wfc3/include/wf3version.h +++ b/pkg/wfc3/include/wf3version.h @@ -4,7 +4,7 @@ /* This string is written to the output primary header as CAL_VER. */ -# define WF3_CAL_VER "3.7.3 (Jan-07-2026)" -# define WF3_CAL_VER_NUM "3.7.3" +# define WF3_CAL_VER "3.8.0 (Mar-12-2026)" +# define WF3_CAL_VER_NUM "3.8.0" #endif /* INCL_WF3VERSION_H */ diff --git a/pkg/wfc3/lib/CMakeLists.txt b/pkg/wfc3/lib/CMakeLists.txt index a59725724..2c9760940 100644 --- a/pkg/wfc3/lib/CMakeLists.txt +++ b/pkg/wfc3/lib/CMakeLists.txt @@ -98,6 +98,7 @@ add_library(${PROJECT_NAME} SHARED wf3ccd/sink.c wf3ccd/wf3ccd.c wf3cte/cte_dobias.c + wf3cte/cte_serial_corr.c wf3cte/getcteflags.c wf3cte/getctepars.c wf3cte/getctesw.c diff --git a/pkg/wfc3/lib/wf3cte/cte_serial_corr.c b/pkg/wfc3/lib/wf3cte/cte_serial_corr.c new file mode 100644 index 000000000..70ed3a965 --- /dev/null +++ b/pkg/wfc3/lib/wf3cte/cte_serial_corr.c @@ -0,0 +1,3032 @@ +/* WFC3 -- Serial CTE (X-CTE) loss correction for UVIS + + P.L. Lim Mar-2026 Adapted for the pipeline from Jay Anderson's + CTE correction code for WFC3 UVIS raw2rzx_wfc3uv.c . + Also see WFC3 ISR 2024-07. No OpenMP support at this time. +*/ +#include "hstio.h" +#include "wf3.h" + +/* This was translated from Fortran that is 1-indexed. + Jay opted to not use the 0 index here, so the + size here is Fortran size + 1. */ +#define XCTE_TRAPS 2103 +#define NUM_TRAP (XCTE_TRAPS + 1) +#define NUM_REG (XCTE_TRAPS + XCTE_TRAPS + 1) /* Includes extra virtual line that comes in from the left side. */ +#define NUM_ELEVEL 611 /* 610 + 1 */ + +/* Local internal functions. */ +static int sub_raz2rzx_wfc3uv(float *, float *,float); +static void zero_reg(float [NUM_REG]); +static void push_reg(float [NUM_REG], float *); +static void read_reg(float [NUM_REG], float *, int [NUM_ELEVEL], float [NUM_TRAP], float [NUM_TRAP], float [NUM_TRAP]); +static void init_e_kn(int [NAMPS][NUM_ELEVEL]); +static void init_ecr_t(float [NUM_TRAP], float [NUM_TRAP], float [NUM_TRAP]); + +/* This was the main routine in raw2rzx_wfc3uv.c + and it is included in cte.h for use in wfc3cte.c . + Jay hardcoded everything and did not rely on PCTETAB + so we also cannot here. + + We do not compute difference image (pixz_d) here + because it is done in wfc3cte.c after Y-CTE anyway. + + raz is both the input and output. When this is called, + raz should have already been bias subtracted and converted + to unit of electrons by applying gain correction. + */ +int sub_xctecor(SingleGroup *raz, double mjd_obs) { + extern int status; + int i, j, k; + float *pixz_i = malloc(RAZ_COLS * RAZ_ROWS * sizeof(float)); + float *pixz_o = malloc(RAZ_COLS * RAZ_ROWS * sizeof(float)); + + /* These values are independent of PCTETAB. */ + float ff_xcte = (mjd_obs - 54962.00) / (60125.0 - 54962.00); + + trlmessage("XCTE_FF: %8.3f\n", ff_xcte); + + /* Make a copy of raz data and use as input. */ + for(j=0; jsci.data, i, j); + } + } + + if (sub_raz2rzx_wfc3uv(pixz_i, pixz_o, ff_xcte)) { + trlerror("X-CTE correction failed! Data is not corrected for X-CTE."); + free(pixz_i); + free(pixz_o); + return (status); + } + + /* Stuff corrected pixels back to raz for Y-CTE next. */ + for(j=0; jsci.data, i, j) = pixz_o[k]; + } + } + + free(pixz_i); + free(pixz_o); + return (status); +} + +/* This is the workhorse function that converts the + uncorrected RAZ image into a corrected RAZ image; + recall that RAZ has the amplifiers CDAB. +*/ +static int sub_raz2rzx_wfc3uv(float *pixz_i, float *pixz_o, float ff_xcte) { + extern int status; + int i, j, k, kk=-1, m; + + /* This array stores the trap profile: The location where there + are 1-electron traps; "n" here refers to the trap number. + Since this was translated from Fortran, 0 index is not used. */ + int e_kn[NAMPS][NUM_ELEVEL]; + + /* "t" in these arrays refers to "trail" */ + float e_t[NUM_TRAP], c_t[NUM_TRAP], r_t[NUM_TRAP]; + + /* To do the XCTE correction, things have to be divided by amp */ + float *amp_i, *amp_o, *amp_d; + + /* This represents the contents of the readout register */ + float reg[NUM_REG]; + + /* Fill array that gives electron-cloud-size that is read. */ + init_e_kn(e_kn); + + /* To activate trap#n for the four different amps, + this lists the trap-release profile (same for all amps). */ + init_ecr_t(e_t, c_t, r_t); + + /* Allocate memory for one amplifier for input, output, and difference, respectively. */ + if ((amp_i = malloc(XCTE_TRAPS * RAZ_ROWS * sizeof(float))) == NULL) { + trlerror("Out of memory for amp_i"); + return (status = OUT_OF_MEMORY); + } + if ((amp_o = malloc(XCTE_TRAPS * RAZ_ROWS * sizeof(float))) == NULL) { + trlerror("Out of memory for amp_o"); + free(amp_i); + return (status = OUT_OF_MEMORY); + } + if ((amp_d = malloc(XCTE_TRAPS * RAZ_ROWS * sizeof(float))) == NULL) { + trlerror("Out of memory for amp_d"); + free(amp_i); + free(amp_o); + return (status = OUT_OF_MEMORY); + } + + /* Go through the four amps one by one */ + for (k=0; k pmax) { + pmax = reg[i]; + } + } + + nmax = 1; + for (n=1; n e_n[n]) { + nmax = n; + } + } + + for (n=nmax; n>0; n--) { + f = 0; + t = 0; + for (i=1; i e_n[n]) { + /* Flush what's left in the trap into this pixel. */ + if (t < NUM_TRAP) { + reg2[i] += f * r_t[t]; + } + /* Reset the trap. */ + t = 0; + /* How much charge is in the trep? (prop to number of i-shifts to reg) + + Jay: 2025 is the number of shifts I calibrated the model to. + Not sure why it wasn’t 2048 (the whole chip), but it’s what I used. + So, f refers to what fraction of the “total” number of serial traps + a particular pixel in a particular column is subject to. + */ + f = i / 2025.0; + /* The first 11 traps count as one, to prevent large shifts in the noise. */ + if (n <= 11) { + f /= 11.0; + } + /* Remove from pixel the charge that gets trapped. */ + rego[i] -= f * c_t[XCTE_TRAPS]; + } + /* The output reg value is the orig value plus what is released. */ + regu[i] += reg2[i] + rego[i]; + } + } + + for (i=1; i @@ -27,7 +30,6 @@ # include # include # include -# include # ifdef _OPENMP # include @@ -80,8 +82,7 @@ float find_raz2rnoival(float *, float *); int WF3cte (char *input, char *output, CCD_Switch *cte_sw, - RefFileInfo *refnames, int printtime, int verbose, int onecpu) { - + RefFileInfo *refnames, int printtime, int verbose, int onecpu) { /* input: filename output: filename @@ -150,8 +151,8 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, float hardset=0.0; /* These are used to find subarrays with physical overscan */ - int sci_bin[2]; /* bin size of science image */ - int sci_corner[2]; /* science image corner location */ + int sci_bin[2]; /* bin size of science image */ + int sci_corner[2]; /* science image corner location */ int ref_bin[2]; int ref_corner[2]; int rsize = 1; /* reference pixel size */ @@ -162,9 +163,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, initHdr(&phdr); initHdr(&scihdr); - float readNoise = 0.0; - - int ret; + float readNoise=0.0; /*check if this is a subarray image. This is necessary because the CTE routine will start with the raw images @@ -194,7 +193,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, /*CONTAIN PARALLEL PROCESSING TO A SINGLE THREAD AS USER OPTION*/ # ifdef _OPENMP trlmessage("Using parallel processing provided by OpenMP inside CTE routine"); - if (onecpu){ + if (onecpu) { omp_set_dynamic(0); max_threads=1; trlmessage("onecpu == TRUE, Using only %i threads/cpu", max_threads); @@ -206,20 +205,20 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, omp_set_num_threads(max_threads); # endif - /* COPY COMMAND-LINE ARGUMENTS INTO WF3. */ WF3Init (&wf3); /*sets default information*/ strcpy (wf3.input, input); strcpy (wf3.output, output); PrBegin ("WFC3CTE"); - if (wf3.printtime) + if (wf3.printtime) { TimeStamp("WFC3CTE Started: ",wf3.rootname); + } /* CHECK WHETHER THE OUTPUT FILE ALREADY EXISTS. */ - if (FileExists (wf3.output)){ + if (FileExists(wf3.output)) { WhichError(status); - return (ERROR_RETURN); + return status; } wf3.pctecorr = cte_sw->pctecorr; @@ -233,44 +232,44 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, PrFileName ("input", wf3.input); PrFileName ("output", wf3.output); - if (wf3.biascorr == COMPLETE){ - trlmessage("BIASCORR complete for input image, CTE can't be performed"); - return(ERROR_RETURN); + if (wf3.biascorr == COMPLETE) { + trlerror("BIASCORR complete for input image, CTE can't be performed"); + return (status=ERROR_RETURN); } - if (wf3.darkcorr == COMPLETE){ - trlmessage("DARKCORR complete for input image, CTE can't be performed"); - return(ERROR_RETURN); + if (wf3.darkcorr == COMPLETE) { + trlerror("DARKCORR complete for input image, CTE can't be performed"); + return (status=ERROR_RETURN); } - if (wf3.blevcorr == COMPLETE){ - trlmessage("BLEVCORR complete for input image, CTE can't be performed"); - return(ERROR_RETURN); + if (wf3.blevcorr == COMPLETE) { + trlerror("BLEVCORR complete for input image, CTE can't be performed"); + return (status=ERROR_RETURN); } /* DETERMINE THE NAMES OF THE TRAILER FILES BASED ON THE INPUT AND OUTPUT FILE NAMES, THEN INITIALIZE THE TRAILER FILE BUFFER WITH THOSE NAMES. */ - if (initCTETrl (input, output)) - return (status); + if (initCTETrl(input, output)) { + return status; + } /* OPEN INPUT IMAGE IN ORDER TO READ ITS PRIMARY HEADER. */ - if (LoadHdr (wf3.input, &phdr) ){ + if (LoadHdr(wf3.input, &phdr)) { WhichError(status); - return (ERROR_RETURN); + return status; } /* GET KEYWORD VALUES FROM PRIMARY HEADER. */ - if (GetKeys (&wf3, &phdr)) { + if (GetKeys(&wf3, &phdr)) { freeHdr (&phdr); - return (status); + return status; } - if (GetCTEFlags (&wf3, &phdr)) { + if (GetCTEFlags(&wf3, &phdr)) { freeHdr(&phdr); - return (status); + return status; } - /*SET UP THE ARRAYS WHICH WILL BE PASSED AROUND*/ initSingleGroup(&raz); allocSingleGroup(&raz, RAZ_COLS, RAZ_ROWS, True); @@ -291,8 +290,8 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, allocSingleGroup(&chg, RAZ_COLS, RAZ_ROWS, True); /*hardset the science arrays*/ - for (i=0;i= P_OVRSCN && finish + V_OVRSCNX2 <= (RAZ_COLS/2) - P_OVRSCN){ - trlmessage("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction\n",start,finish); - return(ERROR_RETURN); + if (start >= P_OVRSCN && finish + V_OVRSCNX2 <= (RAZ_COLS/2) - P_OVRSCN) { + trlerror("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction\n", start, finish); + return (status=ERROR_RETURN); } /*SAVE THE PCTETABLE INFORMATION TO THE HEADER OF THE SCIENCE IMAGE AFTER CHECKING TO SEE IF THE USER HAS SPECIFIED ANY CHANGES TO THE CTE CODE VARIABLES. */ - if (CompareCTEParams(&subcd, &cte_pars)) - return (status); + if (CompareCTEParams(&subcd, &cte_pars)) { + return status; + } /*Put the subarray data into full frame*/ Sub2Full(&wf3, &subcd, &cd, 0, 1, 1); @@ -390,22 +393,20 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, /* Subtract the BIAC file from the subarray before continuing The bias routine will take care of cutting out the correct image location for the subarray.*/ - - if (doCteBias(&wf3,&subcd)){ + if (doCteBias(&wf3,&subcd)) { freeSingleGroup(&subcd); - return(status); + return status; } /*reset the array after bias subtraction*/ Sub2Full(&wf3, &subcd, &cd, 0, 1, 1); - } else { /*chip is 1, ab, sci2*/ start=0; finish=0; initSingleGroup(&subab); getSingleGroup(wf3.input, 1, &subab); - if (hstio_err()){ + if (hstio_err()) { freeSingleGroup(&subab); return (status = OPEN_FAILED); } @@ -416,17 +417,19 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, ab.group_num=2; CreateEmptyChip(&wf3, &ab); - if ( GetCorner(&subab.sci.hdr, rsize, sci_bin, sci_corner)) - return (status); + if (GetCorner(&subab.sci.hdr, rsize, sci_bin, sci_corner)) { + return status; + } - if ( GetCorner(&ab.sci.hdr, rsize, ref_bin, ref_corner)) - return (status); + if (GetCorner(&ab.sci.hdr, rsize, ref_bin, ref_corner)) { + return status; + } start = sci_corner[0] - ref_corner[0]; finish = start + subab.sci.data.nx; - if ( start >= P_OVRSCN && finish + V_OVRSCNX2 <= (RAZ_COLS/2) - P_OVRSCN){ - trlmessage("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction\n",start,finish); - return(ERROR_RETURN); + if (start >= P_OVRSCN && finish + V_OVRSCNX2 <= (RAZ_COLS/2) - P_OVRSCN) { + trlerror("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction\n",start,finish); + return (status=ERROR_RETURN); } /*add subarray to full frame image*/ Sub2Full(&wf3, &subab, &ab, 0, 1, 1); @@ -435,8 +438,9 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, AFTER CHECKING TO SEE IF THE USER HAS SPECIFIED ANY CHANGES TO THE CTE CODE VARIABLES. */ - if (CompareCTEParams(&subab, &cte_pars)) - return (status); + if (CompareCTEParams(&subab, &cte_pars)) { + return status; + } /* now create an empty chip 2*/ initSingleGroup(&cd); @@ -449,9 +453,9 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, /* Subtract the BIAC file from the subarray before continuing*/ subab.group_num=2; - if (doCteBias(&wf3,&subab)){ + if (doCteBias(&wf3,&subab)) { freeSingleGroup(&subab); - return(status); + return status; } /*reset the array after bias subtraction*/ @@ -465,19 +469,19 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, initSingleGroup (&cd); getSingleGroup (wf3.input, 1, &cd); - if (hstio_err()){ + if (hstio_err()) { return (status = OPEN_FAILED); } initSingleGroup (&ab); getSingleGroup (wf3.input, 2, &ab); - if (hstio_err()){ + if (hstio_err()) { return (status = OPEN_FAILED); } /*setup the mask*/ - for(i=0; i< ab.dq.data.nx; i++){ - for(j=0; j< ab.dq.data.ny; j++){ + for(j=0; j< ab.dq.data.ny; j++) { + for(i=0; i< ab.dq.data.nx; i++) { PPix(&ab.dq.data, i, j) = 1; PPix(&cd.dq.data, i, j) = 1; } @@ -487,27 +491,39 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, makeRAZ(&cd,&ab,&raw); /***SUBTRACT THE CTE BIAS FROM BOTH CHIPS IN PLACE***/ - if (doCteBias(&wf3,&cd)){ + if (doCteBias(&wf3,&cd)) { freeSingleGroup(&cd); - return(status); + return status; } - if (doCteBias(&wf3,&ab)){ + if (doCteBias(&wf3,&ab)) { freeSingleGroup(&ab); - return(status); + return status; } /*SAVE THE PCTETABLE INFORMATION TO THE HEADER OF THE SCIENCE IMAGE AFTER CHECKING TO SEE IF THE USER HAS SPECIFIED ANY CHANGES TO THE CTE CODE VARIABLES. */ - if (CompareCTEParams(&cd, &cte_pars)) - return (status); + if (CompareCTEParams(&cd, &cte_pars)) { + return status; + } + } + /* CONVERT TO RAZ, SUBTRACT BIAS AND CORRECT FOR GAIN */ + if (raw2raz(&wf3, &cd, &ab, &raz)) { + return status; } - /*CONVERT TO RAZ, SUBTRACT BIAS AND CORRECT FOR GAIN*/ - if (raw2raz(&wf3, &cd, &ab, &raz)) - return (status); + /* Perform XCTE correction but only for fullframe. + This must be done before YCTE correction. + */ + if (!wf3.subarray) { + trlmessage("XCTE: jumping into the routine..."); + if (sub_xctecor(&raz, wf3.expstart)) { + return status; + } + trlmessage("XCTE: returning from the routine..."); + } SingleGroup fff; initSingleGroup(&fff); @@ -518,13 +534,13 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, cte_ff= (wf3.expstart - cte_pars.cte_date0)/ (cte_pars.cte_date1 - cte_pars.cte_date0); - printf("CTE_FF: %8.3f \n",cte_ff); + trlmessage("YCTE_FF: %8.3f\n", cte_ff); cte_pars.scale_frac=cte_ff; - for(i=0;isubarray){ + if (wf3->subarray) { findPreScanBias(raz, bias_pre, bsig_pre); - for (k=0;kdq.data,i+k*subcol,j)){ + for (k=0;kdq.data,i+k*subcol,j)) { Pix(raz->sci.data,i+k*subcol,j) -= bias_pre[k]; Pix(raz->sci.data,i+k*subcol,j) *= gain; } @@ -740,9 +762,9 @@ int raw2raz(WF3Info *wf3, SingleGroup *cd, SingleGroup *ab, SingleGroup *raz){ } } else { findPostScanBias(raz, bias_post, bsig_post); - for (k=0;ksci.data,i+k*subcol,j) -= bias_post[k]; Pix(raz->sci.data,i+k*subcol,j) *= gain; } @@ -750,7 +772,7 @@ int raw2raz(WF3Info *wf3, SingleGroup *cd, SingleGroup *ab, SingleGroup *raz){ } } - return(status); + return status; } /*calculate the post scan and bias after the biac file has been subtracted @@ -767,7 +789,7 @@ int raw2raz(WF3Info *wf3, SingleGroup *cd, SingleGroup *ab, SingleGroup *raz){ these only exist in full frame images */ -int findPostScanBias(SingleGroup *raz, float *mean, float *sigma){ +int findPostScanBias(SingleGroup *raz, float *mean, float *sigma) { extern int status; int arrsize = 55377; @@ -781,23 +803,23 @@ int findPostScanBias(SingleGroup *raz, float *mean, float *sigma){ float sigreg =7.5; /*sigma clip*/ - int subcol = RAZ_COLS/4; + int subcol = RAZ_COLS/NAMPS; int npix=0; /*track array size for resistant mean*/ /*init plist for full size We'll allocate heap memory for smaller arrays */ - for (i=0;idq.data,i+k*subcol,j)) { plist[npix] = Pix(raz->sci.data,i+k*subcol,j); npix+=1; @@ -805,14 +827,14 @@ int findPostScanBias(SingleGroup *raz, float *mean, float *sigma){ } } } - if (npix > 0 ){ + if (npix > 0 ) { plistSub = (float *) calloc(npix, sizeof(float)); - if (plistSub == NULL){ + if (plistSub == NULL) { trlerror("out of memory for resistmean entrance in findPostScanBias."); free(plistSub); - return (ERROR_RETURN); + return (status=OUT_OF_MEMORY); } - for(i=0; idq.data,i+(k*subcol),j)){ + for (i=5;idq.data,i+(k*subcol),j)) { plist[npix] = Pix(raz->sci.data,i+k*subcol,j); npix+=1; } @@ -871,24 +893,25 @@ int findPreScanBias(SingleGroup *raz, float *mean, float *sigma){ } } - if (0 < npix ){ + if (0 < npix ) { plistSub = (float *) calloc(npix, sizeof(float)); - if (plistSub == NULL){ + if (plistSub == NULL) { trlerror("out of memory for resistmean entrance in findPreScanBias."); free(plistSub); - return (ERROR_RETURN); + return (status=OUT_OF_MEMORY); } - for(i=0; i0) - printf("npix=%i\nmean[%i]=%f\nsigma[%i] = %f\n",npix,k+1,rmean,k+1,rsigma); + if(npix>0) { + trlmessage("npix=%i\nmean[%i]=%f\nsigma[%i] = %f\n",npix,k+1,rmean,k+1,rsigma); + } } return status; } @@ -919,8 +942,7 @@ int findPreScanBias(SingleGroup *raz, float *mean, float *sigma){ the ttrap reference to the image array has to be -1 for C */ - -int sim_colreadout_l(double *pixi, double *pixo, double *pixf, CTEParams *cte){ +int sim_colreadout_l(double *pixi, double *pixo, double *pixf, CTEParams *cte) { extern int status; int j; @@ -955,7 +977,7 @@ int sim_colreadout_l(double *pixi, double *pixo, double *pixf, CTEParams *cte){ /*FIGURE OUT WHICH TRAPS WE DON'T NEED TO WORRY ABOUT IN THIS COLUMN PMAX SHOULD ALWAYS BE POSITIVE HERE */ pmax=10.; - for(j=0; j pmax) pmax=pixo[j]; @@ -963,7 +985,7 @@ int sim_colreadout_l(double *pixi, double *pixo, double *pixf, CTEParams *cte){ /*GO THROUGH THE TRAPS ONE AT A TIME, FROM HIGHEST TO LOWEST Q, AND SEE WHEN THEY GET FILLED AND EMPTIED, ADJUST THE PIXELS ACCORDINGLY*/ - for (w = cte->cte_traps-1; w>=0; w--){ + for (w = cte->cte_traps-1; w>=0; w--) { if ( cte->qlevq_data[w] <= pmax ) { ftrap = 0.0e0; @@ -971,12 +993,12 @@ int sim_colreadout_l(double *pixi, double *pixo, double *pixf, CTEParams *cte){ fcarry = 0.0e0; /*GO UP THE COLUMN PIXEL BY PIXEL*/ - for(j=0; jcte_len) || ( pix_1 >= cte->qlevq_data[w] - 1. ) ){ - if (pixo[j] >= 0 ){ + if ( (ttrap < cte->cte_len) || ( pix_1 >= cte->qlevq_data[w] - 1. ) ) { + if (pixo[j] >= 0 ) { pix_1 = pixo[j] + fcarry; /*shuffle charge in*/ fcarry = pix_1 - floor(pix_1); /*carry the charge remainder*/ pix_1 = floor(pix_1); /*reset pixel*/ @@ -991,14 +1013,14 @@ int sim_colreadout_l(double *pixi, double *pixo, double *pixf, CTEParams *cte){ /*RELEASE THE CHARGE*/ padd_2=0.0; - if (ttrap cte_len){ + if (ttrap cte_len) { ttrap += 1; padd_2 = Pix(rprof->data,w,ttrap-1) *ftrap; } padd_3 = 0.0; prem_3 = 0.0; - if ( pix_1 >= cte->qlevq_data[w]){ + if ( pix_1 >= cte->qlevq_data[w]) { prem_3 = cte->dpdew_data[w] / cte->n_par * pixf[j]; /*dpdew is 1 in file */ if (ttrap < cte->cte_len) padd_3 = Pix(cprof->data,w,ttrap-1)*ftrap; @@ -1013,7 +1035,7 @@ int sim_colreadout_l(double *pixi, double *pixo, double *pixf, CTEParams *cte){ }/*end for w*/ - return(status); + return status; } @@ -1042,7 +1064,6 @@ int initCTETrl (char *input, char *output) { int nsuffix = 1; - /* Start by stripping off suffix from input/output filenames */ if (MkOutName (input, isuffix, trlsuffix, nsuffix, trl_in, CHAR_LINE_LENGTH)) { WhichError (status); @@ -1054,14 +1075,13 @@ int initCTETrl (char *input, char *output) { } /* NOW, CONVERT TRAILER FILENAME EXTENSIONS FROM '.FITS' TO '.TRL' */ - if (MkNewExtn (trl_in, TRL_EXTN) ) { - trlerror("Error with input trailer filename %s", trl_in); WhichError (status); + trlerror("Error with input trailer filename %s", trl_in); } if (MkNewExtn (trl_out, TRL_EXTN) ) { - trlerror("Error with output trailer filename %s", trl_out); WhichError (status); + trlerror("Error with output trailer filename %s", trl_out); } /* If we are working with a RAW file, then see if a TRL file @@ -1080,7 +1100,7 @@ int initCTETrl (char *input, char *output) { ** trailer file into it. */ InitTrlFile (trl_in, trl_out); - return(status); + return status; } @@ -1155,48 +1175,50 @@ int rm_rnZ_colj(double *pixj_chg, /* Bounds of this loop include only the non-zero portion of the data */ for(j=j1; j<=j2; j++) { - /* Compare each pixel to the average of its up/down neighbors */ - dd = pixj_rsz[j]-(pixj_rsz[j+1]+pixj_rsz[j-1])/2.0; - pixj_ffu[j] = 0.0; - - /* If the pixel is within the readnoise window... */ - if (fabs(dd) < RNN) { - /* Cap the adjustments we are willing to make at any time to 20% */ - if (dd > RNMIT/5.0) dd = RNMIT/5.0; - if (dd < -RNMIT/5.0) dd = -RNMIT/5.0; - /* Take half of the maximum adjustment from the current pixel... */ - pixj_rnz[j ] = pixj_rnz[j ] + dd*0.50; - /* ...and give half of the value to the pixel below and above */ - pixj_rnz[j-1] = pixj_rnz[j-1] - dd*0.25; - pixj_rnz[j+1] = pixj_rnz[j+1] - dd*0.25; - /* Flag this pixel to be in the readnoise range so we can track - * the total noise in the nominal pixels - */ - pixj_ffu[j ] = 1.0; - } - } - /* Bounds of this loop include only the non-zero portion of the data */ - for(j=j1; j<=j2; j++) { - pixj_rsz[j] = pixj_chg[j] - pixj_rnz[j]; - } - dtot = 0.; - ntot = 0.; - ftot = 0.; - rtot = 0.; - /* Bounds of this loop include only the non-zero portion of the data */ - for(j=j1; j<=j2; j++) { - ftot = ftot + pixj_ffu[j]; - dtot = dtot + pixj_rnz[j]*pixj_rnz[j]; - dd = pixj_rsz[j]-(pixj_rsz[j+1]+pixj_rsz[j-1])/2.0; - rtot = rtot + dd*dd; - ntot = ntot + 1; - } - RNU = sqrt(dtot/ftot); - - if (RNU > RNMIT) return(0); + /* Compare each pixel to the average of its up/down neighbors */ + dd = pixj_rsz[j]-(pixj_rsz[j+1]+pixj_rsz[j-1])/2.0; + pixj_ffu[j] = 0.0; + + /* If the pixel is within the readnoise window... */ + if (fabs(dd) < RNN) { + /* Cap the adjustments we are willing to make at any time to 20% */ + if (dd > RNMIT/5.0) dd = RNMIT/5.0; + if (dd < -RNMIT/5.0) dd = -RNMIT/5.0; + /* Take half of the maximum adjustment from the current pixel... */ + pixj_rnz[j ] = pixj_rnz[j ] + dd*0.50; + /* ...and give half of the value to the pixel below and above */ + pixj_rnz[j-1] = pixj_rnz[j-1] - dd*0.25; + pixj_rnz[j+1] = pixj_rnz[j+1] - dd*0.25; + /* Flag this pixel to be in the readnoise range so we can track + * the total noise in the nominal pixels + */ + pixj_ffu[j ] = 1.0; + } + } + /* Bounds of this loop include only the non-zero portion of the data */ + for(j=j1; j<=j2; j++) { + pixj_rsz[j] = pixj_chg[j] - pixj_rnz[j]; + } + dtot = 0.; + ntot = 0.; + ftot = 0.; + rtot = 0.; + /* Bounds of this loop include only the non-zero portion of the data */ + for(j=j1; j<=j2; j++) { + ftot = ftot + pixj_ffu[j]; + dtot = dtot + pixj_rnz[j]*pixj_rnz[j]; + dd = pixj_rsz[j]-(pixj_rsz[j+1]+pixj_rsz[j-1])/2.0; + rtot = rtot + dd*dd; + ntot = ntot + 1; + } + RNU = sqrt(dtot/ftot); + + if (RNU > RNMIT) { + return WF3_OK; + } } - return(0); + return WF3_OK; } @@ -1243,8 +1265,8 @@ int sim_colreadout_l_uvis_w(double *pixi, // input column array (JDIM) /* Bounds checking */ if (Ws>WsMAX) { - printf("Ws error\n"); - return(1); + trlerror("Ws error: %d > %d", Ws, WsMAX); + return ERROR_RETURN; } /* Figure out which traps we do not need @@ -1341,7 +1363,7 @@ int sim_colreadout_l_uvis_w(double *pixi, // input column array (JDIM) pixo[j] = pixo[j] + prel_1 + prel_2 - pgrb_3; } } - return(0); + return WF3_OK; } @@ -1396,13 +1418,13 @@ int sub_ctecor_v2c(float *pixz_raz, NITFORs = PCTENFOR; NITPARs = PCTENPAR; - printf(" \n"); - printf(" INSIDE sub_ctecor_v2.f... \n"); - printf(" ---> PCTERNOI: %8.4f \n",PCTERNOI); - printf(" ---> FIX_ROCR: %8.4f \n",FIX_ROCR); - printf(" ---> NITFORs: %5d \n",NITFORs); - printf(" ---> NITPARs: %5d \n",NITPARs); - printf(" \n"); + trlmessage(" \n"); + trlmessage(" INSIDE sub_ctecor_v2 ... \n"); + trlmessage(" ---> PCTERNOI: %8.4f \n",PCTERNOI); + trlmessage(" ---> FIX_ROCR: %8.4f \n",FIX_ROCR); + trlmessage(" ---> NITFORs: %5d \n",NITFORs); + trlmessage(" ---> NITPARs: %5d \n",NITPARs); + trlmessage(" \n"); #pragma omp parallel \ shared(pixz_raz,pixz_fff,pixz_rzc, \ @@ -1415,9 +1437,7 @@ int sub_ctecor_v2c(float *pixz_raz, #pragma omp for - for(i=0;i HIGH_CLIP) RNOIVALo = HIGH_CLIP; - return(RNOIVALo); + return RNOIVALo; } From 7dd8c1a0ef2a4c4b1548cb3fa466990f51737427 Mon Sep 17 00:00:00 2001 From: Pey Lian Lim <2090236+pllim@users.noreply.github.com> Date: Thu, 12 Mar 2026 17:34:30 -0400 Subject: [PATCH 2/5] Clean up TRA outputs --- pkg/wfc3/lib/wf3cte/cte_dobias.c | 2 +- pkg/wfc3/lib/wf3cte/cte_serial_corr.c | 19 +++-- pkg/wfc3/lib/wf3cte/getctepars.c | 103 +++++++++++++------------- pkg/wfc3/lib/wf3cte/wf3cte.c | 72 +++++++++--------- 4 files changed, 99 insertions(+), 97 deletions(-) diff --git a/pkg/wfc3/lib/wf3cte/cte_dobias.c b/pkg/wfc3/lib/wf3cte/cte_dobias.c index dfac7f146..f1d621d54 100644 --- a/pkg/wfc3/lib/wf3cte/cte_dobias.c +++ b/pkg/wfc3/lib/wf3cte/cte_dobias.c @@ -92,7 +92,7 @@ int doCteBias (WF3Info *wf3, SingleGroup *x) { */ if (wf3->verbose){ - trlmessage("Image has starting location of %d,%d in the reference image",x0,y0); + trlmessage(" Image has starting location of %d,%d in the reference image", x0, y0); } /* Subtract the bias image from x. */ diff --git a/pkg/wfc3/lib/wf3cte/cte_serial_corr.c b/pkg/wfc3/lib/wf3cte/cte_serial_corr.c index 70ed3a965..943dbe86b 100644 --- a/pkg/wfc3/lib/wf3cte/cte_serial_corr.c +++ b/pkg/wfc3/lib/wf3cte/cte_serial_corr.c @@ -44,7 +44,7 @@ int sub_xctecor(SingleGroup *raz, double mjd_obs) { /* These values are independent of PCTETAB. */ float ff_xcte = (mjd_obs - 54962.00) / (60125.0 - 54962.00); - trlmessage("XCTE_FF: %8.3f\n", ff_xcte); + trlmessage("XCTE_FF: %8.3f", ff_xcte); /* Make a copy of raz data and use as input. */ for(j=0; jcte_date1=0.0f; pars->cte_traps=0.0f; pars->cte_len=0; - pars->n_forward=0; + pars->n_forward=0; pars->n_par=0; pars->scale_frac=0.0f; /*will be updated during routine run*/ - pars->noise_mit=0; + pars->noise_mit=0; pars->thresh=0.0f; pars->descrip2[0]='\0'; /*static scheduling is faster when there are no dependent loop variables - for dependent variables inside the main loop switch to dynamic so that - the thread wait time can be variable + for dependent variables inside the main loop switch to dynamic so that + the thread wait time can be variable */ for (i=0; iwcol_data[i]=0; + pars->wcol_data[i]=0; pars->qlevq_data[i]=0.0f; pars->dpdew_data[i]=0.0f; } @@ -85,17 +85,17 @@ int GetCTEPars (char *filename, CTEParams *pars) { Filename: wfc3_cte.fits No. Name Type Cards Dimensions Format -0 PRIMARY PrimaryHDU 21 () -1 QPROF BinTableHDU 16 999R x 3C ['i', 'i', 'i'] -2 SCLBYCOL BinTableHDU 20 8412R x 5C ['i', 'e', 'e', 'e', 'e'] -3 RPROF ImageHDU 12 (999, 100) float32 -4 CPROF ImageHDU 12 (999, 100) float32 +0 PRIMARY PrimaryHDU 21 () +1 QPROF BinTableHDU 16 999R x 3C ['i', 'i', 'i'] +2 SCLBYCOL BinTableHDU 20 8412R x 5C ['i', 'e', 'e', 'e', 'e'] +3 RPROF ImageHDU 12 (999, 100) float32 +4 CPROF ImageHDU 12 (999, 100) float32 */ extern int status; /* variable for return status */ int ctraps; /*see if more traps were added to reference file*/ - + /* HSTIO VARIABLES */ Hdr hdr_ptr; @@ -218,12 +218,12 @@ No. Name Type Cards Dimensions Format status = KEYWORD_MISSING; return status; } - + /* trlmessage("FIXROCR: %d",pars->fix_rocr); */ - - + + /* DONE READING STUFF FROM THE PRIMARY HEADER */ freeHdr(&hdr_ptr); @@ -269,9 +269,9 @@ No. Name Type Cards Dimensions Format /* LOOP OVER TABLE ROWS UP TO SIZE TRAPS*/ ctraps=0; /*actual usable traps*/ - + for (j = 0; j < TRAPS; j++) { - + /* GET W FROM THIS ROW */ c_tbegti(tbl_ptr, w_ptr, j+1, &pars->wcol_data[j]); if (c_iraferr()) { @@ -279,7 +279,7 @@ No. Name Type Cards Dimensions Format status = TABLE_ERROR; return status; } - + /* GET QLEVQ FROM THIS ROW */ c_tbegtd(tbl_ptr, qlevq_ptr, j+1, &pars->qlevq_data[j]); if (c_iraferr()) { @@ -287,10 +287,10 @@ No. Name Type Cards Dimensions Format status = TABLE_ERROR; return status; } - + if (pars->qlevq_data[j] < 999999.) ctraps+=1; - + /* GET DPDEW FROM THIS ROW */ c_tbegtd(tbl_ptr, dpdew_ptr, j+1, &pars->dpdew_data[j]); if (c_iraferr()) { @@ -302,7 +302,7 @@ No. Name Type Cards Dimensions Format trlmessage("More TRAPS in reference file than available, update TRAPS: %i -> %i",TRAPS,(int)ctraps); } } - + /*IF CTRAPS EVER OVERFLOWS INT THIS NEEDS TO BE CHANGED*/ pars->cte_traps=(int)ctraps; @@ -310,7 +310,7 @@ No. Name Type Cards Dimensions Format trlmessage("(pctecorr) data check for PCTETAB QPROF, row %i, %i\t%g\t%g\ttraps=%i\n",20, pars->wcol_data[19],pars->qlevq_data[19], pars->dpdew_data[19], pars->cte_traps); */ - + /* CLOSE CTE PARAMETERS FILE FOR EXTENSION 1*/ c_tbtclo(tbl_ptr); @@ -370,7 +370,7 @@ No. Name Type Cards Dimensions Format /* loop over table rows */ for (j = 0; j < RAZ_COLS; j++) { - /* get trap from this row */ + /* get trap from this row */ c_tbegti(tbl_ptr, iz_ptr, j+1, &pars->iz_data[j]); if (c_iraferr()) { @@ -405,40 +405,40 @@ No. Name Type Cards Dimensions Format trlmessage("(pctecorr) data check for PCTETAB SCLBYCOL row %d, %d %g\t%g\t%g\t%g\ntotal traps = %i", j,pars->iz_data[j-1],pars->scale512[j-1],pars->scale1024[j-1],pars->scale1536[j-1],pars->scale2048[j-1],pars->cte_traps); */ - + /* close CTE parameters file for extension 2*/ c_tbtclo(tbl_ptr); /****************************************************************************/ /* extension 3: differential trail profile as image */ - trlmessage("Reading in image from extension 3"); + trlmessage("PCTETAB Reading in image from extension 3"); /* Get the coefficient images from the PCTETAB */ pars->rprof = (FloatHdrData *)calloc(1,sizeof(FloatHdrData)); if (pars->rprof == NULL){ trlerror("Can't allocate memory for RPROF ref data"); - return (status = 1); + return (status = OUT_OF_MEMORY); } initFloatHdrData(pars->rprof); if (getFloatHD (filename, "RPROF", 1, pars->rprof)){ - return (status=1); + return (status = ERROR_RETURN); } /****************************************************************************/ /* ext number 4 : cummulative trail profile as image */ - trlerror("Reading in image from extension 4"); + trlmessage("PCTETAB Reading in image from extension 4"); pars->cprof = (FloatHdrData *)calloc(1,sizeof(FloatHdrData)); if (pars->cprof == NULL){ trlerror("Can't allocate memory for CPROF ref data"); - return (status = 1); + return (status = OUT_OF_MEMORY); } /* Get the coefficient images from the PCTETAB */ initFloatHdrData (pars->cprof); if (getFloatHD (filename, "CPROF", 1, pars->cprof)){ - return (status=1); + return (status = ERROR_RETURN); } return(status); @@ -450,7 +450,7 @@ No. Name Type Cards Dimensions Format are then the values found there will be used instead of the values read from the PCTETAB. If the values are not set they will be populated with the values from the PCTETAB, namely these: - + 'CTE_NAME':'pixelCTE 2012', #name of cte algorithm 'CTE_VER':'1.0' , #version number of algorithm 'CTEDATE0':54962.0, #date of uvis installation in HST in MJD @@ -462,8 +462,8 @@ No. Name Type Cards Dimensions Format 'PCTETRSH':-10.0 ,#over subtraction threshold, always use reference value 'FIXROCR' : 1, #set to 1 for true, fix the readout cr's ***NOTE: The PCTERNOI value is no longer used from the PCTETAB. If the PCTERNOI - keyword value in the raw science image header is non-zero, it will be used - for the CTE computations. Otherwise, the value is computed on-the-fly based + keyword value in the raw science image header is non-zero, it will be used + for the CTE computations. Otherwise, the value is computed on-the-fly based upon the raw image data. (March 2020) */ @@ -476,25 +476,25 @@ int CompareCTEParams(SingleGroup *group, CTEParams *pars) { int n_par; int fix_rocr; int noise_mit; - + cte_len=0; n_forward=0; n_par=0; fix_rocr=0; noise_mit=0; - + /*always put the cte_name and cte_ver from the reference file in*/ if (PutKeyStr(group->globalhdr,"CTE_NAME", pars->cte_name, "CTE algorithm name")){ trlmessage("(pctecorr) Error updating CTE_NAME keyword in image header"); return (status=HEADER_PROBLEM); } - + if(PutKeyStr(group->globalhdr,"CTE_VER", pars->cte_ver, "CTE algorithm version")){ trlmessage("(pctecorr) Error updating CTE_VER keyword in image header"); return (status=HEADER_PROBLEM); } - + if (PutKeyDbl(group->globalhdr, "CTEDATE0", pars->cte_date0,"Date of UVIS installation")) { trlmessage("(pctecorr) Error putting CTEDATE0 keyword in header"); return (status=HEADER_PROBLEM); @@ -510,7 +510,7 @@ int CompareCTEParams(SingleGroup *group, CTEParams *pars) { trlmessage("(pctecorr) Error putting CTEDATE1 keyword in header"); return (status=HEADER_PROBLEM); } - + /*check the PCTENSMD keyword in the header*/ if (GetKeyInt(group->globalhdr, "PCTENSMD", NO_DEFAULT, -999, &noise_mit)){ trlmessage("(pctecorr) Error reading PCTENSMD keyword from header"); @@ -518,8 +518,8 @@ int CompareCTEParams(SingleGroup *group, CTEParams *pars) { } if (noise_mit != pars->noise_mit){ pars->noise_mit = noise_mit; - } - + } + /*check the PCTEDIM keyword in header*/ if (GetKeyInt(group->globalhdr, "PCTETLEN", NO_DEFAULT, -999, &cte_len)) { trlmessage("(pctecorr) Error reading PCTETLEN keyword from header"); @@ -534,13 +534,13 @@ int CompareCTEParams(SingleGroup *group, CTEParams *pars) { return (status=HEADER_PROBLEM); } } - + /* get number of iterations used in forward model */ if (GetKeyInt(group->globalhdr, "PCTENFOR", NO_DEFAULT, -999, &n_forward)) { trlmessage("(pctecorr) Error reading PCTENFOR keyword from header"); return (status=HEADER_PROBLEM); } - + if (n_forward > 1 && n_forward != pars->n_forward){ pars->n_forward = n_forward; } else { @@ -549,15 +549,15 @@ int CompareCTEParams(SingleGroup *group, CTEParams *pars) { return (status=HEADER_PROBLEM); } } - + /* get number of iterations used in parallel transfer*/ if (GetKeyInt(group->globalhdr, "PCTENPAR", NO_DEFAULT, -999, &n_par)) { trlmessage("(pctecorr) Error reading PCTENPAR keyword from header"); return (status=HEADER_PROBLEM); } - - + + if( n_par >1 && n_par != pars->n_par){ pars->n_par = n_par; } else { @@ -566,8 +566,8 @@ int CompareCTEParams(SingleGroup *group, CTEParams *pars) { return (status=HEADER_PROBLEM); } } - - + + /*fix the readout Cr's? */ if (GetKeyInt(group->globalhdr, "FIXROCR", NO_DEFAULT, -999, &fix_rocr)){ trlmessage("(pctecorr) Error reading FIXROCR keyword from header"); @@ -583,4 +583,3 @@ int CompareCTEParams(SingleGroup *group, CTEParams *pars) { } return status; } - diff --git a/pkg/wfc3/lib/wf3cte/wf3cte.c b/pkg/wfc3/lib/wf3cte/wf3cte.c index d6eb89397..b66886db1 100644 --- a/pkg/wfc3/lib/wf3cte/wf3cte.c +++ b/pkg/wfc3/lib/wf3cte/wf3cte.c @@ -366,7 +366,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, start = sci_corner[0] - ref_corner[0]; finish = start + subcd.sci.data.nx; if (start >= P_OVRSCN && finish + V_OVRSCNX2 <= (RAZ_COLS/2) - P_OVRSCN) { - trlerror("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction\n", start, finish); + trlerror("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction", start, finish); return (status=ERROR_RETURN); } @@ -428,7 +428,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, start = sci_corner[0] - ref_corner[0]; finish = start + subab.sci.data.nx; if (start >= P_OVRSCN && finish + V_OVRSCNX2 <= (RAZ_COLS/2) - P_OVRSCN) { - trlerror("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction\n",start,finish); + trlerror("Subarray not taken with physical overscan (%i %i)\nCan't perform CTE correction",start,finish); return (status=ERROR_RETURN); } /*add subarray to full frame image*/ @@ -517,12 +517,14 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, /* Perform XCTE correction but only for fullframe. This must be done before YCTE correction. */ - if (!wf3.subarray) { - trlmessage("XCTE: jumping into the routine..."); + if (wf3.subarray) { + trlmessage("XCTE: Skipped for subarray"); + } else { + trlmessage("XCTE: Jumping into the routine..."); if (sub_xctecor(&raz, wf3.expstart)) { return status; } - trlmessage("XCTE: returning from the routine..."); + trlmessage("XCTE: Returning from the routine..."); } SingleGroup fff; @@ -534,7 +536,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, cte_ff= (wf3.expstart - cte_pars.cte_date0)/ (cte_pars.cte_date1 - cte_pars.cte_date0); - trlmessage("YCTE_FF: %8.3f\n", cte_ff); + trlmessage("YCTE_FF: %8.3f", cte_ff); cte_pars.scale_frac=cte_ff; @@ -552,13 +554,13 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, float FLOAT_RNOIVAL = 0.; float FLOAT_BKGDVAL = 0.; readNoise = wf3.pcternoi; - trlmessage("PCTERNOI: %8.4f (source: primary header of science image)\n\n", readNoise); + trlmessage("PCTERNOI: %8.4f (source: primary header of science image)", readNoise); /* Comparison should be OK - read from FITS header and no computation */ if (readNoise == 0.0) { readNoise = find_raz2rnoival(raz.sci.data.data, &FLOAT_RNOIVAL, &FLOAT_BKGDVAL); - trlmessage("RNOIVAL: %8.4f BKGDVAL: %8.4f\n", FLOAT_RNOIVAL, FLOAT_BKGDVAL); + trlmessage("RNOIVAL: %8.4f BKGDVAL: %8.4f", FLOAT_RNOIVAL, FLOAT_BKGDVAL); trlmessage("PCTERNOI: %8.4f (source: computed on-the-fly from science image)", readNoise); - trlmessage("This computed value supersedes any value obtained from the primary\nheader of the science image.\n\n"); + trlmessage(" This computed value supersedes any value obtained from the primary\n header of the science image."); } /* The PCTERNOI value actually used is written to the PCTERNOI keyword in @@ -568,7 +570,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, /* Invoke the updated CTE correction which does the read noise mitigation in each of the three forward-model iterations. */ - trlmessage("YCTE: jumping into the routine..."); + trlmessage("YCTE: Jumping into the routine..."); if (sub_ctecor_v2c(raz.sci.data.data, fff.sci.data.data, WsMAX, @@ -583,7 +585,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, rzc.sci.data.data)) { return status; } - trlmessage("YCTE: returning from the routine..."); + trlmessage("YCTE: Returning from the routine..."); for (j=0;j0) { - trlmessage("npix=%i\nmean[%i]=%f\nsigma[%i] = %f\n",npix,k+1,rmean,k+1,rsigma); + trlmessage("npix=%i\nmean[%i]=%f\nsigma[%i] = %f", npix, k+1, rmean, k+1, rsigma); } } return status; @@ -1418,13 +1420,11 @@ int sub_ctecor_v2c(float *pixz_raz, NITFORs = PCTENFOR; NITPARs = PCTENPAR; - trlmessage(" \n"); - trlmessage(" INSIDE sub_ctecor_v2 ... \n"); - trlmessage(" ---> PCTERNOI: %8.4f \n",PCTERNOI); - trlmessage(" ---> FIX_ROCR: %8.4f \n",FIX_ROCR); - trlmessage(" ---> NITFORs: %5d \n",NITFORs); - trlmessage(" ---> NITPARs: %5d \n",NITPARs); - trlmessage(" \n"); + trlmessage(" INSIDE sub_ctecor_v2 ..."); + trlmessage(" ---> PCTERNOI: %8.4f", PCTERNOI); + trlmessage(" ---> FIX_ROCR: %8.4f", FIX_ROCR); + trlmessage(" ---> NITFORs: %5d", NITFORs); + trlmessage(" ---> NITPARs: %5d", NITPARs); #pragma omp parallel \ shared(pixz_raz,pixz_fff,pixz_rzc, \ @@ -1520,8 +1520,11 @@ int sub_ctecor_v2c(float *pixz_raz, /* This variable exists for debuggin purposes. */ NDONE++; - /*if (NDONE==(NDONE/100)*100) { trlmessage(" i = %5d %5d %5d \n",i,NCRX,NDONE); }*/ - + /* + if (NDONE % 100 == 0) { + trlmessage(" i = %5d %5d %5d", i, NCRX, NDONE); + } + */ } return status; @@ -1688,13 +1691,10 @@ float find_raz2rnoival(float *raz_cdab, float *FLOAT_RNOIVAL, float *FLOAT_BKGDV } /* For debugging purposes only - trlmessage(" \n"); - trlmessage(" vsum: %12lld \n",vsum); - trlmessage(" nsum: %12lld \n",nsum); - trlmessage(" \n"); - trlmessage(" dbar: %12.2f \n",idmin/2.30*(SPREAD_FOR_HISTO)); - trlmessage(" vbar: %12.2f %12lld %12lld \n",vsum/nsum*(SPREAD_FOR_HISTO),vsum,nsum); - trlmessage(" \n"); + trlmessage(" vsum: %12lld", vsum); + trlmessage(" nsum: %12lld\n", nsum); + trlmessage(" dbar: %12.2f", idmin/2.30*(SPREAD_FOR_HISTO)); + trlmessage(" vbar: %12.2f %12lld %12lld\n", vsum/nsum*(SPREAD_FOR_HISTO), vsum, nsum); */ RNOIVAL = (int)(idmin/2.30/(SPREAD_FOR_HISTO)/sqrt(1+1/8.0)*4+0.5)/4.00; From e0635b1cb08829b92c39ad265c08309127bd6bb1 Mon Sep 17 00:00:00 2001 From: Pey Lian Lim <2090236+pllim@users.noreply.github.com> Date: Thu, 26 Mar 2026 18:08:02 -0400 Subject: [PATCH 3/5] TMP: Write out RZX for INS to compare --- pkg/wfc3/include/cte.h | 1 + pkg/wfc3/lib/razutil.c | 89 +++++++++++++++++++++++++++++++++++- pkg/wfc3/lib/wf3cte/wf3cte.c | 8 ++++ 3 files changed, 96 insertions(+), 2 deletions(-) diff --git a/pkg/wfc3/include/cte.h b/pkg/wfc3/include/cte.h index 1a4e31dee..5c3b68d0b 100644 --- a/pkg/wfc3/include/cte.h +++ b/pkg/wfc3/include/cte.h @@ -98,5 +98,6 @@ int sub_xctecor(SingleGroup *, double); int makeRAZ(SingleGroup *, SingleGroup *, SingleGroup *); int undoRAZ(SingleGroup *, SingleGroup *, SingleGroup *); +int writfits_r4(char *, float *, int, int); #endif /* INCL_CTE_H */ diff --git a/pkg/wfc3/lib/razutil.c b/pkg/wfc3/lib/razutil.c index b3cb32f99..b1f27fe3e 100644 --- a/pkg/wfc3/lib/razutil.c +++ b/pkg/wfc3/lib/razutil.c @@ -6,7 +6,6 @@ /* These routines facilitate moving between the regular WFC3 image structre to the RAZ image structure used in the CTE correction and Sink pixel flagging - convert a raw file to raz file: CDAB longwise amps, save data array for comparison with what jay has during testing @@ -23,7 +22,6 @@ */ - /*convert the sci and dq extensions to the long format*/ int makeRAZ(SingleGroup *cd, SingleGroup *ab, SingleGroup *raz){ @@ -216,3 +214,90 @@ int makeFloatRaz(FloatTwoDArray *x, FloatTwoDArray *raz, int group){ return(status); } + + +/* This is Jay's quick-and-easy FITS writer that can be used to + write out temporary RAZ images for internal testing. + Do not use for production in pipeline. +*/ +int writfits_r4(char *filename, float *data, int nx, int ny) { + + typedef unsigned char Byte; + + char cbuff[2880]; + Byte bbuff[2880]; + Byte bbufr[2880]; + float rbuff[ 720]; + + int i; + int n; + int NRs, NRo; + int nb, NBs; + + Byte btemp4[4]; + Byte *pByte; + + for (i=0;i<2880;i++) { + cbuff[i] = ' '; + } + + /* simplest possible header */ + sprintf(&cbuff[0+ 0*80],"SIMPLE = T"); + sprintf(&cbuff[0+ 1*80],"BITPIX = -32"); + sprintf(&cbuff[0+ 2*80],"NAXIS = 2"); + sprintf(&cbuff[0+ 3*80],"NAXIS1 = %20d",nx); + sprintf(&cbuff[0+ 4*80],"NAXIS2 = %20d",ny); + sprintf(&cbuff[0+ 5*80],"DATATYPE= 'REAL*4' "); + sprintf(&cbuff[0+ 6*80],"COMMENT "); + sprintf(&cbuff[0+ 7*80],"COMMENT Jay's Quickie "); + sprintf(&cbuff[0+ 8*80],"COMMENT C fits writer "); + sprintf(&cbuff[0+ 9*80],"COMMENT writfits_r4() "); + sprintf(&cbuff[0+10*80],"COMMENT "); + sprintf(&cbuff[0+11*80],"END "); + + for (i=0;i<2880;i++) { + if (cbuff[i]==0) { + cbuff[i] = ' '; + } + } + + FILE *fp = fopen(filename,"w"); + fwrite(cbuff,2880,1,fp); + + NRs = nx*ny; + NBs = 1 + ((nx*ny-1)/720); /* the -1 covers the case where the data have an exact multiple + of 2880 bytes; it can happen! */ + + /* distill the data into 2880-byte buffers + we must flip the bytes, since INTEL machines + have the opposite Endian-ness to the FITS + standard + */ + for (nb=0;nb Date: Thu, 4 Jun 2026 17:57:03 -0400 Subject: [PATCH 4/5] Undo row-major loop order --- pkg/wfc3/lib/wf3cte/wf3cte.c | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pkg/wfc3/lib/wf3cte/wf3cte.c b/pkg/wfc3/lib/wf3cte/wf3cte.c index b8dfbceea..2af512fa8 100644 --- a/pkg/wfc3/lib/wf3cte/wf3cte.c +++ b/pkg/wfc3/lib/wf3cte/wf3cte.c @@ -548,8 +548,8 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, cte_pars.scale_frac=cte_ff; - for(j=0;j Date: Fri, 5 Jun 2026 16:27:31 -0400 Subject: [PATCH 5/5] Fix bug, clean more code --- pkg/wfc3/lib/wf3cte/cte_dobias.c | 4 +- pkg/wfc3/lib/wf3cte/cte_serial_corr.c | 23 ++--- pkg/wfc3/lib/wf3cte/wf3cte.c | 116 ++++++++++++++------------ 3 files changed, 72 insertions(+), 71 deletions(-) diff --git a/pkg/wfc3/lib/wf3cte/cte_dobias.c b/pkg/wfc3/lib/wf3cte/cte_dobias.c index f1d621d54..efa86e1d5 100644 --- a/pkg/wfc3/lib/wf3cte/cte_dobias.c +++ b/pkg/wfc3/lib/wf3cte/cte_dobias.c @@ -39,21 +39,19 @@ int doCteBias (WF3Info *wf3, SingleGroup *x) { int scilines; /* number of lines in science image */ int i, j; int update; - int dimx, dimy; + int dimx; int straddle=0; int overstart=0; int offsetx=0; int offsety=0; dimx = x->sci.data.nx; - dimy = x->sci.data.ny; trlmessage("CTE: Subtracting BIACFILE: %s for imset %d",wf3->biac.name, x->group_num); initSingleGroupLine (&y); scilines = x->sci.data.ny; - /* Initialize local variables */ rx = 1; ry = 1; diff --git a/pkg/wfc3/lib/wf3cte/cte_serial_corr.c b/pkg/wfc3/lib/wf3cte/cte_serial_corr.c index 943dbe86b..11f85d55e 100644 --- a/pkg/wfc3/lib/wf3cte/cte_serial_corr.c +++ b/pkg/wfc3/lib/wf3cte/cte_serial_corr.c @@ -37,9 +37,10 @@ static void init_ecr_t(float [NUM_TRAP], float [NUM_TRAP], float [NUM_TRAP]); */ int sub_xctecor(SingleGroup *raz, double mjd_obs) { extern int status; - int i, j, k; - float *pixz_i = malloc(RAZ_COLS * RAZ_ROWS * sizeof(float)); - float *pixz_o = malloc(RAZ_COLS * RAZ_ROWS * sizeof(float)); + int k; + const int raz_totsize = RAZ_COLS * RAZ_ROWS; + float *pixz_i = malloc(raz_totsize * sizeof(float)); + float *pixz_o = malloc(raz_totsize * sizeof(float)); /* These values are independent of PCTETAB. */ float ff_xcte = (mjd_obs - 54962.00) / (60125.0 - 54962.00); @@ -47,11 +48,8 @@ int sub_xctecor(SingleGroup *raz, double mjd_obs) { trlmessage("XCTE_FF: %8.3f", ff_xcte); /* Make a copy of raz data and use as input. */ - for(j=0; jsci.data, i, j); - } + for(k=0; ksci.data.data[k]; } if (sub_raz2rzx_wfc3uv(pixz_i, pixz_o, ff_xcte)) { @@ -62,11 +60,8 @@ int sub_xctecor(SingleGroup *raz, double mjd_obs) { } /* Stuff corrected pixels back to raz for Y-CTE next. */ - for(j=0; jsci.data, i, j) = pixz_o[k]; - } + for(k=0; ksci.data.data[k] = pixz_o[k]; } free(pixz_i); @@ -281,7 +276,7 @@ static void read_reg(float reg[NUM_REG], float *row, int e_n[NUM_ELEVEL], } /* Reset the trap. */ t = 0; - /* How much charge is in the trep? (prop to number of i-shifts to reg) + /* How much charge is in the trap? (prop to number of i-shifts to reg) Jay: 2025 is the number of shifts I calibrated the model to. Not sure why it wasn’t 2048 (the whole chip), but it’s what I used. diff --git a/pkg/wfc3/lib/wf3cte/wf3cte.c b/pkg/wfc3/lib/wf3cte/wf3cte.c index 2af512fa8..ba8704573 100644 --- a/pkg/wfc3/lib/wf3cte/wf3cte.c +++ b/pkg/wfc3/lib/wf3cte/wf3cte.c @@ -138,8 +138,7 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, SingleGroup subcd; /*subarray chip*/ SingleGroup subab; /*subarray chip*/ SingleGroup raz; /* THE LARGE FORMAT COMBINATION OF CDAB*/ - SingleGroup rsz; /* LARGE FORMAT READNOISE CORRECTED IMAGE */ - SingleGroup rsc; /* CTE CORRECTED*/ + SingleGroup razx; SingleGroup rzc; /* FINAL CTE CORRECTED IMAGE */ SingleGroup chg; /* THE CHANGE DUE TO CTE */ SingleGroup raw; /* THE RAW IMAGE IN RAZ FORMAT */ @@ -274,12 +273,6 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, initSingleGroup(&raz); allocSingleGroup(&raz, RAZ_COLS, RAZ_ROWS, True); - initSingleGroup(&rsz); - allocSingleGroup(&rsz, RAZ_COLS, RAZ_ROWS, True); - - initSingleGroup(&rsc); - allocSingleGroup(&rsc, RAZ_COLS, RAZ_ROWS, True); - initSingleGroup(&rzc); allocSingleGroup(&rzc, RAZ_COLS, RAZ_ROWS, True); @@ -294,8 +287,6 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, for (i=0;idata.data, - cte_pars.cprof->data.data, - readNoise, - cte_pars.thresh, - cte_pars.n_forward, - cte_pars.n_par, - rzc.sci.data.data)) { - return status; + if (wf3.subarray) { + i = sub_ctecor_v2c(raz.sci.data.data, + fff.sci.data.data, + WsMAX, + cte_pars.qlevq_data, + cte_pars.dpdew_data, + cte_pars.rprof->data.data, + cte_pars.cprof->data.data, + readNoise, + cte_pars.thresh, + cte_pars.n_forward, + cte_pars.n_par, + rzc.sci.data.data); + } else { + i = sub_ctecor_v2c(razx.sci.data.data, + fff.sci.data.data, + WsMAX, + cte_pars.qlevq_data, + cte_pars.dpdew_data, + cte_pars.rprof->data.data, + cte_pars.cprof->data.data, + readNoise, + cte_pars.thresh, + cte_pars.n_forward, + cte_pars.n_par, + rzc.sci.data.data); + } + if (i) { + return (status = i); } trlmessage("YCTE: Returning from the routine..."); @@ -663,14 +675,14 @@ int WF3cte (char *input, char *output, CCD_Switch *cte_sw, putSingleGroup(output,cd.group_num, &cd,0); putSingleGroup(output,ab.group_num, &ab,0); + + freeSingleGroup(&razx); } /** CLEAN UP ON AISLE 3 **/ freeSingleGroup(&rzc); - freeSingleGroup(&rsc); freeSingleGroup(&chg); freeSingleGroup(&raz); - freeSingleGroup(&rsz); freeSingleGroup(&raw); freeSingleGroup(&cd); freeSingleGroup(&ab); @@ -1552,7 +1564,7 @@ float find_raz2rnoival(float *raz_cdab, float *FLOAT_RNOIVAL, float *FLOAT_BKGDV int ih, iih; long dhist[NUM_BINS], dcum[NUM_BINS], vtot; long vhist[NUM_BINS], vcum[NUM_BINS], dtot; - int ivmin, id1, id2; + int ivmin; int idmin, iv1, iv2; long long vsum; long long nsum; @@ -1570,8 +1582,6 @@ float find_raz2rnoival(float *raz_cdab, float *FLOAT_RNOIVAL, float *FLOAT_BKGDV iv1 = 1; iv2 = 999; - id1 = 1; - id2 = 999; /* * Distill the image variation information and background into quick histograms @@ -1669,8 +1679,6 @@ float find_raz2rnoival(float *raz_cdab, float *FLOAT_RNOIVAL, float *FLOAT_BKGDV for (ih=1; ih<=NUM_BINS-1; ih++) { for (iih=ih+1; iih<=NUM_BINS-1; iih++) { if (dcum[iih-1]-dcum[ih-1] > 0.75*dtot && iih-ih < idmin) { - id1 = ih; - id2 = iih; idmin = iih-ih; }