;
; This script calculates the Stokes number, alpha etc.
;
size=2
if !d.name eq 'PS' then begin
   device,xsize=30,ysize=16,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
;
; 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
;
; Check if there is gravity
;
lgrav=0
grav=0.
if start.gravz gt 0 then begin
    lgrav=1
    grav=start.gravz
endif
;
; Read in reference run without temperature profile
;
lnocoolcase=0
if (lgrav and run.lcooling_general) then begin
    lnocoolcase=1
endif
if lnocoolcase then begin
@distfile
    notempdir='../../gravity_notemp/'
    restore,file=notempdir+notemprun+'/data/N_N0_real.dat'
    Nz_notemp=bin_aver/bin_aver[5]
endif
;
; Finally the data file can be read in
;
restore,file='./data/N_N0_real.dat'
Ntilde=average
;
; Set some defaults
;
S=1000
;
; Find the Stokes number
;
d=start.ap0[0]*2.
nu=run.nu
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
;
; Find alpha
;
if lgrav then begin
    nn=100
    ;
    ; First we must interpolate the density values to the 
    ; points used for the particle density
    ;
    rhointer=fltarr(nn)
    rhointer[0]=rhoaver[0]
    dzaver=zaver[2]-zaver[1]
    for i=1,nn-1 do begin
        zval=zbin[i]
        ind=where(zaver gt zval)
;        print,ind[0]
        dz1=zval-zaver[ind-1]
        dz2=dzaver-dz1        
        rhointer[i]=(dz1*rhoaver[ind]+dz2*rhoaver[ind-1])/dzaver
    end
    ;
    ; Find some constants
    ;
    Wg=grav*tau_p
    Dt=urms/(3*kf)
    Nz=bin_aver/bin_aver[5]
    int=fltarr(nn)
    int[0]=0.
    dz=zbin[2]-zbin[1]
    for i=1,nn-1 do begin
        int[i]=int[i-1]+Wg/Dt*dz
    end
    !p.multi=[0,1,1]
    !y.style=1
    alpha0=1.1
    rhs=exp(-int)
    z0=zbin-zbin[0]
    rhs=exp(-0.5*z0*Wg/Dt)
    lhs=Nz*rhointer^(-alpha0)
    plot_io,zbin,rhs,yr=[min(rhs)*0.5,max(rhs)*2],li=2,ytitle='!8N/N!d0!n!6',xtitle='!8z!6'
    oplot,zbin,lhs*rhs[10]/lhs[10]
    if lnocoolcase then begin
        oplot,zbin,Nz_notemp,col=47
        oplot,zbin,Nz,col=122
    endif
;     dalpha=0.055
;     for ia=1,20 do begin
;         alpha=alpha0+dalpha*ia
;         lhs=Nz*rhointer^(-alpha)
;         oplot,zbin,lhs*rhs[30]/lhs[30],ps=-1,col=122
;     end
    !p.multi=[0,1,1]
endif else begin
    alpha=alog(Ntilde)/alog(rhotilde)
    alpha2=(1-Ntilde)/(1/rhotilde-1)
end
print,'alpha=',alpha
print,'alpha2=',alpha2
;
; Produce output for latex file
;
print,'       Res   kf   urms         nu     Re     g    rho_max      N_max          St      alpha'
print,dim.nx,FORMAT='($,I10,"&")'
print,kf,FORMAT='($,1F4.1,"&")'
print,urms,FORMAT='($,1F6.3,"&")'
print,nu,FORMAT='($,1F10.6,"&")'
print,Re,FORMAT='($,1F6.0,"&")'
print,grav,FORMAT='($,1F6.0,"&")'
print,rhotilde,FORMAT='($,1F10.6,"&")'
print,Ntilde,FORMAT='($,1F10.6,"&")'
print,Stokes,FORMAT='($,1F10.6,"&")'
print,alpha,FORMAT='(1F10.6,"\\")'
;
; Save data to file
;
save,file='./data/ppar_real.dat',kf,Re,Stokes,alpha,grav,alpha2
;
;
;
if !d.name eq 'PS' then begin
    print,'mv idl.ps /home/vsl175/a/nilshau/tex/axel/nathan_igor/nils/fig/al.eps'
endif
;
END