Using the easychem package to approximate the SONORA chemistry grids#
You can find an executable version of this notebook on Google Colab here.
For reference, the SONORA atmosphere and evolution model grids (including chemistry grids) can be found at#
Bobcat: https://zenodo.org/records/5063476 (cloud-free, , Fe/H -0.5 to +0.5; C/O 0.5x to 1.5x solar)
Diamondback: https://zenodo.org/records/12735103 (cloudy, Fe/H -0.5 to +0.5)
Elf Owl: https://zenodo.org/records/10381250 (cloud-free, with diseq chem, Fe/H -1.0 to +1.0; C/O 0.5x to 2.5x solar)
The present notebook approximates equilibrium chemistry of major gas species over a pressure-temperature grid from 75-6000 K and 1E-06 to 1E+04 bar. This 2121-point pressure grid contains all P-T points found in the SONORA models. The first 50 species in this calculation correspond to the species included in the SONORA chemistry outputs. Solar element abundances are taken from Lodders 2020 by default. (NOTE: this notebook uses a truncated species list to approximate the SONORA equilibrium results, and is optimized for chemical equilibrium in a ~solar composition substellar atmosphere. Caution should be taken if extending compositions below -1 dex or beyond +1 dex metallicity. The abundance inputs and species list can be customized as needed.)
We begin by loading any relevant packages: numpy, matplotlib, pandas, and easychem itself.
[ ]:
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
#### Uncomment these lines if you have not yet installed easyCHEM yet
# ! pip install ninja meson-python
# ! pip install easychem --no-build-isolation --no-cache-dir
plt.rcParams['figure.figsize'] = (10, 6)
from numpy import asarray
Chemistry is calculated by creating objects from easychem’s ExoAtmos class.
[19]:
import easychem.easychem as ec
exo = ec.ExoAtmos()
Indicate the path for the thermodynamic data file, which maintains the NASA CEA input data structure. This file may be stored locally or retrieved from a URL.
[20]:
# Indicate path for thermodynamic data file
# Use local file or read in from cvisscher.github.io
# Note: species included in calculation must be present in thermo input file!
# *******************************************************************************************
# set local file
#exo.thermofpath='thermo.inp'
# *******************************************************************************************
# or read file from URL
import requests
# URL for the online thermo input file
url = 'https://raw.githubusercontent.com/cvisscher/cvisscher.github.io/master/thermofile/thermo-sonora-component.inp'
response = requests.get(url)
if response.status_code == 200:
# Write the data to a temp file
with open('temp_file.inp', 'w') as f:
f.write(response.text)
# Set thermo_file_path to this temp file
exo._thermofpath = 'temp_file.inp'
else:
print(f"Failed to retrieve file from {url}. Status code: {response.status_code}")
# *******************************************************************************************
# if desired, show the first lines to confirm:
#f = open(exo._thermofpath, 'r')
#lines = f.readlines()
#f.close()
#for i in range(20):
# print(lines[i][:-1])
# *******************************************************************************************
# Print the whole list of available species in the selected thermodynamic file
exo.print_available_species()
1: H
2: H2
3: H2+
4: H2-
5: H3+
6: He
7: O
8: C
9: N
10: Mg
11: Si
12: Fe
13: S
14: AL
15: Ca
16: Na
17: Ni
18: P
19: K(EC)
20: K
21: Ti
22: CO
23: OH
24: SH
25: N2
26: O2
27: SiO
28: TiO
29: SiS
30: H2O
31: H2O2
32: C2
33: CH
34: CN
35: CS
36: SiC
37: NH
38: N2H2
39: N2H4
40: SiH
41: NO
42: SN
43: SiN
44: S-
45: SH-
46: SO
47: SO2
48: S2
49: C2H
50: HCN
51: C2H2
52: C2H4
53: C2H6
54: C3H8
55: HCHO
56: CH3OH
57: CH2OH
58: CH3O
59: CH4
60: ALH
61: ALO
62: ALOH
63: AL2O
64: ALCL
65: ALCL2
66: ALCL3
67: CaH
68: CaO
69: CaOH
70: MgH
71: MgOH
72: MgO
73: PH3
74: CO2
75: TiO2
76: Si2C
77: SiO2
78: FeO
79: FeS
80: NH2
81: NH3
82: CH2
83: CH3
84: H2S
85: V
86: VO
87: VO2
88: NaCL
89: KCL
90: MgSiO3(c)
91: MgSiO3(L)
92: Mg2SiO4(c)
93: SiC(c)
94: Fe(c)
95: AL2O3(c)
96: Na2S(c)
97: KCL(c)
98: e-
99: H+
100: H-
101: Na+
102: K+
103: K+(EC)
104: KH
105: H2O(c)
106: TiO(c)
107: VO(c)
108: NaH
109: NaOH(CEA)
110: NaOH
111: KOH
112: KOH(EC)
113: Fe(L)
114: Mg2SiO4(L)
115: SiC(L)
116: MgSiO3(L)
117: H2O(L)
118: TiO(L)
119: VO(L)
120: NaAlSi3O8(c)
121: KAlSi3O8(c)
122: MgAl2O4(c)
123: FeO(c)
124: Fe2O3(c)
125: Fe3O4(c)
126: CaMgSi2O6(c)
127: Fe2SiO4(c)
128: P2
129: PS
130: PO
131: PN
132: P4O6(JANAF)
133: P4O6(Gurvich)
134: PH
135: PH2
136: P2H
137: P2H2(cis)
138: P2H2(trans)
139: P4
140: TiO2(c)
141: TiO2(L)
142: H3PO4(L)
143: H3PO4(c)
144: H3PO4
145: HCL
146: SiH4
147: NH4CL(c)
148: SiO2(c)
149: SiO2(L)
150: SiO2(qz)
151: He+
152: Ca+
153: FeH
154: FeOH
155: Fe(OH)2
156: FeCL
157: FeCL2
158: FeCL3
159: Si+
160: Mg+
161: Fe+
162: Ti+
163: V+
164: C-gr
165: COS
166: SiH2
167: SiH3
168: CL
169: CL+
170: CL-
171: CL2
172: H2S2
173: Cr
174: Cr+
175: Cr-
176: CrH
177: CrO
178: CrO2
179: Cr(c)
180: Cr(L)
181: Rb
182: Rb+
183: Rb-
184: RbCL
185: RbH
186: RbO
187: RbOH
188: RbCL(c)
189: atCs
190: Cs+
191: CsH
192: CsCL
193: CsCL(c)
194: Li
195: Li+
196: LiCL
197: LiH
198: LiOH
199: Li2S(c)
200: Diopside(c)
201: Diopside(L)
202: Gehlenite(c)
203: Grossite(c)
204: CaCL
205: CaCL2
206: Ca(OH)2
207: CaS
208: Ca2
209: F
210: F-
211: F2
212: HF
213: LiF
214: LiF(cr)
215: LiCL(cr)
216: LiCL(L)
217: CaF
218: NaF
219: RbF
220: NH4H2PO4(c)
221: SiS2
222: Mg(OH)2
223: MgS
224: C+
225: O+
226: HPO2
227: Mg2
228: MgO(c)
229: MgO(L)
This cell inputs values from published solar abundance databases, including Lodders 2010, Lodders 2020, Asplund+ 2009, and Asplund+ 2020. The last column also allows for the input of custom user-defined elemental abundances.
[21]:
# This cell inputs values from published solar abundance databases (or custom abundances)
# These include (in listed order) Lodders 2010, Lodders 2020, Asplund+ 2009, Asplund+ 2020
abunds_data = ([
#atom,Lo10,Lo20,Asp09,Asp20,custom,
('H',9.106254E-01,9.082387E-01,9.206789E-01,9.232609E-01,9.232609E-01),
('He',8.824980E-02,9.046346E-02,7.836248E-02,7.573985E-02,7.573985E-02),
('Li',1.954856E-09,2.050745E-09,1.033019E-11,8.420239E-12,8.420239E-12),
#('Be',2.151748E-11,2.295826E-11,2.208555E-11,2.214749E-11,2.214749E-11),
#('B',6.609945E-10,6.487419E-10,4.614325E-10,4.627266E-10,4.627266E-10),
('C',2.527952E-04,3.286959E-04,2.478039E-04,2.662713E-04,2.662713E-04),
('N',7.453768E-05,7.893027E-05,6.224553E-05,6.242010E-05,6.242010E-05),
('O',5.520007E-04,5.982842E-04,4.509290E-04,4.521936E-04,4.521936E-04),
('F',2.826806E-08,4.577235E-08,3.342783E-08,2.319126E-08,2.319126E-08),
#('Ne',1.156740E-04,1.575001E-04,7.836248E-05,1.060045E-04,1.060045E-04),
('Na',2.028691E-06,2.083182E-06,1.599957E-06,1.532232E-06,1.532232E-06),
('Mg',3.621406E-05,3.712245E-05,3.665289E-05,3.275853E-05,3.275853E-05),
#('Al',2.974475E-06,2.948892E-06,2.594826E-06,2.484989E-06,2.484989E-06),
('Si',3.515928E-05,3.604122E-05,2.979258E-05,2.987614E-05,2.987614E-05),
('P',2.918220E-07,2.977005E-07,2.366509E-07,2.373146E-07,2.373146E-07),
('S',1.480206E-05,1.575001E-05,1.213691E-05,1.217095E-05,1.217095E-05),
('Cl',1.817735E-07,1.906580E-07,2.911442E-07,1.885057E-07,1.885057E-07),
#('Ar',3.259265E-06,3.521227E-06,2.312641E-06,2.214749E-06,2.214749E-06),
('K',1.321989E-07,1.301448E-07,9.865252E-08,1.084737E-07,1.084737E-07),
#('Ca',2.123621E-06,2.062783E-06,2.014226E-06,1.842148E-06,1.842148E-06),
#('Sc',1.209479E-09,1.214589E-09,1.300493E-09,1.274455E-09,1.274455E-09),
('Ti',8.684343E-08,8.862536E-08,8.205559E-08,8.616372E-08,8.616372E-08),
('V',1.005555E-08,9.911335E-09,7.836248E-09,7.333722E-09,7.333722E-09),
('Cr',4.605866E-07,4.732212E-07,4.018909E-07,3.848792E-07,3.848792E-07),
#('Mn',3.241686E-07,3.276147E-07,2.478039E-07,2.428424E-07,2.428424E-07),
('Fe',2.981507E-05,3.142794E-05,2.911442E-05,2.662713E-05,2.662713E-05),
#('Co',8.262431E-08,8.145315E-08,8.997217E-08,8.041266E-08,8.041266E-08),
#('Ni',1.722805E-06,1.754126E-06,1.527947E-06,1.463270E-06,1.463270E-06),
#('Cu',1.902117E-08,1.928205E-08,1.425963E-08,1.397412E-08,1.397412E-08),
#('Zn',4.570707E-08,4.541193E-08,3.342783E-08,3.352158E-08,3.352158E-08),
#('Ga',1.286830E-09,1.304692E-09,1.009504E-09,9.667728E-10,9.667728E-10),
#('Ge',4.043317E-09,4.324946E-09,4.112522E-09,3.848792E-09,3.848792E-09),
#('As',2.144716E-10,2.187702E-10,1.836996E-10,1.842148E-10,1.842148E-10),
#('Se',2.373252E-09,2.436386E-09,2.014226E-09,2.019875E-09,2.019875E-09),
#('Br',3.762043E-10,4.433070E-10,9.206789E-53,9.232609E-53,9.232609E-53),
#('Kr',1.961888E-09,1.848914E-09,9.206789E-53,9.232609E-53,9.232609E-53),
('Rb',2.542016E-10,2.584155E-10,3.048654E-10,1.928965E-10,1.928965E-10),
#('Sr',8.192113E-10,8.397604E-10,6.825087E-10,6.242010E-10,6.242010E-10),
#('Y',1.627875E-10,1.567793E-10,1.493166E-10,1.497354E-10,1.497354E-10),
#('Zr',3.797202E-10,3.928493E-10,3.500323E-10,3.591902E-10,3.591902E-10),
#('Nb',2.742424E-11,2.811215E-11,2.655267E-11,2.724736E-11,2.724736E-11),
#('Mo',8.965617E-11,9.370717E-11,6.984064E-11,7.003650E-11,7.003650E-11),
#('Ru',6.258352E-11,6.523460E-11,5.177358E-11,5.191877E-11,5.191877E-11),
#('Rh',1.300893E-11,1.218193E-11,7.483559E-12,5.563197E-12,5.563197E-12),
#('Pd',4.781662E-11,4.973688E-11,3.420646E-11,3.430239E-11,3.430239E-11),
#('Ag',1.719289E-11,1.791249E-11,8.018778E-12,8.420239E-12,8.420239E-12),
#('Cd',5.520007E-11,5.694512E-11,4.721806E-11,4.735048E-11,4.735048E-11),
#('In',6.258352E-12,6.451378E-12,5.809091E-12,5.825382E-12,5.825382E-12),
#('Sn',1.265734E-10,1.293880E-10,1.009504E-10,9.667728E-11,9.667728E-11),
#('Sb',1.100486E-11,1.293880E-11,9.421242E-12,9.447664E-12,9.447664E-12),
#('Te',1.648970E-10,1.701145E-10,1.393504E-10,1.397412E-10,1.397412E-10),
#('I',3.867521E-11,5.730554E-11,9.206789E-53,9.232609E-53,9.232609E-53),
#('Xe',1.919697E-10,1.982267E-10,9.206789E-53,9.232609E-53,9.232609E-53),
('Cs',1.304409E-11,1.326317E-11,1.106899E-11,1.110004E-11,1.110004E-11),
#('Ba',1.571620E-10,1.639875E-10,1.393504E-10,1.719192E-10,1.719192E-10),
#('La',1.606779E-11,1.654292E-11,1.159066E-11,1.189390E-11,1.189390E-11),
#('Ce',4.148795E-11,4.180781E-11,3.500323E-11,3.510140E-11,3.510140E-11),
#('Pr',6.047397E-12,6.307213E-12,4.831791E-12,5.191877E-12,5.191877E-12),
#('Nd',3.009635E-11,3.113961E-11,2.421632E-11,2.428424E-11,2.428424E-11),
#('Sm',9.387528E-12,9.839253E-12,8.396691E-12,8.228571E-12,8.228571E-12),
#('Eu',3.515928E-12,3.604122E-12,3.048654E-12,3.057204E-12,3.057204E-12),
#('Gd',1.265734E-11,1.247026E-11,1.081703E-11,1.110004E-11,1.110004E-11),
#('Tb',2.109557E-12,2.252576E-12,1.836996E-12,1.885057E-12,1.885057E-12),
#('Dy',1.420435E-11,1.466878E-11,1.159066E-11,1.162317E-11,1.162317E-11),
#('Ho',3.164335E-12,3.211273E-12,2.780406E-12,2.788203E-12,2.788203E-12),
#('Er',9.211732E-12,9.226552E-12,7.657873E-12,7.858224E-12,7.858224E-12),
#('Tm',1.406371E-12,1.452461E-12,1.159066E-12,1.189390E-12,1.189390E-12),
#('Yb',9.000776E-12,9.082387E-12,6.369542E-12,6.536186E-12,6.536186E-12),
#('Lu',1.336053E-12,1.373170E-12,1.159066E-12,1.162317E-12,1.162317E-12),
#('Hf',5.484848E-12,5.586389E-12,6.517907E-12,6.536186E-12,6.536186E-12),
#('Ta',7.383449E-13,7.748862E-13,6.984064E-13,6.536186E-13,6.536186E-13),
#('W',4.816822E-12,5.189935E-12,6.517907E-12,5.692780E-12,5.692780E-12),
#('Re',2.042754E-12,1.971455E-12,1.675360E-12,1.680059E-12,1.680059E-12),
#('Os',2.383799E-11,2.349887E-11,2.312641E-11,2.066924E-11,2.066924E-11),
#('Ir',2.362704E-11,2.281409E-11,2.208555E-11,1.928965E-11,1.928965E-11),
#('Pt',4.465229E-11,4.469111E-11,3.838028E-11,3.761183E-11,3.761183E-11),
#('Au',6.856060E-12,7.028038E-12,7.657873E-12,7.504546E-12,7.504546E-12),
#('Hg',1.610295E-11,1.355150E-11,1.361784E-11,1.365603E-11,1.365603E-11),
#('Tl',6.398989E-12,6.451378E-12,7.313212E-12,7.679349E-12,7.679349E-12),
#('Pb',1.163772E-10,1.192964E-10,5.177358E-11,8.228571E-11,8.228571E-11),
#('Bi',4.851981E-12,5.081812E-12,4.112522E-12,4.124055E-12,4.124055E-12),
#('Th',1.547008E-12,1.517335E-12,9.640691E-13,9.892919E-13,9.892919E-13),
#('U',6.680264E-14,8.610247E-13,2.655267E-13,2.662713E-13,2.662713E-13)
])
def process_abunds_data(data):
return [list(col) if i == 0 else list(map(float, col))
for i, col in enumerate(zip(*data))]
selected_atoms, Lo10_abunds, Lo20_abunds, Asp09_abunds, Asp20_abunds, custom_abunds = process_abunds_data(abunds_data)
exo.atoms = selected_atoms
natoms = len(exo.atoms)
#print(exo.atoms)
#print(natoms,'atoms selected')
Select the element abundance database from above, and indicate the metallicity and C/O ratio to be applied for the equilibrium calculations.
[41]:
# Select abundances, metallicity, and CO ratio
#******************************************************************
# Lo10_abunds, Lo20_abunds, Asp09_abunds, Asp20_abunds, custom_abunds
#******************************************************************
exo.atomAbunds = np.array(Lo20_abunds)
#******************************************************************
# Indicate metallicity and C/O ratio relative to the "solar" C/O ratio
#******************************************************************
# NOTE: the (truncated) reaction list is optimized for ~solar abundances
# great care should be taking extending calculations Fe/H < -1 or > +1 dex
feh = 0.00 # metallicity value in dex (i.e. +1.0 = 10x solar)
co_factor = 1.0 # C/O ratio **relative to solar** (i.e. 1 = solar ratio)
# note - the filname output will include the *actual* C/O ratio of the run
# for example co_factor = 1.0 will give co0.55 in filename for Lo20 abunds
#******************************************************************
natoms = len(exo.atomAbunds)
# Metallicity adjustment for elements heavier than helium
for i in range(2, natoms):
exo.atomAbunds[i] *= 10.**feh
# Calculate the new C/O ratio
# Note: this keeps C + O = constant to keep total metallicity constant
co_solar = exo.atomAbunds[3]/exo.atomAbunds[5]
co_sum = exo.atomAbunds[3] + exo.atomAbunds[5]
exo.atomAbunds[5] = co_sum/(1+co_factor*co_solar)
exo.atomAbunds[3] = co_solar*co_factor*exo.atomAbunds[5]
# set these values as output_tags for file and plot labels
output_tag_01 = feh
output_tag_02 = co_factor*co_solar
output_tag_02 = round(output_tag_02, 2)
print(natoms,'atoms selected:')
for i in range(len(exo.atoms)):
print(f"{exo.atoms[i]:3}: {exo.atomAbunds[i]:6e}")
20 atoms selected:
H : 9.083916e-01
He : 9.047868e-02
Li : 2.051090e-09
C : 3.287512e-04
N : 7.894355e-05
O : 5.983849e-04
F : 4.578005e-08
Na : 2.083533e-06
Mg : 3.712870e-05
Si : 3.604729e-05
P : 2.977506e-07
S : 1.575266e-05
Cl : 1.906901e-07
K : 1.301667e-07
Ti : 8.864028e-08
V : 9.913003e-09
Cr : 4.733008e-07
Fe : 3.143323e-05
Rb : 2.584590e-10
Cs : 1.326540e-11
Select reactants to be included in calculation
[23]:
# @title
# Select reactants to be included in calculation
exo.reactants=np.array([
# Initial 31 species (up to COS) from *reported output* of original Sonora grids (Summer 2015).
# Li-bearing species and OH, C-gr added with updates following Gharib-Nezhad et al (2021).
# Additional atomics and ions added to output list (summer 2024).
# First 50 species (up to O+) are currently included in output of Sonora chemistry grids.
#
#
# NOTE: For stability and flexibility, this version uses component metal oxides as stoichiometric proxies
# for oxide and silicate condensates in substellar atmopsheres. More precise condensation curves
# may be found by replacing these oxides with the expected species in the condensate sequence.
'e-','H2', 'H', 'H+', 'H-','H2-', 'H2+','H3+',
'He',
'H2O',
'CH4','CO',
'NH3','N2',
'PH3','H2S',
'TiO','VO','Fe','FeH','CrH',
'Na','K','Rb','atCs',
'CO2','HCN','C2H2','C2H4','C2H6','COS',
'SiO','MgH',
'Li','LiOH','LiH','LiCl','Li+','LiF',
'OH','C-gr',
'Mg','Mg+','Si','Fe+','Ti','Ti+','C','O','C+','O+',
# Additional species included in calculation
# MWE list that still approximates Sonora chemistry grids
'He+',
'C2','CH','CN',
'CS','C2H','CH2','CH3','C3H8','HCHO','CH2OH','CH3OH','CH3O',
'N','NH','NH2','NO','N2H2','N2H4',
'O2','H2O2',
'P','PH2', 'P2','PO','PH','P4O6(Gurvich)',
#'P4O6(Gurvich)','P4O6(JANAF)','HPO2','H3PO4','PN','PS',
'S','SH','SN',
'SO', 'S2','SO2','S-','SH-',
'Cr','Cr+','CrO','CrO2',
'FeO','FeOH','FeS','Fe(OH)2','FeCL',
'MgO','MgOH','MgS','Mg(OH)2',
'Si+','SiS','SiH','SiO2','SiH2','SiH3','SiH4',
'Na+', 'NaCL','NaOH','NaH',
'K+','KCl','KH','KOH',
'V', 'V+','VO2',
'TiO2',
'CL-','CL','HCL',#'CL2',
'RbCL','Rb+','RbH','RbO','RbOH','RbF',
'CsCL','Cs+','CsH',
'F','F-','HF','NaF',#'F2',
# Al- and Ca-bearing species (optional)
#'Al','AlH','AlO','AlOH','Al2O','ALCL','ALCL2','ALCL3',
#'Ca','Ca+','CaH','CaO','CaOH',
# Condensates included in the calculation
'NH4H2PO4(c)',
'VO(c)','VO(L)', #proxy
'TiO2(c)','TiO2(L)', #proxy
'MgO(c)','MgO(L)', #proxy
'SiO2(c)','SiO2(L)', #proxy
#'MgSiO3(c)','Mg2SiO4(c)',
'Cr(c)','Cr(L)',
'Fe(c)','Fe(L)',
'H2O(L)','H2O(c)',
'Na2S(c)',
'KCL(c)',
'RbCL(c)',
'CsCL(c)',
'Li2S(c)',
'LiF(cr)'
#'Gehlenite(c)','Grossite(c)',
# 'Al2O3(c)','CaO(c)'
# Additional species as needed:
])
# *****************************************************************************
# list included species by array index (useful for selecting species for plots)
names = asarray(exo.reactants)
number_of_species=len(names)
print(number_of_species,"species included in this calculation")
for i, name in enumerate(names):
print(f"{i}: {name}")
155 species included in this calculation
0: e-
1: H2
2: H
3: H+
4: H-
5: H2-
6: H2+
7: H3+
8: He
9: H2O
10: CH4
11: CO
12: NH3
13: N2
14: PH3
15: H2S
16: TiO
17: VO
18: Fe
19: FeH
20: CrH
21: Na
22: K
23: Rb
24: atCs
25: CO2
26: HCN
27: C2H2
28: C2H4
29: C2H6
30: COS
31: SiO
32: MgH
33: Li
34: LiOH
35: LiH
36: LiCl
37: Li+
38: LiF
39: OH
40: C-gr
41: Mg
42: Mg+
43: Si
44: Fe+
45: Ti
46: Ti+
47: C
48: O
49: C+
50: O+
51: He+
52: C2
53: CH
54: CN
55: CS
56: C2H
57: CH2
58: CH3
59: C3H8
60: HCHO
61: CH2OH
62: CH3OH
63: CH3O
64: N
65: NH
66: NH2
67: NO
68: N2H2
69: N2H4
70: O2
71: H2O2
72: P
73: PH2
74: P2
75: PO
76: PH
77: P4O6(Gurvich)
78: S
79: SH
80: SN
81: SO
82: S2
83: SO2
84: S-
85: SH-
86: Cr
87: Cr+
88: CrO
89: CrO2
90: FeO
91: FeOH
92: FeS
93: Fe(OH)2
94: FeCL
95: MgO
96: MgOH
97: MgS
98: Mg(OH)2
99: Si+
100: SiS
101: SiH
102: SiO2
103: SiH2
104: SiH3
105: SiH4
106: Na+
107: NaCL
108: NaOH
109: NaH
110: K+
111: KCl
112: KH
113: KOH
114: V
115: V+
116: VO2
117: TiO2
118: CL-
119: CL
120: HCL
121: RbCL
122: Rb+
123: RbH
124: RbO
125: RbOH
126: RbF
127: CsCL
128: Cs+
129: CsH
130: F
131: F-
132: HF
133: NaF
134: NH4H2PO4(c)
135: VO(c)
136: VO(L)
137: TiO2(c)
138: TiO2(L)
139: MgO(c)
140: MgO(L)
141: SiO2(c)
142: SiO2(L)
143: Cr(c)
144: Cr(L)
145: Fe(c)
146: Fe(L)
147: H2O(L)
148: H2O(c)
149: Na2S(c)
150: KCL(c)
151: RbCL(c)
152: CsCL(c)
153: Li2S(c)
154: LiF(cr)
Read-in or define the pressure and temperature points for the calculation
[24]:
# @title
# Read-in or define P-T points for calculation
#******************************************************************
# Indicate pressures and temperatures to be used in the calculation
# By default this uses a 2121-point grid containing the SONORA grid
#******************************************************************
p_2121 = [1.000E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04,1.00E-06,3.00E-06,1.00E-05,3.00E-05,1.00E-04,3.00E-04,1.00E-03,3.00E-03,1.00E-02,3.00E-02,1.00E-01,3.00E-01,1.00E+00,3.00E+00,1.00E+01,3.00E+01,1.00E+02,3.00E+02,1.00E+03,3.00E+03,1.00E+04]
t_2121 = [75,75,75,75,75,75,75,75,75,75,75,75,75,75,75,75,75,75,75,75,75,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,110,110,110,110,110,110,110,110,110,110,110,110,110,110,110,110,110,110,110,110,110,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,120,130,130,130,130,130,130,130,130,130,130,130,130,130,130,130,130,130,130,130,130,130,140,140,140,140,140,140,140,140,140,140,140,140,140,140,140,140,140,140,140,140,140,150,150,150,150,150,150,150,150,150,150,150,150,150,150,150,150,150,150,150,150,150,160,160,160,160,160,160,160,160,160,160,160,160,160,160,160,160,160,160,160,160,160,170,170,170,170,170,170,170,170,170,170,170,170,170,170,170,170,170,170,170,170,170,180,180,180,180,180,180,180,180,180,180,180,180,180,180,180,180,180,180,180,180,180,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,190,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200,200,210,210,210,210,210,210,210,210,210,210,210,210,210,210,210,210,210,210,210,210,210,220,220,220,220,220,220,220,220,220,220,220,220,220,220,220,220,220,220,220,220,220,230,230,230,230,230,230,230,230,230,230,230,230,230,230,230,230,230,230,230,230,230,240,240,240,240,240,240,240,240,240,240,240,240,240,240,240,240,240,240,240,240,240,250,250,250,250,250,250,250,250,250,250,250,250,250,250,250,250,250,250,250,250,250,260,260,260,260,260,260,260,260,260,260,260,260,260,260,260,260,260,260,260,260,260,270,270,270,270,270,270,270,270,270,270,270,270,270,270,270,270,270,270,270,270,270,275,275,275,275,275,275,275,275,275,275,275,275,275,275,275,275,275,275,275,275,275,280,280,280,280,280,280,280,280,280,280,280,280,280,280,280,280,280,280,280,280,280,290,290,290,290,290,290,290,290,290,290,290,290,290,290,290,290,290,290,290,290,290,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,300,310,310,310,310,310,310,310,310,310,310,310,310,310,310,310,310,310,310,310,310,310,320,320,320,320,320,320,320,320,320,320,320,320,320,320,320,320,320,320,320,320,320,330,330,330,330,330,330,330,330,330,330,330,330,330,330,330,330,330,330,330,330,330,340,340,340,340,340,340,340,340,340,340,340,340,340,340,340,340,340,340,340,340,340,350,350,350,350,350,350,350,350,350,350,350,350,350,350,350,350,350,350,350,350,350,375,375,375,375,375,375,375,375,375,375,375,375,375,375,375,375,375,375,375,375,375,400,400,400,400,400,400,400,400,400,400,400,400,400,400,400,400,400,400,400,400,400,425,425,425,425,425,425,425,425,425,425,425,425,425,425,425,425,425,425,425,425,425,450,450,450,450,450,450,450,450,450,450,450,450,450,450,450,450,450,450,450,450,450,475,475,475,475,475,475,475,475,475,475,475,475,475,475,475,475,475,475,475,475,475,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,500,525,525,525,525,525,525,525,525,525,525,525,525,525,525,525,525,525,525,525,525,525,550,550,550,550,550,550,550,550,550,550,550,550,550,550,550,550,550,550,550,550,550,575,575,575,575,575,575,575,575,575,575,575,575,575,575,575,575,575,575,575,575,575,600,600,600,600,600,600,600,600,600,600,600,600,600,600,600,600,600,600,600,600,600,650,650,650,650,650,650,650,650,650,650,650,650,650,650,650,650,650,650,650,650,650,700,700,700,700,700,700,700,700,700,700,700,700,700,700,700,700,700,700,700,700,700,725,725,725,725,725,725,725,725,725,725,725,725,725,725,725,725,725,725,725,725,725,750,750,750,750,750,750,750,750,750,750,750,750,750,750,750,750,750,750,750,750,750,800,800,800,800,800,800,800,800,800,800,800,800,800,800,800,800,800,800,800,800,800,850,850,850,850,850,850,850,850,850,850,850,850,850,850,850,850,850,850,850,850,850,900,900,900,900,900,900,900,900,900,900,900,900,900,900,900,900,900,900,900,900,900,950,950,950,950,950,950,950,950,950,950,950,950,950,950,950,950,950,950,950,950,950,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1100,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1200,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1300,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1400,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1500,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1600,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1700,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1800,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,1900,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2000,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2100,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2150,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2200,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2400,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2450,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2500,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2600,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2700,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2800,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,2900,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3000,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3100,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3200,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3250,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3300,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3400,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3600,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3700,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3750,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3800,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,3900,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4000,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4100,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4200,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4300,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4400,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4500,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4600,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4700,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4800,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,4900,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5000,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5100,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5200,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5300,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5400,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5500,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5600,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5700,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5800,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,5900,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000,6000]
p_grid = np.array(p_2121,dtype=np.float64)
t_grid = np.array(t_2121,dtype=np.float64)
# Read-in grid from csv file
#grid=pd.read_csv('pt_2121.csv',dtype=np.float64)
#t_grid = grid['tk'].to_numpy()
#p_grid = grid['pt'].to_numpy()
print(p_grid,t_grid)
[1.e-06 3.e-06 1.e-05 ... 1.e+03 3.e+03 1.e+04] [ 75. 75. 75. ... 6000. 6000. 6000.]
Run the equilibrium model for the selected parameters
[25]:
#******************************************************************
# Run the equilibrium model for the selected parameters
#******************************************************************
%time exo.solve(p_grid, t_grid)
CPU times: user 50.5 s, sys: 1.23 s, total: 51.8 s
Wall time: 57.2 s
Store calculation results in arrays for plotting and printing
[37]:
# Store calculation results in arrays for plotting and printing
names = asarray(exo.reactants) # 1D array of names
data = asarray(exo.reacMols) # 2D array of rows, abundances
Plot the chemistry results in mole fraction contour plots
[44]:
#******************************************************************
# PLOT CHEMISTRY RESULTS IN MOLE FRACTION CONTOUR PLOTS
#******************************************************************
#******************************************************************
# Select the species to include in the output file "values to plot"
#******************************************************************
#specno_values_to_plot = (range(0, 50)) # write first 50 species match SONORA outputs
#specno_values_to_plot = (range(0, number_of_species)) # write ALL species in model
specno_values_to_plot = [10,11] #select to choose specific species from array
#******************************************************************
# Indicate a location to save any figure or file outputs as needed
#output_path = '/content/drive/MyDrive/Colab Notebooks/model-output/'
#******************************************************************
p_log = np.log10(p_grid)
x = np.reshape(t_grid, (101,21))
y = np.reshape(p_log, (101,21))
tmin = 300.
tmax = 4000.
pmin = -3.5
pmax = 3.5
# Iterate over the specified specno values
for specno in specno_values_to_plot:
print(names[specno])
z = data[:,specno]
# Replace zero values with a small value (e.g., 1e-50) to avoid issues with logarithmic scaling
z = np.where(z == 0, 1e-50, z)
z_log = np.log10(z)
z = np.reshape(z_log, (101,21))
contour_levels_major = [-55,-50,-45,-40,-35,-30,-25,-20,-15,-10,-5,0]
contour_levels_minor = [-31,-29,-28,-27,-26,-24,-23,-22,-21,-19,-18,-17,-16,-14,-13,-12,-11,-9,-8,-7,-6,-4,-3,-2,-1]
plt.figure(figsize=(4, 4))
plt.xlabel('Temperature (K)')
plt.ylabel('log(10) Pressure (bar)')
plt.title(names[specno], fontsize=12)
plt.text( 0.98, 0.98,
f"L20 abunds; [Fe/H]={output_tag_01}; [C/O]={output_tag_02}",
fontsize=8, ha='right', va='top',
transform=plt.gca().transAxes,
)
plt.gca().invert_yaxis()
plt.xlim(tmin, tmax)
plt.ylim(pmax, pmin)
contourf = plt.contourf(x, y, z, 45, cmap='viridis', vmin=-40, vmax=0)
contour_plot = plt.contour(x, y, z, levels=contour_levels_major, cmap='gray_r',linewidths=1)
plt.clabel(contour_plot, inline=False, fontsize=8, fmt='%1.1f')
contour_plot = plt.contour(x, y, z, levels=contour_levels_minor, cmap='gray_r', linestyles='dashed',linewidths=1)
# Create a secondary y-axis 10^x scaling
ax2 = plt.gca().twinx()
ax2.set_yscale('log')
ax2.set_ylabel('Pressure (bar)')
y2_min, y2_max = 10**pmin, 10**pmax
ax2.set_ylim(y2_max, y2_min)
# Save each plot (in the specified subdirectory if desired) with a unique filename
filename = f'{specno}_{names[specno]}_plot.pdf'
plt.savefig(f"{filename}")
# Show the plot
plt.show()
plt.close()
CH4
CO
Write grid chemistry model results to output text file
[28]:
# Write grid chemistry model results to output text file
# This writes selected species (first 50) for the SONORA output
import numpy as np
# Load in the data arrays for the output files
p_log = np.log10(p_grid)
x_array = np.asarray(t_grid)
y_array = np.asarray(p_log)
#******************************************************************
# Select the species to include in the output file
#******************************************************************
specno_values_to_plot = (range(0, 50)) # the first 50 species match the SONORA outputs
#specno_values_to_plot = (range(0, number_of_species)) # write all species
#specno_values_to_plot = (0) #select to choose specific species from array
#******************************************************************
# Initialize an empty array to store z_array values
z_array = np.zeros((len(specno_values_to_plot), len(x_array)))
for i, specno in enumerate(specno_values_to_plot):
# print(names[specno])
z = data[:, specno]
# Replace any zero values with very small value (e.g., 1e-50) to avoid issues with log scaling
z = np.where(z == 0, 1e-50, z)
z_array[i, :] = z
# This will place the adopted feh and c/o ratio values into the filename
output_file_name = 'sonora_output_feh{}_co{}.txt'.format(output_tag_01, output_tag_02) # Replace with tags indicating feh and c/o ratio relative to solar
with open(output_file_name, 'w') as file:
# Write and format data output
file.write('{:<6}{:<8}'.format('T(K)', 'P(bar)'))
# Append additional headers for z_array
for specno in specno_values_to_plot:
file.write('{:<13}'.format(' ' + names[specno]))
file.write('\n')
# Write the data for the remaining rows
for i in range(len(x_array)):
file.write('{:6.1f}{:8.3f}'.format(x_array[i], y_array[i]))
for j in range(len(specno_values_to_plot)):
file.write('{:13.4e}'.format(z_array[j, i]))
file.write('\n')