device,decompose=0
loadct,5
;
; Read in data
;
pc_read_ts,obj=ts
pc_read_param,obj=start
pc_read_param,obj=run,/param2
pc_read_kf,kf
@data/index
;
; Find number of printouts
;
dims=size(ts.t)
nt=dims(1)
;
default,lradius,1
print,'lradius=',lradius
if (lradius) then begin
    ad=start.ap0
    adnonzero=where(ad ne 0)
    dims=size(adnonzero)
    nbins=dims[1]
;    print,'nbins=',nbins
endif else begin
;
; Find mass of the different bins
;
    rhods=start.rhods
    deltamd=start.deltamd
    ad1=start.ad1
    md0 = 8*!pi/(3*(1.+deltamd))*ad1^3*rhods
    onethird=1./3.
    mdminus=indgen(ndustspec)*1.0
    mdplus=mdminus
    md=mdminus
    for k=0,ndustspec-1 do begin
        mdminus(k) = md0*deltamd^(k)
        mdplus(k)  = md0*deltamd^(k+1)
        md(k) = 0.5*(mdminus(k)+mdplus(k))
    end
    ad=(0.75*md/(!pi*rhods))^onethird
    print,'md0=',md0
    print,'md=',md
    nbins=ndustspec
end
;print,'ad=',ad
;
; Find the Reynolds number
;
Lx=start.xyz1[0]-start.xyz0[0]
urms=mean(ts.urms(nt/2:nt-1))
nu=run.nu(0)
forcing_scale=Lx/(2*!pi*kf)
Re=urms*forcing_scale/nu
plot,ts.t,ts.urms,ytit='urms [m/s]',xtit='t [s]'
oplot,[ts.t(nt/2),ts.t(nt-1)],[urms,urms],li=2,col=122
tau_integral=forcing_scale/urms
;
; Find Kolmogorov time scale
;
tau_k=forcing_scale/(urms*sqrt(Re))
;print,'tau_k=',tau_k
;
; Find particle time scale
;
density_ratio=1000.
dp=2.*ad
tau_p=density_ratio*dp^2/(18.*nu)
;print,'dp=',dp
;print,'tau_p=',tau_p
;
; Find Stokes number of the particles
;
Stokes=tau_p/tau_k
Stokes_integral=tau_p/tau_integral
;print,'Stokes number=',Stokes
;
; Find Mach number
;
umax=mean(ts.umax(nt/2:nt-1))
Mach_max=umax/start.cs0
Mach_rms=urms/start.cs0
print,' '
print,' '

print,'Mach_max=',Mach_max
print,'Mach_rms=',Mach_rms
print,'Re      =',Re
;
print,' '
print,' '
print,fo='(A7,A19,2A16)','index','radius [micro m]','Stokes_kol','Stokes_int'
print,'----------------------------------------------------------'
for i=0,nbins-1 do begin
    print,fo='(i7,F19.0,2F16.3)',i,ad(i)*1e6,Stokes(i),Stokes_integral(i)
end
;
save,file='data/stokes.dat',ad,Stokes,Stokes_integral

END