OV-GSO-DC / dustem-wrapper_idl

Commit 762cbe4806e4d6984607dd254d012f7a5071d733

Authored by Jean-Philippe Bernard
1 parent 474c4b97
Exists in master

modifed to cope for changing structure of fits file

Showing 2 changed files with 470 additions and 229 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_read_fits_table.pro
  Diff comments   View file @762cbe4
1   -PRO dustem_read_fits_table,filename=filename,help=help,dustem_spectra_st=dustem_spectra_st
  1 +PRO dustem_read_fits_table,filename=filename,dustem_st=dustem_st,help=help
2 2  
3 3 ;+
4 4 ; NAME:
5 5 ; dustem_read_fits_table
6   -;
7 6 ; PURPOSE:
8 7 ; Reads a Dustem fits table into current DustEMWrap system variables
9   -;
10 8 ; CATEGORY:
11 9 ; DustEMWrap, High-Level, Distributed, User convenience
12   -;
13 10 ; CALLING SEQUENCE:
14 11 ; dustem_read_fits_table[,filename=][,/help]
15   -;
16 12 ; INPUTS:
17   -;
18 13 ; OPTIONAL INPUT PARAMETERS:
19 14 ; filename = File name to be read (default='./dustemwrap_results.fits')
20   -;
21 15 ; OUTPUTS:
22 16 ; None
23   -;
24 17 ; OPTIONAL OUTPUT PARAMETERS:
25   -; dustem_spectra_st = emission and extinction spectra (total and for each grain type)
26   -;
  18 +; dustem_st = emission and extinction dustem output structure (total and for each grain type)
27 19 ; ACCEPTED KEY-WORDS:
28 20 ; help = If set, print this help
29   -;
30 21 ; COMMON BLOCKS:
31 22 ; None
32   -;
33 23 ; SIDE EFFECTS:
34 24 ; DustEMWrap system variables are set
35   -;
36 25 ; RESTRICTIONS:
37 26 ; The DustEMWrap IDL code must be installed
38   -;
39 27 ; PROCEDURE:
40 28 ; None
41   -;
42 29 ; EXAMPLES
43 30 ;
44 31 ; MODIFICATION HISTORY:
... ... @@ -61,80 +48,161 @@ ENDIF ELSE BEGIN
61 48 ENDELSE
62 49 message,'Reading '+file,/info
63 50  
64   -unit_input_sed=1
65   -unit_predicted_sed=2
66   -unit_predicted_spectra=3
67   -unit_predicted_extinctions=4
  51 +toto=mrdfits(file,0,main_header)
68 52  
69   -;==== Write the observed SED in extension 1
70   -str_input_SED=mrdfits(file,unit_input_sed,header_input_sed)
71   -;==== add the predicted SED in extension 2
72   -str_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
73   -;==== add the predicted total spectrum in extension 3
74   -;str_predicted_spectrum_tot=mrdfits(file,3,header_predicted_spectrum_tot)
  53 +;stop
  54 +;get extention numbers from main header
  55 +current_ext=1
  56 +ext_names=['']
  57 +exts=[-1]
  58 +finished=0
  59 +WHILE not finished DO BEGIN
  60 + ext_name=sxpar(main_header,'EXT'+strtrim(current_ext,2),count=count)
  61 + IF count NE 0 THEN BEGIN
  62 + ext_names=[ext_names,ext_name]
  63 + exts=[exts,current_ext]
  64 + current_ext=current_ext+1
  65 + ENDIF ELSE BEGIN
  66 + finished=1
  67 + ENDELSE
  68 +ENDWHILE
  69 +ext_names=ext_names[1:*]
  70 +exts=exts[1:*]
  71 +Nextensions=n_elements(exts)
  72 +
  73 +message,'================================================',/info
  74 +message,'Reading '+file,/info
  75 +message,'File extensions detected are: ',/info
  76 +FOR i=0L,Nextensions-1 DO BEGIN
  77 + message,'extension '+strtrim(exts[i],2)+' : '+strtrim(ext_names[i],2),/info
  78 +ENDFOR
  79 +message,'================================================',/info
  80 +
  81 +;extensions_strs=['INPUT_EMISSION_SED','DUSTEM_PREDICTED_SED', $
  82 +; 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED', $
  83 +; 'DUSTEM_PREDICTED_EMISSION_SPECTRA','DUSTEM_PREDICTED_EXTINCTION_SPECTRA', $
  84 +; 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL','DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL']
  85 +
  86 +ind=where(ext_names EQ 'INPUT_EMISSION_SED',count)
  87 +IF count NE 0 THEN BEGIN
  88 + unit_input_emission_sed=exts[ind[0]]
  89 + print,unit_input_emission_sed
  90 +ENDIF ELSE BEGIN
  91 + unit_input_emission_sed=-1
  92 +ENDELSE
  93 +ind=where(ext_names EQ 'DUSTEM_PREDICTED_SED',count)
  94 +IF count NE 0 THEN BEGIN
  95 + unit_dustem_predicted_sed=exts[ind[0]]
  96 + print,unit_dustem_predicted_sed
  97 +ENDIF ELSE BEGIN
  98 + unit_dustem_predicted_sed=-1
  99 +ENDELSE
  100 +
  101 +ind=where(ext_names EQ 'INPUT_EXTINCTION_SED',count)
  102 +IF count NE 0 THEN BEGIN
  103 + unit_input_extinction_sed=exts[ind[0]]
  104 + print,unit_input_extinction_sed
  105 +ENDIF ELSE BEGIN
  106 + unit_input_extinction_sed=-1
  107 +ENDELSE
  108 +ind=where(ext_names EQ 'DUSTEM_PREDICTED_EXTINCTION_SED',count)
  109 +IF count NE 0 THEN BEGIN
  110 + unit_dustem_predicted_extinction_sed=exts[ind[0]]
  111 + print,unit_dustem_predicted_extinction_sed
  112 +ENDIF ELSE BEGIN
  113 + unit_dustem_predicted_extinction_sed=-1 ;That would be abnormal
  114 +ENDELSE
  115 +
  116 +ind=where(ext_names EQ 'DUSTEM_PREDICTED_EMISSION_SPECTRA',count)
  117 +IF count NE 0 THEN BEGIN
  118 + unit_dustem_predicted_emission_spectra=exts[ind[0]]
  119 + print,unit_dustem_predicted_emission_spectra
  120 +ENDIF ELSE BEGIN
  121 + unit_dustem_predicted_emission_spectra=-1 ;That would be abnormal
  122 +ENDELSE
  123 +ind=where(ext_names EQ 'DUSTEM_PREDICTED_EXTINCTION_SPECTRA',count)
  124 +IF count NE 0 THEN BEGIN
  125 + unit_dustem_predicted_extinction_spectra=exts[ind[0]]
  126 + print,unit_dustem_predicted_extinction_spectra
  127 +ENDIF ELSE BEGIN
  128 + unit_dustem_predicted_extinction_spectra=-1 ;That would be abnormal
  129 +ENDELSE
  130 +
  131 +ind=where(ext_names EQ 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL',count)
  132 +IF count NE 0 THEN BEGIN
  133 + unit_dustem_predicted_emission_spectra_pol=exts[ind[0]]
  134 + print,unit_dustem_predicted_emission_spectra_pol
  135 +ENDIF ELSE BEGIN
  136 + unit_dustem_predicted_emission_spectra_pol=-1
  137 +ENDELSE
  138 +ind=where(ext_names EQ 'DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL',count)
  139 +IF count NE 0 THEN BEGIN
  140 + unit_dustem_predicted_extinction_spectra_pol=exts[ind[0]]
  141 + print,unit_dustem_predicted_extinction_spectra_pol
  142 +ENDIF ELSE BEGIN
  143 + unit_dustem_predicted_extinction_spectra_pol=-1
  144 +ENDELSE
  145 +
  146 +;stop
  147 +
  148 +; unit_input_sed=sxpar(main_header,'INPUT_EMISSION_SED')
  149 +; unit_predicted_sed=sxpar(main_header,'DUSTEM_PREDICTED_SED')
  150 +; unit_predicted_spectra=sxpar(main_header,'DUSTEM_PREDICTED_EMISSION')
  151 +; unit_predicted_extinctions=sxpar(main_header,'DUSTEM_PREDICTED_EXTINCTION')
  152 +; unit_predicted_spectra_pol=sxpar(main_header,'DUSTEM_PREDICTED_EMISSION_POL')
  153 +; unit_predicted_extinctions_pol=sxpar(main_header,'DUSTEM_PREDICTED_EXTINCTION_POL')
  154 +
  155 +;==== read the observed SED (if any)
  156 +IF unit_input_emission_sed NE -1 THEN BEGIN
  157 + str_input_SED=mrdfits(file,unit_input_emission_sed,header_input_sed)
  158 +ENDIF ELSE BEGIN
  159 +ENDELSE
  160 +;==== read the predicted SED (if any)
  161 +IF unit_dustem_predicted_sed NE -1 THEN BEGIN
  162 + str_predicted_SED=mrdfits(file,unit_dustem_predicted_sed,header_predicted_SED)
  163 +ENDIF ELSE BEGIN
  164 +ENDELSE
  165 +IF unit_input_extinction_sed NE -1 THEN BEGIN
  166 + str_input_EXT=mrdfits(file,unit_input_extinction_sed,header_input_ext)
  167 +ENDIF ELSE BEGIN
  168 +ENDELSE
  169 +IF unit_dustem_predicted_extinction_sed NE -1 THEN BEGIN
  170 + str_predicted_EXT=mrdfits(file,unit_input_extinction_sed,header_input_ext)
  171 +ENDIF ELSE BEGIN
  172 +ENDELSE
75 173  
76   -used_model=strtrim(sxpar(header_predicted_SED,'MODEL'),2)
77   -used_pol=long(strtrim(sxpar(header_predicted_SED,'POL'),2))
  174 +used_model=strtrim(sxpar(main_header,'MODEL'),2)
  175 +used_pol=long(strtrim(sxpar(main_header,'POL'),2))
  176 +used_version=strtrim(sxpar(main_header,'WRAP_V'),2)
78 177 IF used_model NE 'DEFAULT' THEN BEGIN
79 178 dustem_init,model=used_model,pol=used_pol
80 179 ENDIF ELSE BEGIN
81 180 dustem_init,model=0
82 181 ENDELSE
83 182 !dustem_model=used_model
84   -
85   -;stop
86   -;==== This is to create (*!dustem_data).sed in case it does not exist already
87   -Nsed=n_elements(str_input_SED)
88   -sed=dustem_initialize_internal_sed(Nsed)
89   -(*!dustem_data).sed=ptr_new(sed)
90   -
91   -(*(*!dustem_data).sed).instru_names=dustem_filter2instru(str_input_SED.filter)
92   -(*(*!dustem_data).sed).filt_names=str_input_SED.filter
93   -(*(*!dustem_data).sed).wav=str_input_SED.wavelength
94   -(*(*!dustem_data).sed).values=str_input_SED.I
95   -(*(*!dustem_data).sed).sigma=str_input_SED.varianceII
96   -
97   -IF !run_pol NE 0 THEN BEGIN
98   - ;stop
99   - (*!dustem_data).qsed=ptr_new(sed)
100   - (*!dustem_data).used=ptr_new(sed)
101   - (*(*!dustem_data).qsed).instru_names=dustem_filter2instru(str_input_SED.filter)
102   - (*(*!dustem_data).qsed).filt_names=str_input_SED.filter
103   - (*(*!dustem_data).qsed).wav=str_input_SED.wavelength
104   - (*(*!dustem_data).qsed).values=str_input_SED.Q
105   - (*(*!dustem_data).qsed).sigma=str_input_SED.varianceQQ
106   - (*(*!dustem_data).used).instru_names=dustem_filter2instru(str_input_SED.filter)
107   - (*(*!dustem_data).used).filt_names=str_input_SED.filter
108   - (*(*!dustem_data).used).wav=str_input_SED.wavelength
109   - (*(*!dustem_data).used).values=str_input_SED.U
110   - (*(*!dustem_data).used).sigma=str_input_SED.varianceUU
  183 +IF !dustem_version.version NE used_version THEN BEGIN
  184 + message,'Caution: fits file was not created with the same Dustemwrap version as your current version',/continue
  185 + message,'Fits dustemwrap version is :'+used_version,/continue
  186 + message,'current dustemwrap version is :'+!dustem_version.version,/continue
  187 + stop
111 188 ENDIF
112 189  
113   -;predicted_I_sed=str_predicted_spectrum_tot.I
114   -
115   -best_fit_chi2=sxpar(header_predicted_SED,'CHI2')
116   -best_fit_rchi2=sxpar(header_predicted_SED,'RCHI2')
  190 +best_fit_chi2=sxpar(main_header,'CHI2')
  191 +best_fit_rchi2=sxpar(main_header,'RCHI2')
117 192  
118   -;===== Fill in header with parameter values
119   -;parameters_desc=*(*!dustem_fit).PARAM_DESCS
120   -;parameters_init_values=*(*!dustem_fit).PARAM_INIT_VALUES
121   -;best_parameters=*(*!dustem_fit).CURRENT_PARAM_VALUES
122   -;best_parameters_uncertainties=*(*!dustem_fit).CURRENT_PARAM_ERRORS
123   -;best_fit_chi2=(*!dustem_fit).CHI2
124   -;best_fit_rchi2=(*!dustem_fit).RCHI2
125   -
126   -Nparams=sxpar(header_predicted_SED,'NPARAMS')
  193 +;==== get parameter values
  194 +Nparams=sxpar(main_header,'NPARAMS')
127 195 parameters_desc=strarr(Nparams)
128 196 best_parameters=dblarr(Nparams)
129 197 best_parameters_init_values=dblarr(Nparams)
130 198 best_parameters_uncertainties=dblarr(Nparams)
131 199 parameters_func_values=dblarr(Nparams)
132 200 FOR i=0L,Nparams-1 DO BEGIN
133   - parameters_desc[i]=sxpar(header_predicted_SED,'PARN'+strtrim(i+1,2))
134   - best_parameters[i]=sxpar(header_predicted_SED,'PARV'+strtrim(i+1,2))
135   - best_parameters_uncertainties[i]=sxpar(header_predicted_SED,'PARU'+strtrim(i+1,2))
136   - best_parameters_init_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
137   - parameters_func_values[i]=sxpar(header_predicted_SED,'PARI'+strtrim(i+1,2))
  201 + parameters_desc[i]=sxpar(main_header,'PARN'+strtrim(i+1,2))
  202 + best_parameters[i]=sxpar(main_header,'PARV'+strtrim(i+1,2))
  203 + best_parameters_uncertainties[i]=sxpar(main_header,'PARU'+strtrim(i+1,2))
  204 + best_parameters_init_values[i]=sxpar(main_header,'PARI'+strtrim(i+1,2))
  205 + parameters_func_values[i]=sxpar(main_header,'PARI'+strtrim(i+1,2))
138 206 ENDFOR
139 207  
140 208 (*!dustem_fit).PARAM_DESCS=ptr_new(parameters_desc)
... ... @@ -145,18 +213,90 @@ ENDFOR
145 213 (*!dustem_fit).RCHI2=best_fit_rchi2
146 214 (*!dustem_fit).PARAM_FUNC=ptr_new(parameters_func_values)
147 215  
  216 +
  217 +;stop
  218 +;==== This is to create (*!dustem_data).sed in case it does not exist already
  219 +IF unit_input_emission_sed NE -1 THEN BEGIN
  220 + Nsed=n_elements(str_input_SED)
  221 + sed=dustem_initialize_internal_sed(Nsed)
  222 + (*!dustem_data).sed=ptr_new(sed)
  223 +
  224 + (*(*!dustem_data).sed).instru_names=dustem_filter2instru(str_input_SED.filter)
  225 + (*(*!dustem_data).sed).filt_names=str_input_SED.filter
  226 + (*(*!dustem_data).sed).wav=str_input_SED.wavelength
  227 + (*(*!dustem_data).sed).values=str_input_SED.I
  228 + (*(*!dustem_data).sed).sigma=str_input_SED.varianceII
  229 +
  230 + IF used_pol THEN BEGIN
  231 + ;stop
  232 + (*!dustem_data).qsed=ptr_new(sed)
  233 + (*!dustem_data).used=ptr_new(sed)
  234 + (*(*!dustem_data).qsed).instru_names=dustem_filter2instru(str_input_SED.filter)
  235 + (*(*!dustem_data).qsed).filt_names=str_input_SED.filter
  236 + (*(*!dustem_data).qsed).wav=str_input_SED.wavelength
  237 + (*(*!dustem_data).qsed).values=str_input_SED.Q
  238 + (*(*!dustem_data).qsed).sigma=str_input_SED.varianceQQ
  239 + (*(*!dustem_data).used).instru_names=dustem_filter2instru(str_input_SED.filter)
  240 + (*(*!dustem_data).used).filt_names=str_input_SED.filter
  241 + (*(*!dustem_data).used).wav=str_input_SED.wavelength
  242 + (*(*!dustem_data).used).values=str_input_SED.U
  243 + (*(*!dustem_data).used).sigma=str_input_SED.varianceUU
  244 + ENDIF
  245 +ENDIF
  246 +
  247 +IF unit_input_extinction_sed NE -1 THEN BEGIN
  248 + Next=n_elements(str_input_EXT)
  249 + ext=dustem_initialize_internal_sed(Next)
  250 + (*!dustem_data).ext=ptr_new(ext)
  251 +
  252 + (*(*!dustem_data).ext).instru_names=dustem_filter2instru(str_input_EXT.filter)
  253 + (*(*!dustem_data).ext).filt_names=str_input_EXT.filter
  254 + (*(*!dustem_data).ext).wav=str_input_EXT.wavelength
  255 + (*(*!dustem_data).ext).values=str_input_EXT.I
  256 + (*(*!dustem_data).ext).sigma=str_input_EXT.varianceII
  257 +
  258 + IF used_pol THEN BEGIN
  259 + ;stop
  260 + (*!dustem_data).qext=ptr_new(ext)
  261 + (*!dustem_data).uext=ptr_new(ext)
  262 + (*(*!dustem_data).qext).instru_names=dustem_filter2instru(str_input_EXT.filter)
  263 + (*(*!dustem_data).qext).filt_names=str_input_EXT.filter
  264 + (*(*!dustem_data).qext).wav=str_input_EXT.wavelength
  265 + (*(*!dustem_data).qext).values=str_input_EXT.Q
  266 + (*(*!dustem_data).qext).sigma=str_input_EXT.varianceQQ
  267 + (*(*!dustem_data).uext).instru_names=dustem_filter2instru(str_input_EXT.filter)
  268 + (*(*!dustem_data).uext).filt_names=str_input_EXT.filter
  269 + (*(*!dustem_data).uext).wav=str_input_EXT.wavelength
  270 + (*(*!dustem_data).uext).values=str_input_EXT.U
  271 + (*(*!dustem_data).uext).sigma=str_input_EXT.varianceUU
  272 + ENDIF
  273 +ENDIF
  274 +
148 275 ;=== Below is only to get the same output form as dustem_compute_sed
149 276 ;Note: If the same model is used, the outputs should be the same as what is in the fits file ....
  277 +IF unit_input_emission_sed NE -1 THEN BEGIN
  278 + recomputed_dustem_predicted_emission_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,st=dustem_st)
  279 +ENDIF ELSE BEGIN
  280 + recomputed_dustem_predicted_ext_sed=dustem_compute_ext(*(*!dustem_fit).current_param_values,st=dustem_st)
  281 +ENDELSE
  282 +
150 283 ;stop
151   -dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,out_st=dustem_spectra_st)
  284 +dustem_predicted_em=mrdfits(file,unit_dustem_predicted_emission_spectra,header_predicted_emission)
  285 +dustem_predicted_ext= mrdfits(file,unit_dustem_predicted_extinction_spectra,header_predicted_extinction)
  286 +
  287 +;help,dustem_st.sed,dustem_predicted_em,/str
  288 +;help,dustem_st.ext,dustem_predicted_ext,/str
152 289  
153   -;==== get the predicted emission spectra from extension 4
154   -toto=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
155   -dustem_spectra_st.sed=toto
  290 +dustem_st.sed=dustem_predicted_em
  291 +dustem_st.ext=dustem_predicted_ext
156 292  
157   -;==== get the predicted extinction spectra from extension 5
158   -toto=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)
159   -dustem_spectra_st.ext=toto
  293 +IF used_pol THEN BEGIN
  294 + toto=mrdfits(file,unit_dustem_predicted_emission_spectra_pol,header_predicted_spectra_pol)
  295 + ;help,dustem_st.polsed,toto,/str
  296 + dustem_st.polsed=toto
  297 + toto=mrdfits(file,unit_dustem_predicted_extinction_spectra_pol,header_predicted_extinctions_pol)
  298 + dustem_st.polext=toto
  299 +ENDIF
160 300  
161 301 ;stop
162 302  
... ...
src/idl/dustem_write_fits_table.pro
  Diff comments   View file @762cbe4
... ... @@ -54,62 +54,107 @@ IF keyword_set(help) THEN BEGIN
54 54 goto,the_end
55 55 ENDIF
56 56  
  57 +used_model=!dustem_model
  58 +
  59 +emission_sed_present=ptr_valid((*!dustem_data).sed)
  60 +extinction_sed_present=ptr_valid((*!dustem_data).ext)
  61 +;dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st=st
  62 +;fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7 ;this is correct despite the way it is presented (1.e4/1e.7*1e.20=1.e17 which is the factor to convert from ergs to Mega Janskys)
  63 +;SED_spec = st.sed.em_tot * fact
  64 +
57 65 ;===== run dustem to get predicted spectrum and sed
58   -IF ptr_valid((*!dustem_data).sed) THEN $
59   - dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,out_st=dustem_spectra_st)
60   -
61   -IF !run_pol THEN BEGIN
62   -;This is needed to compute Q,U of both SED and spectra
63   - dustem_predicted_polsed=dustem_compute_stokes(*(*!dustem_fit).current_param_values, $
64   - Q_spec=dustem_predicted_Q_spec,U_spec=dustem_predicted_U_spec,PSI_spec=dustem_predicted_PSI_spec, $
65   - dustem_psi_em=dustem_psi_em)
66   - dustem_predicted_Qsed=dustem_predicted_polsed[0]
67   - dustem_predicted_Used=dustem_predicted_polsed[1]
68   -;Same in extinction, but does not work in polarization extinction was not fitted.
69   -; dustem_predicted_polext=dustem_compute_stokext(*(*!dustem_fit).current_param_values,bidon, $
70   -; dustem_predicted_Qext,dustem_predicted_Uext, $
71   -; dustem_predicted_Qext_spec,dustem_predicted_Uext_spec,dustem_predicted_PSI_ext, $
72   -; dustem_psi_ext,out_st=out_st)
  66 +IF emission_sed_present THEN BEGIN
  67 + dustem_predicted_sed=dustem_compute_sed(*(*!dustem_fit).current_param_values,st=dustem_st)
  68 + IF !run_pol THEN BEGIN
  69 + ;This is needed to compute Q,U of both SED and spectra
  70 + dustem_predicted_polsed=dustem_compute_stokes(*(*!dustem_fit).current_param_values, $
  71 + Q_spec=dustem_predicted_Q_spec,U_spec=dustem_predicted_U_spec,PSI_spec=dustem_predicted_PSI_spec, $
  72 + dustem_psi_em=dustem_psi_em)
  73 + dustem_predicted_Qsed=dustem_predicted_polsed[0]
  74 + dustem_predicted_Used=dustem_predicted_polsed[1]
  75 + ENDIF
  76 +
  77 + ;===== define structures containing results
  78 + one_str_input_SED={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)}
  79 + one_str_predicted_SED={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5)}
  80 + ;one_str_predicted_spectrum_tot={Wavelength:0.,I:0.,Q:0.,U:0.}
  81 + ;one_str_predicted_spectrum_grain1={Wavelength:0.,I:0.,Q:0.,U:0.}
  82 +
  83 + Nsed=n_elements((*(*!dustem_data).sed).filt_names)
  84 + str_input_SED=replicate(one_str_input_SED,Nsed)
  85 + str_input_SED.FILTER=(*(*!dustem_data).sed).filt_names
  86 + str_input_SED.Wavelength=(*(*!dustem_data).sed).wav
  87 + str_input_SED.I=(*(*!dustem_data).sed).values
  88 + str_input_SED.varianceII=(*(*!dustem_data).sed).sigma ;should it be squared ??
  89 +
  90 + N_predicted_sed=n_elements(dustem_predicted_sed)
  91 + str_predicted_SED=replicate(one_str_predicted_SED,N_predicted_sed)
  92 + ;str_predicted_EXT=replicate(one_str_predicted_SED,N_predicted_sed)
  93 + fully_undefined_sed=str_predicted_SED
  94 + str_predicted_SED.FILTER=(*(*!dustem_data).sed).filt_names ;taken from input SED
  95 + str_predicted_SED.Wavelength=(*(*!dustem_data).sed).wav ;taken from input SED
  96 + str_predicted_SED.I=dustem_predicted_sed
  97 + ;stop
  98 + IF !run_pol THEN BEGIN
  99 + str_predicted_SED.Q=dustem_predicted_Qsed
  100 + str_predicted_SED.U=dustem_predicted_Used
  101 + ;str_predicted_EXT.Q=dustem_predicted_Qext
  102 + ;str_predicted_EXT.U=dustem_predicted_Uext
  103 + ENDIF ELSE BEGIN
  104 + str_predicted_SED.Q=fully_undefined_sed.Q
  105 + str_predicted_SED.U=fully_undefined_sed.U
  106 + ;str_predicted_EXT.Q=fully_undefined_sed.Q
  107 + ;str_predicted_EXT.U=fully_undefined_sed.U
  108 + ENDELSE
73 109 ENDIF
74 110  
75   -used_model=!dustem_model
76 111  
77   -Ngrains=n_tags(dustem_spectra_st.sed)-2
78   -
79   -;===== define structures containing results
80   -one_str_input_SED={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)}
81   -one_str_predicted_SED={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5)}
82   -;one_str_predicted_spectrum_tot={Wavelength:0.,I:0.,Q:0.,U:0.}
83   -;one_str_predicted_spectrum_grain1={Wavelength:0.,I:0.,Q:0.,U:0.}
84   -
85   -Nsed=n_elements((*(*!dustem_data).sed).filt_names)
86   -str_input_SED=replicate(one_str_input_SED,Nsed)
87   -str_input_SED.FILTER=(*(*!dustem_data).sed).filt_names
88   -str_input_SED.Wavelength=(*(*!dustem_data).sed).wav
89   -str_input_SED.I=(*(*!dustem_data).sed).values
90   -str_input_SED.varianceII=(*(*!dustem_data).sed).sigma ;should it be squared ??
91   -
92   -N_predicted_sed=n_elements(dustem_predicted_sed)
93   -str_predicted_SED=replicate(one_str_predicted_SED,N_predicted_sed)
94   -str_predicted_EXT=replicate(one_str_predicted_SED,N_predicted_sed)
95   -fully_undefined_sed=str_predicted_SED
96   -str_predicted_SED.FILTER=(*(*!dustem_data).sed).filt_names ;taken from input SED
97   -str_predicted_SED.Wavelength=(*(*!dustem_data).sed).wav ;taken from input SED
98   -str_predicted_SED.I=dustem_predicted_sed
  112 +IF extinction_sed_present THEN BEGIN
  113 + ;stop
  114 + dustem_predicted_ext=dustem_compute_ext(*(*!dustem_fit).current_param_values,st=dustem_st)
  115 + IF !run_pol THEN BEGIN
  116 + ;This is needed to compute Q,U of both SED and spectra
  117 + dustem_predicted_polext=dustem_compute_stokext((*!dustem_fit).current_param_values,st=dustem_st)
  118 + ;dustem_predicted_polext=dustem_compute_polext(*(*!dustem_fit).current_param_values,st=dustem_st, $
  119 + ; POLEXT_spec=POLEXT_spec,SPEXT_spec=SPEXT_spec,dustem_fpolext=dustem_fpolext,dustem_ext=dustem_ext)
  120 + dustem_predicted_Qext=dustem_predicted_polext[0]
  121 + dustem_predicted_Uext=dustem_predicted_polext[1]
  122 + ENDIF
  123 +
  124 + ;===== define structures containing results
  125 + ;one_str_input_EXT=dustem_initialize_internal_sed(1)
  126 + one_str_input_EXT={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)}
  127 + one_str_predicted_EXT={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5)}
  128 + ;one_str_predicted_spectrum_tot={Wavelength:0.,I:0.,Q:0.,U:0.}
  129 + ;one_str_predicted_spectrum_grain1={Wavelength:0.,I:0.,Q:0.,U:0.}
  130 +
  131 + Next=n_elements((*(*!dustem_data).ext).filt_names)
  132 + str_input_EXT=replicate(one_str_input_EXT,Next)
  133 + str_input_EXT.FILTER=(*(*!dustem_data).ext).filt_names
  134 + str_input_EXT.Wavelength=(*(*!dustem_data).ext).wav
  135 + str_input_EXT.I=(*(*!dustem_data).ext).values
  136 + str_input_EXT.varianceII=(*(*!dustem_data).ext).sigma ;should it be squared ??
  137 +
  138 + N_predicted_ext=n_elements(dustem_predicted_polext[0])
  139 + str_predicted_EXT=replicate(one_str_predicted_EXT,N_predicted_EXT)
  140 + fully_undefined_ext=str_predicted_ext
  141 + str_predicted_EXT.FILTER=(*(*!dustem_data).ext).filt_names ;taken from input SED
  142 + str_predicted_EXT.Wavelength=(*(*!dustem_data).ext).wav ;taken from input SED
  143 + str_predicted_EXT.I=dustem_predicted_ext
  144 + IF !run_pol THEN BEGIN
  145 + str_predicted_EXT.Q=dustem_predicted_Qext
  146 + str_predicted_EXT.U=dustem_predicted_Uext
  147 + ENDIF ELSE BEGIN
  148 + str_predicted_EXT.Q=fully_undefined_ext.Q
  149 + str_predicted_EXT.U=fully_undefined_ext.U
  150 + ENDELSE
  151 +ENDIF
  152 +
  153 +Ngrains=n_tags(dustem_st.sed)-2
  154 +
99 155 ;stop
100   -IF !run_pol THEN BEGIN
101   - str_predicted_SED.Q=dustem_predicted_Qsed
102   - str_predicted_SED.U=dustem_predicted_Used
103   - ;str_predicted_EXT.Q=dustem_predicted_Qext
104   - ;str_predicted_EXT.U=dustem_predicted_Uext
105   -ENDIF ELSE BEGIN
106   - str_predicted_SED.Q=fully_undefined_sed.Q
107   - str_predicted_SED.U=fully_undefined_sed.U
108   - ;str_predicted_EXT.Q=fully_undefined_sed.Q
109   - ;str_predicted_EXT.U=fully_undefined_sed.U
110   -ENDELSE
111 156  
112   -N_predicted_spectra=n_elements(dustem_spectra_st.sed)
  157 +;N_predicted_spectra=n_elements(dustem_spectra_st.sed)
113 158  
114 159 ;===== Fill in header with parameter values
115 160 parameters_desc=*(*!dustem_fit).PARAM_DESCS
... ... @@ -122,11 +167,6 @@ parameters_func_values=*(*!dustem_fit).PARAM_FUNC
122 167  
123 168 Nparams=n_elements(parameters_desc)
124 169  
125   -;unit_input_sed=1
126   -;unit_predicted_sed=2
127   -;unit_predicted_spectra=3
128   -;unit_predicted_extinctions=4
129   -
130 170 ;stop
131 171 ;This is what will be saved
132 172 ;help,str_input_SED,/str
... ... @@ -137,143 +177,204 @@ Nparams=n_elements(parameters_desc)
137 177 file='./dustemwrap_tmp.fits'
138 178 ;==== Write the observed SED in extension 1
139 179 unit=1L
140   -mwrfits,str_input_SED,file,/create
141   -unit_input_sed=unit
142   -unit=unit+1
143   -;==== add the predicted SED in extension 2
144   -mwrfits,str_predicted_SED,file
145   -unit_predicted_sed=unit
146   -unit=unit+1
147   -;==== add the predicted per-grain spectrum in extension 3
148   -mwrfits,dustem_spectra_st.sed,file
  180 +;Yes, toto does not exist. This is just to write a dummy input to allow for a general header
  181 +mwrfits,toto,file,/create
  182 +IF emission_sed_present THEN BEGIN
  183 + ;==== add the input emission SED
  184 + mwrfits,str_input_SED,file
  185 + unit_input_sed=unit
  186 + unit=unit+1
  187 + ;==== add the predicted emission SED
  188 + mwrfits,str_predicted_SED,file
  189 + unit_predicted_sed=unit
  190 + unit=unit+1
  191 +ENDIF ELSE BEGIN
  192 + unit_input_sed=-1
  193 + unit_predicted_sed=-1
  194 +ENDELSE
  195 +IF extinction_sed_present THEN BEGIN
  196 + ;stop
  197 + mwrfits,str_input_EXT,file
  198 + unit_input_ext=unit
  199 + unit=unit+1
  200 + ;==== add the predicted emission SED
  201 + mwrfits,str_predicted_EXT,file
  202 + unit_predicted_ext=unit
  203 + unit=unit+1
  204 +ENDIF ELSE BEGIN
  205 + unit_input_ext=-1
  206 + unit_predicted_ext=-1
  207 +ENDELSE
  208 +;==== add the predicted per-grain emission spectrum
  209 +mwrfits,dustem_st.sed,file
149 210 unit_predicted_spectra=unit
150 211 unit=unit+1
151   -;==== add the predicted per-grain extinction in extension 4
152   -mwrfits,dustem_spectra_st.ext,file
  212 +;==== add the predicted per-grain extinction spectrum
  213 +mwrfits,dustem_st.ext,file
153 214 unit_predicted_extinctions=unit
  215 +unit=unit+1
  216 +IF !run_pol EQ 1 THEN BEGIN
  217 + mwrfits,dustem_st.polsed,file
  218 + unit_predicted_spectra_pol=unit
  219 + unit=unit+1
  220 + mwrfits,dustem_st.polext,file
  221 + unit_predicted_extinctions_pol=unit
  222 + unit=unit+1
  223 +ENDIF ELSE BEGIN
  224 + unit_predicted_spectra_pol=-1
  225 + unit_predicted_extinctions_pol=-1
  226 +ENDELSE
154 227  
  228 +extension_numbers=[unit_input_sed,unit_predicted_sed, $
  229 + unit_input_ext,unit_predicted_ext, $
  230 + unit_predicted_spectra,unit_predicted_extinctions, $
  231 + unit_predicted_spectra_pol,unit_predicted_extinctions_pol]
  232 +extensions_strs=['INPUT_EMISSION_SED','DUSTEM_PREDICTED_SED', $
  233 + 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED', $
  234 + 'DUSTEM_PREDICTED_EMISSION_SPECTRA','DUSTEM_PREDICTED_EXTINCTION_SPECTRA', $
  235 + 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL','DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL']
  236 +
  237 +ind=where(extension_numbers NE -1,count)
  238 +extension_numbers=extension_numbers[ind]
  239 +extensions_strs=extensions_strs[ind]
  240 +order=sort(extension_numbers)
  241 +extension_numbers=extension_numbers[order]
  242 +extensions_strs=extensions_strs[order]
  243 +Nextensions=n_elements(extension_numbers)
155 244 ;print,unit_input_sed,unit_predicted_sed,unit_predicted_spectra,unit_predicted_extinctions
156 245 ;stop
157 246  
158   -;===== read back extensions to get proper extension headers
  247 +;stop
  248 +
  249 +;===== read back fits file to get proper extension headers
159 250 file='./dustemwrap_tmp.fits'
160   -sst_input_sed=mrdfits(file,unit_input_sed,header_input_sed)
161   -sst_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
  251 +toto=mrdfits(file,0,hdr) ;this would be to get the main header
  252 +;stop
  253 +IF unit_input_sed NE -1 THEN sst_input_sed=mrdfits(file,unit_input_sed,header_input_sed)
  254 +IF unit_predicted_sed NE -1 THEN sst_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
  255 +IF unit_input_ext NE -1 THEN sst_input_ext=mrdfits(file,unit_input_ext,header_input_ext)
  256 +IF unit_predicted_ext NE -1 THEN sst_predicted_ext=mrdfits(file,unit_predicted_ext,header_predicted_ext)
162 257 ;sst_predicted_spectrum_tot=mrdfits(file,3,header_predicted_spectrum_tot)
163   -dustem_spectra_st.sed=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
164   -dustem_spectra_st.ext=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)
165   -
166   -;===== complement extension header for input SED
167   -sxaddpar,header_input_sed,'TITLE','OBSERVED INPUT SED TO DUSTEMWRAP',before='COMMENT'
168   -sxaddpar,header_input_sed,'TTYPE1','FILTER','Dustemwrap Filter name'
169   -sxaddpar,header_input_sed,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
170   -sxaddpar,header_input_sed,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
171   -sxaddpar,header_input_sed,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
172   -sxaddpar,header_input_sed,'TTYPE5','U','Observed Stokes U Intensity [??]'
173   -sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]'
174   -sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]'
175   -sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]'
  258 +dustem_st.sed=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
  259 +dustem_st.ext=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)
  260 +IF !run_pol EQ 1 THEN BEGIN
  261 + dustem_st.polsed=mrdfits(file,unit_predicted_spectra_pol,header_predicted_spectra)
  262 + dustem_st.polext=mrdfits(file,unit_predicted_extinctions_pol,header_predicted_extinctions)
  263 +ENDIF
  264 +;===== write content of extensions in main header
  265 +FOR i=0L,Nextensions-1 DO BEGIN
  266 + sxaddpar,hdr,'EXT'+strtrim(i+1,2),extensions_strs[i],'fits file extension'
  267 +ENDFOR
176 268  
177   -;===== complement extension header for predicted SED
178   -sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP'
179   -sxaddpar,header_predicted_SED,'TTYPE1','FILTER','Dustemwrap Filter name'
180   -sxaddpar,header_predicted_SED,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
181   -sxaddpar,header_predicted_SED,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
182   -sxaddpar,header_predicted_SED,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
183   -sxaddpar,header_predicted_SED,'TTYPE5','U','Observed Stokes U Intensity [??]'
184   -sxaddpar,header_predicted_SED,'CHI2',best_fit_chi2,'Dustemwrap best fit chi2 [-]'
185   -sxaddpar,header_predicted_SED,'RCHI2',best_fit_rchi2,'Dustemwrap best fit reduced chi2 [-]'
186   -sxaddpar,header_predicted_SED,'NGRAINS',Ngrains,'Number of grain species'
  269 +;====== write dustem wrap parameter best fit values
  270 +sxaddpar,hdr,'COMMENT','======= Dustemwrap dust model used ======',after='COMMENT'
  271 +sxaddpar,hdr,'MODEL',used_model,'Dust model used'
  272 +sxaddpar,hdr,'WRAP_V',!dustem_version.version,'Dustwrap version used'
  273 +sxaddpar,hdr,'COMMENT','===================================================',after='MODEL'
  274 +sxaddpar,hdr,'CHI2',best_fit_chi2,'Dustemwrap best fit chi2 [-]'
  275 +sxaddpar,hdr,'RCHI2',best_fit_rchi2,'Dustemwrap best fit reduced chi2 [-]'
  276 +sxaddpar,hdr,'NGRAINS',Ngrains,'Number of grain species'
187 277 ;====== write dustem wrap parameter names
188   -sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter names ======',after='RCHI2'
189   -sxaddpar,header_predicted_SED,'NPARAMS',Nparams,' Number of DustemWrap free parameters'
  278 +sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter names ======',after='RCHI2'
  279 +sxaddpar,hdr,'NPARAMS',Nparams,' Number of DustemWrap free parameters'
190 280 FOR i=0L,Nparams-1 DO BEGIN
191   - sxaddpar,header_predicted_SED,'PARN'+strtrim(i+1,2),parameters_desc[i],' DustemWrap parameter name'
  281 + sxaddpar,hdr,'PARN'+strtrim(i+1,2),parameters_desc[i],' DustemWrap parameter name'
192 282 ENDFOR
193   -sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARN'+strtrim(i,2)
194   -
195   -;====== write dustem wrap parameter best fit values
196   -sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap dust model used ======',after='COMMENT'
197   -sxaddpar,header_predicted_SED,'MODEL',used_model,'Dust model used'
198   -sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='MODEL'
199   -
200   -sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap dust model used in polarisation ? ======',after='COMMENT'
201   -sxaddpar,header_predicted_SED,'POL',!run_pol,'1 if pol was used, 0 otherwise'
202   -sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='POL'
203   -
  283 +sxaddpar,hdr,'COMMENT','===================================================',after='PARN'+strtrim(i,2)
  284 +sxaddpar,hdr,'COMMENT','======= Dustemwrap dust model used in polarisation ? ======',after='COMMENT'
  285 +sxaddpar,hdr,'POL',!run_pol,'1 if pol was used, 0 otherwise'
  286 +sxaddpar,hdr,'COMMENT','===================================================',after='POL'
204 287 ;====== write dustem wrap parameter best fit values
205   -sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter best fit values ======',after='COMMENT'
  288 +sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter best fit values ======',after='COMMENT'
206 289 FOR i=0L,Nparams-1 DO BEGIN
207   - sxaddpar,header_predicted_SED,'PARV'+strtrim(i+1,2),best_parameters[i],'Param. '+parameters_desc[i]
  290 + sxaddpar,hdr,'PARV'+strtrim(i+1,2),best_parameters[i],'Param. '+parameters_desc[i]
208 291 ENDFOR
209   -sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARV'+strtrim(i,2)
  292 +sxaddpar,hdr,'COMMENT','===================================================',after='PARV'+strtrim(i,2)
210 293  
211 294 ;====== write dustem wrap parameter best fit values uncertainties
212   -sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap best fit parameter uncertainties ======',after='COMMENT'
  295 +sxaddpar,hdr,'COMMENT','======= Dustemwrap best fit parameter uncertainties ======',after='COMMENT'
213 296 FOR i=0L,Nparams-1 DO BEGIN
214   - sxaddpar,header_predicted_SED,'PARU'+strtrim(i+1,2),best_parameters_uncertainties[i],'Param. '+parameters_desc[i]
  297 + sxaddpar,hdr,'PARU'+strtrim(i+1,2),best_parameters_uncertainties[i],'Param. '+parameters_desc[i]
215 298 ENDFOR
216   -sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARU'+strtrim(i,2)
  299 +sxaddpar,hdr,'COMMENT','===================================================',after='PARU'+strtrim(i,2)
217 300  
218 301 ;====== write dustem wrap parameter initial fit values
219   -sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter initial fit values ======',after='COMMENT'
  302 +sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter initial fit values ======',after='COMMENT'
220 303 FOR i=0L,Nparams-1 DO BEGIN
221   - sxaddpar,header_predicted_SED,'PARI'+strtrim(i+1,2),parameters_init_values[i],'Param. '+parameters_desc[i]
  304 + sxaddpar,hdr,'PARI'+strtrim(i+1,2),parameters_init_values[i],'Param. '+parameters_desc[i]
222 305 ENDFOR
223   -sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARI'+strtrim(i,2)
  306 +sxaddpar,hdr,'COMMENT','===================================================',after='PARI'+strtrim(i,2)
224 307  
225 308 ;====== write dustem wrap parameter func values (not sure exactly what this is)
226   -sxaddpar,header_predicted_SED,'COMMENT','======= Dustemwrap parameter func values ======',after='COMMENT'
  309 +sxaddpar,hdr,'COMMENT','======= Dustemwrap parameter func values ======',after='COMMENT'
227 310 FOR i=0L,Nparams-1 DO BEGIN
228   - sxaddpar,header_predicted_SED,'PARF'+strtrim(i+1,2),parameters_func_values[i],'Param. '+parameters_desc[i]
  311 + sxaddpar,hdr,'PARF'+strtrim(i+1,2),parameters_func_values[i],'Param. '+parameters_desc[i]
229 312 ENDFOR
230   -sxaddpar,header_predicted_SED,'COMMENT','===================================================',after='PARF'+strtrim(i,2)
  313 +sxaddpar,hdr,'COMMENT','===================================================',after='PARF'+strtrim(i,2)
  314 +
  315 +;===== complement extension header for input SED
  316 +sxaddpar,header_input_sed,'TITLE','OBSERVED INPUT SED TO DUSTEMWRAP',before='COMMENT'
  317 +sxaddpar,header_input_sed,'TTYPE1','FILTER','Dustemwrap Filter name'
  318 +sxaddpar,header_input_sed,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
  319 +sxaddpar,header_input_sed,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
  320 +sxaddpar,header_input_sed,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
  321 +sxaddpar,header_input_sed,'TTYPE5','U','Observed Stokes U Intensity [??]'
  322 +sxaddpar,header_input_sed,'TTYPE6','VARIANCEII','Observed Total Stokes I Variance [??]'
  323 +sxaddpar,header_input_sed,'TTYPE7','VARIANCEQQ','Observed Stokes Q Variance [??]'
  324 +sxaddpar,header_input_sed,'TTYPE8','VARIANCEUU','Observed Stokes U Variance [??]'
  325 +
  326 +;===== complement extension header for predicted SED
  327 +sxaddpar,header_predicted_SED,'TITLE','OUTPUT BEST FIT SED BY DUSTEMWRAP'
  328 +sxaddpar,header_predicted_SED,'TTYPE1','FILTER','Dustemwrap Filter name'
  329 +sxaddpar,header_predicted_SED,'TTYPE2','WAVELENGTH','Dustemwrap Filter reference wavelength [microns]'
  330 +sxaddpar,header_predicted_SED,'TTYPE3','I','Observed Total Stokes I Intensity [??]'
  331 +sxaddpar,header_predicted_SED,'TTYPE4','Q','Observed Stokes Q Intensity [??]'
  332 +sxaddpar,header_predicted_SED,'TTYPE5','U','Observed Stokes U Intensity [??]'
231 333  
232 334 ;===== complement extension header for predicted spectrum total
233   -;sxaddpar,header_predicted_spectrum_tot,'TITLE','OUTPUT BEST FIT TOTAL SPECTRUM BY DUSTEMWRAP'
234 335 sxaddpar,header_predicted_spectra,'TITLE','OUTPUT BEST FIT PER-GRAIN AND TOTAL EMISSION SPECTRA BY DUSTEMWRAP'
235 336 sxaddpar,header_predicted_extinctions,'TITLE','OUTPUT BEST FIT PER-GRAIN AND TOTAL EXTINCTION BY DUSTEMWRAP'
236 337  
237   -;FOR i=0L,Ngrains-1 DO BEGIN
238   -; hh=*headers_predicted_spectrum_grain[i]
239   -; sxaddpar,hh,'TITLE','OUTPUT BEST FIT SPECTRUM BY DUSTEMWRAP FOR GRAIN '+strtrim(i+1,2)
240   -; *headers_predicted_spectrum_grain[i]=hh
241   -;ENDFOR
242   -
243   -;stop
244   -;hprint,header_input_sed
245   -;hprint,header_predicted_SED
246   -
247   -;write final file
  338 +;============== write final fits file
248 339 IF keyword_set(filename) THEN BEGIN
249 340 file=filename
250 341 ENDIF ELSE BEGIN
251 342 file='./dustemwrap_results.fits'
252 343 ENDELSE
253 344  
254   -;==== Write the observed SED in extension 1
255   -mwrfits,str_input_SED,file,header_input_sed,/create
256   -unit=unit+1
257   -;==== add the predicted SED in extension 2
258   -mwrfits,str_predicted_SED,file,header_predicted_SED
259   -;==== add the predicted total spectrum in extension 3
260   -;mwrfits,str_predicted_spectrum_tot,file,header_predicted_spectrum_tot
261   -;==== add the predicted spectrum for each grain type in extension 4 - 4+Ngrains
262   -;FOR i=0L,Ngrains-1 DO BEGIN
263   -; mwrfits,*str_predicted_spectrum_grains[i],file,*headers_predicted_spectrum_grain[i]
264   -;ENDFOR
265   -;==== add the predicted per-grain spectrum in extension 3
266   -;stop
267   -mwrfits,dustem_spectra_st.sed,file,header_predicted_spectra
268   -;==== add the predicted per-grain extinction in extension 4
269   -mwrfits,dustem_spectra_st.ext,file,header_predicted_extinctions
270   -
  345 +;==== Write the observed SED (if any)
  346 +mwrfits,0,file,hdr,create=1
  347 +IF unit_input_sed NE -1 THEN BEGIN
  348 + mwrfits,str_input_SED,file,header_input_sed
  349 +ENDIF
  350 +;==== add the predicted SED (if any)
  351 +IF unit_input_sed NE -1 THEN BEGIN
  352 + mwrfits,str_predicted_SED,file,header_predicted_SED
  353 +ENDIF
  354 +IF unit_input_ext NE -1 THEN BEGIN
  355 + mwrfits,str_input_EXT,file,header_input_sed
  356 +ENDIF
  357 +;==== add the predicted SED (if any)
  358 +IF unit_predicted_ext NE -1 THEN BEGIN
  359 + mwrfits,str_predicted_EXT,file,header_predicted_SED
  360 +ENDIF
  361 +;==== add the predicted per-grain emission spectra
  362 +mwrfits,dustem_st.sed,file,header_predicted_spectra
  363 +;==== add the predicted per-grain extinction
  364 +mwrfits,dustem_st.ext,file,header_predicted_extinctions
  365 +IF !run_pol EQ 1 THEN BEGIN
  366 + mwrfits,dustem_st.polsed,file
  367 + mwrfits,dustem_st.polext,file
  368 +ENDIF
271 369  
  370 +message,'================================================',/info
272 371 message,'Wrote '+file,/info
273   -
  372 +message,'File extenssions are: ',/info
  373 +FOR i=0L,Nextensions-1 DO BEGIN
  374 + message,'extenssion '+strtrim(extension_numbers[i],2)+' : '+strtrim(extensions_strs[i],2),/info
  375 +ENDFOR
  376 +message,'================================================',/info
274 377 ;stop
275   -;hprint,header_input_sed
276   -;hprint,header_predicted_SED
277 378  
278 379  
279 380 the_end:
... ...