OV-GSO-DC / dustem-wrapper_idl

Commit fcb6eade7b308d098425680732fca2f2af4cc1c5

Authored by Ilyes Choubani
1 parent 5f04fa07
Exists in master

general update - putting some order and cleaning some routines

Showing 7 changed files with 232 additions and 702 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_compute_polsed.pro
  Diff comments   View file @fcb6ead
1   -FUNCTION dustem_compute_polsed,p_dim,st,P_spec,SP_spec,out_st=out_st,_extra=extra
  1 +FUNCTION dustem_compute_polsed,p_dim,st,P_spec,SP_spec,dustem_polfrac,out_st=out_st,_extra=extra
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -134,6 +134,53 @@ IF count_spec NE 0 THEN dustem_polsed(ind_spec)=interpol(P_spec,st.polsed.wav,((
134 134  
135 135 out_st=st
136 136  
  137 +;GENERATING THE INTERPOLATES FOR POLFRAC (dustem_polfrac)
  138 +if !run_lin then begin
  139 +
  140 + if ~isa(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
  141 + if ~isa(dustem_polsed) then dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec) ; this quantity hasn't been computed yet that is why this line is present here.
  142 + If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin
  143 + if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin
  144 +
  145 + dustem_polsed_x = dustem_polsed
  146 + dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
  147 +
  148 + nwaves = n_elements(dustem_polsed)
  149 +
  150 +
  151 + for i=0L,nwaves-1 do begin
  152 +
  153 + j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(i),testwav)
  154 + if testwav ne 0 then dustem_sed_x(i) = dustem_sed(j(0))
  155 +
  156 + endfor
  157 +
  158 + endif ELSE begin
  159 +
  160 + dustem_polsed_x = dustem_sed
  161 + dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
  162 +
  163 + nwaves = n_elements(dustem_sed)
  164 +
  165 + for i=0L,nwaves-1 do begin
  166 + j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(i),testwav)
  167 + if testwav ne 0 then dustem_polsed_x(i) = dustem_polsed(j(0))
  168 + endfor
  169 +
  170 +
  171 + ENDELSE
  172 +
  173 + endif ELSE BEGIN
  174 +
  175 + dustem_polsed_x = dustem_polsed
  176 + dustem_sed_x = dustem_sed
  177 +
  178 + ENDELSE
  179 +
  180 + dustem_polfrac = dustem_polsed_x/dustem_sed_x
  181 +
  182 +ENDIF
  183 +
137 184  
138 185 ;clean pointers
139 186 heap_gc
... ...
src/idl/dustem_compute_stokes.pro
  Diff comments   View file @fcb6ead
1   -FUNCTION dustem_compute_stokes,p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,out_st=out_st,_extra=extra
  1 +FUNCTION dustem_compute_stokes,p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em,out_st=out_st,_extra=extra
2 2  
3 3 ;,SP_spec
4 4  
... ... @@ -155,6 +155,10 @@ ENDIF
155 155  
156 156 out_st=st
157 157  
  158 +;generating interpolates for Psi_em
  159 +dustem_psi_em = 0.5*atan(dustem_used,dustem_qsed)/degtorad
  160 +
  161 +
158 162 ;clean pointers
159 163 heap_gc
160 164  
... ...
src/idl/dustem_fit_sed_orion.pro
  Diff comments   View file @fcb6ead
... ... @@ -245,9 +245,9 @@ st_ie=dustem_set_data(st_sed_fit=fit_struct[ind],st_sed_show=show_st[ind],rchi2_
245 245 ;stop
246 246  
247 247 ;=== RUN fit
248   -tol=1.e-16
249   -xtol=1.e-16
250   -Nitermax=60;0; just an initialization
  248 +tol=1.e-10
  249 +xtol=1.e-10
  250 +Nitermax=3;0; just an initialization
251 251  
252 252 xr=[1.,5e5]
253 253 yr=[5e-8,1.00e6]
... ... @@ -263,17 +263,12 @@ ytit=textoidl(sttr)
263 263  
264 264 ;xtit=''
265 265 t1=systime(0,/sec)
266   -res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,xtit=xtit,ytit=ytit,title=title);
  266 +res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,xtit=xtit,ytit=ytit,title=title,_extra=_extra);
267 267 t2=systime(0,/sec)
268 268  
269   -
270   -
271   -stop
272 269 ;help,(*!dustem_fit),/STRUCTURE
273 270  
274 271  
275   -
276   -
277 272 ;=== SAVE FIT RESULTS
278 273  
279 274 ;NB: THIS LINE IS WRONG BECAUSE I AM NOT USING THE ANGLE CORRELATION CRITERIA
... ... @@ -288,26 +283,33 @@ message,'Saved '+file,/continue
288 283 dustem_restore_system_variables,file
289 284  
290 285 ;=== Plot best fit
  286 +
291 287 win=1 ; who wrote this?
292 288  
293 289 window,win & win=win+1
294 290  
  291 +
295 292 ;=== recover best fit values
296 293 res=*(*!dustem_fit).current_param_values
297   -chi2=(*!dustem_fit).chi2
298   -rchi2=(*!dustem_fit).rchi2
299   -errors=*(*!dustem_fit).current_param_errors
300   -errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
301 294  
302   -;=== Plot best fit
  295 +dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=title
  296 +print, 'DONE!!! IT WORKED'
  297 +
  298 +
  299 +; chi2=(*!dustem_fit).chi2
  300 +; rchi2=(*!dustem_fit).rchi2
  301 +; errors=*(*!dustem_fit).current_param_errors
  302 +; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  303 +
  304 +; ;=== Plot best fit
303 305 ; tit='DUSTEM polarization dust Example (Saved)'
304 306 ; ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
305 307 ; xtit=textoidl('\lambda (\mum)')
306   -
307   -loadct,13
308   -
309   -ps_file=!dustem_wrap_soft_dir+'Docs/'+'Last_dustem_fit.ps'
310   -dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=title,yr=yr,xr=xr,/ysty,/xsty,res=res,chi2=chi2,rchi2=rchi2,/xlog,/ylog,/pol,ps=ps_file,png=png,fpol=fpol
  308 +;
  309 +; loadct,13
  310 +;
  311 +; ps_file=!dustem_wrap_soft_dir+'Docs/'+'Last_dustem_fit.ps'
  312 +; dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=title,yr=yr,xr=xr,/ysty,/xsty,res=res,chi2=chi2,rchi2=rchi2,/xlog,/ylog,/pol,ps=ps_file,png=png,fpol=fpol
311 313  
312 314  
313 315 print,'dustem_mpfit_sed executed in ',t2-t1,' sec'
... ...
src/idl/dustem_fit_sed_readme.pro
  Diff comments   View file @fcb6ead
... ... @@ -121,22 +121,7 @@ CASE use_model OF
121 121 '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
122 122 '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
123 123 'dustem_plugin_continuum_2'] ;Intensity of NIR continuum
124   - ;]
125   -
126   -; pd = [ $
127   -; '(*!dustem_params).gas.G0', $ ;G0
128   -; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
129   -; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
130   -; '(*!dustem_params).grains(2).mdust_o_mh', $ ;Gra
131   -; '(*!dustem_params).grains(3).mdust_o_mh', $ ;Gra
132   -; '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
133   -; 'dustem_create_continuum_2', $
134   -; 'dustem_create_stellar_population_O63',$ ; distance to the O3 stellar population
135   -; 'dustem_create_stellar_population_O65',$ ; number of stars of the O3 stellar population
136   -; 'dustem_create_stellar_population_B43',$ ; distance to the B4 stellar population
137   -; 'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population
138   -; ]
139   -
  124 +
140 125 iv = [1.0,5.4e-4, 5.4e-4,1.8e-4,2.33e-3,8.27e-3,0.001];,10,1.,10.,1.]
141 126  
142 127  
... ... @@ -173,104 +158,8 @@ CASE use_model OF
173 158 7.8e-3 $ ;mass fraction of
174 159 ]
175 160  
176   -;NEW LINES
177   -
178   -;========COMMENTS=======
179   -;amplitudes haven't been used so they're by default set to 1.
180   -;using only the below stellar populations produces amplitudes that exceed 1. I do not understand the physical meaning of this.
181   -
182   -; pd = [ $
183   -; '(*!dustem_params).gas.G0', $ ;G0
184   -; 'dustem_create_continuum_2', $ ;intensity of NIR continuum
185   -; 'dustem_create_stellar_population_O63',$ ; distance to the O3 stellar population
186   -; ;'dustem_create_stellar_population_O64',$ ; amplitude of the O3 stellar population
187   -; 'dustem_create_stellar_population_O65',$ ; number of stars of the O3 stellar population
188   -; 'dustem_create_stellar_population_B43',$ ; distance to the B4 stellar population
189   -; ;'dustem_create_stellar_population_B44',$ ; amplitude of the B4 stellar population
190   -; 'dustem_create_stellar_population_B45',$ ; number of stars of the B4 stellar population
191   -; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
192   -; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
193   -; '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
194   -; '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
195   -; '(*!dustem_params).grains(4).mdust_o_mh' $
196   -; ]
197   -
198   - ; pd = [ $
199   - ; ;'(*!dustem_params).G0', $ ;G0
200   - ; '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
201   - ; '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
202   - ; '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
203   - ; '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
204   - ; '(*!dustem_params).grains(4).mdust_o_mh', $
205   - ; ;'dustem_create_continuum_2', $ ;intensity of NIR continuum
206   - ; 'dustem_create_stellar_population_O63',$ ; distance to the O6 stellar population
207   -
208   - ; ;'dustem_create_stellar_population_O64',$ ; amplitude of the O6 stellar population
209   -
210   - ; ;'dustem_create_stellar_population_O65',$ ; number of stars of the O6 stellar population
211   - ; 'dustem_create_stellar_population_B43'$ ; distance to the B4 stellar population
212   -
213   - ; ;;'dustem_create_stellar_population_B44',$ ; amplitude of the B4 stellar population
214   -
215   - ; ;'dustem_create_stellar_population_B45'$ ; number of stars of the B4 stellar population
216   - ; ]
217   -
218   -
219   - ;initial parameter values for parameters to be fitted
220   - ; iv = [ $
221   - ; 1., $ ;G0
222   - ; 7.8e-4, $ ;mass fraction of PAH0
223   - ; 7.8e-4, $ ;mass fraction of PAH1
224   - ; 1.65e-4, $ ;mass fraction of amCBEx
225   - ; 1.45e-3, $ ;mass fraction of amCBEx
226   - ; 7.8e-3, $ ;mass fraction of big grains?
227   - ; 0.002, $ ;intensity of NIR continuum
228   - ; 10.,$
229   - ; ;;1.,$
230   - ; 1.,$
231   - ; 10.,$
232   - ; ;;1.,$
233   - ; 1.$
234   - ; ]
235   -;
236   -
237   - ; iv = [ $
238   - ; 7.8e-4, $ ;mass fraction of PAH0
239   - ; 7.8e-4, $ ;mass fraction of PAH1
240   - ; 1.65e-4, $ ;mass fraction of amCBEx
241   - ; 1.45e-3, $ ;mass fraction of amCBEx
242   - ; 7.8e-3, $ ;mass fraction of big grains?
243   - ; ;0.002, $ ;intensity of NIR continuum
244   - ; 10.,$
245   - ; ;;1.,$
246   - ; ;1.,$
247   - ; 10.$
248   - ; ;;1.,$
249   - ; ;1.$
250   - ; ]
251   -
252   -
253   -
254   -
255 161 Npar=n_elements(pd)
256   -
257   - ;ulimed=[0 ,0 ,0 ,0 ,0 ,0 , 0 ,0 ,0 ,0 ,0 ,0 ,0]
258   - ;ulimed=[0 ,0 ,0 ,0 ,0 ,0 , 0 ,0 ,0 ,0]
259   -
260   -
261   - ;llimed=[1 , 1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1]
262   - ;llimed=[1 , 1 ,1 ,1 ,1 ,1 ,1 ,1 ,1 ,1]
263   -
264   - ;llims =[0. ,0. ,0. ,0. , 0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0. ,0.]
265   - ;llims =[1. ,0. ,10. ,0. , 0. ,0. ,0. ,0. ,0. ,0. ,0.]
266   -
267   - ;lower limiting the stellar distances to 1pc avoid singularities
268   - ;[1e-9 ,
269   - ;llims =[1e-5 ,0. ,0. ,0. ,0. ,0. ,0. ,1. ,0. , 1. ,0. ]
270   -
271   - ;llims =[0. ,0. ,0. ,0. ,0. ,0. ,1. ,0. , 1. ,0. ]
272   - ;llims =[0. ,0. ,0. ,0. ,0. ,1. , 1. ]
273   -
  162 +
274 163 ulimed=replicate(0,Npar)
275 164 llimed=replicate(1,Npar)
276 165 llims=replicate(0.,Npar)
... ... @@ -343,7 +232,7 @@ CASE use_model OF
343 232 END
344 233 'THEMIS':BEGIN
345 234 pd = [ $
346   - '(*!dustem_params).G0', $ h;G0
  235 + '(*!dustem_params).G0', $ ;G0
347 236 '(*!dustem_params).grains(0).mdust_o_mh',$
348 237 '(*!dustem_params).grains(1).mdust_o_mh',$
349 238 '(*!dustem_params).grains(2).mdust_o_mh', $
... ... @@ -387,13 +276,11 @@ st=dustem_set_data(st_sed_fit=spec);sed=spec)
387 276 ;== SET THE FITTED PARAMETERS
388 277 dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims
389 278  
390   -
391 279 ;dustem_init_fixed_params,fpd,fiv
392 280  
393 281 dustem_init_plugins, pd ;initializing the scopes of the plugins to be read later
394 282  
395 283  
396   -
397 284 ;=== RUN fit
398 285 tol=1.e-16
399 286 ;use_Nitermax=50 ;maximum number of iteration. This is the criterium which will stop the fit procedure , NOTA BENE: 3 iterations aren't enough to fit the entirety of the newly added stellar population parameters.
... ... @@ -428,43 +315,48 @@ file='/tmp/DUSTEM_fit_example.sav'
428 315 dustem_restore_system_variables,file
429 316  
430 317 ;=== Plot best fit
431   -tit='Spectral Energy Distribution (Saved)'
432   -ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
433   -xtit=textoidl('\lambda (\mum)')
434   -errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
435   -chi2=(*!dustem_fit).chi2
436   -rchi2=(*!dustem_fit).rchi2
437   -
438   -window,2
439   -
440   -;=== RESTORE FIT RESULTS
441   -res=*(*!dustem_fit).current_param_values
442   -chi2=(*!dustem_fit).chi2
443   -rchi2=(*!dustem_fit).rchi2
444   -errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
445   -
446   -;=== Plot best fit
447   -loadct,13
448   -IF keyword_set(postcript) THEN BEGIN
449   - set_plot,'PS'
450   -; ps_file=getenv('DUSTEM_SOFT_DIR')+'/Docs/Figures/'+'Last_dustem_fit.ps'
451   - ps_file=!dustem_wrap_soft_dir+'Docs/Figures/'+'Last_dustem_fit.ps'
452   - device,filename=ps_file,/color
453   -ENDIF
454   -
455   -;dustem_sed_plot,(*!dustem_fit_params),ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2
456   -dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos
457   -IF keyword_set(postcript) THEN BEGIN
458   - device,/close
459   - set_plot,'X'
460   - message,'wrote '+ps_file,/info
461   -ENDIF
462   -IF keyword_set(png) THEN BEGIN
463   - ;file_png=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.png'
464   - file_png=dir_png+'Last_dustem_fit_sed_readme.png'
465   - write_png,file_png,tvrd(/true)
466   - message,'Wrote '+file_png,/info
467   -ENDIF
  318 +; tit='Spectral Energy Distribution (Saved)'
  319 +; ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  320 +; xtit=textoidl('\lambda (\mum)')
  321 +; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  322 +; chi2=(*!dustem_fit).chi2
  323 +; rchi2=(*!dustem_fit).rchi2
  324 +;
  325 +; window,2
  326 +;
  327 +; ;=== RESTORE FIT RESULTS
  328 +; res=*(*!dustem_fit).current_param_values
  329 +; chi2=(*!dustem_fit).chi2
  330 +; rchi2=(*!dustem_fit).rchi2
  331 +; errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values)
  332 +;
  333 +; ;=== Plot best fit
  334 +; loadct,13
  335 +; IF keyword_set(postcript) THEN BEGIN
  336 +; set_plot,'PS'
  337 +; ; ps_file=getenv('DUSTEM_SOFT_DIR')+'/Docs/Figures/'+'Last_dustem_fit.ps'
  338 +; ps_file=!dustem_wrap_soft_dir+'Docs/Figures/'+'Last_dustem_fit.ps'
  339 +; device,filename=ps_file,/color
  340 +; ENDIF
  341 +;
  342 +; ;dustem_sed_plot,(*!dustem_fit_params),ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2
  343 +; dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yr,xr=xr,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos
  344 +;
  345 +;
  346 +; IF keyword_set(postcript) THEN BEGIN
  347 +; device,/close
  348 +; set_plot,'X'
  349 +; message,'wrote '+ps_file,/info
  350 +; ENDIF
  351 +; IF keyword_set(png) THEN BEGIN
  352 +; ;file_png=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.png'
  353 +; file_png=dir_png+'Last_dustem_fit_sed_readme.png'
  354 +; write_png,file_png,tvrd(/true)
  355 +; message,'Wrote '+file_png,/info
  356 +; ENDIF
  357 +
  358 +tit='Spectral Energy Distribution (Final run)'
  359 +dustemwrap_plot,res,stp,xr=xr,/xstyle,yr=yr,/ysty,/ylog,/xlog,title=tit
468 360  
469 361 message,'dustem_mpfit_sed executed in '+strtrim(t2-t1,2)+' sec',/info
470 362  
... ...
src/idl/dustem_mpfit_run.pro
  Diff comments   View file @fcb6ead
... ... @@ -52,7 +52,6 @@ IF !run_ionfrac gt 0. THEN toto = dustem_create_ionfrac()
52 52  
53 53 ;RUN THE MODEL AND GET THE SPECTRUM
54 54  
55   -;stop
56 55 dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values),st ;st is an output...
57 56  
58 57 ;stop
... ... @@ -98,17 +97,11 @@ tagnames=tag_names(!dustem_data)
98 97 dustem_out = [0]
99 98 ind_data_tag=dustem_data_tag(count=count_data_tag)
100 99  
101   -;stop
102   -
103   -
104   -;NB: THIS PART CONTAINS THE INTERPOLATES USED FOR THE COMPUTATION OF THE CHI-SQUARED
105   -
106   -
107 100 FOR i_tag=0,count_data_tag-1 DO BEGIN
108 101  
109 102 CASE tagnames(ind_data_tag(i_tag)) OF
110 103  
111   - 'SED' : BEGIN
  104 + 'SED' : BEGIN ;UPDATED. BUT ISSUE WITH THE CONCATENATION TO THE DUSTEM_OUT ARRAY AS OUTLINED BELOW
112 105  
113 106 dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
114 107 chi2_sed = total(((*!dustem_data.sed).values-dustem_sed)^2/(*!dustem_data.sed).sigma^2)
... ... @@ -129,14 +122,14 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
129 122 (*!dustem_fit).chi2=chi2_sed
130 123 (*!dustem_fit).rchi2=rchi2_sed
131 124  
132   - ;I think think this test isn't the correct one because when we want to show the SED data and not fit it we have an issue
133   - ;with the current hiding/showing formalism!!!!
134   - ;figure it out.
  125 + ;I think this test isn't the correct one because when we want to show the SED data and not fit it we have an issue
  126 + ;with the current hiding/showing formalism!!!!
  127 +
135 128 if isa(!dustem_data.sed) then dustem_out = [ dustem_out, dustem_sed ]
136 129  
137 130 END
138 131  
139   - 'EXT' : BEGIN ;STILL HAVEN'T DONE THE EXTINCTION PART
  132 + 'EXT' : BEGIN ;NOT UPDATED
140 133  
141 134 ;ADDING PLUGIN TO SPECTRUM----------------
142 135 IF ptr_valid(*!dustem_plugin) THEN BEGIN
... ... @@ -180,7 +173,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
180 173  
181 174 END
182 175  
183   - 'POLEXT': BEGIN
  176 + 'POLEXT': BEGIN ;NOT UPDATED
184 177  
185 178  
186 179 ;---------------the procedure within this block has not been updated for a vert long time
... ... @@ -271,9 +264,9 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
271 264  
272 265 END
273 266  
274   - 'QSED': BEGIN
  267 + 'QSED': BEGIN ;UPDATED
275 268  
276   - toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec)
  269 + toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em)
277 270  
278 271 ;taking care of the color correction for spass ; SHOULD NOT LEAVE THIS ON
279 272 ;dustem_qsed(-1)=dustem_qsed(-2)
... ... @@ -290,8 +283,8 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
290 283  
291 284 END
292 285  
293   - 'USED': BEGIN
294   - ;THIS IMPLIES THAT QSED HAS BEEN RAN BEFORE
  286 + 'USED': BEGIN ;UPDATED
  287 + ;THIS IMPLIES THAT QSED HAS BEEN RAN BEFORE
295 288 ;toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec)
296 289 chi2_used = total(((*!dustem_data.used).values-dustem_used)^2/(*!dustem_data.used).sigma^2)
297 290 n_used = n_elements((*!dustem_data.used).values)
... ... @@ -303,134 +296,35 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
303 296  
304 297 END
305 298  
306   - 'POLSED': BEGIN
  299 + 'POLSED': BEGIN ;UPDATED
307 300  
308 301 ;PLUGIN IS ADDED INSIDE DUSTEM_COMPUTE_SED
309 302  
310 303 IF !run_pol and !run_lin THEN BEGIN
311 304  
312   - dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
  305 + dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec,dustem_polfrac);dustem_compute_polsed(p_dim,st,cont=cont,freefree=freefree,synchrotron=synchrotron)
313 306  
314 307 ENDIF
315 308 END
316 309  
317   - 'PSI_EM' : BEGIN
  310 + 'PSI_EM' : BEGIN ;This block is present if PSI_EM data is to be shown but not the rest of the stokes parameters ?
318 311  
319 312 IF !run_pol and !run_lin THEN BEGIN
320 313 degtorad = !pi/180
321   - if ~isa(dustem_used) or ~isa(dustem_qsed) then toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec)
322   - dustem_psi_em = 0.5*atan(dustem_used,dustem_qsed)/degtorad
  314 + if ~isa(dustem_used) or ~isa(dustem_qsed) then toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em)
  315 + ;dustem_psi_em = 0.5*atan(dustem_used,dustem_qsed)/degtorad
323 316  
324 317  
325 318 ENDIF
326 319 END
327 320  
328 321  
329   - 'POLFRAC': BEGIN
330   -
331   - if !run_lin then begin
332   -
333   - if ~isa(dustem_sed) then dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
334   - if ~isa(dustem_polsed) then dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec) ; this quantity hasn't been computed yet that is why this line is present here.
335   -
336   -
337   - ;commenting this for now.
338   -
339   -; ind_filt=where((*!dustem_data.polfrac).filt_names NE 'SPECTRUM',count_filt) ;locating the filters of polfrac=smallp in SED xcat
340   -; ind_fspec=where((*!dustem_data.polfrac).filt_names EQ 'SPECTRUM',count_fspec) ;locating the spectrum points in polfrac data. Later testing will have to be on the wavelengths
341   -;
342   -; ;dustem_sed and dustem_polsed might not have the same length thus this block.
343   -;
344   - If n_elements(dustem_sed) ne n_elements(dustem_polsed) then begin
345   - if n_elements(dustem_sed) gt n_elements(dustem_polsed) then begin
346   -
347   - dustem_polsed_x = dustem_polsed
348   - dustem_sed_x = dustem_polsed ;meaning dustem_sed needs to be modified
349   -
350   - nwaves = n_elements(dustem_polsed)
351   - ;nwaves = n_elements((*!dustem_data.polfrac).wav) should give the same result.... TEST this
352   -
353   -
354   - ;NOTA BENE: We do not need to match the filter or the data points. We just need to loop over the wavelengths.
355   - ;This is what we do when we filter the data prior (in the main dustemwrap fitting procedure level).
356   -
357   - ;commenting this for now
358   - ;matching the filter data points
359   -; IF count_filt ne 0 then begin
360   -; for i=0L,count_filt-1 do begin
361   -; j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
362   -; if testfilt ne 0 then dustem_sed_x(ind_filt(i)) = dustem_sed(j(0))
363   -; endfor
364   -; endif
365   -
366   -
367   -; ;matching the spectrum data points
368   -; IF count_fspec ne 0 then begin
369   -; for i=0L,count_fspec-1 do begin
370   -;
371   -; j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
372   -; if testspec ne 0 then dustem_sed_x((ind_fspec(i))) = dustem_sed(j(0))
373   -; endfor
374   -; endif
375   -
376   - ;NEW SECTION
377   - for i=0L,nwaves-1 do begin
378   -
379   - j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(i),testwav)
380   - if testwav ne 0 then dustem_sed_x(i) = dustem_sed(j(0))
381   -
382   - endfor
383   -
384   -
385   -
386   - endif ELSE begin
387   -
388   - dustem_polsed_x = dustem_sed
389   - dustem_sed_x = dustem_sed ;meaning dustem_polsed needs to be modified
390   -
391   - nwaves = n_elements(dustem_sed)
392   -
393   -
394   -
395   -; ;matching the filter data points
396   -; IF count_filt ne 0 then begin
397   -; for i=0L,count_filt-1 do begin
398   -; j=where((*!dustem_data.sed).filt_names EQ (*!dustem_data.polfrac).filt_names(ind_filt(i)),testfilt)
399   -; if testfilt ne 0 then dustem_polsed_x(ind_filt(i)) = dustem_polsed(j(0))
400   -; endfor
401   -; endif
402   -;
403   -
404   -; ;matching the spectrum data points
405   -; IF count_fspec ne 0 then begin
406   -; for i=0L,count_fspec-1 do begin
407   -; j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(ind_fspec(i)),testspec)
408   -; if testspec ne 0 then dustem_polsed_x((ind_fspec(i))) = dustem_polsed(j(0))
409   -; endfor
410   -; endif
411   -
412   -
413   - for i=0L,nwaves-1 do begin
414   - j=where((*!dustem_data.sed).wav EQ (*!dustem_data.polfrac).wav(i),testwav)
415   - if testwav ne 0 then dustem_polsed_x(i) = dustem_polsed(j(0))
416   - endfor
417   -
418   -
419   - ENDELSE
420   - endif ELSE BEGIN
421   - dustem_polsed_x = dustem_polsed
422   - dustem_sed_x = dustem_sed
423   - ENDELSE
424   -
425   - dustem_polfrac = dustem_polsed_x/dustem_sed_x
426   -
427   - ENDIF
428   -
429   - ;stop
430   - END
  322 + 'POLFRAC': BEGIN ;We do not really need this block
  323 +
  324 + END
431 325  
432 326  
433   - 'QEXT': BEGIN ; Q - EXTINCTION CROSS SECTION
  327 + 'QEXT': BEGIN ;NOT UPDATED ; Q - EXTINCTION CROSS SECTION
434 328  
435 329 chi2_qext =total(((*!dustem_data.qext).values-dustem_qext)^2/(*!dustem_data.qext).sigma^2)
436 330 n_qext = n_elements((*!dustem_data.qext).values)
... ... @@ -474,7 +368,7 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
474 368  
475 369 END
476 370 ;
477   - 'UEXT': BEGIN
  371 + 'UEXT': BEGIN ;NOT UPDATED
478 372  
479 373  
480 374 chi2_uext =total(((*!dustem_data.uext).values-dustem_uext)^2/(*!dustem_data.uext).sigma^2)
... ...
src/idl/dustemcgwin_dataset.pro
  Diff comments   View file @fcb6ead
... ... @@ -22,7 +22,6 @@ degtorad = !pi/180 ;factor that the arctan will be devided by in order to have a
22 22  
23 23 clrs_plgns = ['Rosy Brown','Gold','Light Coral','Slate Blue','Dark Khaki','Salmon','Dark Green'] ;FOR NOW
24 24 scopes=tag_names((*!dustem_scope))
25   -
26 25 tgnms_extra = tag_names(_extra)
27 26 ind_xr = where(strmid(tgnms_extra,0,2) eq 'XR')
28 27 ind_yr = where(strmid(tgnms_extra,0,2) eq 'YR')
... ... @@ -98,8 +97,11 @@ if keyword_set(dataset) then begin
98 97 if not keyword_set(refresh) then begin
99 98 ;stop
100 99 xtit=textoidl('\lambda (\mum)')
101   - if !run_pol then xtit = ''
102   - cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,xtickformat='(A1)',color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
  100 + if !run_pol then begin
  101 + xtit = ''
  102 + xtickformat='(A1)'
  103 + endif else xtickformat='(F10.2)' ;MAYBE YOU'LL CHANGE THIS FORMAT TO EXPONENTIAL NOTATION
  104 + cgplot,wavs,wavs/wavs,/xlog,/ys,xs=1,pos=position,noerase=1,xtickformat=xtickformat,color='Black',xr=xr,xtit=xtit,yr=[0.0,2.0],yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.0
103 105 xyouts,pospltxt[0],pospltxt[1],textoidl('norm'),color=0,/normal,charsize=1.1
104 106  
105 107 ;endelse
... ... @@ -129,7 +131,7 @@ if keyword_set(dataset) then begin
129 131  
130 132 xx=((*!dustem_data.sed).wav)[idx_spec]
131 133 yy=dustem_interp[idx_spec]
132   - cgoplot,xx,yy,color='red',pos=position,psym=7,syms=2,noerase=1
  134 + cgoplot,xx,yy,color='Indian Red',pos=position,psym=7,syms=2,noerase=1
133 135 ENDIF
134 136  
135 137  
... ... @@ -315,7 +317,7 @@ if keyword_set(dataset) then begin
315 317  
316 318 xx=((*!dustem_data.polsed).wav)[idx_spec]
317 319 yy=dustem_interp[idx_spec]
318   - cgoplot,xx,yy,color='red',pos=position,psym=7,syms=2,noerase=1
  320 + cgoplot,xx,yy,color='Indian Red',pos=position,psym=7,syms=2,noerase=1
319 321 ENDIF
320 322  
321 323  
... ... @@ -467,7 +469,7 @@ if keyword_set(dataset) then begin
467 469  
468 470 xx=((*!dustem_data.polfrac).wav)[idx_spec]
469 471 yy=dustem_interp[idx_spec]
470   - cgoplot,xx,yy*100,color='red',pos=position,psym=7,syms=2,noerase=1
  472 + cgoplot,xx,yy*100,color='Indian Red',pos=position,psym=7,syms=2,noerase=1
471 473 ENDIF
472 474  
473 475  
... ... @@ -609,7 +611,7 @@ if keyword_set(dataset) then begin
609 611  
610 612 xx=((*!dustem_data.psi_em).wav)[idx_spec]
611 613 yy=dustem_interp[idx_spec]
612   - cgoplot,xx,yy,color='red',pos=position,psym=7,syms=2,noerase=1
  614 + cgoplot,xx,yy,color='Indian Red',pos=position,psym=7,syms=2,noerase=1
613 615 ENDIF
614 616  
615 617  
... ... @@ -844,7 +846,7 @@ if keyword_set(dataset) then begin
844 846 xx=((*!dustem_data.qsed).wav)[idx_spec]
845 847 yy=dustem_interp[idx_spec]
846 848 ;dustem_plot_mlog,xx,yy,ppositions=position,/overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog
847   - dustem_plot_mlog,xx,yy,ppositions=position,/positive_only, /overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
  849 + dustem_plot_mlog,xx,yy,ppositions=position,/positive_only, /overplot, noerase=1, color='Indian Red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
848 850 ENDIF
849 851  
850 852 IF ct_filt NE 0 THEN BEGIN
... ... @@ -895,7 +897,7 @@ if keyword_set(dataset) then begin
895 897 xx=((*!dustem_data.qsed).wav)[idx_spec]
896 898 yy=dustem_interp[idx_spec]
897 899 ;dustem_plot_mlog,xx,yy,ppositions=position,/overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog
898   - dustem_plot_mlog,xx,yy,ppositions=position,/negative_only, /overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
  900 + dustem_plot_mlog,xx,yy,ppositions=position,/negative_only, /overplot, noerase=1, color='Indian Red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
899 901 ENDIF
900 902  
901 903 IF ct_filt NE 0 THEN BEGIN
... ... @@ -1242,7 +1244,7 @@ if keyword_set(dataset) then begin
1242 1244 xx=((*!dustem_data.used).wav)[idx_spec]
1243 1245 yy=dustem_interp[idx_spec]
1244 1246 ;dustem_plot_mlog,xx,yy,ppositions=position,/overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog
1245   - dustem_plot_mlog,xx,yy,ppositions=position,/positive_only, /overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
  1247 + dustem_plot_mlog,xx,yy,ppositions=position,/positive_only, /overplot, noerase=1, color='Indian Red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
1246 1248 ENDIF
1247 1249  
1248 1250 IF ct_filt NE 0 THEN BEGIN
... ... @@ -1290,7 +1292,7 @@ if keyword_set(dataset) then begin
1290 1292 xx=((*!dustem_data.used).wav)[idx_spec]
1291 1293 yy=dustem_interp[idx_spec]
1292 1294 ;dustem_plot_mlog,xx,yy,ppositions=position,/overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog
1293   - dustem_plot_mlog,xx,yy,ppositions=position,/negative_only, /overplot, noerase=1, color='red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
  1295 + dustem_plot_mlog,xx,yy,ppositions=position,/negative_only, /overplot, noerase=1, color='Indian Red',psym=7,syms=2,xr=xr,/xlog,/ylog,ytickformat='(A1)',yr=yr
1294 1296 ENDIF
1295 1297  
1296 1298 IF ct_filt NE 0 THEN BEGIN
... ... @@ -1554,7 +1556,7 @@ if keyword_set(dataset) then begin
1554 1556 parameter_type = dustem_parameter_description2type(parameter_description,string_name=string_name)
1555 1557  
1556 1558  
1557   - IF STRUPCASE(strmid(strtrim(parameter_description,2),0,24)) eq '(*!DUSTEM_PARAMS).GRAINS' and STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O' then begin
  1559 + IF STRUPCASE(strmid(strtrim(parameter_description,2),0,18)) eq '(*!DUSTEM_PARAMS).' and STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O' then begin
1558 1560  
1559 1561 indpop=fix(STRUPCASE(strmid(strtrim(parameter_description,2),12,1,/reverse_offset)))
1560 1562 string_name=(((*!dustem_params). GRAINS).grain_type)[indpop]
... ... @@ -1580,6 +1582,38 @@ if keyword_set(dataset) then begin
1580 1582 k+=1
1581 1583  
1582 1584 ENDIF
  1585 +
  1586 +
  1587 +
  1588 + IF STRUPCASE(strmid(strtrim(parameter_description,2),0,18)) eq '(*!DUSTEM_PARAMS).' and ~(STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O') then begin
  1589 +
  1590 + ;indpop=fix(STRUPCASE(strmid(strtrim(parameter_description,2),12,1,/reverse_offset)))
  1591 + ;string_name=(((*!dustem_params). GRAINS).grain_type)[indpop]
  1592 +
  1593 + ;choosing a step of 0.03 in norm coordinates
  1594 + yypos = 0.95 - k*0.09
  1595 + xxpos = 0.0
  1596 +
  1597 + ;original line
  1598 + ;str=string(string_name+' = ',format=frmt0)+string(res[i],format=frmt1)
  1599 +
  1600 + str = strtrim(string(string_name+' = ',format=frmt0),2)+ string(res[i],format=frmt1) + textoidl(' \pm ') +string(errors(i),format=frmt1)
  1601 + ;will this work?
  1602 + cgtext, xxpos, yypos, str, charsize=-1, width=thiswidth, /normal
  1603 + xxpos = (1-thiswidth)/2
  1604 + str = strtrim(string(string_name+' = ',format=frmt0),2)
  1605 + cgtext, xxpos, yypos, str, charsize=-1, width=thiswidth, /normal
  1606 + xxpos+= thiswidth
  1607 +
  1608 + str= string(res[i],format=frmt1) + textoidl(' \pm ') + string(errors(i),format=frmt1)
  1609 + ;stop
  1610 + xyouts,xxpos,yypos,str,color=0,/normal,charsize=1.
  1611 + k+=1
  1612 +
  1613 + ENDIF
  1614 +
  1615 +
  1616 +
1583 1617  
1584 1618 ENDFOR
1585 1619  
... ... @@ -1597,7 +1631,7 @@ if keyword_set(dataset) then begin
1597 1631 parameter_type = dustem_parameter_description2type(parameter_description,string_name=string_name)
1598 1632  
1599 1633 ;stop
1600   - IF STRUPCASE(strmid(strtrim(parameter_description,2),0,24)) eq '(*!DUSTEM_PARAMS).GRAINS' and STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O' then begin
  1634 + IF STRUPCASE(strmid(strtrim(parameter_description,2),0,18)) eq '(*!DUSTEM_PARAMS).' and STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O' then begin
1601 1635  
1602 1636 indpop=fix(STRUPCASE(strmid(strtrim(parameter_description,2),12,1,/reverse_offset)))
1603 1637 string_name=(((*!dustem_params). GRAINS).grain_type)[indpop]
... ... @@ -1618,6 +1652,29 @@ if keyword_set(dataset) then begin
1618 1652 k+=1
1619 1653 ;stop
1620 1654 ENDIF
  1655 +
  1656 +
  1657 + IF STRUPCASE(strmid(strtrim(parameter_description,2),0,18)) eq '(*!DUSTEM_PARAMS).' and ~(STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O') then begin
  1658 +
  1659 + ;indpop=fix(STRUPCASE(strmid(strtrim(parameter_description,2),12,1,/reverse_offset)))
  1660 + ;string_name=(((*!dustem_params). GRAINS).grain_type)[indpop]
  1661 + ;stop
  1662 + ;choosing a step of 0.03 in norm coordinates
  1663 + yypos = 0.95 - k*0.09
  1664 + xxpos = 0.0
  1665 +
  1666 + str=strtrim(string(string_name+' = ',format=frmt0),2)
  1667 + ;stop
  1668 + str1=str+string(res[i],format=frmt1)+textoidl(' \pm ')+string(errors(i),format=frmt1)
  1669 + ;stop
  1670 + cgtext, xxpos, yypos, str1, charsize=-1, width=thiswidth, /normal
  1671 + xxpos=(1-thiswidth)/2
  1672 + ;stop
  1673 + xyouts,xxpos,yypos,str,color=0,/normal,charsize=1.
  1674 +
  1675 + k+=1
  1676 + ;stop
  1677 + ENDIF
1621 1678  
1622 1679 ENDFOR
1623 1680  
... ...
src/idl/dustemwrap_plot.pro
  Diff comments   View file @fcb6ead
... ... @@ -6,10 +6,17 @@ IF not keyword_set(st) THEN BEGIN
6 6 ;Activation of the plugins is needed because of the plotting of the total emission model that includes them.
7 7 dustem_activate_plugins,p_dim/(*(*!dustem_fit).param_init_values) ; 0/0 division case?
8 8 st=dustem_run(p_dim)
9   - ;should I add something here so that we can run this procedure without the keywords? Probably yes...
  9 + dustem_sed = dustem_compute_sed(p_dim,st,SED_spec)
  10 +
  11 + if !run_pol && !run_lin then begin
  12 +
  13 + dustem_polsed = dustem_compute_polsed(p_dim,st,P_spec,SP_spec,dustem_polfrac)
  14 + toto = dustem_compute_stokes(p_dim,st,dustem_qsed,dustem_used,Q_spec,U_spec,PSI_spec,dustem_psi_em)
  15 +
  16 + endif
  17 +
10 18 ENDIF
11 19  
12   -
13 20 ;if ~windowavailable(cgquery()) then begin
14 21  
15 22 ;LIST OF TESTS THAT NEED TO BE RAN SO THAT THE PLOTTING OF THE MODEL SPECTRA OCCURS
... ... @@ -48,7 +55,6 @@ cmdind_x = 0 ;for extinction
48 55 cmdind_prms = 0 ;for parameters
49 56 cmdind_plgns = 0 ;for plugins
50 57  
51   -
52 58 ;ADDING PLUGIN(S) TO SPECTRUM----------------
53 59  
54 60 ;these following indices will be used to change the keep the same information in the _extra structure
... ... @@ -71,7 +77,7 @@ degtorad = !pi/180 ;factor that the arctan will be devided by in order to have a
71 77  
72 78 ;errors = ((*!dustem_fit).current_param_errors)
73 79 if isa((*!dustem_fit).current_param_errors) then begin
74   - errors = (*(*!dustem_fit).current_param_errors)
  80 + errors = (*(*!dustem_fit).current_param_errors) * (*(*!dustem_fit).param_init_values)
75 81 endif else begin
76 82 errors = (*(*!dustem_fit).param_init_values)*0.+la_undef()
77 83 ENDELSE
... ... @@ -88,14 +94,21 @@ endif else begin
88 94 rchi2 = la_undef();maybe not the best initialization
89 95 endelse
90 96  
91   -if isa((*!dustem_fit).current_param_values) then begin
92   - res = (*(*!dustem_fit).current_param_values)
93   -endif else begin
94   - res = (*(*!dustem_fit).param_init_values)*0.+la_undef();maybe not the best initialization
95   - ;=(*(*!dustem_fit).param_init_values)*0.
96   -endelse
  97 +; if isa((*!dustem_fit).current_param_values) then begin
  98 +; res = (*(*!dustem_fit).current_param_values)
  99 +; endif else begin
  100 +; res = (*(*!dustem_fit).param_init_values)*0.+la_undef();maybe not the best initialization
  101 +; ;=(*(*!dustem_fit).param_init_values)*0.
  102 +; endelse
97 103  
98 104  
  105 +if isa(p_dim) then begin
  106 + res = p_dim
  107 +endif else begin
  108 +
  109 + res = (*(*!dustem_fit).param_init_values)*0.+la_undef();maybe not the best initialization
  110 +
  111 +endelse
99 112  
100 113 ;if iswinsed then begin ;I believe the iswin test should be below the case condition because I am not copying all the first 'HUGE LOOP' but rather replacing the positional arrays.
101 114 if test_m then begin
... ... @@ -511,383 +524,4 @@ ENDIF ELSE BEGIN ;(FIRST RUN)
511 524 ENDELSE
512 525  
513 526  
514   -;stop
515   -
516   -
517   -; ;#COMMENTED BITS OF CODES
518   -;
519   -; IF keyword_set(help) THEN BEGIN
520   -; doc_library,'dustem_plot_fit_sed'
521   -; goto,the_end
522   -; ENDIF
523   -;
524   -; IF keyword_set(ps) THEN BEGIN
525   -; set_plot, 'PS'
526   -; device, filename=ps, /color,set_character_size=[170,250], /encapsulated
527   -; ENDIF ELSE BEGIN
528   -;
529   -; set_plot,'X'
530   -; IF keyword_set(win) then window,win;,xsize=600,ysize=800
531   -; ENDELSE
532   -;
533   -; fact = 1.e4*(*!dustem_HCD)/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
534   -; ;use_col_data_filt=70
535   -; ;use_col_sed_spec=170
536   -; use_col_data_filt='blue'
537   -; ;use_col_sed_spec='red'
538   -; use_col_sed_spec='grey'
539   -; IF not keyword_set(col_sed) THEN BEGIN
540   -; ;use_col_sed_filt=250 ;red
541   -; use_col_sed_filt='red' ;red
542   -; ENDIF ELSE BEGIN
543   -; use_col_sed_filt=col_sed
544   -; ENDELSE
545   -; IF not keyword_set(col_tot) THEN BEGIN
546   -; ;use_col_tot=200
547   -; use_col_tot='black'
548   -; ENDIF ELSE BEGIN
549   -; use_col_tot=col_tot
550   -; ENDELSE
551   -; IF not keyword_set(line_tot) THEN BEGIN
552   -; use_line_tot=0
553   -; ENDIF ELSE BEGIN
554   -; use_line_tot=line_tot
555   -; ENDELSE
556   -;
557   -;
558   -; spec = st.sed.em_tot * fact
559   -; if keyword_set(fpol) then specpol = st.polsed.em_tot * fact
560   -; ;ADDING PLUGIN(S) TO SPECTRUM----------------
561   -; ;if n_tags(!dustem_data.sed) gt 1 then begin
562   -; scopes=tag_names((*!dustem_scope))
563   -; IF scopes[0] NE 'NONE' THEN BEGIN
564   -; ;IF ptr_valid(!dustem_plugin) THEN BEGIN
565   -; FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
566   -; IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_SED') THEN spec+=(*(*!dustem_plugin).(i))[*,0]
567   -; if keyword_set(fpol) then begin
568   -; IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'REPLACE_QSED') THEN specpol=sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)
569   -; endif
570   -; ENDFOR
571   -; FOR i=0L,n_tags(*!dustem_scope)-1 DO BEGIN
572   -; if keyword_set(fpol) then begin
573   -; IF total(strsplit((*(*!dustem_scope).(i)),'+',/extract) EQ 'ADD_QSED') THEN specpol+=sqrt(((*(*!dustem_plugin).(i))[*,1])^2+((*(*!dustem_plugin).(i))[*,2])^2)
574   -; endif
575   -; ENDFOR
576   -; ENDIF
577   -; ;endif
578   -; ;------------------------------------------
579   -;
580   -; ;stop
581   -;
582   -; ;use_cols=[use_col_pah,use_col_vsg,use_col_bg,use_col_cont]
583   -; ;use_lines=[use_line_pah,use_line_vsg,use_line_bg,use_line_cont]
584   -; ;col_off=30
585   -; ;Ngrains=(*!dustem_params).grain.Ngrains
586   -; Ngrains=(*!dustem_params).Ngrains
587   -; ;use_cols=long(findgen(Ngrains)/(Ngrains-1)*(255-col_off)+col_off)
588   -; use_cols=dustem_grains_colors(Ngrains,/cgplot)
589   -; use_lines=replicate(0,Ngrains)
590   -;
591   -; norm = dustem_sed * 0. + 1
592   -;
593   -; ;====== PLOT THE SED
594   -;
595   -; IF keyword_set(title) THEN title = title ELSE title = 'Spectral Energy Distribution (Running)'
596   -;
597   -; IF keyword_set(xr) THEN xr = xr ELSE xr = [1.00E+00,6.00E+04]
598   -;
599   -; IF keyword_set(yr) THEN yr = yr ELSE yr = [1.00E-7,5.00E03]
600   -;
601   -; IF keyword_set(xtit) THEN xtit = xtit ELSE xtit = textoidl('\lambda (\mum)')
602   -;
603   -; IF keyword_set(ytit) THEN ytit = ytit ELSE ytit = textoidl('Brightness/N_H (MJy/sr/H)')
604   -;
605   -;
606   -; ;deffo define a title variable here. The procedures are not communicating with each other.
607   -; ;pos=cgLayout()
608   -; ;###############################
609   -; if keyword_set(fpol) then begin
610   -; if !run_lin then begin
611   -; cgDisplay, 600, 500
612   -; cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit='',/ylog,/xlog,/ys,/xs,position=[0.12,0.25,0.96,0.76],xtickformat='(A1)',_extra=_extra,charsize=1.3,/noerase
613   -; endif
614   -; endif ELSE cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=ytit,tit=title,/ylog,/xlog,/ys,/xs,position=[0.12,0.35,0.96,0.90],xtickformat='(A1)',_extra=_extra,charsize=1.3
615   -; ;###############################
616   -;
617   -;
618   -; ;cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/norm,/nodata,xtit='',ytit=textoidl('Brightness/N_H (MJy/sr/H)'),tit=title,ylog=1,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.12,0.35,0.96,0.90],xtickformat='(A1)'
619   -;
620   -;
621   -; ;stop
622   -; ;cgplot,st.polsed.wav,Q_sed*fact,xtit='',ytit=ytitstq,tit=titstq,ylog=ylog,/xlog,xr=xr,yr=yr,/ys,/xs,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)'
623   -;
624   -; ;this is where you will add the normalized sed
625   -;
626   -;
627   -; ind_filt=where((*!dustem_data.sed).filt_names NE 'SPECTRUM',count_filt)
628   -; ind_spec=where((*!dustem_data.sed).filt_names EQ 'SPECTRUM',count_spec)
629   -; ;=== Plot the data
630   -; ; following lines for fawlty compatibility
631   -; defsysv,'!psym',exists=pexist
632   -; if pexist eq 0 then defsysv,'!psym',0
633   -; IF count_spec NE 0 THEN BEGIN
634   -; ;x&
635   -; xx=((*!dustem_data.sed).wav)[ind_spec]
636   -; yy=((*!dustem_data.sed).values)[ind_spec]/norm[ind_spec]
637   -; rms=3.*((*!dustem_data.sed).sigma)[ind_spec]/2./norm[ind_spec]
638   -; cgoplot,xx,yy,psym ,syms=0.5,color=use_col_sed_spec
639   -; IF not keyword_set(no_spec_error) THEN BEGIN
640   -; cgerrplot,xx,yy-rms,yy+rms,color=use_col_sed_spec
641   -; ;err_bar,((*!dustem_data.sed).wav)(ind_spec),((*!dustem_data.sed).values)(ind_spec)/norm(ind_spec),yrms=3.*((*!dustem_data.sed).sigma)(ind_spec)/2./norm(ind_spec)
642   -; ENDIF
643   -; ENDIF
644   -; IF count_filt NE 0 THEN BEGIN
645   -; ;stop
646   -; xx=((*!dustem_data.sed).wav)[ind_filt]
647   -; yy=((*!dustem_data.sed).values)[ind_filt]/norm[ind_filt]
648   -; rms=3.*((*!dustem_data.sed).sigma)[ind_filt]/2./norm[ind_filt]
649   -; plotsym,0,/fill
650   -; cgoplot,xx,yy,psym=8,color='Dodger Blue';use_col_data_filt
651   -; ;err_bar,((*!dustem_data.sed).wav)(ind_filt),((*!dustem_data.sed).values)(ind_filt)/norm(ind_filt),yrms=3.*((*!dustem_data.sed).sigma)(ind_filt)/2./norm(ind_filt),color=use_col_data_filt
652   -; cgerrplot,xx,yy-rms,yy+rms,color='Dodger Blue';use_col_data_filt
653   -; ENDIF
654   -; ;=== Plot the computed SED
655   -; IF count_filt NE 0 THEN BEGIN
656   -; plotsym,8
657   -; xx=((*!dustem_data.sed).wav)[ind_filt]
658   -; yy=sed[ind_filt]/norm[ind_filt]
659   -; cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2
660   -; ENDIF
661   -; IF count_spec NE 0 THEN BEGIN
662   -; plotsym,0
663   -; xx=((*!dustem_data.sed).wav)[ind_spec]
664   -; yy=sed[ind_spec]/norm[ind_spec]
665   -; cgoplot,xx,yy,color=use_col_sed_filt,psym=6,syms=2
666   -; ENDIF
667   -;
668   -;
669   -; IF !dustem_show_plot EQ 2 THEN BEGIN
670   -; norm = spec
671   -; ENDIF ELSE BEGIN
672   -; norm = spec * 0. + 1
673   -; ENDELSE
674   -; use_cols[1]='Cornflower'
675   -; FOR i=0L,Ngrains-1 DO BEGIN
676   -; cgoplot,st.sed.wav,st.sed.(i+1)*fact/norm,color=use_cols[i],linestyle=use_lines[i]
677   -; ENDFOR
678   -;
679   -; ;stop
680   -;
681   -; ;PLOTTING OF THE PLUGIN(S)--------------- AUTOMATE THIS. QUITE FEASIBLE
682   -; IF tag_exist(*!dustem_scope,'CONTINUUM') THEN BEGIN
683   -; cgoplot,st.sed.wav,(*(*!dustem_plugin).continuum)[*,0],color='Teal',linestyle=3
684   -; ENDIF
685   -; IF tag_exist(*!dustem_scope,'FREEFREE') THEN BEGIN
686   -; cgoplot,st.sed.wav,(*(*!dustem_plugin).freefree)[*,0],color='Dark Red',linestyle=3
687   -; ENDIF
688   -; IF tag_exist(*!dustem_scope,'SYNCHROTRON') THEN BEGIN
689   -; cgoplot,st.sed.wav,(*(*!dustem_plugin).synchrotron)[*,0],color='Crimson',linestyle=3
690   -; ENDIF
691   -;
692   -; IF tag_exist(*!dustem_scope,'MBBDY_ISRF') THEN BEGIN
693   -; cgoplot,st.sed.wav,(*(*!dustem_plugin).mbbdy_isrf)[*,0],color='Gold',linestyle=3
694   -; ENDIF
695   -;
696   -; IF tag_exist(*!dustem_scope,'MBBDY') THEN BEGIN
697   -; cgoplot,st.sed.wav,(*(*!dustem_plugin).mbbdy)[*,0],color='Cornflower',linestyle=3
698   -; ENDIF
699   -;
700   -;
701   -; cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
702   -; if keyword_set(fpol) then cgoplot,st.polsed.wav,specpol,color='Teal',linestyle=4
703   -; ;plot the normalized data as well.
704   -; ;----------------------------------------
705   -;
706   -; ;==== print the legend
707   -; frmt0='(A36)'
708   -; frmt1='(1E10.2)'
709   -; frmt2='(F7.2)'
710   -; use_legend_xpos=0.50
711   -; if keyword_set(fpol) then use_legend_ypos=0.70 ELSE use_legend_ypos=0.84
712   -; use_legend_offset=0.03 ;This is the offset between lines of the legend in normalized units
713   -; legend_charsize=0.86
714   -;
715   -; k=0.
716   -; iscond=0
717   -; n_plgns = n_tags(*!dustem_scope)
718   -; inn=n_plgns/2+1
719   -; IF n_plgns mod 2 ne 0 THEN inn=n_plgns/2+2
720   -;
721   -; IF keyword_set(legend_xpos) THEN use_legend_xpos=legend_xpos
722   -; IF keyword_set(legend_ypos) THEN use_legend_ypos=legend_ypos
723   -; IF keyword_set(legend_offset) THEN use_legend_offset=legend_offset
724   -;
725   -; IF !d.name NE 'PS' THEN cleanplot
726   -;
727   -; tg=0
728   -; prv_str=''
729   -; IF keyword_set(res) THEN BEGIN
730   -; Npar=n_elements(res)
731   -; ;print,'=============='
732   -;
733   -; FOR i=0L,Npar-1 DO BEGIN
734   -; parameter_description=(*(*!dustem_fit).param_descs)[i]
735   -; parameter_type=dustem_parameter_description2type(parameter_description,string_name=string_name)
736   -; str=string(string_name+' = ',format=frmt0)+string(res[i],format=frmt1)
737   -; IF keyword_set(errors) THEN BEGIN
738   -; str=str+textoidl(' \pm ')+string(errors(i),format=frmt1)
739   -; ENDIF
740   -;
741   -; xxpos=use_legend_xpos*0.07
742   -; yypos=use_legend_ypos*0.03
743   -; xyouts,xxpos,yypos,'Model: '+!dustem_model,color=0,/normal,charsize=legend_charsize
744   -;
745   -; IF STRUPCASE(strmid(strtrim(parameter_description,2),0,6)) eq 'DUSTEM' THEN BEGIN
746   -;
747   -; ; xxpos=use_legend_xpos*0.07
748   -; ; yypos=use_legend_ypos*0.03
749   -; ; xyouts,xxpos,yypos,'Model: '+!dustem_model,color=0,/normal,charsize=legend_charsize
750   -;
751   -; xxpos=use_legend_xpos*0.75
752   -; yypos=use_legend_ypos*1.03
753   -; xyouts,xxpos,yypos,'--Plugins--',color=0,/normal,charsize=legend_charsize
754   -;
755   -; ii = strsplit(string_name,'_',count=countx) & ii = ii(countx-1)-1 ; Locating the last underscore to automate the extraction of the plugin's keyword
756   -;
757   -; mm=where(tag_names(*!dustem_scope) eq strupcase(strmid(string_name,7,ii-7)),coun) ; Selecting a plugin through matching the string name of the plugin form the scope system variable with the one read from the parameter description vector
758   -;
759   -; ;tg+=1
760   -; ;if strmid(string_name,7,ii-7) eq prv_str then tg-=1
761   -; prmtg=(*(*!dustem_paramtag).(mm))
762   -; indtg=(strmid(string_name,ii+1)) & indtg=strmid(indtg,0,/reverse_offset)
763   -; ;stop
764   -; indtg=fix(indtg)
765   -; prmtg=prmtg[indtg-1]
766   -;
767   -; str=string(strmid(string_name,7,ii-7)+' ['+strmid(string_name,ii+1)+']: '+prmtg+' = ',format=frmt0)+string(res[i],format=frmt1)
768   -;
769   -; prv_str=strmid(string_name,7,ii-7)
770   -;
771   -;
772   -;
773   -; IF keyword_set(errors) THEN BEGIN
774   -; str=str+textoidl(' \pm ')+string(errors(i),format=frmt1)
775   -; ENDIF
776   -;
777   -; xxpos=use_legend_xpos*0.46
778   -;
779   -; IF n_plgns gt 3 THEN BEGIN ;adapt the display of the plugins if there are more than three (display two columns)
780   -; IF i eq inn+k THEN BEGIN
781   -; k+=inn
782   -; ;iscond=1
783   -; inn+=1
784   -; xxpos=use_legend_xpos*0.06
785   -; ENDIF
786   -; ENDIF
787   -;
788   -; ; IF iscond THEN BEGIN
789   -; ; inn+=1
790   -; ; xxpos=use_legend_xpos*0.06
791   -; ; ENDIF
792   -;
793   -; yypos=use_legend_ypos-(i-k)*use_legend_offset
794   -; xyouts,xxpos,yypos,str,color=0,/normal,charsize=legend_charsize
795   -; endif else begin
796   -; IF STRUPCASE(strmid(strtrim(parameter_description,2),0,24)) eq '(*!DUSTEM_PARAMS).GRAINS' and STRUPCASE(strmid(strtrim(parameter_description,2),3,1,/reverse_offset)) EQ 'O' then begin
797   -; indpop=fix(STRUPCASE(strmid(strtrim(parameter_description,2),12,1,/reverse_offset)))
798   -; ;oo+=1
799   -; ;string_name=(((*!dustem_params).GRAINS).grain_type)[oo-1]
800   -; string_name=(((*!dustem_params).GRAINS).grain_type)[indpop]
801   -; ;stop
802   -; endif
803   -; k+=1
804   -; xxpos=use_legend_xpos*1.5
805   -; yypos=use_legend_ypos*1.03
806   -; xyouts,xxpos,yypos,'--Parameters--',color=0,/normal,charsize=legend_charsize
807   -; yypos=use_legend_ypos-(i)*use_legend_offset
808   -; xxpos=use_legend_xpos*1.14
809   -; str=string(string_name+' = ',format=frmt0)+string(res[i],format=frmt1)
810   -; xyouts,xxpos,yypos,str,color=0,/normal,charsize=legend_charsize
811   -; endelse
812   -; ENDFOR
813   -; ;stop
814   -; ENDIF
815   -; IF keyword_set(chi2) THEN BEGIN
816   -; xxpos=1.29*use_legend_xpos
817   -; yypos=1.13*0.84;use_legend_ypos
818   -; xyouts,xxpos,yypos,string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0,/normal,charsize=legend_charsize
819   -; ENDIF
820   -; IF keyword_set(rchi2) THEN BEGIN
821   -; xxpos=1.29*use_legend_xpos
822   -; yypos=1.13*0.84-use_legend_offset;use_legend_ypos-use_legend_offset
823   -; xyouts,xxpos,yypos,string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0,/normal,charsize=legend_charsize
824   -; ENDIF
825   -;
826   -; xtit=textoidl('\lambda (\mum)')
827   -; ;cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,xr=xr,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.14,0.96,0.35],/noerase,yticks=2,ymino=2,xticklen=0.1
828   -; ;stop
829   -;
830   -; IF keyword_set(_extra) THEN BEGIN
831   -; extra_kept={XR:[0.,0.]}
832   -; extra_tags=tag_names(_extra)
833   -; ind=where(extra_tags EQ 'XR',count)
834   -; IF count NE 0 THEN extra_kept.XR=_extra.(ind[0]) ;ELSE extra_kept=0
835   -; ENDIF
836   -;
837   -; ;stop
838   -; ;,position=[0.12,0.14,0.96,0.35] old position before your polfrac display modification
839   -; if keyword_set(fpol) then begin
840   -; cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,_extra=extra_kept,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.1,0.96,0.25],/noerase,yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3
841   -; endif else cgplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,_extra=extra_kept,/nodata,xtit=xtit,ytit='Normalized',tit='',/xlog,/ys,/xs,yr=[0,2],ylog=0,position=[0.12,0.14,0.96,0.35],/noerase,yticks=2,ymino=2,xticklen=0.1,ytickformat='(F6.2)',charsize=1.3
842   -; ;plot the normalized data as well.
843   -; IF count_spec NE 0 THEN BEGIN
844   -; xx=((*!dustem_data.sed).wav)[ind_spec]
845   -; yy=((*!dustem_data.sed).values/sed)[ind_spec]
846   -; cgoplot,xx,yy,psym=16,symsize=1,thick=2,color=use_col_sed_spec
847   -; trois_sigma=(3.*((*!dustem_data.sed).sigma)/2./sed)[ind_spec]
848   -; cgerrplot,xx,yy-trois_sigma,yy+trois_sigma,color=use_col_sed_spec
849   -; ENDIF
850   -; IF count_filt NE 0 THEN BEGIN
851   -; xx=((*!dustem_data.sed).wav)[ind_filt]
852   -; yy=((*!dustem_data.sed).values/sed)[ind_filt]
853   -; cgoplot,xx,yy,psym=16,symsize=1,thick=2,color='Dodger Blue'
854   -; trois_sigma=(3.*((*!dustem_data.sed).sigma)/2./sed)[ind_filt]
855   -; cgerrplot,xx,yy-trois_sigma,yy+trois_sigma,color='Dodger Blue'
856   -; ENDIF
857   -;
858   -; ;stop
859   -; ;cgoplot,(*!dustem_data.sed).wav,(*!dustem_data.sed).values/sed,psym=16,symsize=1,thick=2,color='Dodger Blue'
860   -; cgoplot,10^!x.crange,replicate(1.,2),color='black',linestyle=0
861   -; ;cgoplot,st.sed.wav,spec/spec,color='black'
862   -; ;cgerrplot,((*!dustem_data.sed).wav),((*!dustem_data.sed).values)/sed-3.*((*!dustem_data.sed).sigma)/2./sed,((*!dustem_data.sed).values)/sed+3.*((*!dustem_data.sed).sigma)/2./sed,color='Dodger Blue'
863   -;
864   -;
865   -; if keyword_set(fpol) then begin
866   -; plotsym,0,/fill
867   -; cgplot, (*!dustem_data.polfrac).wav,(*!dustem_data.polfrac).values*100,/nodata,xtit='',tit=title,ytit='',/xlog,/ylog,/ys,/xs,position=[0.12,0.76,0.96,0.91],charsize=1.3,yr=[1.0E-1,30.0],/noerase,xticklen=0.1,xtickformat='(A1)',_extra=extra_kept;,xr=_extra.xr
868   -; xxpos=use_legend_xpos*0.35
869   -; yypos=use_legend_ypos*1.25
870   -; xyouts,xxpos,yypos,'Polarization fraction (%)',color=cgcolor('purple'),charsize=legend_charsize,/normal
871   -; yypos=use_legend_ypos*1.20
872   -; stringg='Polarization SED '+'('+textoidl('P_{\nu}')+')'
873   -; xyouts,xxpos,yypos,stringg,color=cgcolor('teal'),charsize=legend_charsize,/normal
874   -; pilotfx=fpol(0);((*!dustem_data.polfrac).values)(0);fpol(0);
875   -; cgoplot, (*!dustem_data.polfrac).wav,(*!dustem_data.polfrac).values*100,charsize=1.3,psym=8,syms=1,thick=2,color='Dodger Blue'
876   -; cgoplot, st.polsed.wav,(specpol/spec)*100/pilotfx*((*!dustem_data.polfrac).values)(0),color='purple'
877   -; cgoplot, (*!dustem_data.polfrac).wav,fpol*100/pilotfx*((*!dustem_data.polfrac).values)(0),color='Red',psym=6,syms=1
878   -; ;fpol=abs(fpol)
879   -; ;delp=(*!dustem_data.polfrac).values*100-fpol*100/pilotfx*((*!dustem_data.polfrac).values)(0)
880   -; ;cgoplot, ((*!dustem_data.polfrac).wav),((*!dustem_data.polfrac).values*100+delp),charsize=1.3,color='black',linestyle=2
881   -; ;stop
882   -; endif
883   -;
884   -; IF keyword_set(ps) THEN BEGIN
885   -; device,/close
886   -; set_plot,'X'
887   -; ENDIF
888   -;
889   -; ;stop
890   -;
891   -; the_end:
892   -
893 527 END
... ...