1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
/*****************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*
* Synopsis: cf_reflux input_file output_file aeff_cal1 [aeff_cal2]
*
* Description: Program applies a new flux calibration to an extracted
* spectral file. If two effective-area files are given,
* program interpolates between them based on the MDJ
* of the exposure.
* Program writes a history line to output file header.
*
* Arguments: input_file FUSE extracted spectral file (fcal.fit)
* output_file Copy of input file with modified flux cal.
* aeff_cal1 Effective-area curve
* aeff_cal2 Optional effective-area curve, needed only
* if you wish to interpolate.
*
* Returns: none
*
* History: 06/02/05 1.1 wvd Based on programs cf_uninterp and
* cf_extract_spectra.
*
****************************************************************************/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "calfuse.h"
static char CF_PRGM_ID[] = "cf_reflux";
static char CF_VER_NUM[] = "1.1";
int main(int argc, char *argv[])
{
char aper_act[FLEN_VALUE], datestr[FLEN_VALUE];
int status=0, hdutype;
int aperture, timeref;
long i, frow=1, felement=1, nout;
float exptime, wpc;
float *wave, *flux, *error, *weights, *bkgd;
unsigned char *channel;
fitsfile *infits, *outfits;
/* Enter a timestamp into the log. */
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
/* Initialize error checking */
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
if (argc < 4) {
printf("Usage: cf_reflux input_file output_file aeff_cal1 [aeff_cal2]\n");
return 0;
}
/* Copy input spectral file to output. Close input file. */
FITS_open_file(&infits, argv[1], READONLY, &status);
FITS_create_file(&outfits, argv[2], &status);
FITS_copy_hdu(infits, outfits, 0, &status);
FITS_movrel_hdu(infits, 1, &hdutype, &status);
FITS_copy_hdu(infits, outfits, 0, &status);
FITS_close_file(infits, &status);
/* Modify output file keywords. */
FITS_movabs_hdu(outfits, 1, &hdutype, &status);
FITS_update_key(outfits, TSTRING, "FILENAME", argv[2], NULL, &status);
FITS_update_key(outfits, TSTRING, "AEFF1CAL", argv[3], NULL, &status);
if (argc == 5)
FITS_update_key(outfits, TSTRING, "AEFF2CAL", argv[4], NULL, &status);
else
FITS_update_key(outfits, TSTRING, "AEFF2CAL", argv[3], NULL, &status);
/* Add HISTORY lines to the output file */
fits_get_system_time(datestr, &timeref, &status);
strcat(datestr, " File recalibrated using cf_reflux");
FITS_write_history(outfits, datestr, &status);
/* Through which aperture was this spectrum obtained? */
FITS_read_key(outfits, TSTRING, "APER_ACT", aper_act, NULL, &status);
if (!strcmp(aper_act, "HIRS_LIF")) aperture = 1;
else if (!strcmp(aper_act, "MDRS_LIF")) aperture = 2;
else if (!strcmp(aper_act, "LWRS_LIF")) aperture = 3;
else if (!strcmp(aper_act, "HIRS_SIC")) aperture = 5;
else if (!strcmp(aper_act, "MDRS_SIC")) aperture = 6;
else if (!strcmp(aper_act, "LWRS_SIC")) aperture = 7;
else {
printf("Unable to interpret header keyword APER_ACT = %s\n", aper_act);
FITS_close_file(outfits, &status);
return 0;
}
/* Read WAVE, FLUX, ERROR, WEIGHTS, and BKGD arrays */
FITS_movabs_hdu(outfits, 2, &hdutype, &status);
nout = cf_read_col(outfits, TFLOAT, "WAVE", (void **) &wave);
nout = cf_read_col(outfits, TFLOAT, "FLUX", (void **) &flux);
nout = cf_read_col(outfits, TFLOAT, "ERROR", (void **) &error);
nout = cf_read_col(outfits, TFLOAT, "WEIGHTS", (void **) &weights);
nout = cf_read_col(outfits, TFLOAT, "BKGD", (void **) &bkgd);
channel = (unsigned char *) cf_malloc(sizeof(unsigned char) * nout);
/* Compute inputs to flux-calibration routine */
for (i = 0; i < nout; i++) {
channel[i] = (unsigned char) aperture;
if (flux[i] != 0.) error[i] /= flux[i]; /* Relative error */
flux[i] = weights[i] - bkgd[i]; /* Units are counts */
}
/* Apply flux calibration */
FITS_movabs_hdu(outfits, 1, &hdutype, &status);
FITS_update_key(outfits, TSTRING, "FLUX_COR", "PERFORM", NULL, &status) ;
cf_convert_to_ergs(outfits, nout, flux, flux, channel, wave);
/* Convert flux and error arrays to units of erg/cm2/s/A. */
FITS_movabs_hdu(outfits, 1, &hdutype, &status);
FITS_read_key(outfits, TFLOAT, "EXPTIME", &exptime, NULL, &status);
FITS_read_key(outfits, TFLOAT, "WPC", &wpc, NULL, &status);
if (exptime < 1) {
printf ("EXPTIME = %g. Exiting.\n", exptime);
FITS_close_file(outfits, &status);
return 0;
}
for (i = 0; i < nout; i++) {
flux[i] /= exptime * wpc;
error[i] *= flux[i];
}
/* Write new flux and error arrays to the output file. */
FITS_movabs_hdu(outfits, 2, &hdutype, &status);
FITS_write_col(outfits, TFLOAT, 2, frow, felement, nout, flux, &status);
FITS_write_col(outfits, TFLOAT, 3, frow, felement, nout, error, &status);
/* Close output file. */
FITS_close_file(outfits, &status);
/* Free memory. */
free(wave);
free(flux);
free(error);
free(weights);
free(bkgd);
free(channel);
/* Enter a timestamp into the log. */
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
return 0;
}
|
