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