Blog: Update 11/20: adiabatic_profile.py

#!/usr/bin/env python

#Need to add n_p, the neutral density at planet surface to problem.data

import math

Gamma, M_p, T_p, Lambda_p, rho_p, Lx, Nx, max_L = (None for i in range(8))

# read these from astrobear constant file
G = 6.67e-8
k_B = 1.381e-16
m_H = 1.673e-24
sigma = 6.3e-18
M_J = 1.898e30
R_J = 7.1492e9

# read in dx, Gamma, M_p, Cs_p, R_p
with open('physics.data') as infile:
  for line in infile:
    sline = line.strip().split('!', 1)[0]
    if sline is "":
      continue
    if sline.split()[0] == "gamma":
      Gamma = float(sline.split()[2].replace("d", "e"))
with open('problem.data') as infile:
  for line in infile:
    sline = line.strip().split('!', 1)[0]
    if sline is "":
      continue
    sline = sline.split()
    if sline[0] == "planetMass":
      M_p = float(sline[2].replace("d", "e"))*M_J
    if sline[0] == "planetTemp":
      T_p = float(sline[2].replace("d", "e"))
    if sline[0] == "planetRadius":
      R_p = float(sline[2].replace("d", "e"))*R_J
    if sline[0] == "planetDensity":
      rho_p = float(sline[2].replace("d", "e"))
with open('global.data') as infile:
  for line in infile:
    sline = line.strip().split('!', 1)[0]
    if sline == "":
      continue
    sline = sline.split()
    if sline[0] == "GmX":
      Nx = int(sline[2].split(',')[0].replace("d", "e"))
    if sline[0] == "GxBounds":
      Lx = float(sline[2].split(',')[3].replace("d", "e")) \
          -float(sline[2].split(',')[0].replace("d", "e"))
    if sline[0] == "MaxLevel":
      max_L = int(float(sline[2].replace("d", "e")))

if Gamma < 1.01:
  print ("Warning: default gamma = %f, changing to 5/3 for calculation" % Gamma)
  Gamma = 5./3.

Cs_p2 = (k_B*T_p)/m_H
n_p = rho_p/m_H
Nx=Nx*(2**max_L)
dx = Lx/Nx
print Lx, Nx, dx
Gamma_1 = Gamma - 1.0
r_crit = (Gamma_1*G*M_p)/(Gamma*Cs_p2*R_p)
r_zero = (r_crit)/(r_crit-1.0)
n_rho = 1.0/Gamma_1
n_prs = Gamma/Gamma_1
rScale = rho_p
lScale = R_p
pScale = rho_p*Cs_p2

#atmopshere profile
def atm_profile(r):
  return (r_crit*(1./r - 1./r_zero))**(1./Gamma_1)

#radiative transfer
h = 1.0e-6/3.
tau = 0.0
i=0
if r_zero < 0:
  print "Warning: infinitely extended atmosphere"
  r_start = 3.
else:
  r_start = r_zero

while (tau < 1):
  rad = r_start - i*3.*h
  tau += 3.*h/8.*(sigma*n_p*R_p)* \
         ( atm_profile(rad) + 3.*atm_profile(rad-h) + \
           3.*atm_profile(rad-2.*h) + atm_profile(rad-3.*h) )
  i += 1

r_tau1 = 1
r_ib = r_zero/(1+(r_zero-r_tau1)/r_tau1*math.exp(3*Gamma_1))
r_mask = r_ib - 5*Lx/Nx
r_ob = r_zero/(1+(r_zero-r_tau1)/r_tau1*.1**Gamma_1)

with open('adiabat.data', 'w') as outfile:
  nEntries = int(math.ceil((r_zero-r_mask)/dx))
  outfile.write("%d\n" % (nEntries+1))
  outfile.write("r_mask, r_ib, r_tau1, r_ob, r_crit, r_zero\n")
  outfile.write(str("%e %e %e %e %e %e\n" % (r_mask, r_ib, r_tau1, r_ob, r_crit, r_zero)).replace("e",
"d"))
  outfile.write("r/r_p, rho/rho_p, P/P_p\n")
  for i in range(0, nEntries):
    rad = r_mask + i*dx
    f_rad = r_crit*(1.0/rad-1.0/r_zero)
    rho = f_rad**n_rho
    P = f_rad**n_prs
    outfile.write(str("%e" % rad).replace("e", "d")) #r/r_p
    outfile.write(str(" %e" % rho).replace("e", "d")) #rho/rho_p
    outfile.write(str(" %e" % P).replace("e", "d")) #E/E_p
    outfile.write("\n")
  rad = r_zero
  f_rad = r_crit*(1.0/rad-1.0/r_zero)
  rho = f_rad**n_rho
  P = f_rad**n_prs
  outfile.write(str("%e" % rad).replace("e", "d")) #r/r_p
  outfile.write(str(" %e" % rho).replace("e", "d")) #rho/rho_p
  outfile.write(str(" %e" % P).replace("e", "d")) #E/E_p


print "For local problem suggest setting:"
print "rScale = %s" % (str("%e" % (rScale)).replace("e", "d")) 
print "pScale = %s" % (str("%e" % (pScale)).replace("e", "d"))
print "lScale = %s" % (str("%e" % (lScale)).replace("e", "d"))