import pencil as pc
import pylab as P
import numpy as N
from misc import deriv

par  = pc.read_param(quiet=True)
par2 = pc.read_param(param2=True,quiet=True)
dim  = pc.read_dim()
var  = pc.read_var(ivar=0, run2D=True, quiet=True, magic=['ss','rho','TT'], 
                  trimall=True, param=par)
grid = pc.read_grid(param=par, quiet=True, trim=True)

Lx     = grid.Lx
rhom   = var.rho.mean(axis=1)
tempm  = var.TT.mean(axis=1)
ssm    = var.ss.mean(axis=1)
derssm = deriv(var.z, ssm)
#
length = par.lxyz[2]
grav   = par.gravz
nu     = par2.nu/rhom
kappa  = par2.hcond0
chicp  = kappa/rhom

# Gough's version (ApJ, 206, 1976)
#
#       g d^4    d(s/cp)
# Ra= --------- --------
#       nu chi   dz
#
Ra = grav*length**4 / (nu*chicp)*derssm
Ra_mid=Ra[dim.nz/2]
print 'Rayleigh mid=%.1f'%Ra_mid

#
# theoretical Ra value for a polytrope
#
xi0 = abs(par.cs0**2/(par.gamma*grav*length))
Ra_th = -(1.-(par.gamma-1.)*(par.mpoly0+1.)/par.gamma)*grav*length**3/ \
        ((0.5+xi0*(par.mpoly0+1.))*(nu*chicp)[dim.nz/2])

print '----------------------'
print 'Rayleigh theo mid=%.1f'%Ra_th