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
|
/*****************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*****************************************************************************
*
* Synopsis: cf_rebin_background(infits, aperture, nout, wave_out,
* binx, biny, bnx, bny, bimage, barray)
*
* Arguments: *fitsfile infits : Input IDF file
* int aperture : Aperture ID for the analysis
* long nout : Length of wave_out array
* float wave_out : Wavelength vector
* int binx, biny : Binning factors of the background array
* int bnx, bny : X and Y dimension of the background array
* float bimage : 2-D background image (original)
* float barray : (OUTPUT) Binned background image.
*
* Description:
*
* The input background image is binned in X.
* We need an image binned in wavelength.
* To get it, we use the WAVE_CAL file (shifted to account
* for the FPA position) to derive a wavelength array for the
* background image. For each row in the background image (bimage),
* we build a second image (barray) scaled to the size of the output
* wavelength bins.
*
* History 02/21/03 v1.1 rdr Started work
* 04/01/03 v1.2 rdr Made aeff arrays double precision in
* the cf_read_fluxcal routine to be
* compatible with data in calfiles
* 05/07/03 v1.4 wvd Don't divide by EXPTIME.
* Read WPC from file header.
* Use cf_verbose.
* 06/02/03 v1.5 rdr Correct bug in filling output arrays
* 06/11/03 v1.7 wvd Pass datatype to cf_read_col.
* Change calfusettag.h to calfuse.h
* 08/22/03 v1.8 bjg Move get_extraction_limits to external
* subroutine. Free orphan arrays.
* Change coltype from char to int in
* cf_read_col.
* 08/28/03 v1.9 wvd Debug
* 09/29/03 v1.10 wvd Shift background model to match
* position of FPA.
* 03/16/04 v1.11 wvd Correct error in calculation of dw_ave
* and dweff.
* 04/05/04 v1.12 wvd Modify i/o.
* 04/26/04 v1.13 wvd Do not flux-calibrate 2-D background
* model. Do not extract 1-D background
* spectrum.
* 05/13/04 v1.14 wvd Correct value of CF_PRGM_ID.
* 07/21/04 v1.16 wvd Delete unused variables.
*
*****************************************************************************/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "calfuse.h"
int cf_rebin_background(fitsfile *infits, int aperture, long nout,
float *wave_out, int binx, int biny, int bnx, int bny,
float *bimage, float **barray) {
char CF_PRGM_ID[] = "cf_rebin_background";
char CF_VER_NUM[] = "1.16";
fitsfile *wavefits ;
int status=0;
int dndx, dpix, nyout;
long npix, npix_out, nwave, ndx, ndx1, ndx2;
long i, j, k;
float dwb, dw_out, mf, ctot, dw_ave;
float dpix_fpa;
float *wave, *wavet, *outarr;
char wavefile[FLEN_VALUE], det[FLEN_VALUE] ;
/* Initialize error checking */
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Started execution.");
/* Read wavelength spacing of the output array */
FITS_read_key(infits, TFLOAT, "WPC", &dw_out, NULL, &status) ;
cf_verbose(3, "Bin size of the output wavelength array = %5.3f A",
dw_out) ;
FITS_read_key(infits, TSTRING, "DETECTOR", det, NULL, &status) ;
cf_verbose(3,"Detector = %s, aperture = %d", det, aperture) ;
/* Read wavelength-calibration file */
FITS_read_key(infits, TSTRING, "WAVE_CAL", wavefile, NULL, &status) ;
cf_verbose(3, "Wavelength calibration file = %s",wavefile) ;
FITS_open_file(&wavefits, cf_cal_file(wavefile), READONLY, &status) ;
cf_verbose(3, "HDU to read = %d",aperture+1) ;
FITS_movabs_hdu(wavefits, aperture +1, NULL, &status) ;
nwave=cf_read_col(wavefits,TFLOAT,"WAVELENGTH",(void **) &wave);
cf_verbose(3, "Points read = %d, min & max: %4.2f, %4.2f",
nwave, wave[0], wave[nwave-1]);
/* Read pixel shifts due to FPA positions */
if (aperture < 5)
FITS_read_key(infits, TFLOAT, "FPADXLIF", &dpix_fpa, NULL, &status);
else
FITS_read_key(infits, TFLOAT, "FPADXSIC", &dpix_fpa, NULL, &status);
dpix = cf_nint(dpix_fpa);
/* To correct for the position of the FPA, each photon has been shifted
in X by dpix pixels. Rather than shifting the background to match,
we modify the wavelength array that maps the background onto the
output spectrum. */
cf_verbose(2, "Shifting aperture %d background by %d X pixels "
"to match the FPA shift.", aperture, cf_nint(dpix_fpa));
/* Generate a new wavelength array binned to match the background array */
cf_verbose(3, "X and Y binning factors for background array = %d, %d",
binx, biny) ;
wavet = (float *) cf_calloc(bnx, sizeof(float));
for (i = 0; i < bnx; i++) {
ndx=i * binx + (binx/2) + dpix;
if (ndx > nwave-1) ndx = nwave-1 ;
if (ndx < 0) ndx = 0;
wavet[i] = wave[ndx];
}
/* Error check: average wavelength spacing */
dw_ave = fabs(wavet[bnx-1]-wavet[0])/(bnx-1);
cf_verbose(3,"Average wavelength spacing of background array = %f", dw_ave) ;
if ((aperture < 4 && strncmp(det,"1",1) == 0) ||
(aperture > 4 && strncmp(det,"2",1) == 0)) {
cf_verbose(3, "Maximum tabulated wavelength = %.2f", wavet[bnx-1]);
cf_verbose(3, "Maximum of wave_out = %.2f", wave_out[nout-1] );
}
else {
cf_verbose(3, "Maximum tabulated wavelength = %.2f", wave_out[nout-1] );
cf_verbose(3, "Maximum of wave_out = %.2f", wave_out[nout-1] );
}
/* Define some sizes and allocate space for the output arrays */
nyout = bny * biny ; /* Y dimension of output 2-D bkgd array */
npix_out = nout * nyout ; /* Size of output 2-D bkgd array */
npix = bnx * bny ; /* Size of input 2-D bkgd array */
outarr = (float *) cf_calloc(npix_out, sizeof(float) ); /* 2-D output array */
/* Error check: total counts in the background image */
ctot=0.;
for (i=0; i<npix; i++) ctot+=bimage[i] ;
cf_verbose(3, "Number of counts in the input bkgd image = %f",ctot) ;
/* For each wavelength in the output array, find the nearest point
in the input background array. Scale its value by the relative
sizes of the background and output wavelength bins.
Remember: wavelengths for the Lif channels and the SiC channels run
in opposite directions. For side one, wavelengths increase
with index for Lif but decrease for SiC. The opposite is true
for side two. */
ndx=0 ;
dndx = 1 ;
if ((aperture < 4 && strncmp(det,"2",1) == 0) ||
(aperture > 4 && strncmp(det,"1",1) == 0 ) ) {
ndx = bnx - 1 ;
dndx = -1 ;
}
cf_verbose(3,"Filling output: starting ndx=%d, increment =%d ", ndx, dndx) ;
dwb = dw_ave;
for (i=0; i<nout; i++) {
/* Locate the element of the bkgd array corresponding to this wavelength */
while (wavet[ndx] < wave_out[i] && ndx < bnx && ndx >= 0) ndx += dndx ;
/* compute the relative
wavelength coverage of the background and output arrays */
if (ndx < bnx-2) dwb = fabs((double) (wavet[ndx+1]-wavet[ndx])) ;
mf = dw_out/dwb ;
/* correct the scale factor for binning in y */
mf /= biny ;
/* populate the column of the array */
for (j=0; j< nyout ; j++) {
k = j/biny ;
ndx1 = j*nout + i ;
if (ndx1 > npix_out-1) ndx1 = npix_out-1 ;
ndx2 = k*bnx + ndx ;
if (ndx2 > npix-1) ndx2=npix-1 ;
outarr[ndx1] = mf * bimage[ndx2] ;
}
}
/* Error checking: total counts in the background image */
ctot=0.;
for (i=0; i<npix_out; i++) ctot+=outarr[i] ;
cf_verbose(3, "Number of counts in the output bkgd image = %f", ctot) ;
/* Return output arrays */
*barray = outarr ;
/* Close wavelength calibration file */
FITS_close_file(wavefits, &status) ;
free(wave);
free(wavet);
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished execution.");
return status;
}
|
