;
; This script calculates the Stokes number, alpha etc.
;
size=2
if !d.name eq 'PS' then begin
   device,xsize=20,ysize=30,yoffset=3
   !p.charthick=4 & !p.thick=4 & !x.thick=4 & !y.thick=4
   !p.charsize=size
endif else begin
   !p.charsize=size
   device,decompose=0
end
loadct,5
;
; Read data from files
;
pc_read_param,obj=start
pc_read_param,obj=run,/param2
pc_read_dim,obj=dim
pc_read_kf,kf
restore,file='./data/rho_rho0.dat'
rhotilde=density_ratio
;
; Find how many different particle radii we have
;
npart_radii=0
dims=size(start.ap0)
if (dims[0] gt 0) then begin
    ninit=dims[1]
    for i=0,ninit-1 do begin
        if (start.ap0[i] ne 0) then begin
            npart_radii=npart_radii+1
        endif
    end
endif else begin
    npart_radii=1
end
print,'npart_radii=',npart_radii
;
; Finally the data file can be read in
;
restore,file='./data/N_N0.dat'
npart_radii=n_elements(average)
;
; Set some defaults
;
S=1000
;
; Find the Stokes number
;
d=start.ap0[0:npart_radii-1]*2.
if (run.ivisc[0] eq 'mu-therm') then begin
    cs=run.cs0
    cp=run.cp
    gamma=run.gamma
    temp=cs^2/(cp*(gamma-1))
    nu=run.nu*sqrt(temp)
    print,'nu=',nu
endif else begin
    nu=run.nu
end
tau_p=S*d^2/(18*nu)
tau_f=1./(urms*kf)
Re=urms/(nu*kf)
tau_k=tau_f/sqrt(Re)
Stokes=tau_p/tau_k
print,'Stokes=',Stokes
print,'Ntilde=',Ntilde
;
; Find veights for interpolation
;
Stokes_final=[0.1,1.,10.,14]
dims=size(Stokes_final)
nsf=dims[1]
weight=fltarr(npart_radii,nsf)
for i=0,nsf-1 do begin
    for ipart=0,npart_radii-2 do begin
        if ((Stokes(ipart+1) gt Stokes_final[i]) and (Stokes(ipart) lt Stokes_final[i])) then begin
            delta_stokes=Stokes(ipart+1)-Stokes(ipart)
            weight[ipart,i]=(Stokes(ipart+1)-Stokes_final[i])/delta_stokes
            weight[ipart+1,i]=(Stokes_final[i]-Stokes(ipart))/delta_stokes
            print,i,ipart
        endif
    end
end
;
; Find interpolated values
;
Ntilde_final=fltarr(nsf)
for i=0,nsf-1 do begin
    Ntilde_final[i]=total(weight[*,i]*Ntilde)
end
print,'Ntilde_final=',Ntilde_final
;
; Plot results
;
plot_oi,Stokes,Ntilde,yr=[1,1.6],ystyle=1,ps=-1
oplot,Stokes_final,Ntilde_final,ps=2,color=122



;
; Find width of profile
;
width=fltarr(npart_radii)
for ipart=0,npart_radii-1 do begin
    NN=bin_aver[*,ipart]
    maxNN=max(NN)
    eNN=(max(NN)-min(NN))*0.368+min(NN)
;    print,'maxNN=',maxNN
;    print,'eNN=',eNN
    for k=0,dim.nz-2 do begin
        if ((NN(k+1) gt eNN) and (NN(k) lt eNN)) then begin
            delta_NN=NN(k+1)-NN(k)
            weight_1=(NN(k+1)-eNN)/delta_NN
            weight_2=(eNN-NN(k))/delta_NN
            NNe=weight_1*NN(k)+weight_2*NN(k+1)
            ze =weight_1*zbin(k)+weight_2*zbin(k+1)
;            print,'ez=',ze,NNe
            width[ipart]=-ze
        end
    end
end
;
; Find interpolated values of profile width
;
width_final=fltarr(nsf)
for i=0,nsf-1 do begin
    width_final[i]=total(weight[*,i]*width)
end
print,'width_final=',width_final
;
; Find alpha
;
restore,file='./data/ppar_real.dat'
alpha_final=fltarr(nsf)
for i=0,nsf-1 do begin
    alpha_final[i]=total(weight[*,i]*alpha)
end
print,'alpha_final=',alpha_final



END