dustem_plot_ext.pro
6.96 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
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
PRO dustem_plot_ext,st,ps=ps,UV=UV,win=win,print_Rv=print_Rv,mergePAH=mergePAH, $
ylog=ylog,xlog=xlog,yr=yr,ysty=ysty,xr=xr,xsty=xsty, $
donotclose=donotclose,noerrbars=noerrbars,albedo=albedo,multi=multi, $
_extra=_extra,help=help
;+
; NAME:
; dustem_plot_ext
; PURPOSE:
; Plots a Dustem extinction
; CATEGORY:
; Dustem
; CALLING SEQUENCE:
; dustem_plot_ext,st[,/help]
; INPUTS:
; st = Dustem model output structure
; OPTIONAL INPUT PARAMETERS:
; _extra = extra parameters for the plot routine (currently not passed)
; 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
previous_device=!d.name
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')
;stop
cgplot,xrange,/nodata,xtit=xtit,ytit=ytit,tit=tit,xlog=xlog,ylog=ylog,yr=yr,ys=ys,xr=xr,xsty=xsty
cgoplot,xrange,sca_tot/(sca_tot+abs_tot)
ENDIF ELSE BEGIN
IF multi NE 2 THEN BEGIN
IF keyword_set(ps) THEN BEGIN
previous_device=!d.name
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)
Ngrains=(*!dustem_params).Ngrains
colors=['dark',dustem_grains_colors(Ngrains,ls,/cgplot,ps=ps)] ;dark added because first color taken for total extinction
use_lines=replicate(0,Ngrains)
; stop
ytit=textoidl('\sigma_{ext}/N_H (cm^2/H)')
IF NOT keyword_set(ps) THEN tit='DUSTEMWrapper 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 cgplot,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 cgplot,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 cgplot,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)
IF i EQ 0 THEN cgoplot,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)
cgoplot,xrange,(st.ext.abs_grain(i)+st.ext.sca_grain(i))/norm,color=colors(i+1),psym=psym,line=ls(i+1)
ENDELSE
ENDFOR
;cgoplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
;cgoplot,xrange,(abs_tot+sca_tot+powerlaw)/norm,psym=psym
cgoplot,xrange,(st.ext.ext_tot)/norm,psym=psym
; ; LINES OF CODE TO PLOT THE NEW ADDED PLUGINS TO THE EXTINCTION WILL HAVE TO BE ADDED HERE
;
; cgoplot,xrange,(*(*!dustem_plugin).powerlaw)/norm,psym=psym,color='cadetblue'
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 $
IF keyword_set(noerrbars) THEN BEGIN
cgoplot,xdata,(*!dustem_data.ext).values/norm,psym=4
ENDIF ELSE BEGIN
;else oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
;oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
cgoplot,xdata,(*!dustem_data.ext).values/norm,psym=4
yy=(*!dustem_data.ext).values/norm
rms=3.*(*!dustem_data.ext).sigma/norm/2.
cgerrplot,xdata,yy-rms,yy+rms,color='black'
ENDELSE
;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,previous_device
set_plot,'X'
ENDIF
ENDELSE
the_end:
END