From e69bf24562e7442e7a7e6ecce9e46dc30c6022ef Mon Sep 17 00:00:00 2001 From: Jean-Philippe Bernard Date: Tue, 13 Jul 2021 17:18:07 +0200 Subject: [PATCH] improved --- src/idl/dustem_init.pro | 20 +++++--------------- src/idl/dustem_plot_extinction.pro | 49 +++++++++++++++++++++++-------------------------- src/idl/dustem_plot_nuinu_em.pro | 11 +++++------ src/idl/dustem_read_align.pro | 7 +++++-- src/idl/dustem_read_grain.pro | 3 +++ src/idl/dustem_run_readme.pro | 43 ++++++++++++++++++++++++++----------------- src/idl/dustem_write_grain.pro | 2 ++ src/idl/dustem_write_grain_web3p8.pro | 4 +++- src/idl/dustem_write_qpol.pro | 7 ++++--- 9 files changed, 76 insertions(+), 70 deletions(-) diff --git a/src/idl/dustem_init.pro b/src/idl/dustem_init.pro index 22158ba..97e7954 100644 --- a/src/idl/dustem_init.pro +++ b/src/idl/dustem_init.pro @@ -1,4 +1,4 @@ -PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,mode=mode,help=help,noerase=noerase,st_model=st_model,pol=pol +PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,model=model,help=help,noerase=noerase,st_model=st_model,pol=pol ;+ ; NAME: @@ -8,7 +8,7 @@ PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,mode=mode,help=help,no ; PURPOSE: ; Defines system variables to run the DUSTEM code ; INPUTS: -; mode = Selects one of the dust mixture used by dustem +; 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 @@ -179,8 +179,8 @@ defsysv,'!dustem_show_plot',1 ;0=show no plot dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT',ALIGN:'ALIGN.DAT'} !dustem_inputs=ptr_new(dustem_inputs) ;===This is the Desert option -IF keyword_set(mode) THEN BEGIN - CASE mode OF +IF keyword_set(model) THEN BEGIN + CASE model OF 'COMPIEGNE_ETAL2010':BEGIN (*!dustem_inputs).grain='GRAIN_MC10.DAT' END @@ -213,7 +213,7 @@ IF keyword_set(mode) THEN BEGIN (*!dustem_inputs).align='ALIGN_G17_MODELD.DAT' END ELSE:BEGIN - message,'mode '+mode+' unknown' + message,'model '+model+' unknown' END ENDCASE ENDIF @@ -249,7 +249,6 @@ ENDIF ;=== remove dynamical storage directories ;=== This is needed for iterative use of DUSTEM from IDL if not keyword_set(noerase) then begin -;CASE getenv('DUSTEM_WHICH') OF CASE !dustem_which OF 'VERSTRAETE':BEGIN IF file_test(!dustem_dat,/dir) THEN BEGIN @@ -260,10 +259,6 @@ CASE !dustem_which OF IF file_test(!dustem_res,/dir) THEN BEGIN spawn,'rm -rf '+!dustem_res+'/les_RES' ENDIF -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_DAT' -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_QABS' -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_CAPA' -; spawn,'rm -rf '+getenv('DUSTEM_RES')+'/les_RES' END 'WEB3p8':BEGIN IF file_test(!dustem_dat,/dir) THEN BEGIN @@ -272,10 +267,6 @@ CASE !dustem_which OF spawn,'rm -rf '+!dustem_dat+'/oprop' spawn,'rm -rf '+!dustem_dat+'/out' ENDIF -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/data' -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/hcap' -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/oprop' -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/out' END ELSE: BEGIN END @@ -283,7 +274,6 @@ ENDCASE ;=== create storage directories ;=== This is needed for iterative use of DUSTEM from IDL -;CASE getenv('DUSTEM_WHICH') OF CASE !dustem_which OF 'VERSTRAETE':BEGIN IF file_test(!dustem_dat,/dir) THEN BEGIN diff --git a/src/idl/dustem_plot_extinction.pro b/src/idl/dustem_plot_extinction.pro index 7580312..3a47a9c 100644 --- a/src/idl/dustem_plot_extinction.pro +++ b/src/idl/dustem_plot_extinction.pro @@ -1,4 +1,4 @@ -PRO dustem_plot_extinction,st,st_model,_Extra=extra,help=help +PRO dustem_plot_extinction,st,st_model,_extra=_extra,help=help ;+ ; NAME: @@ -47,45 +47,42 @@ IF keyword_set(help) THEN BEGIN goto,the_end ENDIF -;stop -;Ngrains=n_tags(st.dustem)-2 Ngrains=st_model.Ngrains -colors=dustem_grains_colors(Ngrains) -;colmax=210 & colmin=30 -;colors=findgen(Ngrains)/(1.*Ngrains-1)*(colmax-colmin)+colmin +colors=dustem_grains_colors(Ngrains,/cgplot) -xtit=textoidl('\lambda (\mum)') -;ytit=textoidl('extinction (cm^2/g)') & fact=1. -ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)') -mH=1.66e-24 +mH=1.66e-24 ;g fact=mH*1.e21 -tit='Dustem extinction' ;stop +;st.ext is supposed to be in cm^2/g (of dust) + i=0L -;plot_oo,st.ext.wav,st.ext.abs_grain(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit ffact=fact*st_model.grains[i].MDUST_O_MH -;plot_oo,st.ext.wav,(st.ext.abs_grain(i)+st.ext.sca_grain(i))*ffact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit -cgplot,st.ext.wav,(st.ext.abs_grain[i]+st.ext.sca_grain(i))*ffact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,/xlog,/ylog -abs_tot=st.ext.abs_grain(i)*0. -sca_tot=st.ext.sca_grain(i)*0. +yr=[0,2.5] +xr=[1,10] +;cgplot,1./st.ext.wav,st.ext.EXT_TOT,color='red',xtit=xtit,ytit=ytit,tit=tit,yr=yr,/ysty,xr=xr,/xsty +cgplot,1./st.ext.wav,st.ext.EXT_TOT,color='red',_extra=_extra +abs_tot=st.ext.abs_grain[i]*0. +sca_tot=st.ext.sca_grain[i]*0. FOR i=0L,Ngrains-1 DO BEGIN - ffact=fact*st_model.grains(i).MDUST_O_MH - ;oplot,st.ext.wav,(st.ext.abs_grain(i)+st.ext.sca_grain(i))*ffact,color=colors(i),psym=psym - cgoplot,st.ext.wav,(st.ext.abs_grain[i]+st.ext.sca_grain[i])*ffact,color=colors(i),psym=psym -; oplot,st.ext.wav,st.ext.abs_grain(i)*fact,color=colors(i-1),psym=psym -; oplot,st.ext.wav,st.ext.sca_grain(i)*fact,color=colors(i-1),psym=psym,linestyle=2 + ;ffact=fact*st_model.grains[i].MDUST_O_MH + ffact=1. ;because it looks as if the values are already in 1e-21 cm^2/H + abs=st.ext.abs_grain[i]*ffact + sca=st.ext.sca_grain[i]*ffact + ext=abs+sca + cgoplot,1./(st.ext.wav),abs,color=colors[i],psym=psym,linestyle=1 + cgoplot,1./(st.ext.wav),sca,color=colors[i],psym=psym,linestyle=2 + cgoplot,1./(st.ext.wav),ext,color=colors[i],psym=psym,linestyle=0 abs_tot=abs_tot+st.ext.abs_grain[i]*ffact sca_tot=sca_tot+st.ext.sca_grain[i]*ffact ENDFOR -;oplot,st.ext.wav,(abs_tot+sca_tot),color=255,psym=psym -cgoplot,st.ext.wav,(abs_tot+sca_tot),color=255,psym=psym -;oplot,st.dustem.wav,(abs_tot+sca_tot),color=255,psym=psym -;oplot,st.dustem.wav,abs_tot*fact,color=255,psym=psym -;oplot,st.dustem.wav,sca_tot*fact,color=255,psym=psym,linestyle=2 +total_extinction=abs_tot+sca_tot + +cgoplot,1./(st.ext.wav),total_extinction,color='black',psym=psym the_end: +;stop END diff --git a/src/idl/dustem_plot_nuinu_em.pro b/src/idl/dustem_plot_nuinu_em.pro index 31d2f4f..e751422 100644 --- a/src/idl/dustem_plot_nuinu_em.pro +++ b/src/idl/dustem_plot_nuinu_em.pro @@ -1,4 +1,4 @@ -PRO dustem_plot_nuinu_em,st,_Extra=extra,help=help +PRO dustem_plot_nuinu_em,st,_extra=_extra,help=help ;+ ; NAME: @@ -46,15 +46,14 @@ ENDIF Ngrains=n_tags(st.sed)-2 colors=dustem_grains_colors(Ngrains,/cgplot) -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)') & fact=1./4./!pi*1.e4/1.e7*1.e20 -tit='Dustem Emitted Intensity' +;This is to go to units of W/m^2/sr for N_H=1.e20 H/cm^2 +fact=1./4./!pi*1.e4/1.e7*1.e20 i=1L ;plot_oo,st.sed.wav,st.sed.(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit -cgplot,st.sed.wav,st.sed.(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,/xlog,/ylog +;cgplot,st.sed.wav,st.sed.(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,/xlog,/ylog +cgplot,st.sed.wav,st.sed.(i)*fact,/nodata,_extra=_extra FOR i=1L,Ngrains DO BEGIN ;oplot,st.sed.wav,st.sed.(i)*fact,color=colors(i-1),psym=psym cgoplot,st.sed.wav,st.sed.(i)*fact,color=colors[i-1],psym=psym diff --git a/src/idl/dustem_read_align.pro b/src/idl/dustem_read_align.pro index 9c94dfc..138a7d0 100644 --- a/src/idl/dustem_read_align.pro +++ b/src/idl/dustem_read_align.pro @@ -29,7 +29,10 @@ ENDREP UNTIL first_char NE '#' key_str=str IF stregex(key_str, 'lin', /bool) THEN !run_lin = 1 -IF stregex(key_str, 'univ', /bool) THEN !run_univ = 1 +IF stregex(key_str, 'univ', /bool) THEN BEGIN + !run_univ = 1 + stop +ENDIF IF stregex(key_str, 'circ', /bool) THEN !run_circ = 1 IF stregex(key_str, 'anis', /bool) THEN !run_anis = 1 IF stregex(key_str, 'rrf', /bool) THEN !run_rrf = 1 @@ -84,7 +87,7 @@ FOR i=0L,Ngrains-1 DO BEGIN st[i].plev = strv(ii) & ii=ii+1 ENDIF st[i].aligned = stregex(st_grains(i).type_keywords, 'pol', /bool) - ;stop + stop IF !run_univ eq 1 THEN BEGIN st = replicate(st(i),Ngrains) diff --git a/src/idl/dustem_read_grain.pro b/src/idl/dustem_read_grain.pro index 242ace5..8e3ab58 100644 --- a/src/idl/dustem_read_grain.pro +++ b/src/idl/dustem_read_grain.pro @@ -197,8 +197,11 @@ CASE !dustem_which OF END ENDCASE + the_end: +;stop + RETURN,full_st END diff --git a/src/idl/dustem_run_readme.pro b/src/idl/dustem_run_readme.pro index efd72e8..45e8a69 100644 --- a/src/idl/dustem_run_readme.pro +++ b/src/idl/dustem_run_readme.pro @@ -1,4 +1,4 @@ -PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it +PRO dustem_run_readme,postcript=postcript,model=model,help=help,plot_it=plot_it ;+ ; NAME: @@ -8,7 +8,7 @@ PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it ; CATEGORY: ; Dustem ; CALLING SEQUENCE: -; dustem_run_readme,mode=mode,help=help,/postcript,/help,/plot_it +; dustem_run_readme,model=model,help=help,/postcript,/help,/plot_it ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: @@ -18,7 +18,7 @@ PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it ; OPTIONAL OUTPUT PARAMETERS: ; None ; ACCEPTED KEY-WORDS: -; mode = Selects one of the dust mixture used by dustem +; 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 @@ -35,7 +35,7 @@ PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it ; PROCEDURE: ; None ; EXAMPLES -; dustem_run_readme,mode='COMPIEGNE_ETAL2010' +; 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 @@ -53,14 +53,14 @@ ENDIF loadct,13 -IF keyword_set(mode) THEN BEGIN - use_mode=strupcase(mode) +IF keyword_set(model) THEN BEGIN + use_model=strupcase(model) ENDIF ELSE BEGIN - use_mode='COMPIEGNE_ETAL2010' ;Default is last dustem model + use_model='COMPIEGNE_ETAL2010' ;Default is last dustem model ENDELSE ;== This just initializes -dustem_init,mode=use_mode +dustem_init,model=use_model !dustem_verbose=1 !dustem_show_plot=1 @@ -81,18 +81,18 @@ ENDIF ;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 +;help,st_model,/str ;=== You can here modify the model inputs ;=== e.g. ;stop ;st_model.G0=2. -st_model.grains[0].type_keywords='logn-chrg-spin' -st_model.gas.NH=1.e2 -st_model.gas.Tgas=1.e5 +;===== 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,getenv('DUSTEM_DAT') dustem_write_all,st_model,!dustem_dat ;== This runs the Fortran. The ouput structure contains the results st=dustem_run() @@ -109,12 +109,21 @@ st=dustem_run() ;== The following just plots the resulting emission SED window,0 -;dustem_plot_nuinu_em,st,yr=[1e-28,1.e-22],/ysty,xr=[1,5e3],/xsty -dustem_plot_nuinu_em,st,yr=[1e-11,6.e-7],/ysty,xr=[1,5e3],/xsty +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='Dustem 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 ;== The following just plots the resulting extinction window,1 -;dustem_plot_extinction,st,st_model,yr=[1.e-4,1.e7],/ysty,xr=[1.e-2,1e5],/xsty -dustem_plot_extinction,st,st_model,yr=[1.e-6,1.e0],/ysty,xr=[1.,400],/xsty +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='Dustem extinction '+use_model +dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit the_end: diff --git a/src/idl/dustem_write_grain.pro b/src/idl/dustem_write_grain.pro index b982203..88bb80e 100644 --- a/src/idl/dustem_write_grain.pro +++ b/src/idl/dustem_write_grain.pro @@ -1,5 +1,7 @@ PRO dustem_write_grain,file,st +;stop + OPENW,unit,file,/get_lun Ntypes=n_elements(st) diff --git a/src/idl/dustem_write_grain_web3p8.pro b/src/idl/dustem_write_grain_web3p8.pro index 7d8bd7c..c5e34f1 100644 --- a/src/idl/dustem_write_grain_web3p8.pro +++ b/src/idl/dustem_write_grain_web3p8.pro @@ -3,6 +3,8 @@ PRO dustem_write_grain_web3p8,file,st ;Remark: It is a bit odfd that this is using both st and (*!dustem_params) ;Remark: This should use automated comment detection, as in dustem_write_gas.pro +;stop + Ncomments=100 c=strarr(Ncomments) @@ -39,7 +41,7 @@ printf,unit,(*!dustem_params).G0 ;frmt='(A12,I4,A12,11E14.6)' ;frmt='(A12,I4,A20,11E14.6)' ;frmt='(A20,I4,A20,11E14.6)' -frmt='(A20,I4,A25,11E14.6)' +frmt='(A30,I4,A25,11E14.6)' FOR i=0L,(*!dustem_params).ngrains-1 DO BEGIN printf,unit,st(i).grain_type, $ diff --git a/src/idl/dustem_write_qpol.pro b/src/idl/dustem_write_qpol.pro index 5e6dbf8..99b0888 100644 --- a/src/idl/dustem_write_qpol.pro +++ b/src/idl/dustem_write_qpol.pro @@ -18,12 +18,13 @@ FOR i=0L,Nfiles-1 DO BEGIN ;filename=dir+'Q'+strtrim(i_axis,2)+'_'+st.grain.grains(i).grain_type+Estring filename=dir+'Q'+strtrim(i_axis,2)+'_'+st.grains(i).grain_type+Estring + message,'Will write '+filename,/continue Nsizes=n_elements((*st.qabs(i)).sizes) OPENW,unit,filename,/get_lun printf,unit,Nsizes - printf,unit,(*st.qabs(i)).sizes,format=frmt - NLines=n_elements((*st.qpol(i_axis-1,i)).qabs) - NQabs=n_tags((*st.qpol(i_axis-1,i)).qabs) + printf,unit,(*st.qabs[i]).sizes,format=frmt + NLines=n_elements((*st.qpol[i_axis-1,i]).qabs) + NQabs=n_tags((*st.qpol[i_axis-1,i]).qabs) vars=fltarr(NQabs,Nlines) ;Separator -- libgit2 0.21.2