dustem_plot_fpol.pro
2.07 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
PRO dustem_plot_fpol,dir_in=dir_in,_Extra=extra,help=help,ps=ps,win=win,xr=xr,yr=yr,donotclose=donotclose,polfile=polfile,multi=multi
if not keyword_set(polfile) then polfile = ['blabla','blabla']
if not keyword_set(multi) then multi = 0
nsz_max = 200
; Read DUSTEM_wrap fit files
if not keyword_set(dir_in) then dir_in = !dustem_soft_dir ;!dustem_dat because it was looking for a file only there. My might need to be revised I'm not sure
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
st_grains=dustem_read_grain(dir_in+'data/'+(*!dustem_inputs).grain,/silent)
;st_grains=dustem_read_grain(dir_in+'data/'+'GRAIN.DAT',/silent)
;ngrains = st_grains.ngrains I WILL BE REPLACING THIS LINE
ngrains= ((*!dustem_params).ngrains)
colors=dustem_grains_colors(ngrains,ls,ps=ps)
;j=where(polfile ne '', count)
;if count le 1 then begin
;ls = indgen(ngrains+2)*0
;colors = indgen(ngrains+2)*0 + 255
;if keyword_set(ps) then colors *= 0
;endif
if not keyword_set(ps) then tit='Grain alignment efficiency'
ytit='Grain alignment efficiency'
ytit='f(a)'
xtit=textoidl('Grain size (nm)')
amin = 0.1
amax = 1e8
nrad = 20000
a = amin * exp(alog(amax/amin)/(nrad-1)*indgen(nrad))
if multi eq 0 then cgplot,a,a,/nodata,xr=xr,yr=yr,/xs,/ys,xtit=xtit,ytit=ytit,tit=tit,/xlog
if multi eq 2 then cgplot,a,a,/nodata,xr=xr,yr=yr,/xs,/ys,xtit=xtit,ytit=ytit,tit='',/xlog,position=[0.17,0.14,0.95,0.45],/noerase,ymino=2,xticklen=0.07
st_align=dustem_read_align(dir_in+'data/',st_grains,silent=silent)
for i=0,ngrains-1 do begin
if st_align.grains(i).aligned then begin
f_pol = 0.5*st_align.grains(i).plev * (1 + TANH((alog(a*1e-3/(st_align.grains(i).athresh))) / st_align.grains(i).pstiff ) )
cgoplot,a,f_pol, line=ls(i+1),color='black';,col=colors(i+1) ; this line should be tested with multiple polarizing grain species
endif
endfor
if keyword_set(ps) and not keyword_set(donotclose) then begin
device,/close
set_plot,'X'
endif
END