;
; This script calculates the Stokes number, alpha etc.
;
size=3
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
;
datadir='./forcing_at_k5/gravity/g001_mu/data/'
pc_read_param,obj=start,datadir=datadir
pc_read_param,obj=run,/param2,datadir=datadir
pc_read_dim,obj=dim,datadir=datadir
pc_read_kf,kf,datadir='./forcing_at_k5/gravity/g001_mu'
restore,file='./forcing_at_k5/gravity/g001_mu/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
;
; 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
    gravind=5
    notempdir='./forcing_at_k5/gravity_notemp/g001_mu/'
    restore,file=notempdir+'/data/N_N0.dat'
    Nz_notemp=bin_aver
    for iap=0,npart_radii-1 do begin
        Nz_notemp[iap]=Nz_notemp[iap]/bin_aver[gravind,iap]
    end
endif
;
; Finally the data file can be read in
;
restore,file='./forcing_at_k5/gravity/g001_mu/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.
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=dim.nx
    ;
    ; 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
    ;
    int=fltarr(nn)
    dz=zbin[2]-zbin[1]
    Wg=grav*tau_p
    Dt=urms/(3*kf)
    !p.multi=[0,1,npart_radii]
    !y.style=1
    alpha1=fltarr(npart_radii)
    alpha1[*]=1.0
    alpha1[4]=0.5
    alpha1[3]=0.8
    stokes_tit=strarr(npart_radii)
    stokes_tit[0]=' 0.05'
    stokes_tit[1]=' 0.20'
    stokes_tit[2]=' 0.90'
    stokes_tit[3]=' 3.50'
    stokes_tit[4]='14.00'
    ;
    ymin_rel1=[0.,0.,0.,0.,0.0]
    ;
    ; Loop over all particle sizes
    ;
    for iap=0,npart_radii-1 do begin
        ;
        ; Find some constants
        ;
        Nz=bin_aver[*,iap]/bin_aver[gravind,iap]
        int[0]=0.
        for i=1,nn-1 do begin
            int[i]=int[i-1]+Wg[iap]/Dt*dz
         end        
        alpha0=alpha1[iap]
        alpha=alpha1
        rhs=exp(-int)
        z0=zbin-zbin[0]
        rhs=exp(-z0*Wg[iap]/Dt)
        lhs=Nz*rhointer^(-alpha0)
        ;
        ; Define some help variables
        ;
        result_nocorrected=Nz*rhs[10]/Nz[10]
        result=lhs*rhs[10]/lhs[10]
        nograv_result=Nz_notemp[*,iap]*rhs[10]/Nz_notemp[10,iap]
        ;
        ; Do we want to plot the results relative to analytical
        ; solution for better comparison?
        ;
        relative=0
        if (relative eq 1) then begin
            rel=rhs
            ymin=0.5
            ymax=2.
        endif else begin
            if (relative eq 2) then begin
                rel=nograv_result
                ymin=1-0.3*(iap+1)
                ymax=1+0.6*(iap+1)^0.8
                ymax=max(result_nocorrected[10:*]/rel[10:*])*1.5
            endif else begin
                rel=1.
                ymin=0.
                ymax=max(rhs)*2
            end
         end
        !x.range=[-3,3]
        !x.style=1
        title='St='+Stokes_tit[iap]
        if (relative eq 0) then begin
            plot,zbin,result_nocorrected/rel,yr=[ymin,ymax],li=0,$
              ytitle='!8N/N!d0!n!6',xtitle='!8z!6',tit=title
        endif else begin
           if (relative eq 2) then begin
              plot,zbin,result_nocorrected/rel,yr=[ymin,ymax],li=0,$
                   ytitle='!8N/N!d!6iso!n',xtitle='!8z!6',tit=title
           endif else begin
              plot,zbin,result_nocorrected/rel,yr=[ymin,ymax],li=0,$
                   ytitle='!8N/N!d0!n!6',xtitle='!8z!6',tit=title
           end
        end
        oplot,zbin,nograv_result/rel,li=1
        if (relative ne 0) then begin
           oplot,zbin,rhs/rel,li=2
           oplot,zbin,result/rel,li=3
        end
    end
;     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
if (not lgrav) then begin
    print,'alpha2=',alpha2
end
;
; Produce output for latex file
;
print,'       Res   kf   urms         nu     Re     g    rho_max      N_max           St      alpha'
for i=0,npart_radii-1 do begin
    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[i],FORMAT='($,1F10.6,"&")'
    print,Stokes[i],FORMAT='($,1F11.6,"&")'
    print,alpha[i],FORMAT='(1F10.6,"\\")'
end
;
;
;
if !d.name eq 'PS' then begin
   print,'mv idl.ps ~/tex/axel/nathan_igor/nils/fig/fig4.ps'
endif
;
!x.range=''
;
END