OV-GSO-DC / dustem-wrapper_idl

Commit e69bf24562e7442e7a7e6ecce9e46dc30c6022ef

Authored by Jean-Philippe Bernard
1 parent 84327150
Exists in master

improved

Showing 9 changed files with 76 additions and 70 deletions   Show diff stats
Inline Side-by-side
src/idl/dustem_init.pro
  Diff comments   View file @e69bf24
1   -PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,mode=mode,help=help,noerase=noerase,st_model=st_model,pol=pol
  1 +PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,model=model,help=help,noerase=noerase,st_model=st_model,pol=pol
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -8,7 +8,7 @@ PRO dustem_init,dir=dir,wraptest=wraptest,plot_it=plot_it,mode=mode,help=help,no
8 8 ; PURPOSE:
9 9 ; Defines system variables to run the DUSTEM code
10 10 ; INPUTS:
11   -; mode = Selects one of the dust mixture used by dustem
  11 +; model = Selects one of the dust mixture used by dustem
12 12 ; 'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
13 13 ; 'DBP90' from Desert et al 1990
14 14 ; 'DL01' from Draine & Li 2001
... ... @@ -179,8 +179,8 @@ defsysv,'!dustem_show_plot',1 ;0=show no plot
179 179 dustem_inputs={GRAIN:'GRAIN.DAT',ISRF:'ISRF.DAT',ALIGN:'ALIGN.DAT'}
180 180 !dustem_inputs=ptr_new(dustem_inputs)
181 181 ;===This is the Desert option
182   -IF keyword_set(mode) THEN BEGIN
183   - CASE mode OF
  182 +IF keyword_set(model) THEN BEGIN
  183 + CASE model OF
184 184 'COMPIEGNE_ETAL2010':BEGIN
185 185 (*!dustem_inputs).grain='GRAIN_MC10.DAT'
186 186 END
... ... @@ -213,7 +213,7 @@ IF keyword_set(mode) THEN BEGIN
213 213 (*!dustem_inputs).align='ALIGN_G17_MODELD.DAT'
214 214 END
215 215 ELSE:BEGIN
216   - message,'mode '+mode+' unknown'
  216 + message,'model '+model+' unknown'
217 217 END
218 218 ENDCASE
219 219 ENDIF
... ... @@ -249,7 +249,6 @@ ENDIF
249 249 ;=== remove dynamical storage directories
250 250 ;=== This is needed for iterative use of DUSTEM from IDL
251 251 if not keyword_set(noerase) then begin
252   -;CASE getenv('DUSTEM_WHICH') OF
253 252 CASE !dustem_which OF
254 253 'VERSTRAETE':BEGIN
255 254 IF file_test(!dustem_dat,/dir) THEN BEGIN
... ... @@ -260,10 +259,6 @@ CASE !dustem_which OF
260 259 IF file_test(!dustem_res,/dir) THEN BEGIN
261 260 spawn,'rm -rf '+!dustem_res+'/les_RES'
262 261 ENDIF
263   -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_DAT'
264   -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_QABS'
265   -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/les_CAPA'
266   -; spawn,'rm -rf '+getenv('DUSTEM_RES')+'/les_RES'
267 262 END
268 263 'WEB3p8':BEGIN
269 264 IF file_test(!dustem_dat,/dir) THEN BEGIN
... ... @@ -272,10 +267,6 @@ CASE !dustem_which OF
272 267 spawn,'rm -rf '+!dustem_dat+'/oprop'
273 268 spawn,'rm -rf '+!dustem_dat+'/out'
274 269 ENDIF
275   -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/data'
276   -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/hcap'
277   -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/oprop'
278   -; spawn,'rm -rf '+getenv('DUSTEM_DAT')+'/out'
279 270 END
280 271 ELSE: BEGIN
281 272 END
... ... @@ -283,7 +274,6 @@ ENDCASE
283 274  
284 275 ;=== create storage directories
285 276 ;=== This is needed for iterative use of DUSTEM from IDL
286   -;CASE getenv('DUSTEM_WHICH') OF
287 277 CASE !dustem_which OF
288 278 'VERSTRAETE':BEGIN
289 279 IF file_test(!dustem_dat,/dir) THEN BEGIN
... ...
src/idl/dustem_plot_extinction.pro
  Diff comments   View file @e69bf24
1   -PRO dustem_plot_extinction,st,st_model,_Extra=extra,help=help
  1 +PRO dustem_plot_extinction,st,st_model,_extra=_extra,help=help
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -47,45 +47,42 @@ IF keyword_set(help) THEN BEGIN
47 47 goto,the_end
48 48 ENDIF
49 49  
50   -;stop
51   -;Ngrains=n_tags(st.dustem)-2
52 50 Ngrains=st_model.Ngrains
53 51  
54   -colors=dustem_grains_colors(Ngrains)
55   -;colmax=210 & colmin=30
56   -;colors=findgen(Ngrains)/(1.*Ngrains-1)*(colmax-colmin)+colmin
  52 +colors=dustem_grains_colors(Ngrains,/cgplot)
57 53  
58   -xtit=textoidl('\lambda (\mum)')
59   -;ytit=textoidl('extinction (cm^2/g)') & fact=1.
60   -ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
61   -mH=1.66e-24
  54 +mH=1.66e-24 ;g
62 55 fact=mH*1.e21
63   -tit='Dustem extinction'
64 56  
65 57 ;stop
66 58  
  59 +;st.ext is supposed to be in cm^2/g (of dust)
  60 +
67 61 i=0L
68   -;plot_oo,st.ext.wav,st.ext.abs_grain(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit
69 62 ffact=fact*st_model.grains[i].MDUST_O_MH
70   -;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
71   -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
72   -abs_tot=st.ext.abs_grain(i)*0.
73   -sca_tot=st.ext.sca_grain(i)*0.
  63 +yr=[0,2.5]
  64 +xr=[1,10]
  65 +;cgplot,1./st.ext.wav,st.ext.EXT_TOT,color='red',xtit=xtit,ytit=ytit,tit=tit,yr=yr,/ysty,xr=xr,/xsty
  66 +cgplot,1./st.ext.wav,st.ext.EXT_TOT,color='red',_extra=_extra
  67 +abs_tot=st.ext.abs_grain[i]*0.
  68 +sca_tot=st.ext.sca_grain[i]*0.
74 69 FOR i=0L,Ngrains-1 DO BEGIN
75   - ffact=fact*st_model.grains(i).MDUST_O_MH
76   - ;oplot,st.ext.wav,(st.ext.abs_grain(i)+st.ext.sca_grain(i))*ffact,color=colors(i),psym=psym
77   - cgoplot,st.ext.wav,(st.ext.abs_grain[i]+st.ext.sca_grain[i])*ffact,color=colors(i),psym=psym
78   -; oplot,st.ext.wav,st.ext.abs_grain(i)*fact,color=colors(i-1),psym=psym
79   -; oplot,st.ext.wav,st.ext.sca_grain(i)*fact,color=colors(i-1),psym=psym,linestyle=2
  70 + ;ffact=fact*st_model.grains[i].MDUST_O_MH
  71 + ffact=1. ;because it looks as if the values are already in 1e-21 cm^2/H
  72 + abs=st.ext.abs_grain[i]*ffact
  73 + sca=st.ext.sca_grain[i]*ffact
  74 + ext=abs+sca
  75 + cgoplot,1./(st.ext.wav),abs,color=colors[i],psym=psym,linestyle=1
  76 + cgoplot,1./(st.ext.wav),sca,color=colors[i],psym=psym,linestyle=2
  77 + cgoplot,1./(st.ext.wav),ext,color=colors[i],psym=psym,linestyle=0
80 78 abs_tot=abs_tot+st.ext.abs_grain[i]*ffact
81 79 sca_tot=sca_tot+st.ext.sca_grain[i]*ffact
82 80 ENDFOR
83   -;oplot,st.ext.wav,(abs_tot+sca_tot),color=255,psym=psym
84   -cgoplot,st.ext.wav,(abs_tot+sca_tot),color=255,psym=psym
85   -;oplot,st.dustem.wav,(abs_tot+sca_tot),color=255,psym=psym
86   -;oplot,st.dustem.wav,abs_tot*fact,color=255,psym=psym
87   -;oplot,st.dustem.wav,sca_tot*fact,color=255,psym=psym,linestyle=2
  81 +total_extinction=abs_tot+sca_tot
  82 +
  83 +cgoplot,1./(st.ext.wav),total_extinction,color='black',psym=psym
88 84  
89 85 the_end:
  86 +;stop
90 87  
91 88 END
... ...
src/idl/dustem_plot_nuinu_em.pro
  Diff comments   View file @e69bf24
1   -PRO dustem_plot_nuinu_em,st,_Extra=extra,help=help
  1 +PRO dustem_plot_nuinu_em,st,_extra=_extra,help=help
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -46,15 +46,14 @@ ENDIF
46 46 Ngrains=n_tags(st.sed)-2
47 47 colors=dustem_grains_colors(Ngrains,/cgplot)
48 48  
49   -xtit=textoidl('\lambda (\mum)')
50   -;=== same plot as in Compiegne et al. 2010
51   -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
52   -tit='Dustem Emitted Intensity'
  49 +;This is to go to units of W/m^2/sr for N_H=1.e20 H/cm^2
  50 +fact=1./4./!pi*1.e4/1.e7*1.e20
53 51  
54 52  
55 53 i=1L
56 54 ;plot_oo,st.sed.wav,st.sed.(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit
57   -cgplot,st.sed.wav,st.sed.(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,/xlog,/ylog
  55 +;cgplot,st.sed.wav,st.sed.(i)*fact,/nodata,_Extra=extra,xtit=xtit,ytit=ytit,tit=tit,/xlog,/ylog
  56 +cgplot,st.sed.wav,st.sed.(i)*fact,/nodata,_extra=_extra
58 57 FOR i=1L,Ngrains DO BEGIN
59 58 ;oplot,st.sed.wav,st.sed.(i)*fact,color=colors(i-1),psym=psym
60 59 cgoplot,st.sed.wav,st.sed.(i)*fact,color=colors[i-1],psym=psym
... ...
src/idl/dustem_read_align.pro
  Diff comments   View file @e69bf24
... ... @@ -29,7 +29,10 @@ ENDREP UNTIL first_char NE '#'
29 29 key_str=str
30 30  
31 31 IF stregex(key_str, 'lin', /bool) THEN !run_lin = 1
32   -IF stregex(key_str, 'univ', /bool) THEN !run_univ = 1
  32 +IF stregex(key_str, 'univ', /bool) THEN BEGIN
  33 + !run_univ = 1
  34 + stop
  35 +ENDIF
33 36 IF stregex(key_str, 'circ', /bool) THEN !run_circ = 1
34 37 IF stregex(key_str, 'anis', /bool) THEN !run_anis = 1
35 38 IF stregex(key_str, 'rrf', /bool) THEN !run_rrf = 1
... ... @@ -84,7 +87,7 @@ FOR i=0L,Ngrains-1 DO BEGIN
84 87 st[i].plev = strv(ii) & ii=ii+1
85 88 ENDIF
86 89 st[i].aligned = stregex(st_grains(i).type_keywords, 'pol', /bool)
87   - ;stop
  90 + stop
88 91  
89 92 IF !run_univ eq 1 THEN BEGIN
90 93 st = replicate(st(i),Ngrains)
... ...
src/idl/dustem_read_grain.pro
  Diff comments   View file @e69bf24
... ... @@ -197,8 +197,11 @@ CASE !dustem_which OF
197 197 END
198 198 ENDCASE
199 199  
  200 +
200 201 the_end:
201 202  
  203 +;stop
  204 +
202 205 RETURN,full_st
203 206  
204 207 END
... ...
src/idl/dustem_run_readme.pro
  Diff comments   View file @e69bf24
1   -PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it
  1 +PRO dustem_run_readme,postcript=postcript,model=model,help=help,plot_it=plot_it
2 2  
3 3 ;+
4 4 ; NAME:
... ... @@ -8,7 +8,7 @@ PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it
8 8 ; CATEGORY:
9 9 ; Dustem
10 10 ; CALLING SEQUENCE:
11   -; dustem_run_readme,mode=mode,help=help,/postcript,/help,/plot_it
  11 +; dustem_run_readme,model=model,help=help,/postcript,/help,/plot_it
12 12 ; INPUTS:
13 13 ; None
14 14 ; OPTIONAL INPUT PARAMETERS:
... ... @@ -18,7 +18,7 @@ PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it
18 18 ; OPTIONAL OUTPUT PARAMETERS:
19 19 ; None
20 20 ; ACCEPTED KEY-WORDS:
21   -; mode = Selects one of the dust mixture used by dustem
  21 +; model = Selects one of the dust mixture used by dustem
22 22 ; 'COMPIEGNE_ETAL2010' from Compiegne et al 2010 (default)
23 23 ; 'DBP90' from Desert et al 1990
24 24 ; 'DL01' from Draine & Li 2001
... ... @@ -35,7 +35,7 @@ PRO dustem_run_readme,postcript=postcript,mode=mode,help=help,plot_it=plot_it
35 35 ; PROCEDURE:
36 36 ; None
37 37 ; EXAMPLES
38   -; dustem_run_readme,mode='COMPIEGNE_ETAL2010'
  38 +; dustem_run_readme,model='COMPIEGNE_ETAL2010'
39 39 ; MODIFICATION HISTORY:
40 40 ; Written by J.P. Bernard April 1st 2011
41 41 ; see evolution details on the dustem cvs maintained at CESR
... ... @@ -53,14 +53,14 @@ ENDIF
53 53  
54 54 loadct,13
55 55  
56   -IF keyword_set(mode) THEN BEGIN
57   - use_mode=strupcase(mode)
  56 +IF keyword_set(model) THEN BEGIN
  57 + use_model=strupcase(model)
58 58 ENDIF ELSE BEGIN
59   - use_mode='COMPIEGNE_ETAL2010' ;Default is last dustem model
  59 + use_model='COMPIEGNE_ETAL2010' ;Default is last dustem model
60 60 ENDELSE
61 61  
62 62 ;== This just initializes
63   -dustem_init,mode=use_mode
  63 +dustem_init,model=use_model
64 64 !dustem_verbose=1
65 65 !dustem_show_plot=1
66 66  
... ... @@ -81,18 +81,18 @@ ENDIF
81 81 ;dir_in=!dustem_soft_dir+'/src/'+dir
82 82 ;== The following structure contains the inputs to the model
83 83 st_model=dustem_read_all(dir_in)
84   -help,st_model,/str
  84 +;help,st_model,/str
85 85  
86 86 ;=== You can here modify the model inputs
87 87 ;=== e.g.
88 88 ;stop
89 89 ;st_model.G0=2.
90   -st_model.grains[0].type_keywords='logn-chrg-spin'
91   -st_model.gas.NH=1.e2
92   -st_model.gas.Tgas=1.e5
  90 +;===== This below is to test spinning. Does not work with DBP90 !
  91 +;st_model.grains[0].type_keywords='logn-chrg-spin'
  92 +;st_model.gas.NH=1.e2
  93 +;st_model.gas.Tgas=1.e5
93 94  
94 95 ;== The following line writes inputs to the Fortran
95   -;dustem_write_all,st_model,getenv('DUSTEM_DAT')
96 96 dustem_write_all,st_model,!dustem_dat
97 97 ;== This runs the Fortran. The ouput structure contains the results
98 98 st=dustem_run()
... ... @@ -109,12 +109,21 @@ st=dustem_run()
109 109  
110 110 ;== The following just plots the resulting emission SED
111 111 window,0
112   -;dustem_plot_nuinu_em,st,yr=[1e-28,1.e-22],/ysty,xr=[1,5e3],/xsty
113   -dustem_plot_nuinu_em,st,yr=[1e-11,6.e-7],/ysty,xr=[1,5e3],/xsty
  112 +xtit=textoidl('\lambda (\mum)')
  113 +;=== same plot as in Compiegne et al. 2010
  114 +ytit=textoidl('\nuI_\nu^{em} (W/m^2/sr for N_H=1.e20 H/cm^2)')
  115 +tit='Dustem Emitted Intensity '+use_model
  116 +yr=[1e-11,6.e-7]
  117 +xr=[1,5e3]
  118 +dustem_plot_nuinu_em,st,yr=yr,/ysty,xr=xr,/xsty,/xlog,/ylog,title=tit
114 119 ;== The following just plots the resulting extinction
115 120 window,1
116   -;dustem_plot_extinction,st,st_model,yr=[1.e-4,1.e7],/ysty,xr=[1.e-2,1e5],/xsty
117   -dustem_plot_extinction,st,st_model,yr=[1.e-6,1.e0],/ysty,xr=[1.,400],/xsty
  121 +yr=[0,2.5] ;range of 1/lambda
  122 +xr=[1,10] ;range of sigma
  123 +xtit=textoidl('1/\lambda (\mum^{-1})')
  124 +ytit=textoidl('\sigma_{ext} (1e-21 cm^2/H)')
  125 +tit='Dustem extinction '+use_model
  126 +dustem_plot_extinction,st,st_model,xr=xr,yr=yr,/xsty,/ysty,xtit=xtit,ytit=ytit,title=tit
118 127  
119 128 the_end:
120 129  
... ...
src/idl/dustem_write_grain.pro
  Diff comments   View file @e69bf24
1 1 PRO dustem_write_grain,file,st
2 2  
  3 +;stop
  4 +
3 5 OPENW,unit,file,/get_lun
4 6  
5 7 Ntypes=n_elements(st)
... ...
src/idl/dustem_write_grain_web3p8.pro
  Diff comments   View file @e69bf24
... ... @@ -3,6 +3,8 @@ PRO dustem_write_grain_web3p8,file,st
3 3 ;Remark: It is a bit odfd that this is using both st and (*!dustem_params)
4 4 ;Remark: This should use automated comment detection, as in dustem_write_gas.pro
5 5  
  6 +;stop
  7 +
6 8 Ncomments=100
7 9 c=strarr(Ncomments)
8 10  
... ... @@ -39,7 +41,7 @@ printf,unit,(*!dustem_params).G0
39 41 ;frmt='(A12,I4,A12,11E14.6)'
40 42 ;frmt='(A12,I4,A20,11E14.6)'
41 43 ;frmt='(A20,I4,A20,11E14.6)'
42   -frmt='(A20,I4,A25,11E14.6)'
  44 +frmt='(A30,I4,A25,11E14.6)'
43 45  
44 46 FOR i=0L,(*!dustem_params).ngrains-1 DO BEGIN
45 47 printf,unit,st(i).grain_type, $
... ...
src/idl/dustem_write_qpol.pro
  Diff comments   View file @e69bf24
... ... @@ -18,12 +18,13 @@ FOR i=0L,Nfiles-1 DO BEGIN
18 18  
19 19 ;filename=dir+'Q'+strtrim(i_axis,2)+'_'+st.grain.grains(i).grain_type+Estring
20 20 filename=dir+'Q'+strtrim(i_axis,2)+'_'+st.grains(i).grain_type+Estring
  21 + message,'Will write '+filename,/continue
21 22 Nsizes=n_elements((*st.qabs(i)).sizes)
22 23 OPENW,unit,filename,/get_lun
23 24 printf,unit,Nsizes
24   - printf,unit,(*st.qabs(i)).sizes,format=frmt
25   - NLines=n_elements((*st.qpol(i_axis-1,i)).qabs)
26   - NQabs=n_tags((*st.qpol(i_axis-1,i)).qabs)
  25 + printf,unit,(*st.qabs[i]).sizes,format=frmt
  26 + NLines=n_elements((*st.qpol[i_axis-1,i]).qabs)
  27 + NQabs=n_tags((*st.qpol[i_axis-1,i]).qabs)
27 28 vars=fltarr(NQabs,Nlines)
28 29  
29 30 ;Separator
... ...