From b89f8c1883f144be2f749865de91785820ae4108 Mon Sep 17 00:00:00 2001 From: Ken MacDonald Date: Thu, 11 Apr 2024 14:45:43 -0400 Subject: [PATCH] Updating comments and changing dark current usage due to STCAL PR #254. --- src/stcal/ramp_fitting/src/slope_fitter.c | 67 ++++++++++++++--------- 1 file changed, 41 insertions(+), 26 deletions(-) diff --git a/src/stcal/ramp_fitting/src/slope_fitter.c b/src/stcal/ramp_fitting/src/slope_fitter.c index 484a8f500..c12ca78b3 100644 --- a/src/stcal/ramp_fitting/src/slope_fitter.c +++ b/src/stcal/ramp_fitting/src/slope_fitter.c @@ -46,7 +46,18 @@ typedef enum { /* ========================================================================= */ /* GLOBALS */ /* ------------------------------------------------------------------------- */ + +/* This is mostly used for debugging, but could have other usefulness. */ static npy_intp current_integration; + +/* + * Deals with invalid data. This is one of the ways the python code dealt with + * the limitations of numpy that aren't necessary in this code. The LARGE_VARIANCE + * variable has been removed from use in this code, but due to some strange, + * non-flagged data that still requires the use of LARGE_VARIANCE_THRESHOLD , but + * shouldn't. I think that strange data should have been flagged in previous steps + * but I don't think that has happened. + */ const real_t LARGE_VARIANCE = 1.e8; const real_t LARGE_VARIANCE_THRESHOLD = 1.e6; /* ------------------------------------------------------------------------- */ @@ -55,14 +66,20 @@ const real_t LARGE_VARIANCE_THRESHOLD = 1.e6; /* ========================================================================= */ /* MACROS */ /* ------------------------------------------------------------------------- */ + +/* Formatting to make printing more uniform. */ #define DBL "16.10f" + +/* Is more general and non-type dependent. */ #define BSWAP32(X) ((((X) & 0xff000000) >> 24) | \ (((X) & 0x00ff0000) >> 8) | \ (((X) & 0x0000ff00) << 8) | \ (((X) & 0x000000ff) << 24)) +/* Pointers should be set to NULL once freed. */ #define SET_FREE(X) if (X) {free(X); (X) = NULL;} +/* Ensure all allocated memory gets deallocated properly. */ #define FREE_RAMP_DATA(RD) \ if (RD) { \ clean_ramp_data(rd); \ @@ -70,18 +87,21 @@ const real_t LARGE_VARIANCE_THRESHOLD = 1.e6; (RD) = NULL; \ } +/* Ensure all allocated memory gets deallocated properly. */ #define FREE_PIXEL_RAMP(PR) \ if (PR) { \ clean_pixel_ramp(PR); \ SET_FREE(PR); \ } +/* Ensure all allocated memory gets deallocated properly. */ #define FREE_SEGS_LIST(N, S) \ if (S) { \ clean_segment_list(N, S); \ SET_FREE(S);\ } +/* Complicated dereferencing and casting using a label. */ #define VOID_2_DOUBLE(A) (*((double*)(A))) #define VOID_2_FLOAT(A) (*((float*)(A))) #define VOID_2_REAL(A) (*((real_t*)(A))) @@ -120,7 +140,7 @@ struct ramp_data { ssize_t image_sz; /* The size of an image */ ssize_t ramp_sz; /* The size of a pixel ramp */ - /* Are arrays byteswapped? */ + /* Functions to get the proper data. */ float (*get_data)(PyArrayObject*, npy_intp, npy_intp, npy_intp, npy_intp); float (*get_err)(PyArrayObject*, npy_intp, npy_intp, npy_intp, npy_intp); uint32_t (*get_pixeldq)(PyArrayObject*, npy_intp, npy_intp); @@ -161,10 +181,10 @@ struct ramp_data { * (nints, max_num_segs, nrows, ncols). */ - int save_opt; /* Save optional results value */ - int max_num_segs; - struct simple_ll_node ** segs; - real_t * pedestal; + int save_opt; /* Save optional results value */ + int max_num_segs; /* Max number of segments over all ramps. */ + struct simple_ll_node ** segs; /* The segment list for each ramp. */ + real_t * pedestal; /* The pedestal computed for each ramp. */ /* Meta data */ uint32_t suppress_one_group; /* Boolean to suppress one group */ @@ -178,7 +198,7 @@ struct ramp_data { real_t effintim; /* Effective integration time */ real_t one_group_time; /* Time for ramps with only 0th good group */ weight_t weight; /* The weighting for OLS */ -}; +}; /* END: struct ramp_data */ /* * The ramp fit for a specific pixel. @@ -189,7 +209,7 @@ struct pixel_fit { real_t var_poisson; /* Poisson variance */ real_t var_rnoise; /* Read noise variance */ real_t var_err; /* Total variance */ -}; +}; /* END: struct pixel_fit */ /* * The segment information of an integration ramp is kept track of @@ -215,7 +235,7 @@ struct simple_ll_node { real_t yint; /* Y-intercept */ real_t sigyint; /* Uncertainty in the Y-intercept */ real_t weight; /* Sum of weights */ -}; +}; /* END: struct simple_ll_node */ /* * The list of segments in an integration ramp. The segments form the basis @@ -226,7 +246,7 @@ struct segment_list { struct simple_ll_node * tail; /* The tail node of the list */ ssize_t size; /* The number of segments */ npy_intp max_segment_length; /* The max group length of a segment */ -}; +}; /* END: struct segment_list */ /* * For each integration, count how many groups had certain flags set. @@ -237,7 +257,7 @@ struct integ_gdq_stats { int cnt_dnu_sat; /* SATURATED | DO_NOT_USE count */ int cnt_good; /* GOOD count */ int jump_det; /* Boolean for JUMP_DET */ -}; +}; /* END: struct integ_gdq_stats */ /* * This contains all the information to ramp fit a specific pixel. @@ -252,7 +272,7 @@ struct pixel_ramp { real_t * data; /* The 2-D ramp data (nints, ngroups) */ uint32_t * groupdq; /* The group DQ pixel array */ - uint32_t pixeldq;/* The pixel DQ pixel */ + uint32_t pixeldq; /* The pixel DQ pixel */ real_t gain; /* The pixel gain */ real_t rnoise ; /* The pixel read noise */ real_t zframe; /* The pixel ZEROFRAME */ @@ -260,7 +280,7 @@ struct pixel_ramp { /* Timing bool */ uint8_t * is_zframe; /* Boolean to use ZEROFRAME */ - uint8_t * is_0th; /* Boolean to use ZEROFRAME */ + uint8_t * is_0th; /* XXX Boolean to use ZEROFRAME */ /* C computed values */ real_t median_rate; /* The median rate of the pixel */ @@ -275,11 +295,11 @@ struct pixel_ramp { /* initialize and clean */ struct pixel_fit rate; /* Image information */ struct pixel_fit * rateints; /* Cube information */ -}; +}; /* END: struct pixel_ramp */ struct ols_calcs { real_t sumx, sumxx, sumy, sumxy, sumw; -}; +}; /* END: struct ols_calcs */ struct rate_product { int is_none; @@ -288,7 +308,7 @@ struct rate_product { PyArrayObject * var_poisson; PyArrayObject * var_rnoise; PyArrayObject * var_err; -}; +}; /* END: struct rate_product */ struct rateint_product { int is_none; @@ -297,7 +317,7 @@ struct rateint_product { PyArrayObject * var_poisson; PyArrayObject * var_rnoise; PyArrayObject * var_err; -}; +}; /* END: struct rateint_product */ struct opt_res_product { PyArrayObject * slope; /* Slope of segment */ @@ -313,7 +333,8 @@ struct opt_res_product { PyArrayObject * weights; /* Weights */ PyArrayObject * cr_mag; /* Cosmic ray magnitudes */ -}; +}; /* END: struct opt_res_product */ + /* ------------------------------------------------------------------------- */ /* ========================================================================= */ @@ -2591,9 +2612,7 @@ ramp_fit_pixel_integration_fit_slope_seg_len1( seg->slope = pr->data[idx] / timing; pden = (timing * pr->gain); - /* XXX PR # 254. Is there a JP ticket for this? */ - /* seg->var_p = (pr->median_rate pr->dcurrent) / pden; */ - seg->var_p = pr->median_rate / pden + pr->dcurrent; + seg->var_p = (pr->median_rate pr->dcurrent) / pden; /* Segment read noise variance */ rnum = pr->rnoise / timing; @@ -2638,9 +2657,7 @@ ramp_fit_pixel_integration_fit_slope_seg_len2( /* Segment Poisson variance */ if (pr->median_rate > 0.) { pden = (rd->group_time * pr->gain); - /* XXX PR # 254. Is there a JP ticket for this? */ - /* seg->var_p = (pr->median_rate pr->dcurrent) / pden; */ - seg->var_p = pr->median_rate / pden + pr->dcurrent; + seg->var_p = (pr->median_rate pr->dcurrent) / pden; } else { seg->var_p = pr->dcurrent; } @@ -2801,9 +2818,7 @@ ramp_fit_pixel_integration_fit_slope_seg_default_weighted_seg( /* Segment Poisson variance */ if (pr->median_rate > 0.) { pden = (rd->group_time * pr->gain * (seglen - 1.)); - /* XXX PR # 254. Is there a JP ticket for this? */ - /* seg->var_p = (pr->median_rate pr->dcurrent) / pden; */ - seg->var_p = pr->median_rate / pden + pr->dcurrent; + seg->var_p = (pr->median_rate pr->dcurrent) / pden; } else { seg->var_p = pr->dcurrent; }