dustem_plot_ext.pro
5.49 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
PRO dustem_plot_ext,st,ps=ps,_Extra=extra,help=help,UV=UV,win=win,print_Rv=print_Rv,mergePAH=mergePAH,ylog=ylog,xlog=xlog,donotclose=donotclose,yr=yr,ysty=ysty,xr=xr,xsty=xsty,noerrbars=noerrbars,albedo=albedo,multi=multi
;+
; NAME:
; dustem_plot_ext
; PURPOSE:
; Plots a Dustem extinction
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; dustem_plot_ext,st,,/help][,_extra=]
; INPUTS:
; st = Dustem model output structure
; OPTIONAL INPUT PARAMETERS:
; _extra = extra parameters for the plot routine
; OUTPUTS:
; None
; OPTIONAL OUTPUT PARAMETERS:
; None
; ACCEPTED KEY-WORDS:
; help = If set, print this help
; COMMON BLOCKS:
; None
; SIDE EFFECTS:
; model is plotted
; RESTRICTIONS:
; The dustem idl wrapper must be installed
; PROCEDURE:
; None
; EXAMPLES
; dustem_init,/wrap
; dir_in=getenv('DUSTEM_SOFT_DIR')+'/src/dustem3.8_web/'
; dir=getenv('DUSTEM_RES')
; st=dustem_read_all_res(dir_in)
; loadct,13
; dustem_plot_ext,st,yr=[1.e-6,1.e0],/ysty,xr=[1.,400],/xsty
; MODIFICATION HISTORY:
;
; Modified by VG with flags /UV,/IR
;
; Written by J.-Ph. Bernard
; 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_plot_ext'
goto,the_end
ENDIF
if not keyword_set(multi) then multi=0
defsysv,'!dustem_data',exists=data_exists
if keyword_set(print_Rv) then begin
Vband = 0.55
Bband = 0.44
tauV = interpol(st.ext.ext_tot,st.ext.wav,Vband)
tauB = interpol(st.ext.ext_tot,st.ext.wav,Bband)
print
print,"E(B-V) model/NH",1.086*(tauB-tauV)/1d21
print,"Av model",1.086*tauV
print,"Rv model = ",tauV/(tauB-tauV)
print,"NH/Av model = ",1e21/(tauv*1.086)
print
if data_exists then begin
tauV = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Vband)
tauB = interpol((*!dustem_data.ext).values,(*!dustem_data.ext).wav,Bband)
print,"E(B-V) data/NH",1.086*(tauB-tauV)/1d21
print,"Av data",1.086*tauV
print,"Rv data = ",tauv/(tauB-tauV)
print,"NH/Av data = ",1e21/(tauv*1.086)
print
endif
endif else if keyword_set(albedo) then begin
if keyword_set(ps) then begin
set_plot, 'PS'
device, filename=ps, /color, /encapsulated
endif else begin
set_plot,'X'
if keyword_set(win) then window,win
endelse
Ngrains=n_tags(st.sed)-2
abs_tot=st.ext.abs_grain(0)*0.
sca_tot=st.ext.sca_grain(0)*0.
FOR i=0L,Ngrains-1 DO BEGIN
abs_tot=abs_tot+st.ext.abs_grain(i)
sca_tot=sca_tot+st.ext.sca_grain(i)
ENDFOR
xrange=1./st.ext.wav
xr=[0,10]
yr=[0,1]
xtit=textoidl('1/\lambda (\mum^{-1})')
ytit=textoidl('albedo')
plot,xrange,/nodata,xtit=xtit,ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,ys=ys,xr=xr,xsty=xsty
oplot,xrange,sca_tot/(sca_tot+abs_tot)
endif else begin
if multi ne 2 then begin
if keyword_set(ps) then begin
set_plot, 'PS'
device, filename=ps, /color, /encapsulated
endif else begin
set_plot,'X'
if keyword_set(win) then window,win
endelse
endif
Ngrains=n_tags(st.sed)-2
colors=dustem_grains_colors(Ngrains,ls,ps=ps)
; stop
ytit=textoidl('\sigma_{ext}/N_H (cm^2/H)')
if not keyword_set(ps) then tit='Dustem extinction'
abs_tot=st.ext.abs_grain(0)*0.
sca_tot=st.ext.sca_grain(0)*0.
FOR i=0L,Ngrains-1 DO BEGIN
abs_tot=abs_tot+st.ext.abs_grain(i)
sca_tot=sca_tot+st.ext.sca_grain(i)
ENDFOR
if !dustem_show_plot EQ 2 or multi eq 2 then begin
norm = abs_tot+sca_tot
yr=[0.,2.]
ylog = 0
ys = 1
endif else begin
norm = abs_tot * 0. + 1.d21
if keyword_set(UV) then ylog = 0 else ylog = 1
endelse
;VG
if keyword_set(UV) then begin
xtit=textoidl('1/\lambda (\mum^{-1})')
xrange=1./st.ext.wav
xlog = 0
endif else begin
xtit=textoidl('\lambda (\mum)')
xrange=st.ext.wav
xlog = 1
endelse
if multi eq 0 then plot,xrange,/nodata,xtit=xtit,ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,ys=ys,xr=xr,xsty=xsty
if multi eq 1 then plot,xrange,/nodata,xtit='',ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,/ys,xr=xr,xsty=xsty,position=[0.17,0.35,0.95,0.95],xtickformat='(A1)',yticks=3,yminor=10
if multi eq 2 then plot,xrange,/nodata,xtit=xtit,xlog=xlog,ylog=ylog,ys=ys,xsty=xsty,ytit='Normalized',tit='',xr=xr,yr=[0,2],position=[0.17,0.14,0.95,0.35],/noerase,yticks=2,yminor=5,xticklen=0.1
;VG
;VG
; Plot model
FOR i=0L,Ngrains-1 DO BEGIN
if i le 1 and keyword_set(mergePAH) then begin
if i eq 0 then oplot,xrange,(st.ext.abs_grain(0)+st.ext.sca_grain(0)+st.ext.abs_grain(1)+st.ext.sca_grain(1))/norm,color=colors(1),psym=psym,line=ls(1)
endif else begin
oplot,xrange,(st.ext.abs_grain(i)+st.ext.sca_grain(i))/norm,color=colors(i+1),psym=psym,line=ls(i+1)
endelse
ENDFOR
oplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
if !dustem_show_plot EQ 2 or multi eq 2 then begin
norm = interpol(st.ext.ext_tot,st.ext.wav,(*!dustem_data.ext).wav)
endif else begin
norm = st.ext.ext_tot * 0. + 1.d21
endelse
; Plot data
if data_exists then begin
if ptr_valid(!dustem_data.ext) then begin
if keyword_set(UV) then xdata = 1./(*!dustem_data.ext).wav else xdata = (*!dustem_data.ext).wav
if keyword_set(noerrbars) then oplot,xdata,(*!dustem_data.ext).values/norm,psym=4 $
else oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
endif
endif
if keyword_set(ps) and not keyword_set(donotclose) then begin
device,/close
set_plot,'X'
endif
endelse
the_end:
END