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=2 ; this should be always 2 because of the bug in "particles_stalker.f90"
deltaivar=1
;
;read stalk file-----------------------------------------------------------
pc_read_pstalk,obj=stalk,quiet=quiet
ap=stalk.ap
npswarm=stalk.npswarm
t=stalk.t
save,file='stalk.sav',ap,npswarm,t
ivar1=n_elements(stalk.t)-1
pdfi_n,stalk.ap[*,0],stalk.npswarm[*,0],xx,yy,n=nbins-1,/log
norma=total(yy)/n0
radii=fltarr(nbins,ivar1+1)
radii_norm=radii
radii_norm[*,0]=xx/microm
prob=radii 
prob_norm=radii
prob_norm[*,0]=yy/norma
for ivar=ivar0,ivar1 do begin
    print,'t=',stalk.t[ivar]
    pdfi_n,stalk.ap[*,ivar],stalk.npswarm[*,ivar],xx,yy,n=nbins-1,/log
    radii[*,ivar-1]=xx
    prob[*,ivar-1]=yy
    radii_norm[*,ivar-1]=xx/microm
    ;prob_norm[*,ivar]=smooth(yy,10)/norma
    prob_norm[*,ivar-1]=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