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