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
|
/*******************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*******************************************************************************
*
* Usage: cf_xcorr exposure_list exp_shift_list
*
* Description: Determine pixel shifts between exposures. The
* reference spectrum is the one having the longest exposure
* time. The input list (ASCII) is the set of exposures to process,
* the output list (ASCII) has 3 columns: exposure name, shift and
* sigma_shift. Sigma_shift = -1 corresponds to non-detected shift.
* Sigma_shift = 1 indicates a valid shift value.
* The reference exposure always has shift and sigma_shift of 0.
*
* History: 04/07/05 tc v1.0 First release
* 04/15/05 tc v1.1 Use cf_read_col
* 04/27/05 wvd v1.2 Use variance arrays properly.
* Use QUALITY array to identify
* bad data.
* Change spec array to type INT.
* Use chi2, not reduced chi2,
* to set error bars.
* 05/20/05 wvd v1.3 For continuum spectra, use
* single region between 1045 and
* 1070 A. For emission-line
* targets, use O VI lines.
* 06/03/05 wvd v1.4 Delete unused variables.
* 07/21/05 wvd v1.5 Adopt the sign convention for
* spectral shifts used by
* FUSE_REGISTER.
* 08/01/05 wvd v1.6 Store EXPNIGHT of reference
* spectrum in ref_expnight.
* 09/13/05 wvd v1.7 If called with no arguments,
* return calling info, not error.
* 04/27/06 wvd v1.8 Scale each spectrum to match
* mean of reference spectrum.
* 05/16/06 wvd v1.9 Use Lyman beta to align
* background exposures.
* Use O VI and C II to align WD's.
* 03/28/08 wvd v1.10 Write value of NORM to output.
* For HIST data only: if XCORRR
* fails, but NORM > 0.5, then
* set shift and sigma to 0.
* 08/15/08 wvd v1.11 If EXP_STAT = 2, treat as
* background observation.
*
******************************************************************************/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "calfuse.h"
/* Calfuse variables */
static char CF_PRGM_ID[] = "cf_xcorr";
static char CF_VER_NUM[] = "1.11";
/* Parameters */
#define MAXSHIFT 20
#define N_SHIFT 41
int main(int argc, char *argv[])
{
/* Variables */
char *spec_lis, *out_name, src_type[FLEN_VALUE];
char spec_name[FLEN_FILENAME], ref_spec_name[FLEN_FILENAME];
char program[FLEN_VALUE];
char instmode[FLEN_CARD] ;
double *ref_flux, *ref_error, *flux, *error;
double scale, sum, variance, norm, ref_tot, total;
float chi_square[N_SHIFT], chi_square_min, chi_square_max;
float exptime, exptime_max;
float *wave, w0, wpc, wmin, wdelta;
int ly_beta=FALSE, exp_stat, objclass, shift, sigma;
int ind_chi2_min, N, N_tmp, nchi, n_spec, ind_ref_spec=0;
long expnight, ref_expnight, i, k, n_pix, ref_start, start;
FILE *pt_file, *pt_file_out;
fitsfile *in_fits;
int status = 0, hdutype = 0;
/***********************************
** Enter a timestamp into the log **
***********************************/
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Begin Processing");
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
/***********************************************************
** Check for proper command-line usage and read arguments **
***********************************************************/
if (argc != 3) {
printf("Usage: cf_xcorr exposure_list exp_shift_list\n");
return(0);
}
spec_lis = argv[1];
out_name = argv[2];
/********************************************
** Open input list and create output file **
********************************************/
if ((pt_file = fopen(spec_lis, "r")) == NULL)
{
cf_if_error("Unable to open file %s", spec_lis);
}
pt_file_out = fopen(out_name, "w");
/********************************************************************
** Set the reference spectrum to the one which has maximum EXPTIME **
********************************************************************/
n_spec = 0;
exptime_max = -1;
while (fscanf(pt_file, "%80s", spec_name) != EOF)
{
FITS_open_file(&in_fits, spec_name, READONLY, &status);
FITS_read_key(in_fits, TFLOAT, "EXPTIME", &exptime, NULL, &status);
FITS_read_key(in_fits, TLONG, "EXPNIGHT", &expnight, NULL, &status);
FITS_close_file(in_fits, &status);
if (exptime > exptime_max)
{
exptime_max = exptime;
ind_ref_spec = n_spec;
strcpy(ref_spec_name, spec_name);
}
n_spec++;
}
if (n_spec == 0) cf_if_error("At least one spectrum is required");
if (n_spec == 1)
{
fprintf(pt_file_out, "%s %2d %2d %5.0f %5ld %c%c%c 1.0\n", ref_spec_name, 0, 0,
exptime_max, expnight, ref_spec_name[8], ref_spec_name[9], ref_spec_name[10]);
fclose(pt_file_out);
return(0);
}
printf("n_spec: %d, max exptime: %g, ind_ref: %d, ref_name: %s\n",
n_spec, exptime_max, ind_ref_spec, ref_spec_name);
/* Read source type and set XCORR limits accordingly. */
FITS_open_file(&in_fits, ref_spec_name, READONLY, &status);
FITS_read_key(in_fits, TINT, "EXP_STAT", &exp_stat, NULL, &status);
FITS_read_key(in_fits, TSTRING, "PRGRM_ID", &program, NULL, &status);
FITS_read_key(in_fits, TINT, "OBJCLASS", &objclass, NULL, &status);
FITS_read_key(in_fits, TSTRING, "SRC_TYPE", &src_type, NULL, &status);
FITS_read_key(in_fits, TSTRING, "INSTMODE", instmode, NULL, &status) ;
FITS_read_key(in_fits, TFLOAT, "W0", &w0, NULL, &status);
FITS_read_key(in_fits, TFLOAT, "WPC", &wpc, NULL, &status);
if (src_type[1] == 'E') { /* Emission-line source */
/* Use Lyman beta to align background observations. */
if (exp_stat == (int) TEMPORAL_LIMB ||
objclass == 1 || objclass == 7 || objclass == 90) {
printf("Assuming background observation.\n");
wmin = 1024.;
wdelta = 6.;
ly_beta = TRUE;
} else { /* Use O VI emission for everything else. */
printf("SRC_TYPE = %s. Emission-line target.\n", src_type);
wmin = 1030.;
wdelta = 9.;
}
}
else if (objclass == 17 || objclass == 29 || objclass == 37) {
/* Use O VI and C II to align white dwarf spectra. */
printf("OBJCLASS = %d. Assuming nearby white dwarf.\n", objclass);
wmin = 1030.;
wdelta = 9.;
}
else { /* Use 1045-1070 A region for all other continuum sources. */
printf("SRC_TYPE = %s. Assuming continuum target.\n", src_type);
wmin = 1045.;
wdelta = 25.;
}
start = cf_nlong((wmin - w0) / wpc);
ref_start = start + MAXSHIFT;
n_pix = cf_nlong(wdelta / wpc) - (N_SHIFT - 1);
/* Read wave, flux, and error arrays from reference exposure. */
FITS_read_key(in_fits, TLONG, "EXPNIGHT" , &ref_expnight, NULL, &status);
FITS_movabs_hdu(in_fits, 2, &hdutype, &status);
if (hdutype != BINARY_TBL) cf_if_error("FITS files must contain BINARY TABLE in HDU #2");
N = cf_read_col(in_fits, TFLOAT, "WAVE", (void **) &wave);
N = cf_read_col(in_fits, TDOUBLE, "FLUX", (void **) &ref_flux);
N = cf_read_col(in_fits, TDOUBLE, "ERROR", (void **) &ref_error);
FITS_close_file(in_fits, &status);
/* Compute total flux of reference spectrum in region of interest. */
ref_tot = 0.;
for (i = 0; i < n_pix; i++) ref_tot += ref_flux[start+i];
/*****************************************************************
** Compute the shift of each spectrum relative to the reference **
******************************************************************/
n_spec = 0;
rewind(pt_file);
while (fscanf(pt_file, "%80s", spec_name) != EOF)
{
/* Do nothing if the current spectrum is the reference. */
if (n_spec++ == ind_ref_spec) {
fprintf(pt_file_out, "%s %3d %2d %5.0f %5ld %c%c%c 1.0\n",
ref_spec_name, 0, 0, exptime_max, ref_expnight,
ref_spec_name[8], ref_spec_name[9], ref_spec_name[10]);
continue;
}
/* Read the next spectrum */
FITS_open_file(&in_fits, spec_name, READONLY, &status);
FITS_read_key(in_fits, TFLOAT, "EXPTIME", &exptime, NULL, &status);
FITS_read_key(in_fits, TLONG, "EXPNIGHT", &expnight, NULL, &status);
FITS_movabs_hdu(in_fits, 2, &hdutype, &status);
if (hdutype != BINARY_TBL) cf_if_error("FITS files must contain BINARY TABLE in HDU #2");
N_tmp = cf_read_col(in_fits, TDOUBLE, "FLUX", (void **) &flux);
N_tmp = cf_read_col(in_fits, TDOUBLE, "ERROR", (void **) &error);
if (N_tmp != N) cf_if_error("Tables must have the same number of elements");
FITS_close_file(in_fits, &status);
/* If EXPTIME < 1, move on to the next data file. */
if (exptime < 1) {
fprintf(pt_file_out, "%s %3d %2d %5.0f %5ld %c%c%c 0.0\n",
spec_name, 0, -1, exptime, expnight,
spec_name[8], spec_name[9], spec_name[10]);
continue;
}
/* Compute total flux of data spectrum in region of interest. */
norm = 1.;
if (!ly_beta) { /* Don't rescale if aligning on airglow. */
total = 0.;
for (i = 0; i < n_pix; i++) total += flux[start+i];
norm = ref_tot / total;
}
/* Compute chi-squared for shifts between +/- MAXSHIFT pixels. */
for (k = 0; k < N_SHIFT; k++) {
nchi = 0;
sum = 0;
for (i = 0; i < n_pix; i++) {
variance = ref_error[ref_start + i] * ref_error[ref_start + i] +
norm * norm * error[start + i + k] * error[start + i + k];
if (variance > 0) {
sum += (ref_flux[ref_start + i] - norm * flux[start + i + k]) *
(ref_flux[ref_start + i] - norm * flux[start + i + k]) / variance;
nchi++;
}
}
scale = (double) nchi / n_pix;
chi_square[k] = sum / scale;
}
/* Compute min and max values of chi-squared. */
ind_chi2_min = 0;
chi_square_min = 1E5;
chi_square_max = -1E5;
for (k = 0; k < N_SHIFT; k++) {
if (chi_square[k] < chi_square_min)
{
chi_square_min = chi_square[k];
ind_chi2_min = k;
}
else if (chi_square[k] > chi_square_max)
chi_square_max = chi_square[k];
}
shift = MAXSHIFT - ind_chi2_min;
/* If chi-squared changes by less than 20% over range of shifts, bail out.
* Otherwise, return shift corresponding to lowest value of chi-square. */
if (chi_square_max / chi_square_min < 1.2)
sigma = -1;
else
sigma = 1;
/* If we're in HIST mode and sigma = -1 and NORM < 2,
* then set shift = sigma = 0. */
if (!(strncasecmp(instmode, "HIST", 4)) && sigma == -1 && norm < 2.0)
shift = sigma = 0;
/* printf ("shift = %d\tchi_square_max / chi_square_min = %g\n",
ind_chi2_min - MAXSHIFT, chi_square_max / chi_square_min); */
fprintf(pt_file_out, "%s %3d %2d %5.0f %5ld %c%c%c %f\n",
spec_name, shift, sigma, exptime, expnight,
spec_name[8], spec_name[9], spec_name[10], norm);
free(flux);
free(error);
}
free(wave);
free(ref_flux);
free(ref_error);
fclose(pt_file);
fclose(pt_file_out);
return(0);
}
|
