if !d.name eq 'PS' then begin
  device,filename='idl.ps',xsize=18,ysize=14,yoffset=3,/color
  !p.charthick=3 & !p.thick=3 & !x.thick=3 & !y.thick=3
end
!p.charsize=2
!x.margin=[10,3.5]
!y.margin=[3.2,.5]
!p.multi=[0,1,1]
;Define parameters
nbins=50 ;choose number of bins for histogram
onethird=1.0/3.0
microm=1e-6
;
;----Obtain the number density of the particles
pc_read_pdim,obj=pdim
npar=pdim.npar
pc_read_dim,obj=dim
nx=dim.nx
ny=dim.ny
nz=dim.nz
pc_read_param,obj=start
np_swarm0=start.np_swarm0
n0=npar*np_swarm0/(nx*ny*nz)
print,'n0=',n0
;
;number density------------------------------------------------------------
ivar0=0
deltaivar=1
;ivar1=118
ivar1=5
restore,'./pvar_file/pvar0.sav'
pdfi_n,a,npswarm,xx,yy,n=nbins-1,/log
norma=total(yy)/n0
radii=fltarr(nbins,ivar1+1)
radii_norm=radii
prob=radii 
prob_norm=radii
for ivar=ivar0,ivar1 do begin
    restore,'./pvar_file/'+'pvar'+str(ivar)+'.sav'
    print,'t=',t
    pdfi_n,a,npswarm,xx,yy,n=nbins-1,/log
    radii[*,ivar]=xx
    prob[*,ivar]=yy
    radii_norm[*,ivar]=xx/microm
    ;prob_norm[*,ivar]=smooth(yy,10)/norma
    prob_norm[*,ivar]=yy/norma
endfor
save,file='./data/radii_norm.sav',radii_norm
save,file='./data/prob_norm.sav',prob_norm
;plot---------------------------------------
xr=[10,5000]
yr=[1e-2,1e8]
!x.title='!8a!6 [!7l!6m]'
!y.title='!8n!dk!n!6 [m!u-3!n]'
!X.TickFormat='exponent'
!Y.TickFormat='exponent'
plot_oo,radii_norm[*,0],prob_norm[*,0],xr=xr,yr=yr,xstyle=1
w=0.2
for ivar=ivar0+1,ivar1,deltaivar do begin
    wait,w
    oplot,radii_norm[*,ivar],prob_norm[*,ivar]
endfor
END