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
|
/*****************************************************************************
* Johns Hopkins University
* Center For Astrophysical Sciences
* FUSE
*****************************************************************************
*
* Synopsis: cf_rebin_probability_array(infits, extended, aperture, nout,
* wave_out, pny, pycent, parray)
*
* Arguments: *fitsfile infits : Input IDF file
* int extended : TRUE if target is extended
* int aperture : Aperture ID for the analysis
* long nout : number of output wavelength bins
* float wave_out : output wavelength vector
* int pny : y dimension of probability array
* float pycent : y position of the peak of the
* probability array
* float parray : probability array
*
* History: 02/20/03 v1.1 rdr Started work
* 04/03/03 v1.2 rdr Normalize so that integral along y
* at each sample point is 1
* 06/02/03 v1.3 rdr Correct error in filling output array
* 09/30/03 v1.5 wvd Use cf_read_col to read WAVE_CAL.
* 03/22/04 v1.6 wvd Change pycent from int to float.
* 03/23/04 v1.7 wvd Shift weights array to account for
* FPA shifts.
* 06/14/04 v1.8 wvd Because the weights array is defined
* wrt the FARF, we need not correct
* for FPA shifts here.
* 05/17/05 v1.9 wvd Don't normalize regions where the
* total probability is less than 0.1.
* Set them to zero.
* 06/15/05 v1.10 wvd BUG FIX: Program always read the
* probability arrays for point-source
* targets from the WGTS_CAL file. Now
* it distinguishes between point-source
* and extended targets.
*
*****************************************************************************/
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include "calfuse.h"
int cf_rebin_probability_array(fitsfile *infits, int extended, int aperture,
long nout, float *wave_out, int *pny, float *pycent, float **parray)
{
char CF_PRGM_ID[] = "cf_rebin_probability_array";
char CF_VER_NUM[] = "1.10";
fitsfile *wgtsfits, *wavefits ;
int status=0, fpixel=1, nullval=0, anynull=0 ;
int hdu, nxwt, nywt, dndx, wtbin;
long npix, npix_out, nwave, ndx, ndx1, ndx2, i, j ;
float *wave, *wavet, *wgts_arr, *ynorm, *outarr ;
float total, wgt_cent ;
char wgtsfile[FLEN_VALUE]={'\0'}, wavefile[FLEN_VALUE]={'\0'} ;
char buffer[FLEN_CARD], det[FLEN_CARD] ;
/* Initialize error checking */
cf_error_init(CF_PRGM_ID, CF_VER_NUM, stderr);
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Started execution.");
/* Read calibration file names from the IDF file header. */
FITS_movabs_hdu(infits, 1, NULL, &status) ;
FITS_read_key(infits, TSTRING, "WGTS_CAL", wgtsfile, NULL, &status) ;
cf_verbose(3,"weights cal file = %s ",wgtsfile) ;
FITS_read_key(infits, TSTRING, "WAVE_CAL", wavefile, NULL, &status) ;
cf_verbose(3, "wavelength cal file = %s ",wavefile) ;
FITS_read_key(infits, TSTRING, "DETECTOR", det, NULL, &status) ;
cf_verbose(3,"Detector = %s ",det) ;
/* Open the weights file and read the array */
FITS_open_file(&wgtsfits, cf_cal_file(wgtsfile), READONLY, &status) ;
hdu = extended * 8 + aperture + 1;
cf_verbose(3, "WGTS_CAL: extended = %d, aperture = %d, hdu = %d",
extended, aperture, hdu);
FITS_movabs_hdu(wgtsfits, hdu, NULL, &status) ;
FITS_read_key(wgtsfits, TINT, "NAXIS1", &nxwt, buffer, &status) ;
FITS_read_key(wgtsfits, TINT, "NAXIS2", &nywt, buffer, &status) ;
FITS_read_key(wgtsfits, TFLOAT, "WGT_CENT", &wgt_cent, buffer, &status) ;
npix=nxwt * nywt ;
cf_verbose(3, "prob array: nx=%d, ny=%d, wgt_cent=%f ",nxwt, nywt, wgt_cent);
wgts_arr = (float *) calloc(npix, sizeof(float)) ;
FITS_read_img(wgtsfits, TFLOAT, fpixel, npix, &nullval, wgts_arr,
&anynull, &status) ;
FITS_close_file(wgtsfits, &status) ;
/* Read the wavelength array */
FITS_open_file(&wavefits, cf_cal_file(wavefile), READONLY, &status) ;
FITS_movabs_hdu(wavefits, aperture +1, NULL, &status) ;
nwave=cf_read_col(wavefits,TFLOAT,"WAVELENGTH", (void **) &wave);
FITS_close_file(wavefits, &status) ;
/* Determine the binning factor for the weights and generate a new
wavelength array that corresponds to the binning factor used. */
wtbin = NXMAX / nxwt ;
cf_verbose(3, "Binning factor for weights array = %d ", wtbin) ;
wavet = (float *) cf_calloc(nxwt, sizeof(float) ) ;
for (i=0; i<nxwt; i++) {
ndx=i * wtbin + (wtbin/2);
if (ndx > nwave-1) ndx = nwave-1 ;
if (ndx < 0) ndx = 0 ;
wavet[i] = wave[ndx] ;
}
/* Allocate space for the output array */
npix_out = nout * nywt ;
if ((aperture < 4 && strncmp(det,"1",1) == 1) ||
(aperture > 4 && strncmp(det,"2",1) == 1))
cf_verbose(3, "maximum tabulated wavelength = %g, "
"maximum of wave_out=%g", wavet[nxwt-1],wave_out[nout-1] );
else
cf_verbose(3, "maximum tabulated wavelength = %g, "
"maximum of wave_out=%g", wavet[0],wave_out[nout-1] );
outarr = (float *) cf_calloc(npix_out, sizeof(float)) ;
/* Fill in the output array - use the weight profile which is closest
to the individual wavelengths in the output wavelength array.
Remember : the Lif channels and the SiC channels have wavelengths
that run in opposite directions, e.g. for side 1 wavelengths
increase with index for Lif but decrease for SiC. The opposite is
true for side 2. */
ndx=0 ;
dndx = 1 ;
if ((aperture < 4 && strncmp(det,"2",1) == 0) ||
(aperture > 4 && strncmp(det,"1",1) == 0 ) ) {
ndx = nxwt - 1 ;
dndx = -1 ;
}
cf_verbose(3, "filling output: starting index = %d, increment = %d ",
ndx, dndx) ;
for (i=0; i<nout; i++) {
/* locate the wavelength corresponding to this entry */
while (wavet[ndx] < wave_out[i] && ndx <= nxwt-1 && ndx >= 0 )
ndx += dndx ;
/* populate the column of the array */
for (j=0; j<nywt ; j++) {
ndx1 = j*nout + i ;
if (ndx1 > npix_out-1) ndx1 = npix_out-1 ;
ndx2 = j*nxwt + ndx ;
if (ndx2 > npix-1) ndx2=npix-1 ;
outarr[ndx1] = wgts_arr[ndx2] ;
}
}
/* Adjust the probabilities so that the integral along y at every
wavelength position is 1 */
/* First sum the points along y (ynorm) at each point */
ynorm= (float *) cf_calloc(nout, sizeof(float)) ;
for (i=0; i<nout; i++) {
total=0. ;
for (j=0; j<nywt; j++) total += outarr[j*nout+i] ;
ynorm[i]=total ;
}
/* Normalize all points by dividing by the total */
for (i=0; i<nout; i++) {
if (ynorm[i] > 0.1)
for (j=0; j<nywt; j++) outarr[j*nout+i] /= ynorm[i] ;
else
for (j=0; j<nywt; j++) outarr[j*nout+i] = 0. ;
}
/* Redirect the output pointers. */
*pny = nywt ;
*pycent = wgt_cent ;
*parray = outarr ;
cf_timestamp(CF_PRGM_ID, CF_VER_NUM, "Finished execution.");
return status;
}
|
