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_init_support.c | 1344 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1344 insertions(+) create mode 100644 src/libcf/cf_init_support.c (limited to 'src/libcf/cf_init_support.c') diff --git a/src/libcf/cf_init_support.c b/src/libcf/cf_init_support.c new file mode 100644 index 0000000..68db0f6 --- /dev/null +++ b/src/libcf/cf_init_support.c @@ -0,0 +1,1344 @@ +/***************************************************************************** + * Johns Hopkins University + * Center For Astrophysical Sciences + * FUSE + ***************************************************************************** + * + * Description: Routines used by cf_ttag_init and cf_hist_init. + * + * History: 07/15/03 v1.1 wvd Move routines from cf_ttag_init + * 07/24/03 1.2 wvd If tsunrise and tsunset get huge, + * store them as floats, not integers. + * 07/30/03 1.5 wvd Omit TZERO and TSCALE for these + * timeline table entries: + * HIGH_VOLTAGE + * LIF_CNT_RATE + * SIC_CNT_RATE + * FEC_CNT_RATE + * AIC_CNT_RATE + * BKGD_CNT_RATE + * This will effectively convert them + * to arrays of shorts. + * Add 0.5 to all count rates in + * anticipation of truncation. + * 08/22/03 1.7 wvd Populate MIN_LIMB keyword. + * 09/03/03 1.8 wvd Implement cf_verbose throughout. + * Add EXP_JITR and EXPNIGHT to IDF header. + * 09/22/03 1.9 wvd Correct comment field for YCENT1 + * 10/02/03 1.10 wvd Modify timeline table: + * - no time gaps + * - set TEMPORAL_OPUS for times outside + * of OPUS good-time intervals. + * Delete cf_get_times. + * Use cf_read_col in cf_get_gti + * and cf_get_geocorona. + * Don't modify gtis[0] for HIST data. + * 10/13/03 1.11 wvd Add 1s to end of timeline table for + * TTAG data. + * 10/15/03 1.12 wvd Swap voltages for detectors 2A and 2B + * if housekeeping file was written + * before April of 2003. + * 10/22/03 1.13 wvd Compute EXPNIGHT for raw data file + * and write to output file header. + * Add YQUALLIF & YQUALSIC keywords. + * 11/25/03 1.14 wvd Raise upper limits on LIF_CNT_RATE and + * SIC_CNT_RATE to 32000. Raise limits + * on FEC_CNT_RATE and AIC_CNT_RATE + * to 64000. Use TZERO = 32768 for + * FEC_CNT_RATE and AIC_CNT_RATE. + * They must be read as floats. + * 01/15/04 1.15 bjg Now performs cnt0 -= 16777216 when + * cnt0 < cnt before computing + * cnt_rate for the first cnt and + * cnt0 in function + * fill_cntrate_array + * cf_timeline 1.6 -> 1.7 + * free the arrays only at the end + * version 1.6 tries to read tflag + * array after it has been freed + * 02/06/04 1.16 wvd Modify fill_cntrate_array() so + * that the output array no longer runs + * 16 seconds behind the HSKP data. + * 02/18/04 1.17 wvd Voltages for detectors 2A and 2B + * are still swapped. Set HKERR_ENDS + * to 20100000. + * 03/04/04 1.18 bjg Read begin and end counters as float + * Set countrate to zero when cntb is more + * than 16777216. + * 03/19/04 bjg FITS data types (TFLOAT, TSTRING, TBYTE) + * now match type of C containers (float, + * char *, char) when reading and writing + * in FITS header. + * 04/06/04 1.19 bjg Include ctype.h + * Functions now return EXIT_SUCCESS + * Remove unused variables + * Remove static keyword in struct key + * definition + * If cnt_rate>32000 set cnt_rate=0 + * 05/07/04 1.20 wvd Delete HKERR_ENDS. If HSKPVERS keyword + * is not present in HSKP file, test to + * determine whether high-voltage arrays + * for 2A and 2B are swapped. + * 05/27/04 1.21 bjg Test counter keywords. If bad, use + * cnt_rate=NEVENTS/exptime or 0 + * 06/14/04 1.22 wvd Add BRIT_OBJ keyword to header. + * Rewrite fill_hv_array() to make it more + * robust. + * If HV values in header are good, use + * them to determine whether HV + * arrays for 2A and 2B are swapped. + * 10/17/04 bjg Replaced lots of + * "while (val +#include +#include +#include +#include +#include +#include "calfuse.h" +#include "sgp4.h" + + +struct key { + char keyword[FLEN_KEYWORD]; + float value; + }; + + +/*********************************************************************** + * + * procedure to get the GTI values from the input file + * + **********************************************************************/ + +int cf_get_gti(fitsfile *infits, double **gtis, double **gtie) { + + char instmode[FLEN_VALUE]; + int ngti, status=0; + float ttperiod; + + cf_verbose(3, "Entering cf_get_gti."); + + /* Read header keywords. */ + FITS_movabs_hdu(infits, 1, NULL, &status); + FITS_read_key(infits, TSTRING, "INSTMODE", instmode, NULL, &status); + FITS_read_key(infits, TFLOAT, "TTPERIOD", &ttperiod, NULL, &status); + if (ttperiod < 1E-6) ttperiod = 1.0; + + /* Move to the second extension and read the GTI's */ + FITS_movabs_hdu(infits, 3, NULL, &status); + ngti = cf_read_col(infits, TDOUBLE, "START", (void **) gtis); + ngti = cf_read_col(infits, TDOUBLE, "STOP", (void **) gtie); + + /* + * If gtis[0] = 0.000, then it is probably wrong. + * Set it equal to the fractional part of gtie[0]. + * In HIST mode, gtis[0] is always 0., so don't change it. + * This trick only works if TTPERIOD = 1.0. + */ + if (**gtis < FRAME_TOLERANCE && !strncmp(instmode, "TTAG", 4) + && fabs(ttperiod - 1.0) < FRAME_TOLERANCE) + **gtis = fmod(*gtie[0], 1.); + + cf_verbose(3, "Exiting cf_get_gti."); + return ngti; +} + +/*********************************************************************** + * + * procedure to get a list of geocoronal line locations from the AIRG + * calibration file + * + ***********************************************************************/ + +int +cf_get_geocorona(fitsfile *outfits, short **gxmin, short **gxmax, + short **gymin, short **gymax) +{ + char airgfile[FLEN_VALUE]; + int nrow, status=0; + fitsfile *airgfits; + + cf_verbose(3, "Entering cf_get_geocorona."); + + /* Read the name of the AIRG calibration file to use and open it */ + FITS_movabs_hdu(outfits, 1, NULL, &status); + FITS_read_key(outfits, TSTRING, "AIRG_CAL", airgfile, NULL, &status); + FITS_open_file(&airgfits, cf_cal_file(airgfile), READONLY, &status); + cf_verbose(3, "AIRG_CAL = %s", airgfile); + + /* + * Move to the first extension and read the positions of the + * geocoronal lines. + */ + FITS_movabs_hdu(airgfits, 2, NULL, &status); + nrow = cf_read_col(airgfits, TSHORT, "XMIN", (void **) gxmin); + nrow = cf_read_col(airgfits, TSHORT, "XMAX", (void **) gxmax); + nrow = cf_read_col(airgfits, TSHORT, "YMIN", (void **) gymin); + nrow = cf_read_col(airgfits, TSHORT, "YMAX", (void **) gymax); + + FITS_close_file(airgfits, &status); + cf_verbose(3, "Exiting cf_get_geocorona."); + return nrow; +} + + +/*********************************************************************** + * + * Procedures to compute times of most recent sunrise and sunset + * + ***********************************************************************/ + +SGP4 set_orbit_parms(fitsfile *); + +static double +sunrise(double mjd_start, SGP4 sgp4) +{ + double ang_sep, mjd, ptime, pos[3], vel[3]; + int i, dn_flag, dn_flag_last; + + dn_flag_last = 1; /* Assume that we begin in night. */ + + for (i=1; i>-8000; i--) { + ptime = (double) i; + mjd = mjd_start + ptime/86400.0; + + /* Get the state vector at time mjd */ + SGP4_getStateVector(sgp4, mjd, pos, vel); + SGP4_precess(pos, mjd, MJD2000); + + dn_flag = eclipse(pos, mjd, &ang_sep); /* 0=day, 1=night */ + + /* We're going backward in time, so sunrise is day into night. */ + if ((dn_flag == 1) && (dn_flag_last == 0)) break; + + dn_flag_last=dn_flag; + } + /* Historically, we've defined sunrise as the time increment after + the change, so add 1 to the number we've just calculated. */ + return (ptime + 1.); +} + + +static double +sunset(double mjd_start, SGP4 sgp4) +{ + double ang_sep, mjd, ptime, pos[3], vel[3]; + int i, dn_flag, dn_flag_last; + + dn_flag_last = 0; /* Assume that we begin in day. */ + + for (i=1; i>-8000; i--) { + ptime = (double) i; + mjd = mjd_start + ptime/86400.0; + + /* Get the state vector at time mjd */ + SGP4_getStateVector(sgp4, mjd, pos, vel); + SGP4_precess(pos, mjd, MJD2000); + + dn_flag = eclipse(pos, mjd, &ang_sep); /* 0=day, 1=night */ + + /* We're going backward in time, so sunset is night into day. */ + if ((dn_flag == 0) && (dn_flag_last == 1)) break; + + dn_flag_last = dn_flag; + } + /* Historically, we've defined sunset as the time increment after + the change, so add 1 to the number we've just calculated. */ + return (ptime + 1.); +} + +/****************************************************************************/ + +static void +fill_cntrate_array(long nsam_hk, long *hk_colval, double *time_hk, + long ntime, double *ttime, float *tarray) +{ + /* + * Procedure to fill an array with count rate housekeeping (HK) + * information, which is tabulated on an approx 16s time interval. + * + * nsam_hk = number of time samples in the housekeeping file + * hk_colval = array containing the housekeeping data + * time_hk = time from the exposure start for each second in the + * houskeeping file + * ntime = number of tabulated rows in the timeline array + * ttime = tabulated times in the timeline array + * tarray = timeline array to be filled, one sample per second + * starting at the start of the exposure + */ + + double hk_time, hk_time0; + long hk_cnt, hk_cnt0; + long i; /* index through timeline table */ + long k; /* index through housekeeping array */ + float cnt_rate = 0.; + + i = k = 0; /* Initialize indices */ + + /* + * Find the first non-negative value in the HK array + */ + while(k < nsam_hk && hk_colval[k] < 0) k++; + hk_cnt0 = hk_colval[k]; + hk_time0 = time_hk[k]; + + /* + * Begin big loop through the timeline table + */ + while (i < ntime) { + + /* Find the next non-negative value in the HK array. */ + k++; + while(k < nsam_hk && hk_colval[k] < 0) k++; + + /* If we've run out of HK entries, + fill in the rest of the timeline and quit. */ + if (k == nsam_hk) { + for ( ; i < ntime; i++) + tarray[i] = cnt_rate; + break; + } + + /* Otherwise, calculate the count rate between + times hk_time0 and hk_time. */ + hk_cnt = hk_colval[k]; + hk_time = time_hk[k]; + if (hk_cnt < hk_cnt0) + hk_cnt0 -= 16777216; + cnt_rate = (hk_cnt - hk_cnt0) / (hk_time - hk_time0); + + /* Populate the timeline table for times less than hk_time */ + while (i < ntime && cf_nlong(ttime[i]) < hk_time) + tarray[i++] = cnt_rate; + + /* Save latest HK time and count values. */ + hk_cnt0 = hk_cnt; + hk_time0 = hk_time; + } +} + + +/****************************************************************************/ + +static void +hv_from_header(fitsfile *outfits, char *det, char *side, long ntime, short *hv) +{ + + char hv_key[FLEN_CARD]; + int det_hv, hv_flag, status=0; + long i; + + fits_read_key(outfits, TINT, "HV_FLAG", &hv_flag, NULL, &status); + if (hv_flag == -1) { + det_hv = -999; + cf_verbose(1, "Bad HV keywords: setting HV to -999 in timeline table."); + } + else { + sprintf(hv_key, "DET%.1sHV%.1sH", det, side); + FITS_read_key(outfits, TINT, hv_key, &det_hv, NULL, &status); + } + if (hv_flag == 0) + cf_verbose(1, "Applying max voltage value to whole exposure."); + + cf_verbose(2, "High voltage = %d", det_hv); + for (i=0; i MAX_EXPTIME) { + TIMES_TOO_BIG = TRUE; + cf_if_warning("Exposure appears to span %0.0f ks. Truncating " + "timeline table.", (gtie[ngti-1] - gtis[0])/1000.); + } + + /* For HIST data, use the good-time intervals + to set the times in the timeline table. */ + if (!strncmp(instmode, "HIST", 4)) { + if (TIMES_TOO_BIG) { + ntime = 0; + for (i = 0; i < ngti; i++) + ntime += (long) (gtie[i] - gtis[i] + 1.5); + ptime = (double *) cf_malloc(sizeof(double) * ntime); + k = 0; + for (i = 0; i < ngti; i++) { + ntime = (long) (gtie[i] - gtis[i] + 1.5); + for (j = 0; j < ntime; j++) + ptime[k++] = gtis[i] + j; + } + ntime = k; + } + else { + ntime = (long) (gtie[ngti-1] - gtis[0] + 1.5); + ptime = (double *) cf_malloc(sizeof(double) * ntime); + for (i = 0; i < ntime; i++) + ptime[i] = gtis[0] + i; + } + } + + /* For TTAG data, use the photon times themselves. */ + else { + FITS_movabs_hdu(outfits, 2, NULL, &status); + nevents = cf_read_col(outfits, TFLOAT, "TIME", (void **) &time); + FITS_movabs_hdu(outfits, 1, NULL, &status); + n_bad_times = 0; + if (TIMES_TOO_BIG) { + last_time = 0; + i = nevents - 1; + while (time[i] > MAX_EXPTIME && i > 0) { + if (fabs(time[i] - last_time) > FRAME_TOLERANCE && n_bad_times < 1024) { + bad_times[n_bad_times++] = time[i]; + last_time = time[i]; + } + i--; + } + last_time = time[i] + 1.; + } + else + last_time = time[nevents-1] + 1.; + first_time = time[0]; + if (first_time < fmod(last_time, 1.)) + first_time = fmod(last_time, 1.) - 1.; + ntime = (long) (last_time - first_time + 1.5) + n_bad_times; + ptime = (double *) cf_malloc(sizeof(double) * ntime); + for (i = 0; i < ntime - n_bad_times; i++) + ptime[i] = first_time + i; + for (j = n_bad_times-1; i < ntime; i++, j--) + ptime[i] = bad_times[j]; + free(time); + } + cf_verbose(2, "Timeline table will contain %ld entries", ntime); + + /* Now generate all of the other arrays in the timeline table. */ + + lon = (double *) cf_calloc(ntime, sizeof(double)); + lat = (double *) cf_calloc(ntime, sizeof(double)); + limbang = (double *) cf_calloc(ntime, sizeof(double)); + vorb = (double *) cf_calloc(ntime, sizeof(double)); + tsunrise = (double *) cf_calloc(ntime, sizeof(double)); + tsunset = (double *) cf_calloc(ntime, sizeof(double)); + hv = (short *) cf_calloc(ntime, sizeof(short)); + lifcr = (float *) cf_calloc(ntime, sizeof(float)); + siccr = (float *) cf_calloc(ntime, sizeof(float)); + feccr = (float *) cf_calloc(ntime, sizeof(float)); + aiccr = (float *) cf_calloc(ntime, sizeof(float)); + bkgdcr = (float *) cf_calloc(ntime, sizeof(float)); + tflag = (unsigned char *) cf_calloc(ntime, sizeof(unsigned char)); + ycentl = (float *)cf_calloc(ntime, sizeof(float)); + ycents = (float *)cf_calloc(ntime, sizeof(float)); + + for (i=0; i limbang[i]) min_limb = limbang[i]; + + dn_flag = eclipse(pos, mjd, &ang_sep); /* 0=day, 1=night */ + + /* The first time through, and after any breaks in the timeline, + * compute times of most recent sunrise and sunset. */ + if ((i == 0) || (ptime[i] - ptime[i-1] - 1. > FRAME_TOLERANCE)) { + lastsunrise = ptime[i] + sunrise(mjd, sgp4); + lastsunset = ptime[i] + sunset(mjd, sgp4); + if (lastsunset > lastsunrise) dn_flag_last = 1 ; + else dn_flag_last = 0 ; + /* The first time through, compute the orbital period and + * write it to file header. */ + if (i == 0) { + mjd += (lastsunset - ptime[i] + 7000.)/86400.0; + period = 7000. + sunset(mjd, sgp4); + FITS_movabs_hdu(outfits, 1, NULL, &status); + FITS_update_key(outfits, TDOUBLE, "ORBPERID", &period, + "[s] Estimate of orbital period", &status); + } + } + /* If we've moved from day into night (or night into day), update + last sunrise/sunset times. */ + else if ((dn_flag == 0) && (dn_flag_last == 1)) + lastsunrise = ptime[i]; + else if ((dn_flag == 1) && (dn_flag_last == 0)) + lastsunset = ptime[i]; + + /* If this is daytime, set the day bit in TFLAG array. */ + if (dn_flag == 0) + tflag[i] |= TEMPORAL_DAY; + + lon[i] = geo_lon; + lat[i] = geo_lat; + vorb[i] = dx; + tsunrise[i] = ptime[i] - lastsunrise; + tsunset[i] = ptime[i] - lastsunset; + dn_flag_last = dn_flag; + } + + + /* + * Set the OPUS flag for times outside of the OPUS good-time + * intervals. Since there are no photons associated with the + * last second of a GTI, we make sure that it is flagged as + * bad. The first second of a good-time interval is good. + */ + j = 0; + for (i = 0; i < ntime; i++) { + while (j < ngti-1 && ptime[i] > gtie[j] - 0.5) + j++; + if (ptime[i] < gtis[j] - FRAME_TOLERANCE || + ptime[i] > gtie[j] + FRAME_TOLERANCE) + tflag[i] |= TEMPORAL_OPUS; + } + tflag[ntime-1] |= TEMPORAL_OPUS; + + /* + * No exposure should be longer than 55 ks. For HIST data, we'll just + * put up with it, but for TTAG data, we'll flag as bad any timeline + * table entries with absurd arrival times. We use the OPUS flag for + * this. What we're really saying is that we don't believe that the + * times associated with these photons are correct. + */ + for (i = 0; i < ntime; i++) + if (ptime[i] > MAX_EXPTIME) + tflag[i] |= TEMPORAL_OPUS; + + /* Fill housekeeping columns (count rates and high voltage). */ + FITS_movabs_hdu(outfits, 1, NULL, &status); + FITS_read_key(outfits, TSTRING, "HKEXISTS", hkexists, NULL, &status); + FITS_read_key(outfits, TSTRING, "DETECTOR", det, NULL, &status); + strncpy(side, det+1, 1); + + /* Populate MIN_LIMB keyword. */ + FITS_update_key(outfits, TDOUBLE, "MIN_LIMB", &min_limb, NULL, &status); + + /* If no housekeeping file exists, use header info to fill timeline. */ + loop: + if(!strncasecmp(hkexists, "N", 1)) { + char cntb_key[FLEN_CARD], cnte_key[FLEN_CARD]; + float cntb, cnte; + float cnt_rate, cnt_rate2, eng_time, exp_time; + double ctimeb, ctimee; + + cf_verbose(1, "No housekeeping file: filling timeline with header info."); + + /* move to the top level header to read the relevant keywords */ + FITS_movabs_hdu(outfits, 1, NULL, &status); + FITS_read_key(outfits, TFLOAT, "EXPTIME", &exp_time, NULL, &status); + FITS_read_key(outfits, TLONG, "NEVENTS", &nevts, NULL, &status); + + /* calculate engineering count rates */ + FITS_read_key(outfits, TDOUBLE, "CTIME_B", &ctimeb, NULL, &status); + FITS_read_key(outfits, TDOUBLE, "CTIME_E", &ctimee, NULL, &status); + eng_time = (ctimee-ctimeb) * 86400.; + + /* If eng_time is unreasonable, use the EXPTIME keyword */ + if (eng_time < 1) { + cf_if_warning("Engineering snapshot time less than or equal to zero."); + if ((eng_time = (mjd_end-mjd_start) * 86400.) < MAX_EXPTIME) + cf_if_warning("Estimating LiF, SiC, FEC and AIC count rates from EXPSTART and EXPEND."); + else { + eng_time = exp_time; + cf_if_warning("Estimating LiF, SiC, FEC and AIC count rates from EXPTIME."); + } + } + + /* SiC count rate timeline */ + sprintf(cntb_key, "C%.2sSIC_B", det); + FITS_read_key(outfits, TFLOAT, cntb_key, &cntb, NULL, &status); + sprintf(cnte_key, "C%.2sSIC_E", det); + FITS_read_key(outfits, TFLOAT, cnte_key, &cnte, NULL, &status); + + if ((cnte<0) || (cntb<0) || (cntb>cnte) || eng_time < 1. || + (cnte==0xEB90EB90) || (cnte==0xDEADDEAD) || + (cntb==0xEB90EB90) || (cntb==0xDEADDEAD) || + ((cnt_rate = (cnte - cntb) / eng_time) > 32000)) { + cf_if_warning("Bad SiC counter. SiC count rate will be set to zero."); + cnt_rate=0; + } + + cf_verbose(2, "SiC count rate = %d", (int) cnt_rate); + for (i=0; icnte) || eng_time < 1. || + (cnte==0xEB90EB90) || (cnte==0xDEADDEAD) || + (cntb==0xEB90EB90) || (cntb==0xDEADDEAD) || + ((cnt_rate = (cnte - cntb) / eng_time) > 32000)) { + cf_if_warning("Bad LiF counter. LiF count rate will be set to zero."); + cnt_rate=0; + } + + cf_verbose(2, "LiF count rate = %d", (int) cnt_rate); + for (i=0; icnte) || eng_time<1 || + (cnte==0xEB90EB90) || (cnte==0xDEADDEAD) || + (cntb==0xEB90EB90) || (cntb==0xDEADDEAD) || + ((cnt_rate = (cnte - cntb) / eng_time) > 64000) || + (cnt_rate < 0.8*nevts/exp_time)) { + cf_if_warning("Bad FEC counter"); + cf_if_warning("-- FEC count rate will be computed from NEVENTS and EXPTIME."); + cf_if_warning("-- Electronic deadtime correction will be underestimated."); + cf_if_warning("-- Y stretch will be underestimated."); + cnt_rate=nevts/exp_time; + } + + cf_verbose(2, "FEC count rate = %d", (int) cnt_rate); + for (i=0; icnte) || (eng_time<1) || + (cnte==0xEB90EB90) || (cnte==0xDEADDEAD) || + (cntb==0xEB90EB90) || (cntb==0xDEADDEAD) || + ((cnt_rate2 = (cnte - cntb)/ eng_time) > 64000)) || + ((!(strncasecmp(dets[i],det,2))) && (cnt_rate2<0.8*nevts/exp_time))) { + cf_if_warning("Bad AIC counter %s",dets[i]); + cf_if_warning("-- AIC count rate will be computed from NEVENTS and EXPTIME."); + cf_if_warning("-- IDS deadtime correction will be underestimated."); + cnt_rate=nevts/exp_time; + break; + } + cnt_rate+=cnt_rate2; + } + + cf_verbose(2, "AIC count rate = %d", (int) cnt_rate); + for (i=0; i vmax) + vmax = hk_colval[j]; + if (vmax <= 0) { + cf_if_warning("Data missing from housekeeping column. " + "Will treat file as missing."); + strncpy(hkexists, "N", 1); + status=0; + goto loop; } + else fill_cntrate_array(nsam_hk, hk_colval, time_hk, ntime, + ptime, siccr);} + else { + cf_if_warning("Data column missing from housekeeping file. " + "Will treat file as missing."); + strncpy(hkexists, "N", 1); + status=0; + goto loop; + } + cf_verbose(3,"SiC count-rate info read from housekeeping file."); + + /* LiF count rate timeline */ + sprintf(cntb_key, "I_DET%.1sCLIF%.1s", det, side); + FITS_get_colnum(hskpfits, TRUE, cntb_key, &hk_colnum, &status); + FITS_read_col(hskpfits, TLONG, hk_colnum, 1L, 1L, nsam_hk, &intnull, + hk_colval, &anynull, &status); + fill_cntrate_array(nsam_hk, hk_colval, time_hk, ntime, ptime, lifcr); + cf_verbose(3,"LiF count-rate info read from housekeeping file."); + + /* FEC count rate timeline */ + sprintf(cntb_key, "I_DET%.1sCFE%.1s", det, side); + FITS_get_colnum(hskpfits, TRUE, cntb_key, &hk_colnum, &status); + FITS_read_col(hskpfits, TLONG, hk_colnum, 1L, 1L, nsam_hk, &intnull, + hk_colval, &anynull, &status); + fill_cntrate_array(nsam_hk, hk_colval, time_hk, ntime, ptime, feccr); + cf_verbose(3,"FEC count-rate info read from housekeeping file."); + + /* AIC count rate timeline */ + aiccr2 = (float *) cf_calloc(ntime, sizeof(float)); + for (j=0;j -1) { + FITS_read_key(outfits, TINT, "DET2HVAH", &hdr_max2a, NULL, &status); + cf_verbose(3, " Max for 2A: header says %d, housekeeping file says %d.", + hdr_max2a, max2a); + if (abs(max2a - hdr_max2a) < 5) + cf_verbose(3, " HV arrays are OK. No need to swap."); + else + swap_HV = TRUE; + } + + /* + * If HV values in file header are bad, compare with expected + * values from VOLT_CAL file. + */ + else { + FITS_open_file(&voltfits, cf_cal_file(volt_cal), READONLY, &status); + n = 0L; + do { + n++; + sprintf(mjd_str, "MJD%ld", n); + FITS_read_key(voltfits, TDOUBLE, mjd_str, &mjd_volt, NULL, &status); + } while (mjd_start > mjd_volt); + n--; + + sprintf(full_str, "FULL%ld", n); + sprintf(saa_str, "SAA%ld", n); + FITS_read_key(voltfits, TINT, full_str, &full, NULL, &status); + FITS_read_key(voltfits, TINT, saa_str, &saa, NULL, &status); + FITS_close_file(voltfits, &status); + + /* If max for 2A matches expected full or SAA voltage -- and + we're on side 2A -- then we're OK. Otherwise, must swap. */ + cf_verbose(3, " 2A max = %d; Expected values for %s: full = %d, SAA = %d", + max2a, det, full, saa); + if ((abs(max2a - full) < 10 || abs(max2a - saa) < 5) && + (*side == 'A')) + cf_verbose(3, " HV arrays are OK. No need to swap."); + else swap_HV = TRUE; + } + } + + /* High voltage timeline */ + if (swap_HV) { + if (*side == 'A') sprintf(cntb_key, "I_DET2HVBIASBST"); + else sprintf(cntb_key, "I_DET2HVBIASAST"); + cf_verbose(3, " Swapping HV values for detectors 2A and 2B"); + } + else sprintf(cntb_key, "I_DET%.1sHVBIAS%.1sST", det, side); + FITS_get_colnum(hskpfits, TRUE, cntb_key, &hk_colnum, &status); + FITS_read_col(hskpfits, TLONG, hk_colnum, 1L, 1L, nsam_hk, &intnull, + hk_colval, &anynull, &status); + if (!hv_from_hskp(nsam_hk, hk_colval, time_hk, ntime, ptime, hv)) + cf_verbose(3,"HV info read from housekeeping file."); + else { + cf_verbose(1,"Housekeeping file corrupted: reading HV info from file header."); + hv_from_header(outfits, det, side, ntime, hv); + } + + free(time_hk); + free(hk_mjd); + free(hk_colval); + FITS_close_file(hskpfits, &status); + } + + /* + * Set TFORM, TSCALE, and TZERO values for output table. + */ + sprintf(fmt_byte, "%ldB", ntime); + sprintf(fmt_float, "%ldE", ntime); + sprintf(fmt_short, "%ldI", ntime); + + /* First set default values. */ + tform[0] = fmt_float; + tform[1] = fmt_byte; + for (i=2; i max_lif) { + lifcr[i] = 0.; + biglif = TRUE; + } + if (siccr[i] > max_sic) { + siccr[i] = 0.; + bigsic = TRUE; + } + if (feccr[i] > max_fec) { + feccr[i] = 0.; + bigfec = TRUE; + } + if (aiccr[i] > max_aic) { + aiccr[i] = 0.; + bigaic = TRUE; + } + } + + if (biglif) cf_if_warning("LIF_CNT_RATE out of bounds. Setting bad values to zero."); + if (bigsic) cf_if_warning("SIC_CNT_RATE out of bounds. Setting bad values to zero."); + if (bigfec) cf_if_warning("FEC_CNT_RATE out of bounds. Setting bad values to zero."); + if (bigaic) cf_if_warning("AIC_CNT_RATE out of bounds. Setting bad values to zero."); + + /* + * Write the timeline table to the output file. + */ + FITS_movabs_hdu(outfits, 3, NULL, &status); + FITS_create_hdu(outfits, &status); + + FITS_create_tbl(outfits, BINARY_TBL, 1L, tfields, ttype, + tform, tunit, extname, &status); + + /* Write TSCALE and TZERO entries to header */ + for (i=2; i