OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_read_fits_table.pro 4.96 KB
Edit Raw Normal View History
a1839067   Jean-Philippe Bernard   First commit 
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