diff options
Diffstat (limited to 'src/libcf/cf_identify_channel.c')
| -rw-r--r-- | src/libcf/cf_identify_channel.c | 201 |
1 files changed, 201 insertions, 0 deletions
diff --git a/src/libcf/cf_identify_channel.c b/src/libcf/cf_identify_channel.c new file mode 100644 index 0000000..7d6a8e3 --- /dev/null +++ b/src/libcf/cf_identify_channel.c @@ -0,0 +1,201 @@ +/***************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ***************************************************************************** + * + * Synopsis: cf_identify_channel(fitsfile *header, long nevents, + * float *xfarf, float *yfarf, unsigned char *chan, + * unsigned char *locflags, int pad, int final_call) + * + * Description: Assigns an aperture to each photon in the IDF. + * + * Change of logic in v1.11: This routine may be called + * multiple times. All photons are examined and assigned to + * a channel unless final_call = TRUE. On the final call, only + * photons with channel > 0 are considered. Because only photons + * assigned to a channel are motion corrected, the spectrum can + * move on top of previously-rejected photons. We don't want + * to include them in the final spectrum, so leave their channel + * values unchanged. + * + * Change of logic in v1.16: We've decided to use extended- + * source apertures for the non-target channels. These + * apertures can overlap, so we assign disputed photons to + * their most likely source. In decreasing order, these are + * the target aperture, the larger non-target aperture, and + * the smaller non-target aperture. + * + * Arguments: fitsfile *header Pointer to the location of the FITS + * header of the Intermediate Data File + * long nevents Number of photons in the file + * float *xfarf X positions for each photon event + * float *yfarf Y positions for each photon event + * unsigned char *chan Pointer to pointer to the array + * containing the channel information + * unsigned char *locflags Location flags for each photon + * int pad Switch saying whether to expand the + * width of the tabulated extraction + * window (yes if > 0). The number refers + * to the amount of padding (e.g. pad=5 + * will add 5 pixels to each side of the + * extraction window. + * int final_call This subroutine may be called + * many times. The last time, set + * final_call = TRUE. + * + * Return: 0 on success + * + * History: 08/09/02 1.0 RDR Begin work + * 09/04/02 1.1 rdr Change program so that it just fills + * the channel array rather than creating + * and filling it. Decreased the size of + * the mask. Use the value of pad to + * specify degree of window padding. + * Check for extended source. + * 12/11/02 1.2 RDR Adjust aperture limits by measured + * y centroid (if available) + * 02/11/03 1.3 rdr changed test for ycent values + * 02/12/03 1.6 wvd Don't crash if YCENT# is undefined + * 03/25/03 1.7 rdr Ignore pinhole aperture + * 03/31/03 1.10 wvd Ignore regions where YMAX = YMIN + 1. + * They are not valid channel regions. + * 04/07/03 1.11 wvd Shift channel boundaries to match FPA + * positions. Implement use of cf_verbose. + * Update CHID_COR if final_call = TRUE. + * 05/20/03 1.12 rdr Added call to cf_proc_check + * 07/23/03 1.14 rdr Correct procedure which reads CHID + * calibration file + * 08/21/03 1.15 wvd Change channel to unsigned char. + * 08/25/03 1.16 wvd Use extended windows for non-target + * apertures. Replace aperture mask + * with a grid of aperture limits. + * Use cf_nint() to round floats. + * 08/25/03 1.17 wvd Change coltype from string to int + * in cf_read_col. + * 09/08/03 1.18 wvd Default value is channel[i] = 0. + * Identify low and high values of x + * by setting y limits to -NYMAX. + * 09/18/03 1.19 wvd Consider only photons falling on the + * active area of the detector. + * 11/12/03 1.20 wvd Don't bother with non-target channels + * in HIST mode. + * 06/14/04 1.21 wvd Don't bother shifting extraction + * windows to correct for FPA position. + * If known, shift is not needed. + * 06/17/04 1.22 wvd Because the extraction windows don't + * quite match the WGTS arrays, we pad + * windows by XPAD pixels in X. + * 07/21/04 1.23 wvd Change i from long to int. + * 03/15/05 1.24 wvd Call cf_extraction_limits to get limits + * of extraction windows, rather than + * reading CHID_CAL file directly. + * Do not pad windows in X. + * 03/16/05 1.25 wvd Pass srctype to cf_extraction_limits. + * + ****************************************************************************/ + +#include <string.h> +#include "calfuse.h" + +int +cf_identify_channel(fitsfile* header, long npts, float* x, float* y, + unsigned char* channel, unsigned char* locflags, int pad, + int final_call) +{ + char CF_PRGM_ID[] = "cf_identify_channel"; + char CF_VER_NUM[] = "1.25"; + + char instmode[FLEN_VALUE], zero=0; + int errflg=0, extended, n_chan=6, status=0; + int active_ap[2], chan, srctype[8]; + int i; /* Always steps through targ_ap */ + int *targ_ap, + lwrs[]={3,7,2,6,1,5}, /* LWRS, MDRS, HIRS (LiF & SiC) */ + mdrs[]={2,6,3,7,1,5}, /* MDRS, LWRS, HIRS (LiF & SiC) */ + hirs[]={1,5,3,7,2,6}; /* HIRS, LWRS, MDRS (LiF & SiC) */ + long j, /* Always steps through extraction window arrays */ + k; /* Always steps through photon list */ + short xint, yint, xmin, xmax, *ymin[6], *ymax[6]; + + /* Initialize error checking */ + cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); + + /* Enter a time stamp into the log */ + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing"); + + /* Check whether routine is appropriate for this data file. */ + if ((errflg = cf_proc_check(header, CF_PRGM_ID))) return errflg; + + /* If data were taken in HIST mode, ignore non-target channels. */ + FITS_read_key(header, TSTRING, "INSTMODE", instmode, NULL, &status); + if (!strncmp(instmode, "HIST", 4)) n_chan = 2; + + /* Determine source type and active apertures */ + extended = cf_source_aper(header, active_ap); + + /* Set the order in which we will assign photons to channels. */ + switch (active_ap[0]) { + case 1: targ_ap = hirs; + break; + case 2: targ_ap = mdrs; + break; + default: targ_ap = lwrs; + break; + } + + /* For target apertures, srctype reflects target. + Use extended apertures for the rest. */ + for (i = 0; i < 2; i++) srctype[i] = extended; + for ( ; i < n_chan; i++) srctype[i] = 1; + + /* While we're at it, here's one more bit of i/o. */ + if (pad > 0) + cf_verbose(3, "Padding extraction windows by %d Y pixels", pad); + + /* + * Loop through all channels, read extraction limits, + * and apply any requested padding. + */ + for (i = 0; i < n_chan; i++) { + chan = targ_ap[i]; + + /* Read the extraction limits for this aperture. */ + (void) cf_extraction_limits(header, chan, srctype[i], + (ymin+i), (ymax+i), &xmin, &xmax); + + /* Add any padding requested by the user. */ + if (pad > 0) { + for (j = xmin; j <= xmax; j++) { + ymin[i][j] -= pad; + ymax[i][j] += pad; + } + } + } + + /* Loop through photon list. Assign a channel ID to each. */ + cf_verbose(3, "Entering loop to assign channel ID's."); + for (k = 0; k < npts; k++) { + if (final_call && channel[k] == zero) continue; + channel[k] = zero; + if (!(locflags[k] & LOCATION_SHLD)) { + xint = (short) cf_nint(x[k]); + yint = (short) cf_nint(y[k]); + if (xint < xmin || xint > xmax) continue; + for (i = 0; i < n_chan; i++) { + if (yint >= ymin[i][xint] && yint <= ymax[i][xint]) { + channel[k] = (char) targ_ap[i]; + break; + } + } + } + } + + /* If final_call = TRUE, update IDF header keyword for this subroutine. */ + if (final_call) + cf_proc_update(header, CF_PRGM_ID, "COMPLETE"); + + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done Processing"); + + return status; +} |