dustem_read_fits_table.pro
4.96 KB
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
PRO dustem_read_fits_table,filename=filename,help=help,dustem_spectra_st=dustem_spectra_st
;+
; NAME:
; dustem_read_fits_table
; PURPOSE:
; Reads a Dustem fits table into current Dustem system variables
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; dustem_read_fits_table[,filename=][,/help]
; INPUTS:
; OPTIONAL INPUT PARAMETERS:
; filename = File name (default='/tmp/dustem_wrap_ouput_final.fits')
; OUTPUTS:
; None
; OPTIONAL OUTPUT PARAMETERS:
; dustem_spectra_st = emission and extinction spectra (total and for each grain type)
; ACCEPTED KEY-WORDS:
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; dustem system variables are set
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
;
; MODIFICATION HISTORY:
; Written by J.-Ph. Bernard
; see evolution details on the dustem cvs maintained at IRAP
; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems.
;-
; !DUSTEM_DAT = '/tmp/ramdisk8MB/'
; !DUSTEM_DATA = -> <Anonymous> Array[1]
; !DUSTEM_DO_CC = 1
; !DUSTEM_F90_EXEC = '$HOME/Soft_Librairies/dustem4.3_wk/src/dustem'
; !DUSTEM_FILTERS = <PtrHeapVar570>
; !DUSTEM_FIT = <PtrHeapVar1>
; !DUSTEM_F_HI = 1.00000
; !DUSTEM_IDL_CONTINUUM = 0.00000
; !DUSTEM_IDL_FREEFREE = 0.00000
; !DUSTEM_IDL_SYNCHROTRON = 0.00000
; !DUSTEM_INPUTS = <PtrHeapVar2>
; !DUSTEM_INSTRUMENT_DESCRIPTION = -> <Anonymous> Array[120]
; !DUSTEM_NEVER_DO_CC = 0
; !DUSTEM_PARAMS = <PtrHeapVar15>
; !DUSTEM_PARINFO = <PtrHeapVar943>
; !DUSTEM_PREVIOUS_CC = <PtrHeapVar3346>
; !DUSTEM_RES = '/tmp/ramdisk8MB/'
; !DUSTEM_SHOW_PLOT = 1
; !DUSTEM_SOFT_DIR = '$HOME/Soft_Librairies/dustem4.3_wk/'
; !DUSTEM_VERBOSE = 1
; !DUSTEM_WHICH = 'WEB3p8'
; !DUSTEM_WRAP_SOFT_DIR = '$HOME/Soft_Librairies/dustem-wrapper_idl/'
; !FIT_RCHI2_WEIGHT = -> <Anonymous> Array[1]
; !RUN_ANIS = 0.00000
; !RUN_CIRC = 0.00000
; !RUN_IONFRAC = 0.00000
; !RUN_LIN = 1.00000
; !RUN_POL = 1
; !RUN_RRF = 0.00000
; !RUN_TLS = 0.00000
; !RUN_UNIV = 0.00000
IF keyword_set(help) THEN BEGIN
doc_library,'dustem_read_fits_table'
goto,the_end
ENDIF
;write final file
IF keyword_set(filename) THEN BEGIN
file=filename
ENDIF ELSE BEGIN
file='/tmp/dustem_wrap_ouput_final.fits'
ENDELSE
message,'Reading '+file,/info
unit_input_sed=1
unit_predicted_sed=2
unit_predicted_spectra=3
unit_predicted_extinctions=4
;==== Write the observed SED in extension 1
str_input_SED=mrdfits(file,unit_input_sed,header_input_sed)
;==== add the predicted SED in extension 2
str_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
;==== add the predicted total spectrum in extension 3
;str_predicted_spectrum_tot=mrdfits(file,3,header_predicted_spectrum_tot)
Nsed=n_elements(str_input_SED)
(*!dustem_data.sed).instru_names=dustem_filter2instru(str_input_SED.filter)
(*!dustem_data.sed).filt_names=str_input_SED.filter
(*!dustem_data.sed).wav=str_input_SED.wavelength
(*!dustem_data.sed).values=str_input_SED.I
(*!dustem_data.sed).sigma=str_input_SED.varianceII
IF !run_pol NE 0 THEN BEGIN
ENDIF
;predicted_I_sed=str_predicted_spectrum_tot.I
best_fit_chi2=sxpar(header_predicted_SED,'CHI2')
best_fit_rchi2=sxpar(header_predicted_SED,'RCHI2')
;===== Fill in header with parameter values
parameters_desc=*(*!dustem_fit).PARAM_DESCS
parameters_init_values=*(*!dustem_fit).PARAM_INIT_VALUES
best_parameters=*(*!dustem_fit).CURRENT_PARAM_VALUES
best_parameters_uncertainties=*(*!dustem_fit).CURRENT_PARAM_ERRORS
best_fit_chi2=(*!dustem_fit).CHI2
best_fit_rchi2=(*!dustem_fit).RCHI2
Nparams=sxpar(header_predicted_SED,'NPARAMS')
parameters_desc=strarr(Nparams)
best_parameters=dblarr(Nparams)
best_parameters_init_values=dblarr(Nparams)
best_parameters_uncertainties=dblarr(Nparams)
FOR i=0L,Nparams-1 DO BEGIN
parameters_desc[i]=sxpar(header_predicted_SED,'PARN'+strtrim(i+1,2))
best_parameters[i]=sxpar(header_predicted_SED,'PARV'+strtrim(i+1,2))
best_parameters_uncertainties[i]=sxpar(header_predicted_SED,'PARU'+strtrim(i+1,2))
best_parameters_init_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
ENDFOR
*(*!dustem_fit).PARAM_DESCS=parameters_desc
*(*!dustem_fit).CURRENT_PARAM_VALUES=best_parameters
*(*!dustem_fit).CURRENT_PARAM_ERRORS=best_parameters_uncertainties
*(*!dustem_fit).PARAM_INIT_VALUES=best_parameters_init_values
;=== Below is only to get the same output form as dustem_compute_sed
;Note: If the same model is used, the outputs should be the same as what is in the fits file ....
dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,out_st=dustem_spectra_st)
;==== get the predicted emission spectra from extension 4
toto=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
dustem_spectra_st.sed=toto
;==== get the predicted extinction spectra from extension 5
toto=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)
dustem_spectra_st.ext=toto
;stop
the_end:
END