;;----define the gamma distribution
;function meanT, shape, slope, N0, r
;  return, N0*r^shape*exp(-slope*r)/gamma(shape+1)
;end


;restore,'vpz.sav'
;restore,'time.sav'
;
microm=1.e6
;
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
a0=start.a0_initdist
n0=npar*np_swarm0/(nx*ny*nz)
ap0=start.ap0[0]
ap1=start.ap1
;
pc_read_pstalk,obj=stalk
vpz=stalk.vpz
;
;tau0_inverse=n0*(ap0+ap1)^2*abs(vpz[0,1]-vpz[1,1])
tau0_inverse=n0*(ap0+ap1)^2*abs(vpz[0,2]-vpz[1,2])
tau0=1./tau0_inverse
;
print, 'tau0_inverse=', tau0_inverse
print, 'tau0=', tau0
;
END