Commit 17f02a5ae0a601f25b0d5361c52335e78b67c3bb

Authored by Jean-Philippe Bernard
1 parent 244eeac8
Exists in master

updated

src/idl/dustem_read_fits_table.pro
... ... @@ -90,8 +90,8 @@ FOR i=0L,Nextensions-1 DO BEGIN
90 90 ENDFOR
91 91 message,'================================================',/info
92 92  
93   -;extensions_strs=['INPUT_EMISSION_SED','DUSTEM_PREDICTED_SED', $
94   -; 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED', $
  93 +;extensions_strs=['INPUT_EMISSION_SED','DUSTEM_PREDICTED_SED','DUSTEM_SHOWN_SED', $
  94 +; 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED','DUSTEM_SHOWN_EXTINCTION_SED', $
95 95 ; 'DUSTEM_PREDICTED_EMISSION_SPECTRA','DUSTEM_PREDICTED_EXTINCTION_SPECTRA', $
96 96 ; 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL','DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL']
97 97  
... ... @@ -109,6 +109,13 @@ IF count NE 0 THEN BEGIN
109 109 ENDIF ELSE BEGIN
110 110 unit_dustem_predicted_sed=-1
111 111 ENDELSE
  112 +ind=where(ext_names EQ 'DUSTEM_SHOWN_SED',count)
  113 +IF count NE 0 THEN BEGIN
  114 + unit_dustem_shown_sed=exts[ind[0]]
  115 + ;print,unit_dustem_predicted_sed
  116 +ENDIF ELSE BEGIN
  117 + unit_dustem_shown_sed=-1
  118 +ENDELSE
112 119  
113 120 ind=where(ext_names EQ 'INPUT_EXTINCTION_SED',count)
114 121 IF count NE 0 THEN BEGIN
... ... @@ -124,6 +131,13 @@ IF count NE 0 THEN BEGIN
124 131 ENDIF ELSE BEGIN
125 132 unit_dustem_predicted_extinction_sed=-1 ;That would be abnormal
126 133 ENDELSE
  134 +ind=where(ext_names EQ 'DUSTEM_SHOWN_EXTINCTION_SED',count)
  135 +IF count NE 0 THEN BEGIN
  136 + unit_dustem_shown_extinction_sed=exts[ind[0]]
  137 + ;print,unit_dustem_predicted_extinction_sed
  138 +ENDIF ELSE BEGIN
  139 + unit_dustem_shown_extinction_sed=-1 ;That would be abnormal
  140 +ENDELSE
127 141  
128 142 ind=where(ext_names EQ 'DUSTEM_PREDICTED_EMISSION_SPECTRA',count)
129 143 IF count NE 0 THEN BEGIN
... ... @@ -174,6 +188,12 @@ IF unit_dustem_predicted_sed NE -1 THEN BEGIN
174 188 str_predicted_SED=mrdfits(file,unit_dustem_predicted_sed,header_predicted_SED)
175 189 ENDIF ELSE BEGIN
176 190 ENDELSE
  191 +;==== read the shown SED (if any)
  192 +;stop
  193 +IF unit_dustem_shown_sed NE -1 THEN BEGIN
  194 + str_shown_SED=mrdfits(file,unit_dustem_shown_sed,header_shown_SED)
  195 +ENDIF ELSE BEGIN
  196 +ENDELSE
177 197 IF unit_input_extinction_sed NE -1 THEN BEGIN
178 198 str_input_EXT=mrdfits(file,unit_input_extinction_sed,header_input_ext)
179 199 ENDIF ELSE BEGIN
... ... @@ -182,6 +202,10 @@ IF unit_dustem_predicted_extinction_sed NE -1 THEN BEGIN
182 202 str_predicted_EXT=mrdfits(file,unit_input_extinction_sed,header_input_ext)
183 203 ENDIF ELSE BEGIN
184 204 ENDELSE
  205 +IF unit_dustem_shown_extinction_sed NE -1 THEN BEGIN
  206 + str_shown_EXT=mrdfits(file,unit_shown_extinction_sed,header_shown_ext)
  207 +ENDIF ELSE BEGIN
  208 +ENDELSE
185 209  
186 210 used_model=strtrim(sxpar(main_header,'MODEL'),2)
187 211 used_pol=long(strtrim(sxpar(main_header,'POL'),2))
... ... @@ -226,62 +250,51 @@ ENDFOR
226 250 (*!dustem_fit).PARAM_FUNC=ptr_new(parameters_func_values)
227 251  
228 252  
229   -;stop
230 253 ;==== This is to create (*!dustem_data).sed in case it does not exist already
231 254 IF unit_input_emission_sed NE -1 THEN BEGIN
232 255 Nsed=n_elements(str_input_SED)
233 256 sed=dustem_initialize_internal_sed(Nsed)
234   - (*!dustem_data).sed=ptr_new(sed)
235   -
236   - (*(*!dustem_data).sed).instru_names=dustem_filter2instru(str_input_SED.filter)
237   - (*(*!dustem_data).sed).filt_names=str_input_SED.filter
238   - (*(*!dustem_data).sed).wav=str_input_SED.wavelength
239   - (*(*!dustem_data).sed).values=str_input_SED.STOKESI
240   - (*(*!dustem_data).sed).sigma=la_power(str_input_SED.varianceII,0.5)
241   -
242   - IF used_pol THEN BEGIN
243   - ;stop
244   - (*!dustem_data).qsed=ptr_new(sed)
245   - (*!dustem_data).used=ptr_new(sed)
246   - (*(*!dustem_data).qsed).instru_names=dustem_filter2instru(str_input_SED.filter)
247   - (*(*!dustem_data).qsed).filt_names=str_input_SED.filter
248   - (*(*!dustem_data).qsed).wav=str_input_SED.wavelength
249   - (*(*!dustem_data).qsed).values=str_input_SED.STOKESQ
250   - (*(*!dustem_data).qsed).sigma=la_power(str_input_SED.varianceQQ,0.5)
251   - (*(*!dustem_data).used).instru_names=dustem_filter2instru(str_input_SED.filter)
252   - (*(*!dustem_data).used).filt_names=str_input_SED.filter
253   - (*(*!dustem_data).used).wav=str_input_SED.wavelength
254   - (*(*!dustem_data).used).values=str_input_SED.STOKESU
255   - (*(*!dustem_data).used).sigma=la_power(str_input_SED.varianceUU,0.5)
256   - ENDIF
  257 + dustem_fillup_systvar_from_fits,*!dustem_data,sed,str_input_SED,used_pol
  258 +ENDIF
  259 +
  260 +;help,str_shown_SED
  261 +;stop
  262 +
  263 +;==== This is to create (*!dustem_show).sed in case it does not exist already
  264 +IF unit_dustem_shown_sed NE -1 THEN BEGIN
  265 + Nsed=n_elements(str_shown_SED)
  266 + sed=dustem_initialize_internal_sed(Nsed)
  267 + dustem_fillup_systvar_from_fits,*!dustem_show,sed,str_shown_SED,used_pol
257 268 ENDIF
258 269  
259 270 IF unit_input_extinction_sed NE -1 THEN BEGIN
260 271 Next=n_elements(str_input_EXT)
261 272 ext=dustem_initialize_internal_sed(Next)
262   - (*!dustem_data).ext=ptr_new(ext)
263   -
264   - (*(*!dustem_data).ext).instru_names=dustem_filter2instru(str_input_EXT.filter)
265   - (*(*!dustem_data).ext).filt_names=str_input_EXT.filter
266   - (*(*!dustem_data).ext).wav=str_input_EXT.wavelength
267   - (*(*!dustem_data).ext).values=str_input_EXT.I
268   - (*(*!dustem_data).ext).sigma=la_power(str_input_EXT.varianceII,0.5)
269   -
270   - IF used_pol THEN BEGIN
271   - ;stop
272   - (*!dustem_data).qext=ptr_new(ext)
273   - (*!dustem_data).uext=ptr_new(ext)
274   - (*(*!dustem_data).qext).instru_names=dustem_filter2instru(str_input_EXT.filter)
275   - (*(*!dustem_data).qext).filt_names=str_input_EXT.filter
276   - (*(*!dustem_data).qext).wav=str_input_EXT.wavelength
277   - (*(*!dustem_data).qext).values=str_input_EXT.Q
278   - (*(*!dustem_data).qext).sigma=la_power(str_input_EXT.varianceQQ,0.5)
279   - (*(*!dustem_data).uext).instru_names=dustem_filter2instru(str_input_EXT.filter)
280   - (*(*!dustem_data).uext).filt_names=str_input_EXT.filter
281   - (*(*!dustem_data).uext).wav=str_input_EXT.wavelength
282   - (*(*!dustem_data).uext).values=str_input_EXT.U
283   - (*(*!dustem_data).uext).sigma=la_power(str_input_EXT.varianceUU,0.5)
284   - ENDIF
  273 + dustem_fillup_systvar_from_fits,*!dustem_data,ext,str_input_EXT,used_pol
  274 +
  275 + ; (*!dustem_data).ext=ptr_new(ext)
  276 +
  277 + ; (*(*!dustem_data).ext).instru_names=dustem_filter2instru(str_input_EXT.filter)
  278 + ; (*(*!dustem_data).ext).filt_names=str_input_EXT.filter
  279 + ; (*(*!dustem_data).ext).wav=str_input_EXT.wavelength
  280 + ; (*(*!dustem_data).ext).values=str_input_EXT.I
  281 + ; (*(*!dustem_data).ext).sigma=la_power(str_input_EXT.varianceII,0.5)
  282 +
  283 + ; IF used_pol THEN BEGIN
  284 + ; ;stop
  285 + ; (*!dustem_data).qext=ptr_new(ext)
  286 + ; (*!dustem_data).uext=ptr_new(ext)
  287 + ; (*(*!dustem_data).qext).instru_names=dustem_filter2instru(str_input_EXT.filter)
  288 + ; (*(*!dustem_data).qext).filt_names=str_input_EXT.filter
  289 + ; (*(*!dustem_data).qext).wav=str_input_EXT.wavelength
  290 + ; (*(*!dustem_data).qext).values=str_input_EXT.Q
  291 + ; (*(*!dustem_data).qext).sigma=la_power(str_input_EXT.varianceQQ,0.5)
  292 + ; (*(*!dustem_data).uext).instru_names=dustem_filter2instru(str_input_EXT.filter)
  293 + ; (*(*!dustem_data).uext).filt_names=str_input_EXT.filter
  294 + ; (*(*!dustem_data).uext).wav=str_input_EXT.wavelength
  295 + ; (*(*!dustem_data).uext).values=str_input_EXT.U
  296 + ; (*(*!dustem_data).uext).sigma=la_power(str_input_EXT.varianceUU,0.5)
  297 + ; ENDIF
285 298 ENDIF
286 299  
287 300 ;=== Below is only to get the same output form as dustem_compute_sed
... ...
src/idl/dustem_write_fits_table.pro
... ... @@ -33,15 +33,12 @@ PRO dustem_write_fits_table,filename=filename,help=help
33 33 ;
34 34 ; SIDE EFFECTS:
35 35 ; Output FITS file is written to disk.
36   -;
  36 +; All dustem system variables are reinitialized
37 37 ; RESTRICTIONS:
38 38 ; The DustEMWrap IDL code must be installed
39   -;
40 39 ; PROCEDURE:
41 40 ; None
42   -;
43 41 ; EXAMPLES
44   -;
45 42 ; MODIFICATION HISTORY:
46 43 ; Written by J.-Ph. Bernard July 2022
47 44 ; Evolution details on the DustEMWrap gitlab.
... ... @@ -53,6 +50,8 @@ IF keyword_set(help) THEN BEGIN
53 50 goto,the_end
54 51 ENDIF
55 52  
  53 +;stop
  54 +
56 55 used_model=!dustem_model
57 56  
58 57 emission_sed_present=ptr_valid((*!dustem_data).sed)
... ... @@ -72,83 +71,10 @@ IF emission_sed_present THEN BEGIN
72 71 dustem_predicted_Qsed=dustem_predicted_polsed[0]
73 72 dustem_predicted_Used=dustem_predicted_polsed[1]
74 73 ENDIF
75   -
76   - ;===== define structures containing results
77   - one_str_input_SED={FILTER:'',Wavelength:la_undef(4),STOKESI:la_undef(5),STOKESQ:la_undef(5),STOKESU:la_undef(5),varianceII:la_undef(5),varianceQQ:la_undef(5),varianceUU:la_undef(5)}
78   - one_str_predicted_SED={FILTER:'',Wavelength:la_undef(4),STOKESI:la_undef(5),STOKESQ:la_undef(5),STOKESU:la_undef(5)}
79   - ;one_str_predicted_spectrum_tot={Wavelength:0.,I:0.,Q:0.,U:0.}
80   - ;one_str_predicted_spectrum_grain1={Wavelength:0.,I:0.,Q:0.,U:0.}
81   -
82   - Nsed=n_elements((*(*!dustem_data).sed).filt_names)
83   - str_input_SED=replicate(one_str_input_SED,Nsed)
84   - str_input_SED.FILTER=(*(*!dustem_data).sed).filt_names
85   - str_input_SED.Wavelength=(*(*!dustem_data).sed).wav
86   - str_input_SED.STOKESI=(*(*!dustem_data).sed).values
87   - str_input_SED.varianceII=la_power((*(*!dustem_data).sed).sigma,2) ;This is the variance
88   -
89   - N_predicted_sed=n_elements(dustem_predicted_sed)
90   - str_predicted_SED=replicate(one_str_predicted_SED,N_predicted_sed)
91   - ;str_predicted_EXT=replicate(one_str_predicted_SED,N_predicted_sed)
92   - fully_undefined_sed=str_predicted_SED
93   - str_predicted_SED.FILTER=(*(*!dustem_data).sed).filt_names ;taken from input SED
94   - str_predicted_SED.Wavelength=(*(*!dustem_data).sed).wav ;taken from input SED
95   - str_predicted_SED.STOKESI=dustem_predicted_sed
96   - ;stop
97   - IF !run_pol THEN BEGIN
98   - qsed=aos2soa(*(*!dustem_data).qsed)
99   - ind=where(qsed.filt_names NE 'SPECTRUM',count,complement=inds,Ncomplement=counts)
100   - IF count NE 0 THEN BEGIN
101   - qsed=qsed[ind]
102   - FOR i=0L,n_elements(qsed.wav)-1 DO BEGIN
103   - ind=where(str_input_SED.FILTER EQ qsed[i].filt_names,count)
104   - IF count NE 0 THEN BEGIN
105   - str_input_SED[ind[0]].STOKESQ=qsed[i].values
106   - str_input_SED[ind[0]].varianceQQ=la_power(qsed[i].sigma,2.)
107   - ENDIF
108   - str_predicted_SED[ind[0]].STOKESQ=dustem_predicted_Qsed[i]
109   - ENDFOR
110   - ENDIF
111   - qsed=aos2soa(*(*!dustem_data).qsed)
112   - IF counts NE 0 THEN BEGIN
113   - qsed=qsed[inds]
114   - ind2=where(str_predicted_SED.filter EQ 'SPECTRUM',count2)
115   - str_predicted_SED[ind2].STOKESQ=interpol(qsed.values,qsed.wav,str_predicted_SED[ind2].Wavelength)
116   - str_predicted_SED[ind2].varianceQQ=interpol(la_power(qsed.sigma,2.),qsed.wav,str_predicted_SED[ind2].Wavelength)
117   - ENDIF
118   - used=aos2soa(*(*!dustem_data).used)
119   - ind=where(used.filt_names NE 'SPECTRUM',count,complement=inds,Ncomplement=counts)
120   - IF count NE 0 THEN BEGIN
121   - used=used[ind]
122   - FOR i=0L,n_elements(used.wav)-1 DO BEGIN
123   - ind=where(str_input_SED.FILTER EQ used[i].filt_names,count)
124   - IF count NE 0 THEN BEGIN
125   - str_input_SED[ind[0]].STOKESU=used[i].values
126   - str_input_SED[ind[0]].varianceUU=la_power(used[i].sigma,2.)
127   - ENDIF
128   - str_predicted_SED[ind[0]].STOKESU=dustem_predicted_Used[i]
129   - ENDFOR
130   - ENDIF
131   - used=aos2soa(*(*!dustem_data).used)
132   - IF counts NE 0 THEN BEGIN
133   - used=used[inds]
134   - ind2=where(str_predicted_SED.filter EQ 'SPECTRUM',count2)
135   - str_predicted_SED[ind2].STOKESU=interpol(used.values,used.wav,str_predicted_SED[ind2].Wavelength)
136   - str_predicted_SED[ind2].varianceUU=interpol(la_power(used.sigma,2.),used.wav,str_predicted_SED[ind2].Wavelength)
137   - ENDIF
138   - ;stop
139   - ;str_predicted_SED.Q=dustem_predicted_Qsed
140   - ;str_predicted_SED.U=dustem_predicted_Used
141   - ;str_predicted_SED.Q=interpol(dustem_predicted_Qsed,(*(*!dustem_data).qsed).wav,(*(*!dustem_data).sed).wav)
142   - ;str_predicted_SED.U=interpol(dustem_predicted_Used,(*(*!dustem_data).used).wav,(*(*!dustem_data).sed).wav)
143   - ;ENDFOR
144   - ENDIF ELSE BEGIN
145   - str_predicted_SED.STOKESQ=fully_undefined_sed.STOKESQ
146   - str_predicted_SED.STOKESU=fully_undefined_sed.STOKESU
147   - ENDELSE
  74 + str_predicted_SED=dustem_make_fits_predicted_sed(*!dustem_data,dustem_predicted_sed,str_input_sed=str_input_sed)
  75 + str_predicted_shown_SED=dustem_make_fits_predicted_sed(*!dustem_show,dustem_predicted_sed,str_input_sed=str_shown_sed)
148 76 ENDIF
149 77  
150   -;stop
151   -
152 78 IF extinction_sed_present THEN BEGIN
153 79 ;stop
154 80 dustem_predicted_ext=dustem_compute_ext(*(*!dustem_fit).current_param_values,st=dustem_st)
... ... @@ -158,38 +84,8 @@ IF extinction_sed_present THEN BEGIN
158 84 dustem_predicted_Qext=dustem_predicted_polext[0]
159 85 dustem_predicted_Uext=dustem_predicted_polext[1]
160 86 ENDIF
161   -
162   - ;===== define structures containing results
163   - 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)}
164   - one_str_predicted_EXT={FILTER:'',Wavelength:la_undef(4),I:la_undef(5),Q:la_undef(5),U:la_undef(5)}
165   -
166   - Next=n_elements((*(*!dustem_data).ext).filt_names)
167   - str_input_EXT=replicate(one_str_input_EXT,Next)
168   - str_input_EXT.FILTER=(*(*!dustem_data).ext).filt_names
169   - str_input_EXT.Wavelength=(*(*!dustem_data).ext).wav
170   - str_input_EXT.I=(*(*!dustem_data).ext).values
171   - str_input_EXT.varianceII=la_power((*(*!dustem_data).ext).sigma,2.)
172   -
173   - ;N_predicted_ext=n_elements(dustem_predicted_polext[0])
174   - N_predicted_ext=n_elements(dustem_predicted_ext)
175   - str_predicted_EXT=replicate(one_str_predicted_EXT,N_predicted_EXT)
176   - fully_undefined_ext=str_predicted_ext
177   - IF N_predicted_ext NE Next THEN BEGIN ;This is for cases when the extinction data to fit did not have the dimension of computed extinction curve
178   - message,'I have a problem with dimensions ...',/continue
179   - stop
180   - ENDIF ELSE BEGIN
181   - str_predicted_EXT.FILTER=(*(*!dustem_data).ext).filt_names ;taken from input SED
182   - str_predicted_EXT.Wavelength=(*(*!dustem_data).ext).wav ;taken from input SED
183   - str_predicted_EXT.I=dustem_predicted_ext
184   - ENDELSE
185   - ;stop
186   - IF !run_pol THEN BEGIN
187   - str_predicted_EXT.Q=interpol(dustem_predicted_Qext,(*(*!dustem_data).qext).wav,(*(*!dustem_data).ext).wav)
188   - str_predicted_EXT.U=interpol(dustem_predicted_Uext,(*(*!dustem_data).uext).wav,(*(*!dustem_data).ext).wav)
189   - ENDIF ELSE BEGIN
190   - str_predicted_EXT.Q=fully_undefined_ext.Q
191   - str_predicted_EXT.U=fully_undefined_ext.U
192   - ENDELSE
  87 + str_predicted_EXT=dustem_make_fits_predicted_ext(*!dustem_data,dustem_predicted_ext,str_input_ext=str_input_ext)
  88 + str_predicted_shown_EXT=dustem_make_fits_predicted_ext(*!dustem_show,dustem_predicted_ext,str_input_ext=str_shown_ext)
193 89 ENDIF
194 90  
195 91 Ngrains=n_tags(dustem_st.sed)-2
... ... @@ -219,22 +115,32 @@ IF emission_sed_present THEN BEGIN
219 115 mwrfits,str_predicted_SED,file
220 116 unit_predicted_sed=unit
221 117 unit=unit+1
  118 + ;==== add the shown emission SED
  119 + mwrfits,str_shown_SED,file
  120 + unit_shown_sed=unit
  121 + unit=unit+1
222 122 ENDIF ELSE BEGIN
223 123 unit_input_sed=-1
224 124 unit_predicted_sed=-1
  125 + unit_shown_sed=-1
225 126 ENDELSE
226 127 IF extinction_sed_present THEN BEGIN
227 128 ;stop
228 129 mwrfits,str_input_EXT,file
229 130 unit_input_ext=unit
230 131 unit=unit+1
231   - ;==== add the predicted emission SED
  132 + ;==== add the predicted extinction SED
232 133 mwrfits,str_predicted_EXT,file
233 134 unit_predicted_ext=unit
234 135 unit=unit+1
  136 + ;==== add the shown extinction SED
  137 + mwrfits,str_shown_EXT,file
  138 + unit_shown_ext=unit
  139 + unit=unit+1
235 140 ENDIF ELSE BEGIN
236 141 unit_input_ext=-1
237 142 unit_predicted_ext=-1
  143 + unit_shown_ext=-1
238 144 ENDELSE
239 145 ;==== add the predicted per-grain emission spectrum
240 146 mwrfits,dustem_st.sed,file
... ... @@ -256,12 +162,12 @@ ENDIF ELSE BEGIN
256 162 unit_predicted_extinctions_pol=-1
257 163 ENDELSE
258 164  
259   -extension_numbers=[unit_input_sed,unit_predicted_sed, $
260   - unit_input_ext,unit_predicted_ext, $
  165 +extension_numbers=[unit_input_sed,unit_predicted_sed,unit_shown_sed, $
  166 + unit_input_ext,unit_predicted_ext,unit_shown_ext, $
261 167 unit_predicted_spectra,unit_predicted_extinctions, $
262 168 unit_predicted_spectra_pol,unit_predicted_extinctions_pol]
263   -extensions_strs=['INPUT_EMISSION_SED','DUSTEM_PREDICTED_SED', $
264   - 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED', $
  169 +extensions_strs=['INPUT_EMISSION_SED','DUSTEM_PREDICTED_SED','DUSTEM_SHOWN_SED', $
  170 + 'INPUT_EXTINCTION_SED','DUSTEM_PREDICTED_EXTINCTION_SED','DUSTEM_SHOWN_EXTINCTION_SED', $
265 171 'DUSTEM_PREDICTED_EMISSION_SPECTRA','DUSTEM_PREDICTED_EXTINCTION_SPECTRA', $
266 172 'DUSTEM_PREDICTED_EMISSION_SPECTRA_POL','DUSTEM_PREDICTED_EXTINCTION_SPECTRA_POL']
267 173  
... ... @@ -278,8 +184,10 @@ file='./dustemwrap_tmp.fits'
278 184 toto=mrdfits(file,0,hdr) ;this would be to get the main header
279 185 IF unit_input_sed NE -1 THEN sst_input_sed=mrdfits(file,unit_input_sed,header_input_sed)
280 186 IF unit_predicted_sed NE -1 THEN sst_predicted_SED=mrdfits(file,unit_predicted_sed,header_predicted_SED)
  187 +IF unit_shown_sed NE -1 THEN sst_shown_SED=mrdfits(file,unit_shown_sed,header_shown_SED)
281 188 IF unit_input_ext NE -1 THEN sst_input_ext=mrdfits(file,unit_input_ext,header_input_ext)
282 189 IF unit_predicted_ext NE -1 THEN sst_predicted_ext=mrdfits(file,unit_predicted_ext,header_predicted_ext)
  190 +IF unit_shown_ext NE -1 THEN sst_shown_ext=mrdfits(file,unit_shown_ext,header_shown_ext)
283 191 dustem_st.sed=mrdfits(file,unit_predicted_spectra,header_predicted_spectra)
284 192 dustem_st.ext=mrdfits(file,unit_predicted_extinctions,header_predicted_extinctions)
285 193 IF !run_pol EQ 1 THEN BEGIN
... ... @@ -370,18 +278,27 @@ ENDELSE
370 278 ;==== Write the observed SED (if any)
371 279 mwrfits,0,file,hdr,create=1
372 280 IF unit_input_sed NE -1 THEN BEGIN
373   - mwrfits,str_input_SED,file,header_input_sed
  281 + mwrfits,str_input_SED,file,header_input_SED
374 282 ENDIF
375 283 ;==== add the predicted SED (if any)
376   -IF unit_input_sed NE -1 THEN BEGIN
  284 +IF unit_predicted_sed NE -1 THEN BEGIN
377 285 mwrfits,str_predicted_SED,file,header_predicted_SED
378 286 ENDIF
  287 +;==== add the shown SED (if any)
  288 +IF unit_shown_sed NE -1 THEN BEGIN
  289 + mwrfits,str_shown_SED,file,header_shown_SED
  290 +ENDIF
  291 +;==== Write the observed extinction SED (if any)
379 292 IF unit_input_ext NE -1 THEN BEGIN
380   - mwrfits,str_input_EXT,file,header_input_sed
  293 + mwrfits,str_input_EXT,file,header_input_EXT
381 294 ENDIF
382 295 ;==== add the predicted SED (if any)
383 296 IF unit_predicted_ext NE -1 THEN BEGIN
384   - mwrfits,str_predicted_EXT,file,header_predicted_SED
  297 + mwrfits,str_predicted_EXT,file,header_predicted_EXT
  298 +ENDIF
  299 +;==== add the show extinction SED (if any)
  300 +IF unit_shown_ext NE -1 THEN BEGIN
  301 + mwrfits,str_shown_EXT,file,header_shown_EXT
385 302 ENDIF
386 303 ;==== add the predicted per-grain emission spectra
387 304 mwrfits,dustem_st.sed,file,header_predicted_spectra
... ...