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
|
pro fuse_sic_only, input
;+
; NAME:
; FUSE_EXTRACT
;
;*PURPOSE:
; To extract full-resolution spectra from the primary LiF and SiC
; apertures of an intermediate data file.
;
; THIS VERSION EXTRACTS ONLY A SIC SPECTRUM.
;
;*CATEGORY:
; INPUT/OUTPUT
;
;*CALLING SEQUENCE:
; FUSE_EXTRACT, input
;
;*INPUT:
; input: intermediate data file
;
;*OUTPUT:
; LiF and SiC spectral files in v2.x format.
;
;
;*PROCEDURES USED:
; FXADDPAR, FXHMAKE, FXBADDCOL, FXBCREATE, FXBWRITE, FXBFINISH, MRDFITS,
; SXPAR
;
;*HISTORY:
; Written by: V. Dixon March 2003
; 03/23/2003 V. Dixon Correct X values for Doppler shift.
; 03/24/2003 V. Dixon Fix bug in call to where().
; 06/27/2003 R. Robinson Update to new version of IDF
; - added fscale keyword to mrdfits call
; - divide flux by exposure time
; 06/08/2004 V. Dixon Use EXPTIME rather than RAWTIME.
; Correct calculation of COUNTS
; 07/01/2004 V. Dixon Don't read headers of extensions.
; Employ faster scheme for calculating
; Doppler shift.
; 07/02/2004 V. Dixon Don't compute dl per pixel. It suffers
; from truncation errors. Instead, use
; value from file header.
; 07/15/2004 V. Dixon Extract only SiC spectrum.
;-
;-----------------------------------------------------------------------------
;
; print calling sequence if no parameters supplied
;
if n_params(0) lt 1 then begin
print,'CALLING SEQUENCE: fuse_sic_only, input'
return
end
;
; Physical constants
;
C = 2.99792458e5 ; km/s
;
; Open input file and read keywords from header.
;
dummy = mrdfits(input, 0, hdr, status=status)
rootname = SXPAR(hdr, 'ROOTNAME')
vhelio = SXPAR(hdr, 'V_HELIO')
detector = SXPAR(hdr, 'DETECTOR')
aperture = SXPAR(hdr, 'APERTURE')
wavecal = SXPAR(hdr, 'WAVE_CAL')
exptime = SXPAR(hdr, 'EXPTIME')
ap_array = ['HIRS', 'MDRS', 'LWRS', 'PINH']
ap_code = ['3','2','4','1']
for k = 0, 3 do if (strcmp(aperture, ap_array[k], 1)) then ap = k
channel = ap+1
GET_DATE, FITS_date
;
; Read photon information in the first extension of the IDF.
;
a = mrdfits(input, 1, status=status, /fscale)
;
; Read timeline table from the third extension of the IDF.
;
b = mrdfits(input, 3, status=status, /fscale)
nseconds = n_elements(b.time)
if (nint(b.time[nseconds-1] - b.time[0]) gt nseconds) then $
print, 'WARNING: '+input+' timeline table not continuous.'
; for k = channel, channel+4, 4 do begin
k = channel + 4
print, 'Analyzing extension ', k
;
; Read corresponding wavelength array from WAVE_CAL file.
;
wcalfile=strtrim('/data1/fuse/calfuse/v3.0/calfiles/'+wavecal,2)
print, 'Reading wavelength data from ', wcalfile
w=mrdfits(wcalfile,k,w_hdr)
wave = w.wavelength
wpc = SXPAR(w_hdr, 'DISPAPIX')
;
; Select photons from aperture of interest.
;
g = where (a.channel eq k and (a.timeflgs and not 1B) eq 0 and $
(a.loc_flgs and 16B) eq 0, n)
print, 'Number of photons in spectrum = ', n
x = a.x[g]
;
; Correct X values for Doppler shift.
;
print, 'Applying Doppler shift'
t = nint(a.time[g] - b.time[0])
ovc = (b.ORBITAL_VEL + vhelio)/C
for i = 0L, n-1L do x[i] = x[i] + ovc[t[i]] * a.lambda[g[i]]/wpc
;
; Extract spectrum and populate output arrays.
;
print, 'Extracting spectrum'
x = nint(x)
flux = fltarr(16384)
error = fltarr(16384)
cntserr = fltarr(16384)
for i = 0L, n-1L do flux[x[i]] = flux[x[i]] + a.ergcm2[g[i]]
flux = flux / exptime / abs(wpc)
counts = histogram(x, min=0, max=16383)
g = where (counts gt 0)
cntserr[g] = sqrt(counts[g])
error[g] = flux[g]/cntserr[g]
quality = indgen(16384)
;
; Create primary HDU of output file
;
if (k lt 5) then ch = 'lif' else ch = 'sic'
filename = strcompress(rootname+strlowcase(detector)+ch+ap_code[ap]+'ttagfcal.fit',/remove_all)
FXADDPAR, hdr, 'DATE', FITS_date
FXADDPAR, hdr, 'FILENAME', filename
FXADDPAR, hdr, 'FILETYPE', 'CALIBRATED EXTRACTED SPECTRUM'
FXADDPAR, hdr, 'APER_ACT', strupcase(strcompress(ap_array[ap]+'_'+ch,/remove_all))
FXHMAKE, hdr, /EXTEND
FXWRITE, filename, hdr
FXBHMAKE, HEADER, 1, 'FUSE 1D Spectrum'
FXBADDCOL, iwave, HEADER, wave, 'WAVE', TUNIT='ANGSTROMS'
FXBADDCOL, iflux, HEADER, flux, 'FLUX', TUNIT='ERG/CM2/S/A'
FXBADDCOL, ierror, HEADER, error, 'ERROR', TUNIT='ERG/CM2/S/A'
FXBADDCOL, iquality, HEADER, quality, 'QUALITY', TUNIT='UNITLESS'
FXBADDCOL, icounts, HEADER, counts, 'COUNTS', TUNIT='COUNTS'
FXBADDCOL, icntserr, HEADER, cntserr, 'CNTSERR', TUNIT='COUNTS'
FXBCREATE, out, filename, HEADER
FXBWRITE, out, wave, iwave, 1
FXBWRITE, out, flux, iflux, 1
FXBWRITE, out, error, ierror, 1
FXBWRITE, out, quality, iquality, 1
FXBWRITE, out, counts, icounts, 1
FXBWRITE, out, cntserr, icntserr, 1
FXBFINISH, out
;endfor
close, /all
return
end
|
