PRO dustem_run_readme,model=model,help=help,plot_it=plot_it,png=png,postcript=postcript ;+ ; NAME: ; dustem_run_readme ; PURPOSE: ; This is an example of how to run dustem SEDs with the dustem wrapper. ; CATEGORY: ; Dustem ; CALLING SEQUENCE: ; dustem_run_readme,model=model,help=help,/postcript,/help,/plot_it ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: ; None ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: ; model = Selects one of the dust mixture used by dustem ; 'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default) ; 'DBP90' from Desert et al 1990 ; 'DL01' from Draine & Li 2001 ; 'DL07' from Draine & Li 2007 ; 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 fortran code must be installed ; The dustem idl wrapper must be installed ; PROCEDURE: ; None ; EXAMPLES ; dustem_run_readme,model='COMPIEGNE_ETAL2010' ; MODIFICATION HISTORY: ; Written by J.P. Bernard April 1st 2011 ; see evolution details on the dustem cvs maintained at CESR ; Contact J.-Ph. Bernard (Jean-Philippe.Bernard@cesr.fr) in case of problems. ;- ;This Readme shows how to run dustem from the wrapper ;Obviously, the dustem package must have been installed succesfully (see ;dustem_cvs_readme.txt for install instructions). IF keyword_set(help) THEN BEGIN doc_library,'dustem_run_readme' goto,the_end ENDIF loadct,13 IF keyword_set(model) THEN BEGIN use_model=strupcase(model) ENDIF ELSE BEGIN use_model='COMPIEGNE_ETAL2010' ;Default is last dustem model ENDELSE ;== This just initializes dustem_init,model=use_model !dustem_verbose=1 !dustem_show_plot=1 IF keyword_set(png) THEN BEGIN dir_png=!dustem_dat+'/Figures/' force_mkdir,dir_png ENDIF IF keyword_set(postcript) THEN BEGIN dir_ps=!dustem_dat+'/Figures/' force_mkdir,dir_ps ENDIF IF not keyword_set(dir) THEN BEGIN ; CASE getenv('DUSTEM_WHICH') OF CASE !dustem_which OF 'DESERT':dir='DESERT_POST_ISO_MORE/' 'COMPIEGNE': dir='MC_DAT/' 'VERSTRAETE': dir='d_3.5/' 'VERSTRAETE': dir='LV_DAT/' ;'WEB3p8': dir_in='dustem3.8_web/' ; this is the original 'WEB3p8': dir_in=!dustem_soft_dir ; I changed this ELSE: dir='les_DAT/' ENDCASE ENDIF ;== This is the input directory selected for Dustem ;dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/'+dir ;dir_in=!dustem_soft_dir+'/src/'+dir ;== The following structure contains the inputs to the model st_model=dustem_read_all(dir_in) ;help,st_model,/str ;=== You can here modify the model inputs ;=== e.g. ;stop ;st_model.G0=2. ;===== This below is to test spinning. Does not work with DBP90 ! ;st_model.grains[0].type_keywords='logn-chrg-spin' ;st_model.gas.NH=1.e2 ;st_model.gas.Tgas=1.e5 ;== The following line writes inputs to the Fortran dustem_write_all,st_model,!dustem_dat ;== This runs the Fortran. The ouput structure contains the results st=dustem_run() ;== Optionaly, this creates an SED to be overploted ;stop ;filters=['IRAS1','IRAS4'] ;Nfilt=n_elements(filters) ;sed=dustem_initialize_sed(Nfilt) ;sed.filter=filters ;sed.wave=dustem_filter2wav(filters) ;sed.instru=dustem_filter2instru(filters) ;stop ;stt=dustem_set_data(sed=sed) ;== The following plots the resulting emission SED win=0L window,win & win=win+2 xtit=textoidl('\lambda (\mum)') ;=== same plot as in Compiegne et al. 2010 ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)') tit='DUSTEMWrapper Emitted Intensity '+use_model yr=[1e-11,6.e-7] xr=[1,5e3] dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit IF keyword_set(postcript) THEN BEGIN file_ps=dir_ps+'Last_dustem_run_readme_nuinuem.ps' previous_device=!d.name set_plot,'PS' device,filename=file_ps dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit device,/close set_plot,previous_device message,'Wrote '+file_ps,/continue ENDIF ;== The following plots the resulting extinction window,win & win=win+1 yr=[0,2.5] ;range of 1/lambda xr=[1,10] ;range of sigma xtit=textoidl('1/\lambda (\mum^{-1})') ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)') tit='DUSTEMWrapper extinction '+use_model dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit IF keyword_set(postcript) THEN BEGIN file_ps=dir_ps+'Last_dustem_run_readme_extinction.ps' previous_device=!d.name set_plot,'PS' device,filename=file_ps dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit device,/close set_plot,previous_device message,'Wrote '+file_ps,/continue ENDIF IF keyword_set(png) THEN BEGIN win=10L window,win ;& win=win+2 file_png=dir_png+'Last_dustem_run_readme_nuinuem.png' xtit=textoidl('\lambda (\mum)') ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)') tit='DUSTEMWrapper Emitted Intensity '+use_model yr=[1e-11,6.e-7] xr=[1,5e3] dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit,xtit=xtit,ytit=ytit write_png,file_png,tvrd(/true) message,'Wrote '+file_png,/info ;===== file_png=dir_png+'Last_dustem_run_readme_extinction.png' yr=[0,2.5] ;range of 1/lambda xr=[1,10] ;range of sigma xtit=textoidl('1/\lambda (\mum^{-1})') ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)') tit='DUSTEMWrapper extinction '+use_model dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit write_png,file_png ,tvrd(/true),/verbose message,'Wrote '+file_png,/info ENDIF the_end: END