<jupyter_start><jupyter_text>Setup: ModulesBe sure to run this code block first.  Imports all of the modules necessary for the rest of the notebook.  Note that the ALS module also has scipy as a dependency, even though we do not need to import it within this notebook.This example notebook was developed and tested using the following packages/versions.  Other versions may also work.     - python (3.8.8) - numpy (1.20.1) - pandas (1.2.4) - scipy (1.6.2) - matplotlib (3.3.4) - ipython (7.22.0) - ipympl (0.7.0) - ALS (1.2.0) - StaticCell (1.0.0)<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
import ALS
import StaticCell as SC<jupyter_output><empty_output><jupyter_text>Setup: User Model Step 1: Define the user model. This is the step that requires the most coding on behalf of the user.  I prefer to code the user model in a separate text editor (Sublime Text), save it as a module, and then import it.  The 'ex_model_1.py' file contains the example model_H2O2_depletion with comments providing an explanation how to structure the user model.  It is recommended to follow the template provided there, although any function that has the correct arguments / returns (see below) will work.**Important:** See 'create_model_code.py' for an optional utility function that converts models written in an A+B==>C format (ex: Kintecus) to code representing a system of differential rate equations.  The code generated by this utility function can be used when creating the user model function to reduce typing errors / save time.  See the comments at the top of 'create_model_code.py' for additional details / documentation.*function* **user_model**(*t, model_params*)**Parameters:**>**t : *ndarray***>>Time axis (ms) over which to integrate the model.  You may assume that the points are evenly spaced in ascending order.>**model_params : *dict***>>Keys (strings) are the names of the parameters used by the model; values (floats) are the parameter values.  Any parameters that will be fit or included in a monte carlo simulation of systematic error should be included.>>Only one parameter is required: 'X0' is the key and the initial radical concentration immediately after photolysis is its value.**Returns:** It is only required to return species for which there is observable data to fit; returning other species from the model is optional, but could be useful when plotting model output.  The keys of *m* and *c* must be the same.>**m : *dict***>>Keys (strings) are the names of species returned by the model; values (floats) are their masses (amu).>**c : *dict***>>Keys (strings) are the names of species returned by the model; values (ndarray) are concentrations (molc/cm3) corresponding to times in *t*.<jupyter_code>from ex_model_2 import model_H2O2_depletion<jupyter_output><empty_output><jupyter_text>Step 2: Instantiate the StaticCell object.*class* **StaticCell**(*user_model*)**Parameters:**>**user_model : *function***>>The user model function defined in Step 1.<jupyter_code>model_StaticCell = SC.StaticCell(model_H2O2_depletion)<jupyter_output><empty_output><jupyter_text>Setup: Model Step 1: Calculations for fixed parameters.This is a workspace for any fixed model parameters that require calculation.  Feel free to modify as necessary; there is no recommended format as this will be model-specific.  Variables computed here are used in later steps.**Important:** If converting from a Kintecus model, code for the rate constants can be automatically generated using the 'create_model_code.py' utility function.  See the comments at the top of that file for documentation / an example.  The example corresponds to the running example in this notebook.  Using this utility function may reduce typing errors / save time.<jupyter_code># Fix initial H2O2 concentration and H2O2 scale factor - ignore uncertainties for now
c_H2O2_0 = 5.25e14
c_H2O2_0_err = 0

T = 298 # K
P = 30  # Torr
M = (P*133.3224)/(1.38e-23*T)/1e6 # molc/cm3

# Helper function to compute bimolecular rate constants
def calc_k_bi(A, E_R):
    k_T = A*np.exp(-E_R/T)
    return k_T

# Helper function to compute rate constant uncertainty factor (JPL format)
def calc_f(f_298, g):
    f_T = f_298*np.exp(np.abs(g*(1/T - 1/298.0)))
    return f_T

# Calculate bimolecular rate constants and their uncertainties (cm3/molc/s)

# OH + H2O2 --> HO2 + H2O
# JPL 2015, temperature independent over 200 - 300 K
k1 = 1.8e-12
k1_err = k1*(calc_f(1.15, 45) - 1)

# OH + HO2 --> H2O + O2
# JPL 2015, T = 252-420 K
k2 = calc_k_bi(4.8e-11, -250)
k2_err = k2*(calc_f(1.15, 50) - 1)

# HO2 + HO2 --> H2O2 + O2
# JPL 2015, T = 222-1120 K, M is air
k3a = calc_k_bi(3e-13, -460)
k3a_err = k3a*(calc_f(1.15, 100) - 1)
k3b = calc_k_bi(2.1e-33*M, -920)
k3b_err = k3b*(calc_f(1.2, 200) - 1)
k3 = k3a + k3b
k3_err = (k3a_err**2 + k3b_err**2)**0.5<jupyter_output><empty_output><jupyter_text>Step 2: Organize model parameters into a DataFrame.The *df_model_params* DataFrame will be passed to the StaticCell methods.  Rows correspond to each model parameter; the indices are parameter names and must match each parameter in the user model.  The columns are detailed below.  It is recommended to use the code template in the cell below.**Important:** If converting from a Kintecus model, code for the rate constants can be automatically generated using the 'create_model_code.py' utility function.  See the comments at the top of that file for documentation / an example.  The example corresponds to the running example in this notebook.  Using this utility function may reduce typing errors / save time.**Columns:**>**val : *float***>>Value of the parameter.>**err : *float***>>Absolute uncertainty in the parameter (1 standard error).  Used to vary the value of the parameter during monte carlo simulations of systematic model uncertainty.  Ignored if set to 0 or if *fit* is True. **Currently ignored by StaticCell**>**fit : *bool***>>If True, then this parameter will be optimized during a fit.  If False, then the parameter will be fixed during a fit. **Currently ignored by StaticCell.**<jupyter_code>df_model_params = {}
df_model_params['k_OH_wall'] =  {'val':15,       'err':0,            'fit':True }
df_model_params['k_HO2_wall'] = {'val':4,        'err':0,            'fit':True }
df_model_params['k1'] =         {'val':k1,       'err':k1_err,       'fit':False}
df_model_params['k2'] =         {'val':k2,       'err':k2_err,       'fit':False}
df_model_params['k3'] =         {'val':k3,       'err':k3_err,       'fit':False}
df_model_params = pd.DataFrame.from_dict(df_model_params, orient='index')

print('Inputted Model Params DataFrame:')
display(df_model_params)<jupyter_output>Inputted Model Params DataFrame:<jupyter_text>Step 3: Set Initial ConcentrationsThis is where you set your pre-photolysis concentrations. *initial_concentrations* is a Python dictionary containing the pre-photolysis concentrations for all species in the model (including products).<jupyter_code>initial_concentrations = {
    'H2O2': c_H2O2_0,
    'OH':   0,
    'HO2':  0,
    'H2O':  0,
    'O2' :  0,
}<jupyter_output><empty_output><jupyter_text>Step 4: Photolysis ParametersSet up your photolysis parameters in the next block.**df_photolysis_params : *DataFrame***>Each row is a species to be photolyzed, with the following columns:>**xsn : *float***>>The cross section of the species, in cm^2.>**products : *list of strings***>>List of products formed. Must correspond in order with *qyields*.>**qyields : *list of floats***>>List of quantum yields. Must correspond in order with *products*.<jupyter_code>fluence = 3.2e16# in photons cm^-2

df_photolysis_params = {} ## Cross sections in cm^2, 
df_photolysis_params['H2O2'] = {'xsn': 8.3e-20, 'products': ['OH'], 'qyields': [2.0]}
df_photolysis_params = pd.DataFrame.from_dict(df_photolysis_params, orient='index')

print('Inputted Photolysis Params DataFrame:')
display(df_photolysis_params)<jupyter_output>Inputted Photolysis Params DataFrame:<jupyter_text>Method: Plot the modelIntegrates the user model with the parameters specified in *df_model_params* (no fitting).  The output is plotted (and optionally saved) in concentration units.  All species returned by the user model are plotted, regardless of the *fit* field in *df_data*.See the file 'ex_init_model_2.csv' for an example of how the output is formatted if *save_fn* is specified.*method* **StaticCell.StaticCell.plot_model**(*t_start, t_end, tbin, df_model_params, , delta_xtick=20.0, save_fn=None*)**Parameters:**>**t_prephoto : *float***>>Amount of time to plot pre-photolysis.  Must be an integer multiple of *tbin* \* 0.02 ms and cannot be less than -20 ms.>**t_react : *float***>>How long the reaction goes before the the next photolysis cycle.  Must be an integer multiple of *tbin* \* 0.02 ms and must be greater than *t_start*.>**tbin : *int***>>The time axis step size will be *tbin* \* 0.02 ms.>**df_model_params : *DataFrame***>>Contains the model parameters.  See setup above for formatting.>**initial_concentrations : *dictionary***>>Contains the pre-photolysis concentrations of all species, including products and intermediates.>**df_photolysis_params : *DataFrame***>>Contains the photolysis parameters.  See setup above for formatting.>**fluence : *float***>>Laser fluence in photons cm^-2>**photolysis_cycles : *int***>>Number of photolysis pulses.>**delta_xtick : *float, optional***>>Tick marks and labels for the time axis include zero and are spaced by *delta_xtick* (ms).>**save_fn : *str, optional***>>The points in the plots are saved to *save_fn* if this parameter is specified.  First column is the time axis (ms) and the remaining columns contain the concentrations (molc/cm3) for each species.<jupyter_code>%matplotlib widget
model_StaticCell.plot_model(-20, 60, 10, df_model_params, initial_concentrations, df_photolysis_params, fluence=3.2e16, photolysis_cycles=3, delta_xtick=40.0, save_fn='ex_init_model_2.csv')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10323482428115016, 'fraction_numerical': 0.028454472843450478, 'mean_word_length': 4.951871657754011, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 101, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '445076', 'license': 'Apache-2.0', 'n_tokens_mistral': 3269, 'n_tokens_neox': 3038, 'n_words': 1191, 'path': 'ex_notebook_2.ipynb', 'repo_name': 'bapride11/ALS-Kinetic-Model', 'size': '10016'}