From d54fe7c1f704a63824c5bfa0ece65245572e9b27 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 4 Mar 2015 21:21:30 -0500 Subject: Initial commit --- src/libcf/cf_grating_motion.c | 432 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 432 insertions(+) create mode 100644 src/libcf/cf_grating_motion.c (limited to 'src/libcf/cf_grating_motion.c') diff --git a/src/libcf/cf_grating_motion.c b/src/libcf/cf_grating_motion.c new file mode 100644 index 0000000..fe50d63 --- /dev/null +++ b/src/libcf/cf_grating_motion.c @@ -0,0 +1,432 @@ +/**************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ***************************************************************************** + * + * Synopsis: cf_grating_motion(header, nevents, time, x, y, channel, + * ntimes, ttime, tsunrise) + * + * Description: Reads spectral shift caused by thermal grating motions + * and corrects the X and Y locations of each photon event. + * + * fitsfile *header Pointer to the FITS header of the + * Intermediate Data File + * long nevents Number of photons in the file + * int *time Time (in seconds) since exposure start + * float *x Position of the photon (in pixels) + * float *y Position of the photon (in pixels) + * unsigned char channel Channel ID of each photon + * long ntimes Number of entries in timeline table + * float ttime Time array of timeline table + * short tsunrise Time since last sunrise + * + * Returns: 0 on success + * + * History: 09/05/02 1.0 RDR Begin work; adapted from cf_make_shift + * 04/21/03 1.3 wvd Use tsunrise array from timeline table. + * Do not assume that ttime is continuous. + * 05/20/03 1.4 rdr Added call to cf_proc_check + * 08/21/03 1.5 wvd Change channel to unsigned char. + * 06/25/04 1.6 wvd Estimate time between sunrises using + * orbit period in file header. + * 03/22/05 1.7 wvd Change tsunrise from float to short. + * Read orbital period from file header. + * 04/15/05 1.8 wvd Close GRAT_CAL file before exiting. + * If ORBPERID keyword is not found, + * assume that orbit period is 6000 s. + * 12/30/05 1.9 wvd Current algorithm corrects for orbital + * drift of spectrum, but wavelength zero + * point still varies with beta angle and + * observation date. New subroutine + * corrects LiF 1 data for these effects. + * 07/06/06 2.0 wvd Rewrite to use Dave Sahnow's new + * grating-correction file (version 003). + * 07/12/06 2.1 wvd Modify to read an arbitrary number of + * extensions, each with phase information + * in file header keywords. + * 10/20/06 2.2 wvd Grating motion also depends on APER_PA. + * Now we read it from the file header + * and select appropriate coefficients. + * If no correction is defined for a + * particular set of parameters, set + * 11/28/06 2.3 wvd Grating motion includes a long-term, + * non-periodic drift. Read it from the + * first four extensions of GRAT_CAL file. + * 04/07/07 2.4 wvd Clean up compiler warnings. + * + ****************************************************************************/ + +#include +#include +#include +#include +#include "calfuse.h" + +static int +read_nonperiodic_shifts(fitsfile *header, fitsfile *shiftfits, + long ntimes, float *dxlif, float *dylif, float *dxsic, float *dysic) +{ + char channel[FLEN_VALUE], comment[FLEN_CARD], detector[FLEN_VALUE]; + int nrows, tndx; + int status=0, hdutype=0, hdunum=0; + long i; + float *mjd, *xshift, *yshift; + double expstart; + + /* Read header keywords from data file. */ + fits_read_key(header, TDOUBLE, "EXPSTART", &expstart, NULL, &status); + FITS_read_key(header, TSTRING, "DETECTOR", detector, NULL, &status); + + /* First do the LiF channels */ + if (!strncmp(detector,"1",1)) { + hdunum=2; + } else if (!strncmp(detector,"2",1)) { + hdunum=4; + } else { + cf_if_error("Cannot parse DETECTOR keyword in data file."); + } + cf_verbose(3, "Reading HDU %d of GRAT_CAL file.", hdunum); + + FITS_movabs_hdu(shiftfits, hdunum, &hdutype, &status); + FITS_read_key(shiftfits, TSTRING, "CHANNEL", channel, NULL, &status); + sprintf(comment, "LiF%c", detector[0]); + if (strncmp(channel, comment, 4)) + cf_if_error("Cannot find %s correction info in extension %d", + comment, hdunum); + + /* Read X and Y corrections as a function of time. */ + nrows = cf_read_col(shiftfits, TFLOAT, "MJD", (void *) &mjd); + nrows = cf_read_col(shiftfits, TFLOAT, "XSHIFT", (void *) &xshift); + nrows = cf_read_col(shiftfits, TFLOAT, "YSHIFT", (void *) &yshift); + + /* Select appropriate MJD for this exposure. */ + i = 0; + while (expstart > mjd[i] && i < nrows) i++; + tndx = i-1; + if (tndx < 0) tndx = 0; + cf_verbose(3,"EXPSTART = %g", expstart); + cf_verbose(3,"Using LiF corrections for MJD >= %g", mjd[tndx]); + if (tndx < nrows-1) + cf_verbose(3," and MJD < %g", mjd[tndx+1]); + cf_verbose(3,"xshift = %f, yshift = %f", xshift[tndx], yshift[tndx]); + + /* Copy X and Y corrections to output arrays. */ + for (i = 0; i < ntimes; i++) { + dxlif[i] = xshift[tndx]; + dylif[i] = yshift[tndx]; + } + + /* Now do the SiC channels */ + hdunum += 1; + cf_verbose(3, "Reading HDU %d of GRAT_CAL file.", hdunum); + FITS_movabs_hdu(shiftfits, hdunum, &hdutype, &status); + FITS_read_key(shiftfits, TSTRING, "CHANNEL", channel, NULL, &status); + sprintf(comment, "SiC%c", detector[0]); + if (strncmp(channel, comment, 4)) + cf_if_error("Cannot find %s correction info in extension %d", + comment, hdunum); + + /* Read X and Y corrections as a function of time. */ + nrows = cf_read_col(shiftfits, TFLOAT, "MJD", (void *) &mjd); + nrows = cf_read_col(shiftfits, TFLOAT, "XSHIFT", (void *) &xshift); + nrows = cf_read_col(shiftfits, TFLOAT, "YSHIFT", (void *) &yshift); + + /* Select appropriate MJD for this exposure. */ + i = 0; + while (expstart > mjd[i] && i < nrows) i++; + tndx = i-1; + if (tndx < 0) tndx = 0; + cf_verbose(3,"Using SiC corrections for MJD >= %g", mjd[tndx]); + if (tndx < nrows-1) + cf_verbose(3," and MJD < %g", mjd[tndx+1]); + cf_verbose(3,"xshift = %f, yshift = %f", xshift[tndx], yshift[tndx]); + + /* Copy X and Y corrections to output arrays. */ + for (i = 0; i < ntimes; i++) { + dxsic[i] = xshift[tndx]; + dysic[i] = yshift[tndx]; + } + + return 0; +} + + +static int +read_periodic_shifts(fitsfile *header, fitsfile *shiftfits, int ncorrection, + double *lif_mjd0, double *lif_period, double *sic_mjd0, double *sic_period, + double *lif_x_coef, double *lif_y_coef, double *sic_x_coef, double *sic_y_coef) +{ + char detector[FLEN_CARD], channel[FLEN_CARD], comment[FLEN_CARD]; + int status=0, hdutype=0, hdunum=0, anynul=0; + int i, nregions, region, typecode; + long max_coeff, ncoeff, width; + float aper_pa, *aperpalo, *aperpahi; + float *betalo, *betahi, *polelo, *polehi; + double beta, pole, sunang; + + /* Initialize variables. */ + *lif_mjd0 = *lif_period = *sic_mjd0 = *sic_period = 0.; + for (i = 0; i < 11; i++) { + lif_x_coef[i] = 0.; + lif_y_coef[i] = 0.; + sic_x_coef[i] = 0.; + sic_y_coef[i] = 0.; + } + + FITS_read_key(header, TFLOAT, "APER_PA", &aper_pa, NULL, &status); + FITS_read_key(header, TDOUBLE, "SUNANGLE", &sunang, NULL, &status); + FITS_read_key(header, TDOUBLE, "POLEANGL", &pole, NULL, &status); + FITS_read_key(header, TSTRING, "DETECTOR", detector, NULL, &status); + while (999.9 > aper_pa && aper_pa > 360.) aper_pa -= 360.; + while (aper_pa < 0.) aper_pa += 360.; + beta = 180.0 - sunang; + + cf_verbose(3, "S/C orientation: beta = %5.1f, pole = %5.1f, aper_pa = %5.1f", + beta, pole, aper_pa); + + /* First do the LiF channels */ + if (!strncmp(detector,"1",1)) { + hdunum=2; + } else if (!strncmp(detector,"2",1)) { + hdunum=4; + } else { + cf_if_error("Cannot parse DETECTOR keyword in read_shift_file"); + } + hdunum += ncorrection * 4; + cf_verbose(3, "Reading HDU %d of GRAT_CAL file.", hdunum); + + FITS_movabs_hdu(shiftfits, hdunum, &hdutype, &status); + FITS_read_key(shiftfits, TSTRING, "CHANNEL", channel, NULL, &status); + sprintf(comment, "LiF%c", detector[0]); + if (strncmp(channel, comment, 4)) + cf_if_error("Cannot find %s correction info in extension %d", + comment, hdunum); + + if (ncorrection > 1) { + FITS_read_key(shiftfits, TDOUBLE, "PERIOD", lif_period, NULL, &status); + FITS_read_key(shiftfits, TDOUBLE, "MJD_ZERO", lif_mjd0, NULL, &status); + cf_verbose(3, "%s: period = %g, mjd_zero = %g", + channel, *lif_period, *lif_mjd0); + } + + nregions = cf_read_col(shiftfits, TFLOAT, "BETALO", (void *) &betalo); + nregions = cf_read_col(shiftfits, TFLOAT, "BETAHI", (void *) &betahi); + nregions = cf_read_col(shiftfits, TFLOAT, "POLELO", (void *) &polelo); + nregions = cf_read_col(shiftfits, TFLOAT, "POLEHI", (void *) &polehi); + nregions = cf_read_col(shiftfits, TFLOAT, "APERPALO", (void *) &aperpalo); + nregions = cf_read_col(shiftfits, TFLOAT, "APERPAHI", (void *) &aperpahi); + FITS_get_coltype(shiftfits, 7, &typecode, &ncoeff, &width, &status); + max_coeff = ncoeff; + + for (i = 0; i < nregions; i++) + if (betalo[i] <= beta && beta < betahi[i] && + polelo[i] <= pole && pole < polehi[i] && + aperpalo[i] <= aper_pa && aper_pa < aperpahi[i]) break; + + if (i == nregions) + cf_verbose(3, "LiF correction not defined for beta = %g, pole = %g, " + "aper_pa = %g", beta, pole, aper_pa); + else { + cf_verbose(3, "LIF REGION %2d: beta = %g - %g, pole = %g - %g,", + i, betalo[i], betahi[i], polelo[i], polehi[i]); + cf_verbose(3, " aper_pa = %g - %g, %d coefficients", + aperpalo[i], aperpahi[i], ncoeff); + + region = i+1; /* CFITSIO counts from 1, not 0. */ + FITS_read_col(shiftfits, TDOUBLE, 7, region, 1, ncoeff, NULL, lif_x_coef, + &anynul, &status); + FITS_read_col(shiftfits, TDOUBLE, 8, region, 1, ncoeff, NULL, lif_y_coef, + &anynul, &status); + } + + /* Now the SiC channels */ + hdunum += 1; + cf_verbose(3, "Reading HDU %d of GRAT_CAL file.", hdunum); + FITS_movabs_hdu(shiftfits, hdunum, &hdutype, &status); + FITS_read_key(shiftfits, TSTRING, "CHANNEL", channel, NULL, &status); + sprintf(comment, "SiC%c", detector[0]); + if (strncmp(channel, comment, 4)) + cf_if_error("Cannot find %s correction info in extension %d", + comment, hdunum); + + if (ncorrection > 1) { + FITS_read_key(shiftfits, TDOUBLE, "PERIOD", sic_period, NULL, &status); + FITS_read_key(shiftfits, TDOUBLE, "MJD_ZERO", sic_mjd0, NULL, &status); + cf_verbose(3, "%s: period = %g, mjd_zero = %g", + channel, *sic_period, *sic_mjd0); + } + + nregions = cf_read_col(shiftfits, TFLOAT, "BETALO", (void *) &betalo); + nregions = cf_read_col(shiftfits, TFLOAT, "BETAHI", (void *) &betahi); + nregions = cf_read_col(shiftfits, TFLOAT, "POLELO", (void *) &polelo); + nregions = cf_read_col(shiftfits, TFLOAT, "POLEHI", (void *) &polehi); + nregions = cf_read_col(shiftfits, TFLOAT, "APERPALO", (void *) &aperpalo); + nregions = cf_read_col(shiftfits, TFLOAT, "APERPAHI", (void *) &aperpahi); + FITS_get_coltype(shiftfits, 7, &typecode, &ncoeff, &width, &status); + if (max_coeff < ncoeff) max_coeff = ncoeff; + + for (i = 0; i < nregions; i++) + if (betalo[i] <= beta && beta < betahi[i] && + polelo[i] <= pole && pole < polehi[i] && + aperpalo[i] <= aper_pa && aper_pa < aperpahi[i]) break; + + if (i == nregions) + cf_verbose(3, "SiC correction not defined for beta = %g, pole = %g, " + "aper_pa = %g", beta, pole, aper_pa); + else { + cf_verbose(3, "SIC REGION %2d: beta = %g - %g, pole = %g - %g,", + i, betalo[i], betahi[i], polelo[i], polehi[i]); + cf_verbose(3, " aper_pa = %g - %g, %d coefficients", + aperpalo[i], aperpahi[i], ncoeff); + + region = i+1; /* CFITSIO counts from 1, not 0. */ + FITS_read_col(shiftfits, TDOUBLE, 7, region, 1, ncoeff, NULL, sic_x_coef, + &anynul, &status); + FITS_read_col(shiftfits, TDOUBLE, 8, region, 1, ncoeff, NULL, sic_y_coef, + &anynul, &status); + } + + cf_verbose(2, "Fourier coefficients of grating-motion correction:"); + cf_verbose(2," LiF X LiF Y SiC X SiC Y"); + for (i=0; i 0 && channel[j] < 4) { + x[j] += dxlif[k]; + y[j] += dylif[k]; + } + else if (channel[j] > 4) { + x[j] += dxsic[k]; + y[j] += dysic[k]; + } + } + + /* Release memory. */ + free(dxlif); + free(dylif); + free(dxsic); + free(dysic); + + cf_proc_update(header, CF_PRGM_ID, "COMPLETE"); + cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished execution."); + + return status; +} -- cgit