dustem_plot_sdist.pro
2.26 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
PRO dustem_plot_sdist,dir_out=dir_out,_Extra=extra,help=help,ps=ps,win=win,xr=xr,yr=yr,mergePAH=mergePAH,donotclose=donotclose,multi=multi
nsz_max = 200
if not keyword_set(dir_out) then dir_out=!dustem_dat
if not keyword_set(multi) then multi = 0
; Size distribution
file=dir_out+'out/'+'SDIST.RES'
OPENR, unit, file, /get_lun
tt = '#'
WHILE STRPOS(tt,'#') EQ 0 do begin
READF, unit, tt
ENDWHILE
ntype=uint(tt)
nsize=intarr(ntype+1)
ax = dblarr(ntype,nsz_max)
ava = dblarr(ntype,nsz_max)
FOR i=0,ntype-1 do begin
READF,unit,tt
sp=strsplit(tt,/extract)
if nsize(i) eq 0 then begin
nsize(i)=uint(sp(1))
endif
for is=0,nsize(i)-1 do begin
READF, unit, tx, ty
ax(i,is) = tx
ava(i,is) = ty
endfor
ENDFOR
close,unit
free_lun,unit
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
if not keyword_set(ps) then tit='Dust size distributions'
yscl = 1.D29
yr /= yscl
yscl = 1.
;xtit='a (nm)'
xtit=textoidl('a (\mu m)')
norm_size = 1.d7
norm_size = 1.d4
xr /=1d3
;ytit='Mass log( 10!u29!n n!dH!u-1!n a dm/da ) '
ytit=textoidl('4\pi/3 a^3dn/d ln a (cm^3/H)')
if multi eq 0 then plot_oo, xr, /nodata,xr=xr,/xs,xtit=xtit,/ys,yr=yr,ytit=ytit,tit=''
if multi eq 1 then plot_oo, xr, /nodata,xr=xr,/xs,xtit='',/ys,yr=yr,ytit=ytit,tit='',position=[0.17,0.45,0.95,0.95],xtickformat='(A1)'
colors=dustem_grains_colors(ntype,ls,ps=ps)
;stop
FOR i=0,ntype-1 DO begin
;if i le 1 and keyword_set(mergePAH) then begin
; if i eq 0 then oplot, ax(0,0:nsize(0)-1)*1.e7, alog10(ava(0,0:nsize(0)-1)*yscl)+alog10(ava(1,0:nsize(1)-1)*yscl), line=ls(1),color = colors(1)
;endif else begin
; oplot, ax(i,0:nsize(i)-1)*1.e7, alog10(ava(i,0:nsize(i)-1)*yscl), line=ls(i+1),color = colors(i+1)
;endelse
if i le 1 and keyword_set(mergePAH) then begin
if i eq 0 then oplot, ax(0,0:nsize(0)-1)*norm_size, ava(0,0:nsize(0)-1)+ava(1,0:nsize(1)-1), line=ls(1),color = colors(1)
endif else begin
oplot, ax(i,0:nsize(i)-1)*norm_size, ava(i,0:nsize(i)-1), line=ls(i+1),color = colors(i+1)
endelse
END
if keyword_set(ps) and not keyword_set(donotclose) then begin
device,/close
set_plot,'X'
endif
END