OV-GSO-DC / dustem-wrapper_idl

Blame view

src/idl/dustem_plot_ext.pro 6.88 KB
Edit Raw Normal View History
427f1205   Jean-Michel Glorian   version 4.2 merged 
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
 
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')
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
101
102
103
104
105
 
	stop
	;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)
	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)
36e8b879   Jean-Michel Glorian   update from deborah 
106
427f1205   Jean-Michel Glorian   version 4.2 merged 
107
108
109
110
111
112
113
114
115
116
117
118
119
 
	
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
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
120
121
122
123
124
125
 
	;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)]
	use_lines=replicate(0,Ngrains)
36e8b879   Jean-Michel Glorian   update from deborah 
126
 
;	stop
427f1205   Jean-Michel Glorian   version 4.2 merged 
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
 
	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
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
157
158
159
160
161
162
 
	;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
	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
427f1205   Jean-Michel Glorian   version 4.2 merged 
163
164
165
166
167
168
 
	;VG
	
	;VG
	; Plot model
	FOR i=0L,Ngrains-1 DO BEGIN
		if i le 1 and keyword_set(mergePAH) then begin
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
169
170
 
  			;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)
427f1205   Jean-Michel Glorian   version 4.2 merged 
171
 
		endif else begin
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
172
173
 
  			;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)
427f1205   Jean-Michel Glorian   version 4.2 merged 
174
175
 
		endelse
	ENDFOR
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
176
177
 
	;oplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
	cgoplot,xrange,(abs_tot+sca_tot)/norm,psym=psym
427f1205   Jean-Michel Glorian   version 4.2 merged 
178
179
180
181
182
183
184
185
186
187
188
 
	
	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
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
189
190
191
 
      			;if keyword_set(noerrbars) then oplot,xdata,(*!dustem_data.ext).values/norm,psym=4 $
      			if keyword_set(noerrbars) then cgoplot,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
427f1205   Jean-Michel Glorian   version 4.2 merged 
192
 
        		else oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
cca12c00   Jean-Philippe Bernard   impreoved for ext... 
193
 
        		;oploterror,xdata,(*!dustem_data.ext).values/norm,xdata*0.,(*!dustem_data.ext).sigma/norm,psym=4
427f1205   Jean-Michel Glorian   version 4.2 merged 
194
195
196
197
198
199
200
201
202
203
204
205
206
 
   		endif
	endif
	
	if keyword_set(ps) and not keyword_set(donotclose) then begin
		device,/close
		set_plot,'X'
	endif

endelse

the_end:

END