Chemical equilibrium composition of WASP-43b atmosphere#

 1% -------------------------------------------------------------------------
 2% EXAMPLE: EXOPLANET WASP-43b - METALLICITY 1
 3%
 4% Compute equilibrium vertical composition with a metallicity 1 of WASP-43b
 5%   
 6% URL RESULTS TEA:
 7% https://github.com/dzesmin/RRC-BlecicEtal-2015a-ApJS-TEA/tree/master/Fig6/WASP43b-solar
 8%
 9% @author: Alberto Cuadra Lara
10%                 
11% Last update April 02 2024
12% -------------------------------------------------------------------------
13
14% Import packages
15import combustiontoolbox.databases.SolarAbundances
16import combustiontoolbox.databases.NasaDatabase
17import combustiontoolbox.core.*
18import combustiontoolbox.equilibrium.*
19import combustiontoolbox.utils.display.*
20
21% Definitions
22listSpecies = {'C2H2_acetylene', 'C2H4', 'C', 'CH4', 'CO2', 'CO', 'H2', 'H2O', 'H2S', 'H', 'HCN', 'He', 'HS_M', 'N2', 'N', 'NH3', 'O', 'S'};
23species = {'H', 'He', 'C', 'N', 'O', 'S'};
24metallicity = 1;
25
26% Get initial composition from solar abundances
27DB_solar = SolarAbundances();
28moles = DB_solar.abundances2moles(species, metallicity);
29
30% Get Nasa database
31DB = NasaDatabase();
32
33% Define chemical system
34system = ChemicalSystem(DB, listSpecies);
35
36% Initialize mixture
37mix = Mixture(system);
38
39% Define chemical state
40set(mix, species, moles);
41
42% Define properties
43mixArray = setProperties(mix, 'temperature', linspace(100, 4000, 300), 'pressure', logspace(-5, 2, 300));
44
45% Initialize solver
46solver = EquilibriumSolver('problemType', 'TP');
47
48% Solve problem
49solver.solveArray(mixArray);
50
51% Plot molar fractions
52plotComposition(mixArray(1), mixArray, 'Xi', 'p', 'mintol', 1e-14, 'ydir', 'reverse', 'xscale', 'log');