Commit 9be94157d7de369f5991d52158041ab140c1b9d3

Authored by Jean-Philippe Bernard
1 parent 39a73c0d
Exists in master

improved

src/idl/dustem_fit_sed_readme.pro
... ... @@ -165,17 +165,15 @@ CASE !dustem_which OF
165 165 ; dustem_init_fixed_params,fpd,fiv
166 166 END
167 167 'COMPIEGNE_ETAL2010':BEGIN
168   -; PUTTING G0 last DOES NOT WORK BECAUSE OF ORDER ISSUE
169 168 pd = [ $
170   - '(*!dustem_params).G0', $ ;G0
171   - '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
  169 + '(*!dustem_params).gas.G0', $ ;G0
172 170 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction
  171 + '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction
173 172 '(*!dustem_params).grains(2).mdust_o_mh', $ ;amCBEx
174 173 '(*!dustem_params).grains(3).mdust_o_mh', $ ;amCBEx
175 174 '(*!dustem_params).grains(4).mdust_o_mh', $ ;aSil
176 175 'dustem_create_continuum_2'] ;Intensity of NIR continuum
177 176 iv = [1.0, 7.8e-4, 7.8e-4,1.65e-4,1.45e-3,7.8e-3,0.001]
178   -; iv = [7.8e-4, 7.7e-4,1.65e-4,1.45e-3,7.8e-3,0.001]
179 177 Npar=n_elements(pd)
180 178 ulimed=replicate(0,Npar)
181 179 llimed=replicate(1,Npar)
... ... @@ -243,10 +241,23 @@ tol=1.e-14
243 241 use_Nitermax=2 ;maximum number of iteration. This is the criterium which will stop the fit procedure
244 242 IF keyword_set(itermax) THEN use_Nitermax=itermax
245 243  
246   -loadct,13
247   -!y.range=[1e-4,10] ;This is to ajust plot range from outside the routine
  244 +;loadct,13
  245 +;!y.range=[1e-4,10] ;This is to ajust plot range from outside the routine
  246 +;t1=systime(0,/sec)
  247 +;res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol,/xlog,/ylog)
  248 +;t2=systime(0,/sec)
  249 +
248 250 t1=systime(0,/sec)
249   -res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol,/xlog,/ylog)
  251 +
  252 +yr=[1e-4,10]
  253 +xr=[1,1.e5]
  254 +legend_xpos=0.6
  255 +legend_ypos=0.8
  256 +tit='DUSTEMWrap Intensity SED (Running)'
  257 +ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
  258 +xtit=textoidl('\lambda (\mum)')
  259 +res=dustem_mpfit_data(tol=tol,Nitermax=use_Nitermax,gtol=gtol,/xlog,/ylog,xr=xr,yr=yr,xtit=xtit,ytit=ytit,title=tit,legend_xpos=legend_xpos,legend_ypos=legend_ypos,errors=errors,chi2=chi2,rchi2=rchi2)
  260 +
250 261 t2=systime(0,/sec)
251 262  
252 263 ;=== SAVE FIT RESULTS
... ... @@ -264,8 +275,6 @@ file='/tmp/DUSTEM_fit_example.sav'
264 275 dustem_restore_system_variables,file
265 276  
266 277 ;=== Plot best fit
267   -yr=[1e-4,10]
268   -xr=[1,2000]
269 278 tit='DUSTEM Intensity SED fit Example (Saved)'
270 279 ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2')
271 280 xtit=textoidl('\lambda (\mum)')
... ... @@ -289,8 +298,6 @@ IF keyword_set(postcript) THEN BEGIN
289 298 ps_file=!dustem_wrap_soft_dir+'/Docs/Figures/'+'Last_dustem_fit.ps'
290 299 device,filename=ps_file,/color
291 300 ENDIF
292   -legend_xpos=2.
293   -legend_ypos=1.
294 301  
295 302 ;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
296 303 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
... ...
src/idl/dustem_mpfit_run.pro
... ... @@ -40,10 +40,8 @@ FUNCTION dustem_mpfit_run,x,p_min,_EXTRA=extra
40 40  
41 41 ;note that x is not used.
42 42  
43   -;COMMON MPFIT_ERROR, ERROR_CODE
44   -
45 43 IF keyword_set(help) THEN BEGIN
46   - doc_library,'dustem_mpfit_run_vg'
  44 + doc_library,'dustem_mpfit_run'
47 45 sed=0.
48 46 goto,the_end
49 47 ENDIF
... ... @@ -126,13 +124,14 @@ FOR i_tag=0,count_data_tag-1 DO BEGIN
126 124 / total(dustem_sed^2/(*!dustem_data.sed).sigma^2)
127 125  
128 126 print,"Best SED chi2 obtained if f_HI = ",alpha, " or if we multiply model by ",beta
  127 + (*!dustem_fit).chi2=chi2_sed
  128 + (*!dustem_fit).rchi2=rchi2_sed
129 129  
130   - ;dustem_plot_nuinu_em,st,dustem_sed,cont,/print_radiance,p_dim=p_dim
131   -
132 130 ; Plot if needed
133 131 IF !dustem_show_plot NE 0 THEN BEGIN
134 132 window,0
135 133 dustem_plot_fit_sed,st,dustem_sed,cont,freefree,synchrotron,_extra=extra,res=p_min*(*(*!dustem_fit).param_init_values),chi2=(*!dustem_fit).chi2,rchi2=(*!dustem_fit).rchi2
  134 + ;stop
136 135 ENDIF
137 136  
138 137 dustem_out = [ dustem_out, dustem_sed ]
... ...
src/idl/dustem_plot_fit_sed.pro
... ... @@ -2,7 +2,7 @@ PRO dustem_plot_fit_sed,st,sed,cont,freefree,synchrotron, $
2 2 res=res,errors=errors,chi2=chi2,rchi2=rchi2, $
3 3 no_spec_error=no_spec_error, $
4 4 help=help,win=win,ps=ps, $
5   - legend_xpos=legend_xpos,legend_ypos=legend_ypos, $
  5 + legend_xpos=legend_xpos,legend_ypos=legend_ypos,legend_offset=legend_offset, $
6 6 _extra=_extra
7 7  
8 8 ;+
... ... @@ -67,7 +67,8 @@ fact = 1.e4*1.E20/(4.*!pi)/(3.e8/1.e-6/st.sed.wav)*1.e20/1.e7
67 67 ;use_col_data_filt=70
68 68 ;use_col_sed_spec=170
69 69 use_col_data_filt='blue'
70   -use_col_sed_spec='red'
  70 +;use_col_sed_spec='red'
  71 +use_col_sed_spec='grey'
71 72 IF not keyword_set(col_sed) THEN BEGIN
72 73 ;use_col_sed_filt=250 ;red
73 74 use_col_sed_filt='red' ;red
... ... @@ -136,7 +137,6 @@ col_off=30
136 137 Ngrains=(*!dustem_params).Ngrains
137 138 ;use_cols=long(findgen(Ngrains)/(Ngrains-1)*(255-col_off)+col_off)
138 139 use_cols=dustem_grains_colors(Ngrains,/cgplot)
139   -
140 140 use_lines=replicate(0,Ngrains)
141 141  
142 142 norm = sed * 0. + 1
... ... @@ -168,9 +168,9 @@ IF count_spec NE 0 THEN BEGIN
168 168 xx=((*!dustem_data.sed).wav)[ind_spec]
169 169 yy=((*!dustem_data.sed).values)[ind_spec]/norm[ind_spec]
170 170 rms=3.*((*!dustem_data.sed).sigma)[ind_spec]/2./norm[ind_spec]
171   - cgoplot,xx,yy,psym=8,_extra=extra,syms=0.5
  171 + cgoplot,xx,yy,psym=8,_extra=extra,syms=0.5,color=use_col_sed_spec
172 172 IF not keyword_set(no_spec_error) THEN BEGIN
173   - cgerrplot,xx,yy-rms,yy+rms
  173 + cgerrplot,xx,yy-rms,yy+rms,color=use_col_sed_spec
174 174 ;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)
175 175 ENDIF
176 176 ENDIF
... ... @@ -182,7 +182,7 @@ IF count_filt NE 0 THEN BEGIN
182 182 plotsym,0,/fill
183 183 cgoplot,xx,yy,psym=8,_extra=extra,color=use_col_data_filt
184 184 ;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
185   - cgerrplot,xx,yy-rms,yy+rms
  185 + cgerrplot,xx,yy-rms,yy+rms,color=use_col_data_filt
186 186 ENDIF
187 187 ;=== Plot the computed SED
188 188 IF count_filt NE 0 THEN BEGIN
... ... @@ -220,7 +220,7 @@ ENDIF
220 220 cgoplot,st.sed.wav,spec/norm,color=use_col_tot,linestyle=use_line_tot
221 221  
222 222 ;==== print the legend
223   -frmt0='(A25)'
  223 +frmt0='(A26)'
224 224 frmt1='(1E10.2)'
225 225 frmt2='(F7.2)'
226 226 ;WE have to find a way to pass this through keywords ...
... ... @@ -230,10 +230,11 @@ frmt2='(F7.2)'
230 230 ;use_legend_ypos=2.5
231 231 use_legend_xpos=0.8
232 232 use_legend_ypos=0.8
233   -offset=0.05
  233 +use_legend_offset=0.03 ;This is the offset between lines of the legend in normalized units
234 234 legend_charsize=1.0
235 235 IF keyword_set(legend_xpos) THEN use_legend_xpos=legend_xpos
236 236 IF keyword_set(legend_ypos) THEN use_legend_ypos=legend_ypos
  237 +IF keyword_set(legend_offset) THEN use_legend_offset=legend_offset
237 238  
238 239 IF !d.name NE 'PS' THEN cleanplot
239 240 IF keyword_set(res) THEN BEGIN
... ... @@ -241,39 +242,46 @@ IF keyword_set(res) THEN BEGIN
241 242 FOR i=0L,Npar-1 DO BEGIN
242 243 ;stop
243 244 sstr=(str_sep((*(*!dustem_fit).param_descs)[i],'(*!dustem_params).'))
  245 + ;hprint,sstr
  246 + ;stop
244 247 ;sstr=strsplit((*(*!dustem_fit).param_descs)[i],'(*!dustem_params).',/extract,count=count,/preserve_null)
245 248 IF n_elements(sstr) EQ 2 AND sstr[0] EQ '' THEN BEGIN ;case of a parameter description starting with (*!dustem_params).
246 249 ssstr=sstr[n_elements(sstr)-1]
247 250 str=string(ssstr+'=',format=frmt0)+string(res(i),format=frmt1)
248   - ;stop
249 251 IF keyword_set(errors) THEN BEGIN
250 252 str=str+textoidl('\pm')+string(errors(i),format=frmt1)
251 253 ENDIF
252   - yypos=use_legend_ypos-(i+1)*offset
253   - ;stop
  254 + yypos=use_legend_ypos-(i+1)*use_legend_offset
  255 + xyouts,use_legend_xpos,yypos,str,color=0,/normal,charsize=legend_charsize
  256 + ENDIF
  257 + IF n_elements(sstr) EQ 1 THEN BEGIN ;pluggins parameters
  258 + ssstr=sstr[n_elements(sstr)-1]
  259 + str=string(ssstr+'=',format=frmt0)+string(res(i),format=frmt1)
  260 + IF keyword_set(errors) THEN BEGIN
  261 + str=str+textoidl('\pm')+string(errors(i),format=frmt1)
  262 + ENDIF
  263 + yypos=use_legend_ypos-(i+1)*use_legend_offset
254 264 xyouts,use_legend_xpos,yypos,str,color=0,/normal,charsize=legend_charsize
255 265 ENDIF
256 266 sstr=(str_sep((*(*!dustem_fit).param_descs)[i],'!dustem_')) ;case of a parameter description starting with !dustem_.
257 267 IF n_elements(sstr) EQ 2 AND sstr[0] EQ '' THEN BEGIN ;case of a parameter description starting with (*!dustem_params).
258 268 ssstr=sstr[n_elements(sstr)-1]
259 269 str=string(ssstr+'=',format=frmt0)+string(res(i),format=frmt1)
260   - ;stop
261 270 IF keyword_set(errors) THEN BEGIN
262 271 str=str+textoidl('\pm')+string(errors(i),format=frmt1)
263 272 ENDIF
264   - yypos=use_legend_ypos-(i+1)*offset
265   - ;stop
  273 + yypos=use_legend_ypos-(i+1)*use_legend_offset
266 274 xyouts,use_legend_xpos,yypos,str,color=0,/normal,charsize=legend_charsize
267 275 ENDIF
268 276 ;stop
269 277 ENDFOR
270 278 ENDIF
271 279 IF keyword_set(chi2) THEN BEGIN
272   - yypos=use_legend_ypos-(Npar+1)*offset
  280 + yypos=use_legend_ypos-(Npar+1)*use_legend_offset
273 281 xyouts,use_legend_xpos,yypos,string('chi2=',format=frmt0)+string(chi2,format=frmt2),color=0,/normal,charsize=legend_charsize
274 282 ENDIF
275 283 IF keyword_set(rchi2) THEN BEGIN
276   - yypos=use_legend_ypos-(Npar+2)*offset
  284 + yypos=use_legend_ypos-(Npar+2)*use_legend_offset
277 285 xyouts,use_legend_xpos,yypos,string('red. chi2=',format=frmt0)+string(rchi2,format=frmt2),color=0,/normal,charsize=legend_charsize
278 286 ENDIF
279 287  
... ...