dustem_read_fits_table.pro
4.42 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
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 July 2022
; Evolution details on the DustEMWrap gitlab.
; See http://dustemwrap.irap.omp.eu/ for FAQ and help.
;-
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='./dustemwrap_results.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)
used_model=strtrim(sxpar(header_predicted_SED,'MODEL'),2)
IF used_model NE 'DEFAULT' THEN BEGIN
dustem_init,model=used_model
ENDIF ELSE BEGIN
dustem_init,model=0
ENDELSE
!dustem_model=used_model
;stop
;==== This is to create (*!dustem_data).sed in case it does not exist already
Nsed=n_elements(str_input_SED)
sed=dustem_initialize_internal_sed(Nsed)
(*!dustem_data).sed=ptr_new(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)
parameters_func_values=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))
parameters_func_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
ENDFOR
(*!dustem_fit).PARAM_DESCS=ptr_new(parameters_desc)
(*!dustem_fit).PARAM_INIT_VALUES=ptr_new(best_parameters_init_values)
(*!dustem_fit).CURRENT_PARAM_VALUES=ptr_new(best_parameters)
(*!dustem_fit).CURRENT_PARAM_ERRORS=ptr_new(best_parameters_uncertainties)
(*!dustem_fit).CHI2=best_fit_chi2
(*!dustem_fit).RCHI2=best_fit_rchi2
(*!dustem_fit).PARAM_FUNC=ptr_new(parameters_func_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