pc_read_pvar,obj=pvar0,varfile='PVAR0',/quiet
pc_read_pvar,obj=pvar,/quiet
pc_read_param,obj=start
pc_read_pdim,obj=pdim
pc_read_dim,obj=dim

ap0=start.ap0[0]
ap1=start.ap0[1]

nbin0=indgen(3)*0.0
nbin0[0]=total(pvar0.npswarm(where(pvar0.a lt ap0*1.0001)))/dim.nx
nbin0[1]=total(pvar0.npswarm(where((pvar0.a lt ap1*1.0001) and $
                                 (pvar0.a gt ap1*0.9999))))/dim.nx
;nbin0[2]=total(pvar0.npswarm(where(pvar0.a gt ap1*1.0001)))/dim.nx

print,'nbin0=',nbin0

nbin=indgen(3)*0.0
nbin[0]=total(pvar.npswarm(where(pvar.a lt ap0*1.0001)))/dim.nx
nbin[1]=total(pvar.npswarm(where((pvar.a lt ap1*1.0001) and $
                                 (pvar.a gt ap1*0.9999))))/dim.nx
nbin[2]=total(pvar.npswarm(where(pvar.a gt ap1*1.0001)))/dim.nx

print,'nbin=',nbin
dt=pvar.t

dnidt=(nbin(0)-nbin0(0))/dt

print,'dnidt=',dnidt
;
; Find terminal velocity of the two initial particle sizes
;
g=9.81
S=1000.
d=2*start.ap0
nu=1e-5
v=d
for j=0,1 do begin
   tau0=S*d(j)^2/(18*nu)
   v(j)=g*tau0
   for i=1,20 do begin
      Re=d(j)*v(j)/nu
      tau=tau0/(1+0.15*Re^0.687)
      v(j)=g*tau
   end
end
print,'v=',v
;
; Test
;
deltav=v(1)-v(0)
kernel=!pi*deltav*(start.ap0(0)+start.ap0(1))^2
dnidt_test=-kernel*nbin0(0)*nbin0(1)
print,'kernel=',kernel
print,'dnidt_test=',dnidt_test

END