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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
|
/*****************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*****************************************************************************
*
* Synopsis: cf_remove_motions options intermediate_file
*
* Description: Associates a photon with an aperture and corrects for grating,
* mirror, spacecraft (jitter) and focal plane assembly (FPA)
* motion. The position of the Y centroid is also determined
* both as a function of time (before correction for s/c motion)
* and for the entire exposure (after motion correction).
* The resulting coordinates represent the data obtained by a
* motionless instrument.
*
* By default, all corrections are performed. Command line
* options allow one or more corrections to be omitted.
*
* Arguments: input_file Intermediate data file
*
* Calibration files:
*
* Returns: 0 on successful completion
*
* Calls: cf_find_spectra, cf_calculate_ycent_motion,
* cf_grating_motion, cf_fpa_position, cf_mirror_motion,
* cf_satellite_jitter, cf_calculate_y_centroid
*
* History: 09/30/02 peb 1.1 Begin work
* 12/04/02 wvd 1.5 Install Y centroid routines.
* Modify call to cf_satellite_jitter
* 12/12/02 wvd 1.6 Change infits -> header throughout
* 03/04/03 peb 1.10 Minor code changes to remove gcc -Wall
* messages: unused j in main.
* Added unistd.h, so getopt is portable
* 03/05/03 wvd 1.11 Added cf_target_count_rate
* 03/10/03 peb 1.12 Added verbose_level, unistd.h for
* getopt portability, and changed
* locflags to unsigned char *.
* 03/13/03 rdr 1.13 corrected small error in tabulating
* count rate_sic
* 04/17/03 wvd 1.17 Add cf_find_spectra for initial
* assignment of photons to
* target apertures. Add last_call to
* cf_identify_channel, last_call and
* weight to cf_calculate_y_centroid
* 04/21/03 wvd 1.18 Pass tsunrise to cf_grating_motion,
* tsunset to cf_mirror_motion,
* channel to cf_satellite_jitter.
* 06/02/03 wvd 1.20 Don't modify count-rate arrays
* for HIST data.
* 06/11/03 wvd 1.21 Pass datatype to cf_read_col and
* cf_write_col.
* 08/01/03 wvd 1.23 Add cf_error_init after calls to
* subroutines.
* 08/06/03 wvd 1.24 Change channel to unsigned char.
* Remove GTI's from cf_satellite_jitter.
* 08/25/03 wvd 1.25 Change coltype from string to int in
* cf_read_col and cf_write_col.
* 09/18/03 wvd 1.26 Pass locflags to cf_identify_channel
* 10/21/03 wvd 1.27 For HIST data, don't write blank
* LiF and SiC count-rate arrays to IDF.
* 10/27/03 wvd 1.28 Pass pha array to cf_find_spectra.
* Delete last_call from argument list
* of cf_calculate_y_centroid.
* 10/31/03 wvd 1.29 Treat channel as unsigned char
* throughout.
* 04/06/04 bjg 1.30 Remove unused variables
* 05/03/04 wvd 1.31 Improve documentation.
* Move first call of cf_identify_channel
* from cf_find_spectra to this routine.
* 06/02/04 wvd 1.32 Copy/move cf_set_photon_flags(),
* cf_set_good_time_intervals(), and
* cf_modify_hist_times() from
* cf_screen_photons. Run them only if
* cf_satellite_jitter rejects some time
* AND if they were already run before.
* Pass timeflags to cf_find_spectra()
* and cf_calculate_y_centroid().
* Get exp_jitr from cf_satellite_jitter.
* 07/21/04 wvd 1.33 Add include file <string.h>
* 03/10/05 wvd 1.34 Change use of variable pad. Now add
* 10 to user_pad on first call to
* cf_identify_channel, 0 on second call.
* 03/22/05 wvd 1.35 Change TIME_SUNRISE and TIME_SUNSET
* from floats to shorts.
* 05/20/05 wvd 1.36 Clean up i/o.
* 09/09/05 wvd 1.37 Add -j to list of allowed options.
* 11/29/05 wvd 1.38 Separate cf_satellite_jitter into two
* routines. Screening is performed
* earlier, so we need not copy flags.
* 08/30/06 wvd 1.39 Add -a to list of allowed options.
* 12/29/06 wvd 1.40 Pass weights array and initial LiF
* and SiC counts arrays to
* cf_target_count_rate. Make arrays
* rate_lif and rate_sic floats. If
* they are too large to store as shorts,
* store as unsigned shorts.
* 02/02/07 wvd 1.41 If LiF or SiC count rate is too large
* for an unsigned int, set to zero.
*
****************************************************************************/
#include <unistd.h>
#include <stdlib.h>
#include <string.h>
#include "calfuse.h"
static char CF_PRGM_ID[] = "cf_remove_motions";
static char CF_VER_NUM[] = "1.41";
int
main(int argc, char *argv[])
{
unsigned char *channel=NULL, *timeflags=NULL, *locflags=NULL,
*statflag=NULL;
int grating_motion=1, fpa_position=1, mirror_motion=1, sat_jitter=1;
int last_call, set_times=FALSE, set_gtis=FALSE, airglow_centroid=FALSE;
int user_pad=0, status=0, set_rate=FALSE, optc;
long i, nevents, ntimes;
float *time=NULL, *ttime=NULL, *xfarf=NULL, *yfarf=NULL, *x, *y;
float *weight=NULL, *ycent_lif, *ycent_sic;
float max_rate_lif=0, max_rate_sic=0, *rate_lif=NULL, *rate_sic=NULL;
float tscale=1, tzero=32768;
short *tsunrise=NULL, *tsunset=NULL;
GTI gti;
fitsfile *header;
char opts[] = "hagfmjp:v:";
char usage[] =
"Usage:\n"
" cf_remove_motions [-hagfmj] [-p pad] [-v level] idf_file\n";
char option[] =
"Options:\n"
" -h: this help message\n"
" -v: verbose level (default is 1; 0 is silent)\n"
" -a: use airglow lines to determine spectral centroid\n"
" -p: aperture padding in pixels (default is 10)\n"
" -g: no grating motion correction\n"
" -f: no fpa position correction\n"
" -m: no mirror motion correction\n"
" -j: no satellite jitter correction\n";
verbose_level = 1;
while ((optc = getopt(argc, argv, opts)) != -1) {
switch(optc) {
case 'h':
printf("%s\n%s", usage, option);
return 0;
case 'v':
verbose_level = atoi(optarg);
break;
case 'p':
user_pad = atoi(optarg);
break;
case 'g':
grating_motion = 0;
break;
case 'f':
fpa_position = 0;
break;
case 'a':
airglow_centroid = TRUE;
break;
case 'm':
mirror_motion = 0;
break;
case 'j':
sat_jitter = 0;
break;
}
}
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
if (argc <= optind) {
printf("%s", usage);
cf_if_error("Incorrect number of command-line arguments");
}
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
FITS_open_file(&header, argv[optind], READWRITE, &status);
FITS_movabs_hdu(header, 2, NULL, &status);
nevents = cf_read_col(header, TFLOAT, "TIME", (void **) &time);
nevents = cf_read_col(header, TFLOAT, "WEIGHT", (void **) &weight);
nevents = cf_read_col(header, TFLOAT, "XFARF", (void **) &xfarf);
nevents = cf_read_col(header, TFLOAT, "YFARF", (void **) &yfarf);
nevents = cf_read_col(header, TBYTE, "TIMEFLGS", (void **) &timeflags);
nevents = cf_read_col(header, TBYTE, "LOC_FLGS", (void **) &locflags);
channel = cf_calloc(nevents, sizeof(unsigned char));
FITS_movabs_hdu(header, 4, NULL, &status) ;
ntimes = cf_read_col(header, TFLOAT, "TIME", (void **) &ttime) ;
ntimes = cf_read_col(header, TBYTE, "STATUS_FLAGS", (void **) &statflag);
ntimes = cf_read_col(header, TSHORT, "TIME_SUNRISE", (void **) &tsunrise) ;
ntimes = cf_read_col(header, TSHORT, "TIME_SUNSET", (void **) &tsunset) ;
ntimes = cf_read_col(header, TFLOAT, "LIF_CNT_RATE", (void **) &rate_lif) ;
ntimes = cf_read_col(header, TFLOAT, "SIC_CNT_RATE", (void **) &rate_sic) ;
ycent_lif=(float *) cf_calloc(ntimes, sizeof(float)) ;
ycent_sic=(float *) cf_calloc(ntimes, sizeof(float)) ;
FITS_movabs_hdu(header, 1, NULL, &status);
x = xfarf;
y = yfarf;
/* Find all six spectra and compute centroids. */
cf_find_spectra(header, nevents, weight, x, y, channel, timeflags,
locflags, airglow_centroid);
/* Initial assignment of channel numbers to photons */
cf_identify_channel(header, nevents, xfarf, yfarf, channel, locflags,
user_pad+10, last_call=FALSE);
/* Track Y centroids of moving spectra. */
cf_calculate_ycent_motion(header, nevents, time, y, channel, locflags,
ntimes, ttime, ycent_lif, ycent_sic);
/* Correct for detector motion during exposure. */
if (grating_motion)
cf_grating_motion(header, nevents, time, x, y, channel,
ntimes, ttime, tsunrise);
/* Shift spectra in X to account for position of the FPA. */
if (fpa_position)
cf_fpa_position(header, nevents, x, channel);
/* Correct for mirror motion during exposure. */
if (mirror_motion)
cf_mirror_motion(header, nevents, time, x, y, channel,
ntimes, ttime, tsunset);
/* Correct for spacecraft motion during exposure. */
if (sat_jitter)
cf_satellite_jitter(header, nevents, time, x, y, channel,
ntimes, ttime, statflag);
/* Compute centroids of motion-corrected spectra. */
cf_calculate_y_centroid(header, nevents, weight, x, y, channel,
timeflags, locflags);
/* Final assignment of channel numbers to photons */
cf_identify_channel(header, nevents, x, y, channel, locflags,
user_pad, last_call=TRUE);
/* Compute count rates for TTAG target spectra. */
set_rate = !(cf_target_count_rate(header, nevents, time, weight,
channel, locflags, ntimes, ttime, rate_lif, rate_sic));
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
cf_verbose(3, "Writing X, Y, CHANNEL to IDF");
FITS_movabs_hdu(header, 2, NULL, &status);
if (set_times)
cf_write_col(header, TFLOAT, "TIME", (void *) time, nevents);
cf_write_col(header, TFLOAT, "X", (void *) x, nevents);
cf_write_col(header, TFLOAT, "Y", (void *) y, nevents);
cf_write_col(header, TBYTE, "CHANNEL", (void *) channel, nevents);
if (set_gtis) {
cf_verbose(3, "Writing good-time intervals to IDF");
FITS_movabs_hdu(header, 3, NULL, &status);
cf_write_col(header, TDOUBLE, "START", (void *) gti.start, gti.ntimes);
cf_write_col(header, TDOUBLE, "STOP", (void *) gti.stop, gti.ntimes);
free(gti.stop);
free(gti.start);
}
FITS_movabs_hdu(header, 4, NULL, &status);
cf_verbose(3, "Writing timeline information to IDF");
if (set_rate) {
for (i = 0; i < ntimes; i++) {
if (rate_lif[i] > 65535) rate_lif[i] = 0;
if (rate_sic[i] > 65535) rate_sic[i] = 0;
if (max_rate_lif < rate_lif[i]) max_rate_lif = rate_lif[i];
if (max_rate_sic < rate_sic[i]) max_rate_sic = rate_sic[i];
}
cf_verbose(2, "max_rate_lif = %.0f", max_rate_lif);
cf_verbose(2, "max_rate_sic = %.0f", max_rate_sic);
if (max_rate_lif > 32767) {
FITS_update_key(header, TFLOAT, "TSCAL10", &tscale, NULL, &status);
FITS_update_key(header, TFLOAT, "TZERO10", &tzero, NULL, &status);
fits_set_tscale(header, 10, (double)tscale, (double)tzero, &status);
}
if (max_rate_sic > 32767) {
FITS_update_key(header, TFLOAT, "TSCAL11", &tscale, NULL, &status);
FITS_update_key(header, TFLOAT, "TZERO11", &tzero, NULL, &status);
fits_set_tscale(header, 11, (double)tscale, (double)tzero, &status);
}
FITS_write_col(header, TFLOAT, 10, 1, 1, ntimes, rate_lif, &status);
FITS_write_col(header, TFLOAT, 11, 1, 1, ntimes, rate_sic, &status);
}
cf_write_col(header, TFLOAT, "YCENT_LIF", (void *) ycent_lif, ntimes);
cf_write_col(header, TFLOAT, "YCENT_SIC", (void *) ycent_sic, ntimes);
cf_verbose(3, "IDF updates complete");
FITS_close_file(header, &status);
free(time);
free(weight);
free(xfarf);
free(yfarf);
free(locflags);
free(timeflags);
free(channel);
free(ttime);
free(rate_lif);
free(rate_sic);
free(ycent_lif);
free(ycent_sic);
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Done processing");
return status;
}
|
