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]
TVLCT, 53, 156, 83, 100
;Define parameters
nbins=50 ;choose number of bins for histogram
factor=15.*1.e-5
var1=10
microm=1.e-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)
a0=start.a0_initdist
fac=a0/n0
print,'n0=',n0
;
proc=[10, 1000, 2000, 4000]
colors=[0, 122, 55, 100]
;
!x.title='!8r!6 [!7l!6m]'
!y.title='!8f!6 !8r!6!dini!n/!8n!6!d0!n'
yr=[1.e-6,1.e2]
xr=[3.,1.e2]
plot_oo,proc,/nodata,yr=yr,xr=xr
;
for i=0,n_elements(proc)-1 do begin
		iproc=proc[i]
		pc_read_pvar,obj=pvar,ivar=var1,proc=iproc
    save, file='pvar_ipro'+str(iproc),pvar
		pdfi_n,pvar.a,pvar.npswarm,xx,yy,n=nbins-1,/log
		;oplot,xx,yy,li=1,col=colors[i]
		norma=total(yy)/n0
		dlnad=alog(xx[3])-alog(xx[2])
		oplot,xx/microm,fac*yy/norma/(dlnad*xx),col=colors[i]
		oplot,xx/microm,fac*yy/norma/(dlnad*xx),col=colors[i],ps=1
		print,'max and minimum radius in the selected subvolume:',minmax(pvar.a)
endfor
;
;
; The global size spectra at t=10s (ivar=10)
run='SW512condens0_coag_grav0_turb_shima_f2em2_ngp_a24_L025'
restore,'../'+run+'/data/radii_norm.sav'
restore,'../'+run+'/data/prob_norm.sav'
radii_norm1=radii_norm
prob_norm1=prob_norm
dlnad=alog(radii_norm1[3,1])-alog(radii_norm1[2,1])
xx0=radii_norm1[*,1]
yy0=(fac/microm)*prob_norm1[*,1]/(dlnad*xx0)
oplot,xx0,yy0,li=1
;
print,'mv idl.ps ~/tex/nils/particles_coag/fig/f_local_subvolume_comp.ps'
END