;
; 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
;
; 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='../../gravity_notemp/'
    restore,file=notempdir+notemprun+'/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='./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
;
; 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
    ;
    ; 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(-0.5*z0*Wg[iap]/Dt)
;        rhs=exp(-0.4*z0*Wg[iap]/Dt)
        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=2
        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=min(rhs)*0.5
                ymax=max(rhs)*2
                ymin=0.01^(iap+1)

            end
        end
        title='St='+str(Stokes[iap])+'  !4a!6='+str(alpha0)
        if (relative eq 0) then begin
            plot_io,zbin,rhs/rel,yr=[ymin,ymax],li=2,$
              ytitle='!8N/N!d0!n!6',xtitle='!8z!6',tit=title
        endif else begin
            plot,zbin,rhs/rel,yr=[ymin,ymax],li=2,$
              ytitle='!8N/N!d0!n!6',xtitle='!8z!6',tit=title
        end
        oplot,zbin,nograv_result/rel
        oplot,zbin,result_nocorrected/rel,li=0,col=122
        oplot,zbin,result/rel,li=3,col=47
    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
;
; Save data to file
;
if (lgrav) then begin
    save,file='./data/ppar_real.dat',kf,Re,Stokes,alpha,grav
endif else begin
    save,file='./data/ppar_real.dat',kf,Re,Stokes,alpha,grav,alpha2
end
;
;
;
if !d.name eq 'PS' then begin
    if (lgrav) then begin
        print,'mv idl.ps /home/vsl175/a/nilshau/tex/axel/nathan_igor/nils/fig/al_grav.ps'
    endif else begin
        print,'mv idl.ps /home/vsl175/a/nilshau/tex/axel/nathan_igor/nils/fig/al.ps'
    end
endif
;
END