dustem_read_fits_table.pro
3.8 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
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)
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