Skip to content

Commit

Permalink
Updating comments and changing dark current usage due to STCAL PR spa…
Browse files Browse the repository at this point in the history 
…cetelescope#254.
  • Loading branch information
@kmacdonald-stsci
kmacdonald-stsci committed Apr 11, 2024
1 parent 6068890 commit b89f8c1
Showing 1 changed file with 41 additions and 26 deletions.
67 changes: 41 additions & 26 deletions src/stcal/ramp_fitting/src/slope_fitter.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
/* ------------------------------------------------------------------------- */
Expand All @@ -55,33 +66,42 @@ 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); \
free(RD); \
(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)))
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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 */
Expand All @@ -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.
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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.
Expand All @@ -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.
Expand All @@ -252,15 +272,15 @@ 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 */
real_t dcurrent; /* The pixel average dark current */

/* 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 */
Expand All @@ -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;
Expand All @@ -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;
Expand All @@ -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 */
Expand All @@ -313,7 +333,8 @@ struct opt_res_product {
PyArrayObject * weights; /* Weights */

PyArrayObject * cr_mag; /* Cosmic ray magnitudes */
};
}; /* END: struct opt_res_product */

/* ------------------------------------------------------------------------- */

/* ========================================================================= */
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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;
}
Expand Down Expand Up @@ -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;
}
Expand Down

0 comments on commit b89f8c1

Please sign in to comment.