PRO dustem_fit_sed_extsed_readme,postscript=postscript,model=model,help=help,itermax=itermax ;This Readme describes how to fit total intensity and polarization SEDs jointly ;Obviously, the dustem package must have been installed succesfully ;see DustemWrap documentation for install instructions. ;+ ; NAME: ; dustem_fit_sed_polsed_readme ; PURPOSE: ; This is an example of how to fit total intensity and polarization SEDs jointly ; It is meant to be an example to follow when writing your own ; programs using the dustem IDL wrapper. ; CATEGORY: ; Dustem ; CALLING SEQUENCE: ; dustem_fit_sed_polsed_readme,postcript=postcript,mode=mode,help=help ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: ; None ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: ; postcript = if set plot is done in DUSTEM/Docs/Figures/Last_dustem_fit.ps ; help = If set print this help ; COMMON BLOCKS: ; None ; SIDE EFFECTS: ; None ; RESTRICTIONS: ; The dustem idl wrapper must be installed ; PROCEDURE: ; None ; EXAMPLES ; dustem_fit_sed_extsed_readme,/postscript ; MODIFICATION HISTORY: ; Written by J.P. Bernard June 29th 2011 ; see evolution details on the dustem cvs maintained at CESR ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems. ;- IF keyword_set(help) THEN BEGIN doc_library,'dustem_fit_sed_extsed_readme' goto,the_end ENDIF ;=== initialise dustem ;dustem_init,mode='COMPIEGNE_ETAL2010',/pol use_model='COMPIEGNE_ETAL2010' IF keyword_set(model) THEN use_model=model dustem_init,mode=use_model !dustem_verbose=1 !dustem_show_plot=1 IF keyword_set(postscript) THEN BEGIN ps_dir=!dustem_dat+'/Docs/Figures/' force_mkdir,ps_dir ENDIF ;stop ;==== read emission SED dir=!dustem_wrap_soft_dir+'/Data/SEDs/' file=dir+'Gal_composite_spectrum.xcat' ;file=dir+'Mathis1990_DISM.xcat' ;file=dir+'Gal_composite_spectrum_wAKARI.xcat' ;stop spec=read_xcat(file,/silent) ind=where(spec.error EQ 0.,count) IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec ind=where(spec.instru EQ 'FIRAS',count) IF count NE 0 THEN spec(ind).error=0.2*spec(ind).spec ;==== read extinction SED dir=!dustem_wrap_soft_dir+'/Data/EXTs/' file=dir+'Mathis1990_DISM.xcat' extsed=read_xcat(file,/silent) ind=where(extsed.error EQ 0.,count) IF count NE 0 THEN extsed[ind].error=0.2*spec(ind).spec st=dustem_set_data(ext=extsed,sed=spec,rchi2_weight=rchi2_weight,f_HI=f_HI) ;=== Set which parameters you want to fit pd = [ $ 'dustem_plugin_continuum_2',$ '(*!dustem_params).gas.G0', $ ;G0 '(*!dustem_params).grains(0).mdust_o_mh',$ ;PAH0 mass fraction '(*!dustem_params).grains(1).mdust_o_mh',$ ;PAH1 mass fraction '(*!dustem_params).grains(2).mdust_o_mh' $ ;amCBEx ] p_truth=[0.002,1.,7.8000E-04,7.8000E-04,1.6500E-04] ;=== set initial value of parameters you want to fit ;iv=p_truth ;exact solution iv=p_truth+[0.01,0.,3.e-5,1.e-5,-2.e-5] ;shifted from solution ulimed=[0, 0 ,0 , 0 ,0 ] llimed=[1 ,1 ,1 , 1 ,1 ] llims =[0. ,0., 0. ,0. ,0. ] ; dustem_init_plugins, pd ;will have to be uncommented for this procedure to work. ;== SET THE FITTED PARAMETERS (now done in dustem_init_parinfo) dustem_init_parinfo,pd,iv,up_limited=ulimed,lo_limited=llimed,up_limits=ulims,lo_limits=llims ;=== RUN fit tol=1.e-12 xtol=1.e-12 Nitermax=3 ;maximum number of iteration. This is the criterium which will stop the fit procedure IF keyword_set(itermax) THEN Nitermax=itermax xrange=[1.,2.e4] yrange=[1e-4,1.e3] loadct,13 ;!y.range=[1e-8,100] ;This is to ajust plot range from outside the routine t1=systime(0,/sec) yr=[1e-4,10] xr=[1,1.e5] legend_xpos=0.6 legend_ypos=0.8 tit='DUSTEMWrapper Intensity SED (Running)' ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') res=dustem_mpfit_data(tol=tol,xtol=xtol,Nitermax=Nitermax,xrange=xrange,/xstyle,yrange=yrange,/ysty,/ylog,/xlog,xtit=xtit,ytit=ytit,legend_xpos=legend_xpos,legend_ypos=legend_ypos,title=tit) t2=systime(0,/sec) print,(res-p_truth)/p_truth*100. ;=== SAVE FIT RESULTS file_out='/tmp/DUSTEM_sedextsed_fit_example.sav' dustem_save_system_variables,file_out message,'Saved '+file_out,/continue ;====================================== ;====You can exit IDL here and re-enter ;====================================== file='/tmp/DUSTEM_sedextsed_fit_example.sav' dustem_restore_system_variables,file ;=== Plot best fit intensity SED win=0L window,win & win=win+1 cleanplot xrange=[1.,2.e4] yrange=[1e-4,1.e3] res=*(*!dustem_fit).current_param_values chi2=(*!dustem_fit).chi2 rchi2=(*!dustem_fit).rchi2 errors=(*(*!dustem_fit).current_param_errors)*(*(*!dustem_fit).param_init_values) tit='DUSTEMWrapper polarization dust Example (Saved)' ytit=textoidl('I_\nu (MJy/sr) for N_H=10^{20} H/cm^2') xtit=textoidl('\lambda (\mum)') dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos loadct,13 IF keyword_set(postscript) THEN BEGIN previous_device=!d.name set_plot,'PS' !p.charsize=0.8 ps_file=ps_dir+'Last_dustem_sed_fit.ps' device,filename=ps_file,/color dustem_sed_plot,*(*!dustem_fit).current_param_values,ytit=ytit,xtit=xtit,title=tit,yr=yrange,xr=xrange,/ysty,/xsty,res=res,errors=errors,chi2=chi2,rchi2=rchi2,/xlog,/ylog,legend_xpos=legend_xpos,legend_ypos=legend_ypos device,/close set_plot,previous_device message,'wrote '+ps_file,/info ENDIF ;stop ;=== Plot best fit Extinction SED loadct,13 window,win & win=win+1 cleanplot xr=[0.3,40] yr=[1.e-26,1.e-21] dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty IF keyword_set(postscript) THEN BEGIN previous_device=!d.name set_plot,'PS' !p.charsize=0.8 ps_file=ps_dir+'Last_dustem_ext_ir_fit.ps' device,filename=ps_file,/color dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty device,/close set_plot,previous_device message,'wrote '+ps_file,/info ENDIF ;stop loadct,13 window,win & win=win+1 cleanplot xr=[0.2,10] yr=[0,3e-21] dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty,/UV IF keyword_set(postscript) THEN BEGIN previous_device=!d.name set_plot,'PS' !p.charsize=0.8 ps_file=ps_dir+'Last_dustem_ext_uv_fit.ps' device,filename=ps_file,/color dustem_extsed_plot,*(*!dustem_fit).current_param_values,pol=pol,xr=xr,yr=yr,/xsty,/ysty,/UV device,/close set_plot,previous_device message,'wrote '+ps_file,/info ENDIF print,'dustem_mpfit_sed executed in ',t2-t1,' sec' stop the_end: END