diff options
Diffstat (limited to 'src/libcf/cf_astigmatism.c')
| -rw-r--r-- | src/libcf/cf_astigmatism.c | 198 |
1 files changed, 198 insertions, 0 deletions
diff --git a/src/libcf/cf_astigmatism.c b/src/libcf/cf_astigmatism.c new file mode 100644 index 0000000..e02f361 --- /dev/null +++ b/src/libcf/cf_astigmatism.c @@ -0,0 +1,198 @@ +/***************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ***************************************************************************** + * + * Description: Modifies X array to correct for instrumental astigmatism. + * + * Since we have no astigmatism correction for extended + * sources, this program does not attempt to apply one. + * + * In keeping with the CalFUSE standard, the astigmatism + * correction is added to the measured X position of each photon. + * Versions of the ASTG_CAL file prior to v009 have the opposite + * sign and are incompatible with this version of the program. + * + * Returns: 0 upon successful completion. + * + * History: 12/02/02 1.1 jch Adapt from cf_astig.c + * 12/20/02 1.2 wvd Install + * 02/11/03 1.3 wvd Apply astig correction only if + * YCENT# > O for target aperture. + * Add ASTG_KEY keyword. Comment field: + * Performed? 1=LiF, 2=SiC, 3=Both + * 02/13/03 1.4 wvd Change all indices to type long. + * 02/24/03 1.5 peb Changed include file to calfusettag.h + * 03/04/03 1.6 peb Removed unused variables and math.h. + * Corrected sprintf:185 argument error. + * 03/11/03 1.7 wvd Change channel to type char + * 03/12/03 1.8 wvd If wavelength undefined, set + * channel[k] = 0 + * 04/04/03 1.9 wvd Test for overflow of astig array. + * 04/07/03 1.10 wvd Delete test for quality of YCENT. + * Delete ASTG_KEY keyword. + * Don't round yoff to nearest integer. + * 04/09/03 1.11 wvd Write name of ASTG_CAL file to trailer. + * 05/20/03 1.12 rdr Added call to cf_proc_check + * 09/16/03 1.13 wvd Change sign of ASTG_CAL file so that + * correction is additive. + * 09/22/03 1.14 wvd Initialize overflow counter. + * Print out number of overflows. + * 10/17/03 1.15 wvd Discard photons that fall outside of + * the astigmatism-correction window. + * 10/31/03 1.16 wvd Change channel to unsigned char. + * 11/11/03 1.17 wvd Interpolate between tabulated + * wavelength values. Use cf_read_col + * to read wavecal files. + * 02/17/04 1.18 wvd Replace cf_nint() with cf_nlong() + * 03/02/04 1.19 wvd Change name of assign_wavelengths + * to cf_x2lambda and make it + * generally accessible. + * 07/16/04 1.20 wvd Must change ASTG_COR to COMPLETE or + * SKIPPED by hand. + * 02/18/05 1.21 wvd Add some diagnostic output. + * 11/30/05 1.22 wvd Don't complain when photon events + * fall outside of astig images. + * 05/15/06 1.23 wvd Divide cf_astigmatism_and_dispersion + * into two separate routines. + * + ****************************************************************************/ + +#include <string.h> +#include <stdio.h> +#include <stdlib.h> +#include <math.h> +#include "calfuse.h" + +static char CF_PRGM_ID[] = "cf_astigmatism"; +static char CF_VER_NUM[] = "1.23"; + + +/* + * Loop over the apertures, read appropriate astigmatism correction, + * and add to the input X shifts of pixels within extraction window. + */ +static int +add_astig_shifts(fitsfile *infits, long nevents, float *x, float *y, + unsigned char *channel) +{ + char astig_file[FLEN_VALUE], keyword[FLEN_VALUE]; + int status=0; + int ap, active_ap[2], ast_cen, extno, i, overflow; + long k, ii, jj, kk, npix, nx, ny; + float centroid, yoff; + float *astig=NULL; + fitsfile *astigfits; + + /* Read target apertures from the file header */ + (void) cf_source_aper (infits, active_ap); + + /* Open the astigmatism file */ + FITS_read_key(infits, TSTRING, "ASTG_CAL", astig_file, NULL, &status); + FITS_open_file(&astigfits, cf_cal_file(astig_file), READONLY, &status); + + /* Since we presently have no astig correction for extended + * sources, we apply correction only to target aperture. + */ + for (i = 0; i < 2; i++) { + ap = active_ap[i]; + + /* Read the centroid of the aperture from the input file. */ + sprintf(keyword, "YCENT%1d", ap); + FITS_read_key(infits, TFLOAT, keyword, ¢roid, NULL, &status); + cf_verbose(3, "Target aperture: %d\t Y centroid: %f", ap, centroid); + + /* Read astigmatism correction for this aperture */ + extno = ap+1; + cf_verbose(3, "Reading extension %d of %s", extno, astig_file); + FITS_movabs_hdu(astigfits, extno, NULL, &status); + FITS_read_key(astigfits, TINT, "SLIT_CEN", &ast_cen, NULL, &status); + FITS_read_key(astigfits, TLONG, "NAXIS1", &nx, NULL, &status); + FITS_read_key(astigfits, TLONG, "NAXIS2", &ny, NULL, &status); + + npix = nx * ny; + astig = (float *) cf_malloc(sizeof(float) * npix); + FITS_read_img(astigfits, TFLOAT, 1L, npix, 0, astig, NULL, &status); + + /* Determine offset between spectrum and astig correction file */ + yoff = centroid - ast_cen; + cf_verbose(3, "Offset between spectrum and astig correction: %f", yoff); + + /* Go through the event list and find the astigmatism correction + * appropriate for the photon's (x, y) position. + */ + overflow = 0; + for (k = 0; k < nevents; k++) { + if (channel[k] == ap) { + if (((ii = cf_nlong(x[k])) < 0 || ii >= nx) || + ((jj = cf_nlong(y[k] - yoff)) < 0 || jj >= ny) || + ((kk = jj * nx + ii) >= npix)) { + channel[k] = 0; + overflow++; + cf_verbose(3, "Cannot correct pixel %05ld (%f, %f)", + k, x[k], y[k]); + } + else x[k] += astig[kk]; + } + } + /* If overflow flag is set, there is a problem either with + the ASTG_CAL file or the input photon list. + if (overflow) + cf_if_warning("%d events fell outside of aperture %d ASTG_CAL window", + overflow, ap); + Comment this out, as it is always triggered by cf_bad_pixels. - wvd (11/30/05) */ + + /* Space for astig array is allocated in each loop. */ + free(astig); + cf_verbose(3, "End of loop for aperture %d", ap); + } + + FITS_close_file(astigfits, &status); + return status; +} + + +/* + * If APERTURE = RFPT or target is not a point source, exit without + * applying an astigmatism correction. + */ +int cf_astigmatism(fitsfile *infits, long nevents, float *x, + float *y, unsigned char *channel) +{ + char aperture[FLEN_VALUE]; + char source_type[FLEN_VALUE]; + int errflg=0, fileok=TRUE, status=0; + + /* Enter a timestamp into the log. */ + cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr); + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing"); + + if ((errflg = cf_proc_check(infits, CF_PRGM_ID))) return errflg; + + FITS_read_key(infits, TSTRING, "APERTURE", aperture, NULL, &status); + if (!strncmp(aperture, "RFPT", 4)) { + cf_verbose(1, "Aperture is RFPT. No astig correction applied."); + fileok = FALSE; + } + + /* Check source type */ + FITS_read_key(infits, TSTRING, "SRC_TYPE", source_type, NULL, &status); + if (strncmp(source_type, "P", 1)) { + cf_verbose(1, "Source type is %s. No astig correction applied.", + source_type); + fileok = FALSE; + } + + if (fileok) { + add_astig_shifts(infits, nevents, x, y, channel); + cf_proc_update(infits, CF_PRGM_ID, "COMPLETE"); + } + else + cf_proc_update(infits, CF_PRGM_ID, "SKIPPED"); + + /* Update processing flags. */ + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done Processing"); + + return fileok; +} |