Source code for engine.ComputeZ

from scipy.io.idl import readsav
import os
from bokeh.plotting import figure
from bokeh.io import output_file, show
from bokeh.layouts import column
import numpy as np
from .elements import ELEMENTS as ele
try:
    from sympy.mpmath import *
except:
    from mpmath import *

'''
molDict = {
    "mols": {
        "H2O":1e-3, 
        "H2": 99.0, 
        "CH4": 1e-3, 
        "CO":1e-3, 
        "C2H2":1e-3, 
        "HCN": 1e-3}, 
    "T":1500.0, 
    "g":, 
    "Rp":, 
    "R*":, 
    "P": 1e-4, 
    "P0": 10.0, 
    "taueq": 0.56
}
'''
k = 1.380658e-16 #cgs 
rsun = 6.96e10 #cm 
barye = 1000000.0 #1 bar to 1e6 barye
amu =1.66e-24 #grams
color = ["red", "blue", "green", "purple", "black", "yellow", "orange", "pink","cyan","brown"]
__XSEC__ = os.path.join(os.path.dirname(__file__), "reference/CrossSections/")

[docs]def readXSecs(molDict): MOL = molDict["mols"].keys() files = os.listdir(__XSEC__) output_file("cross.html") TOOLS = "pan,wheel_zoom,box_zoom,resize,reset,save" xsec = {} waves = {} #compute concentration of background gas molDict = addBackground(molDict) #create dictionary of wavelengths and cross sections for nmol in range(0,len(MOL)): if MOL[nmol] == "H2": for i in files: if i.find("H2H2") != -1: readthis = i elif MOL[nmol] == "He": for i in files: if i.find("H2He") != -1: readthis = i else: for i in files: if i.find(MOL[nmol]+'_') != -1: readthis = i cross = readsav(os.path.join(__XSEC__,readthis)) Tind = find_nearest(cross.T, molDict["T"]) Pind = find_nearest(cross.P, molDict["P"]) actualT = cross.T[Tind] actualP = cross.P[Pind] xsec[MOL[nmol]] = cross.xsecarr[Pind, Tind] waves[MOL[nmol]] = cross.wno return {"waves":waves, "xsec":xsec}
[docs]def computeAlpha(molDict, plot = False): plot1 = False plot2 = False MOL = molDict["mols"].keys() read = readXSecs(molDict) waves = read["waves"] xsec = read["xsec"] TOOLS = "pan,wheel_zoom,box_zoom,resize,reset,save" if plot: plot2 = figure(plot_width=800, plot_height=200, tools=TOOLS, responsive =False, y_axis_type="log", x_axis_label='Wavelength [um]', x_axis_type="log", y_axis_label='Weighted Cross Section', y_range = [1e-29, 1e-17]) #compute mean molecular weight and beta mu = computeMu(molDict["mols"])*amu beta = molDict["P0"]*barye / molDict["taueq"] * np.sqrt(2.0*np.pi*molDict["Rp"]*rsun) H = k*molDict["T"]/mu/molDict["g"] sumprod = computeZSum(xsec, molDict["mols"], waves, plot2,plot) Z = H * np.log( 1.0 / np.sqrt(k*molDict["T"]*mu*molDict["g"]) * beta * sumprod) C = 2.0*(molDict["Rp"]*rsun)/(molDict["R*"]*rsun)**2.0 alpha = C*Z if plot: plot1 = figure(plot_width=800, plot_height=200, tools=TOOLS, responsive =False, #y_axis_type="log", x_axis_label='Wavelength [um]', x_axis_type="log", y_axis_label='Alpha Lambda', y_range = [min(alpha[alpha> 0.0]), max(alpha[alpha> 0.0])]) plot1.line(1e4/waves[MOL[0]][alpha> 0.0], alpha[alpha> 0.0], alpha = 0.5, line_width = 3) #"plot1":plot1,"plot2":plot2, return {"w":1e4/waves[MOL[0]][alpha> 0.0], "Z":Z, "alpha":alpha[alpha> 0.0], "H":H, "mu":mu, "xsec":xsec,"mols":molDict["mols"]}
[docs]def find_nearest(array,value): #small program to find the a temp and pressure match in the cross section data idx = (np.abs(array-value)).argmin() return idx
[docs]def computeZSum(xsec, squig, waves, plot2,plot): sumprod = 0.0 num = -1 for i in squig.keys(): if i=="H2": sumprod += squig[i]*xsec[i]*squig["H2"] num +=1 if plot: plot2.line(1e4/waves[i],squig[i]*xsec[i]*squig["H2"], color = color[num], legend = i) elif i=="He": sumprod += squig[i]*xsec[i]*squig["H2"] num +=1 if plot: plot2.line(1e4/waves[i],squig[i]*xsec[i]*squig["H2"], color = color[num], legend = i) else: sumprod += squig[i]*xsec[i] num +=1 if plot: plot2.line(1e4/waves[i],squig[i]*xsec[i], color = color[num], legend = i) return sumprod
[docs]def computeMu(squig): """ Compute mean molecular weight of an atmosphere Input: Dict of molecules with mixing ratios i.e. {"H2O": 0.98, "H2":0.2} Output: mean molecular weight, float """ mu = 0 for i in squig.keys(): totmass = 0.0 for j in range(0, len(i)): try: check = float(i[j]) except: if i[j].isupper(): elem = ele[i[j]] try: num = float(i[j+1]) except: if i[j].islower(): elem = ele[i[j-1:j+1]] num = 1.0 else: num = 1.0 totmass += elem.mass * num mu += squig[i] * totmass return mu
[docs]def addBackground(molDict): #compute concentration of background gas if (molDict["mols"]["H2"] == "bkg") & (molDict["mols"]["He"] == "bkg"): molDict["mols"]["H2"] = 0.0 molDict["mols"]["He"] = 0.0 fH2He = 1.0 - sum(molDict["mols"].values()) try: check1 = np.sqrt(fH2He) check2 = 1.0/fH2He except: raise Exception("Mixing Ratios Exceed total sum = 1, no background gas") frac = molDict["fracH2He"] molDict["mols"]["H2"] = fH2He/(1.+frac) molDict["mols"]["He"] = frac*molDict["mols"]["H2"] elif (molDict["mols"]["H2"] == "bkg") & (molDict["mols"]["He"] != "bkg"): molDict["mols"]["H2"] = 0.0 molDict["mols"]["H2"] = 1.0 - sum(molDict["mols"].values()) elif (molDict["mols"]["H2"] != "bkg") & (molDict["mols"]["He"] == "bkg"): molDict["mols"]["He"] = 0.0 molDict["mols"]["He"] = 1.0 - sum(molDict["mols"].values()) else: for name in molDict["mols"].keys(): if molDict["mols"][name] == "bkg": molDict["mols"][name] = 0.0 molDict["mols"][name] = 1.0 - sum(molDict["mols"].values()) return molDict
####### Should probably eventually separate into two classes starting here. #######
[docs]def computeJac(molDict): """ Function to construct Jacobian of equation alpha http://arxiv.org/pdf/1511.09443v2.pdf eqn. 4 Inputs: T, mu, g, beta: float xsec: dict of cross sections squig: dict of mixing ratios """ squig = molDict["mols"] T = molDict["T"] g = molDict["g"] read = readXSecs(molDict) xsec = read["xsec"] #xsec= { # "H2O":[1.0,2.0], # "H2": [1.0,2.0], # "He": [1.0,2.0], # "CH4": [1.0,2.0], # "CO":[1.0,2.0], # "C2H2":[1.0,2.0], # "HCN": [1.0,2.0]} #compute mean molecular weight and beta mu = computeMu(squig)*amu #mu = 1.0 beta = molDict["P0"]*barye / molDict["taueq"] * np.sqrt(2.0*np.pi*molDict["Rp"]*rsun) flag = [0]*int(len(squig.keys())+3.0) lenW = len(xsec["CO"]) J = np.zeros((lenW,len(flag))) mp.dps = 15; mp.pretty = True for i in range(0,lenW): for j in range(0,len(flag)): flag[j] = 1 aCO = xsec["CO"][i] aH2H2= xsec["H2"][i] aH2O= xsec["H2O"][i] aCH4= xsec["CH4"][i] aHCN= xsec["HCN"][i] aC2H2= xsec["C2H2"][i] aH2He= xsec["He"][i] val = diff(lambda Td,mud,gd, xCO, xH2,xH2O,xCH4,xHCN,xC2H2,xHe: k*Td/mud/gd * log(1.0 / sqrt(k*Td*mud*gd) * beta * (xCO*aCO+xH2*xH2*aH2H2+xH2O*aH2O+xCH4*aCH4+ xHCN*aHCN+xC2H2*aC2H2+xH2*xHe*aH2He)), (T,mu,g,squig["CO"],squig["H2"],squig["H2O"],squig["CH4"],squig["HCN"],squig["C2H2"],squig["He"]), (flag[0],flag[1],flag[2],flag[3],flag[4],flag[5],flag[6],flag[7],flag[8],flag[9])) #since this does not stable enough to handle mu derivative .. calculate analytically instead if j == 1: val = -k*T/mu**2.0/g*np.log(beta/np.sqrt(k*T*mu*g)*(squig["CO"]*aCO+squig["H2"]*squig["H2"] *aH2H2+squig["H2O"]*aH2O+squig["CH4"]*aCH4+ squig["HCN"]*aHCN+squig["C2H2"]*aC2H2+squig["H2"]*squig["He"]*aH2He)) - ( k*T/2.0/mu**2.0/g) flag[j] = 0 #print(i, j, val, type(val)) J[i,j] = float(val) return J