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
|
/**************************************************************************
* Johns Hopkins University
* Center for Astrophysical Sciences
* FUSE
**************************************************************************
*
* Synopsis: cf_calculate_y_centroid(infits, nevents, weight, x, y,
* channel, timeflags, locflags)
*
* Description: Determines the y centroid of the target and airglow spectra
* in each aperture. The value written to the header depends
* on the quality flag for each channel in the IDF header.
* If QUALITY = HIGH, the target centroid is used; if MEDIUM,
* the airglow centroid is used; if LOW, the default centroid
* (from the CHID_CAL file) is used.
*
* Arguments: fitsfile infits Pointer to the input Intermediate Data File
* long nevents Number of photon events in the file
* float weight Scale factor for each photon
* float *x, *y X and Y positions of the photon
* unsigned char *channel Aperture associated with each photon
* unsigned char *timeflags Time flags - used to identify photons
* arriving during bursts, etc.
* unsigned char *locflags Location flags - used to identify photons
* in geocoronal line regions
*
* Calibration files required:
*
* Returns: 0 on success
*
* HISTORY: 10/16/02 v1.1 RDR started work
* 12/02/02 v1.2 RDR cleaned up variables
* 12/12/02 v1.3 wvd added include files
* 12/26/02 v1.4 RDR corrected error in calculating centroid
* 02/24/03 v1.5 peb changed include file to calfitsio.h
* 03/11/03 v1.6 wvd Changed screen to unsigned char
* 03/25/03 v1.7 rdr ignore pinhole aperture
* 04/03/03 v1.9 wvd Simplify calculation of centroid.
* Include airglow in centroid
* calculation of non-target apertures.
* Replace printf() with cf_verbose().
* 04/17/03 v1.10 wvd Fix bug in calculation of non-target
* centroids.
* 04/18/03 v1.11 wvd Call cf_proc_update only if final_call
* is TRUE. Scale YCENT by photon weight
* for use with HIST data. Changed
* screen to locflags.
* 05/20/03 v1.12 rdr Added call to cf_proc_check
* 05/22/03 v1.14 wvd cf_error_init to stderr
* 08/21/03 v1.16 wvd Change channel to unsigned char
* 09/02/03 v1.17 bjg Now read user-specified centroid
* information from PARM_CAL and
* expected centroids from CHID_CAL.
* 09/02/03 v1.18 wvd Tidy up code.
* 09/17/03 v1.19 wvd Return -1 if a target centroid
* is too far from expected value.
* 10/28/03 v1.20 wvd Rewrite program with new logic.
* Delete use of final_call.
* 11/12/03 v1.21 wvd If no photons in channel, don't
* crash, just use default centroid.
* If quality = LOW, ygeo = ycent.
* 04/09/04 v1.22 bjg Include string.h and stdio.h
* 06/04/04 v1.23 wvd Use photon time flags to exclude
* bursts, etc. from calculation.
* 08/18/04 v1.24 wvd Test all spectra to see whether
* centroid is out of bounds.
* 03/11/05 v1.25 wvd Tabulate centroids as doubles,
* then convert to floats.
* Write target centroid values
* to trailer file if verbose >= 2.
* 04/15/05 v1.26 wvd Clean up i/o.
* 04/26/05 v1.27 wvd Don't populate YGEO keywords, as
* they are wrong for point sources.
* Consider only events with good
* LOCATION flags.
* Use the largest aperture with a valid
* YGEO to determine offset from tabulated
* position. Use this shift to calculate
* ygeo for the smaller apertures.
* 05/18/05 v1.28 wvd Test both SPEX_LIF and SPEX_SIC.
* Don't override user-specified centroid,
* even if it's far from expected value.
* Fix bug in tabulation of expected
* Y centroids.
* 02/01/06 v1.29 wvd If ycent is too far from default value,
* use default value, not airglow centroid.
* Allow computation of yshift from
* airglow photons in target aperture.
* 11/25/08 1.30 wvd In HIST mode, there are no spectra in
* other apertures to cause confusion,
* so we need not require that the target
* centroid lie within X pixels of either
* the airglow or tabulated centroid.
*
**************************************************************************/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include "calfuse.h"
int cf_calculate_y_centroid(fitsfile *infits, long nevents, float *weight,
float *x, float *y, unsigned char *channel, unsigned char *timeflags,
unsigned char *locflags)
{
char CF_PRGM_ID[] = "cf_calculate_y_centroid";
char CF_VER_NUM[] = "1.30";
char aperture[FLEN_VALUE], instmode[FLEN_VALUE];
char ycentname[FLEN_KEYWORD], file_name[FLEN_VALUE];
char key[FLEN_KEYWORD], quality[8][FLEN_VALUE];
double dncent, dngeo, dycent, dygeo;
float ycent, ygeo, yshift, ycent_tab[8];
int chan_num, errflg=0, status=0, active_ap[2];
int got_shift, i, hdu, target_ap, user;
int spex_sic, spex_lif, emax_sic, emax_lif, emax, extended;
long n;
fitsfile *parm_cal, *chid_cal;
/* Initialize error checking */
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
/* Enter a time stamp into the log */
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
/* Check whether routine is appropriate for this data file. */
if ( (errflg = cf_proc_check(infits, CF_PRGM_ID)) ) return errflg;
/* Read the source aperture from the file header. */
FITS_read_key(infits, TSTRING, "APERTURE", aperture, NULL, &status);
extended=cf_source_aper(infits, active_ap) ;
/* Were data taken in HIST or TTAG mode? */
FITS_read_key(infits, TSTRING, "INSTMODE", instmode, NULL, &status);
/* Read centroid quality keywords from file header. */
for (i = 1; i < 8; i++) {
if (i == 4) continue;
sprintf(key, "YQUAL%1d", i);
FITS_read_key(infits, TSTRING, key, &quality[i][0], NULL, &status);
}
/* Read extraction parameters from PARM_CAL file. */
FITS_read_key(infits, TSTRING, "PARM_CAL", file_name, NULL, &status);
FITS_open_file(&parm_cal, cf_parm_file(file_name), READONLY, &status);
FITS_read_key(parm_cal,TINT,"SPEX_SIC",&spex_sic,NULL,&status);
FITS_read_key(parm_cal,TINT,"SPEX_LIF",&spex_lif,NULL,&status);
FITS_read_key(parm_cal,TINT,"EMAX_SIC",&emax_sic,NULL,&status);
FITS_read_key(parm_cal,TINT,"EMAX_LIF",&emax_lif,NULL,&status);
FITS_close_file(parm_cal,&status);
/* Read expected centroids from CHID_CAL file. Use extended-aperture
values except for apertures containing point-source targets. */
FITS_read_key(infits, TSTRING, "CHID_CAL", file_name, NULL, &status);
FITS_open_file(&chid_cal, cf_cal_file(file_name), READONLY, &status);
for (i = 1; i < 8; i++) {
if (i==4) continue;
if ((i == active_ap[0] || i == active_ap[1]) && !extended) hdu=i+1;
else hdu=i+9;
FITS_movabs_hdu(chid_cal, hdu, NULL, &status);
FITS_read_key(chid_cal,TFLOAT,"CENTROID",ycent_tab+i,NULL,&status);
}
FITS_close_file(chid_cal,&status);
/* Determine the centroids for each channel */
yshift = 0;
got_shift = FALSE;
for (chan_num = 7; chan_num > 0; chan_num--) {
if (chan_num == 4) {
yshift = 0;
got_shift = FALSE;
continue;
}
dngeo = dncent = 0.;
dygeo = dycent = 0.;
sprintf (ycentname, "YCENT%d", chan_num);
/* Is this a target aperture? */
if (chan_num == active_ap[0] || chan_num == active_ap[1])
target_ap = TRUE;
else
target_ap = FALSE;
/* If not, and we're in histogram mode, skip to the next channel. */
if (!target_ap && !strncmp(instmode, "HIST", 4)) continue;
/* Calculate separate target and geocoronal centroids. */
for (n = 0; n < nevents; n++) {
if (channel[n] == chan_num && !(timeflags[n] & ~TEMPORAL_DAY) &&
!(locflags[n] & ~LOCATION_AIR)) {
if(locflags[n] & LOCATION_AIR) { /* Airglow photon */
dygeo += y[n] * weight[n];
dngeo += weight[n];
}
else { /* Target photon */
dycent += y[n] * weight[n];
dncent += weight[n];
}
}
}
/* If we can't find any photons but expect to, set quality to LOW. */
if ((quality[chan_num][0] == 'H' && dncent < 1) ||
(quality[chan_num][0] == 'M' && dngeo < 1)) {
quality[chan_num][0] = 'L';
sprintf(key, "YQUAL%1d", chan_num);
FITS_update_key(infits, TSTRING, key, "LOW", NULL, &status) ;
cf_verbose(1, "No photon events in channel %d", chan_num);
}
/* Compute centroids of target and airglow spectra */
ycent = ygeo = 0;
if (dncent > 0) ycent = dycent / dncent;
if (dngeo > 0) ygeo = dygeo / dngeo;
/* If yshift is already known, use it to calculate ygeo.
* If not, and ygeo is valid for this aperture, compute yshift.
*/
if (got_shift)
ygeo = yshift + ycent_tab[chan_num];
else if (dngeo > 0) {
yshift = ygeo - ycent_tab[chan_num];
got_shift = TRUE;
}
/* Now decide which values to write to the header. */
/* If we are in the target aperture and the user has specified
* YCENT for this channel, use it.
*/
if (target_ap && ((chan_num<5 && spex_lif>=0 && spex_lif<1024) ||
(chan_num>4 && spex_sic>=0 && spex_sic<1024))) {
if (chan_num<5) ycent=spex_lif;
else ycent=spex_sic;
cf_verbose(1, "Channel %d: User-specified Y centroid = %g",
chan_num, ycent);
user = TRUE;
}
else
user = FALSE;
/* If the centroid quality is HIGH, use the calculated target
* value, regardless of whether we are in a target aperture.
* (This happens by default, so we don't have to do anything.)
*/
/* If the centroid quality is MED, use the airglow centroid. */
if (quality[chan_num][0] == 'M')
ycent = ygeo;
/* If the quality is LOW, use the default value. */
else if (quality[chan_num][0] == 'L')
ycent = ycent_tab[chan_num];
/*
* Test whether YCENT is too far from the expected value.
* If so, use the tabulated value and set the quality to LOW.
* Don't test HIST data. (wvd, 11/25/2008)
*/
if (chan_num<5) emax=emax_lif; else emax=emax_sic;
if (fabs(ycent-ycent_tab[chan_num]) > emax && !user && !strncmp(instmode, "TTAG", 4)) {
if(target_ap) {
cf_verbose(1, "Channel %d: computed centroid (%.1f) is too far "
"from expected value.", chan_num, ycent);
cf_verbose(1, "Channel %d: using default centroid of %.1f",
chan_num, ycent_tab[chan_num]);
}
else {
cf_verbose(2, "Channel %d: using default centroid of %.1f",
chan_num, ycent_tab[chan_num]);
}
ycent = ycent_tab[chan_num];
sprintf(key, "YQUAL%1d", chan_num);
FITS_update_key(infits, TSTRING, key, "LOW", NULL, &status) ;
}
/* Write centroid of target spectrum to the header */
FITS_update_key(infits, TFLOAT, ycentname, &ycent, NULL, &status) ;
cf_verbose(2, "For channel %d, Y centroid = %.1f", chan_num, ycent) ;
} /* End of loop over channels. */
cf_proc_update(infits, CF_PRGM_ID, "COMPLETE");
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished processing");
return status;
}
|
