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_dat

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
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 plot,a,a,/nodata,xr=xr,yr=yr,/xs,/ys,xtit=xtit,ytit=ytit,tit='',/xlog
if multi eq 2 then plot,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 ) )
		oplot,a,f_pol, line=ls(i+1),col=colors(i+1)
	endif
endfor

if keyword_set(ps) and not keyword_set(donotclose) then begin
        device,/close
        set_plot,'X'
endif

END