ex_notebook_2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup: Modules\n",
    "\n",
    "Be 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.\n",
    "\n",
    "This example notebook was developed and tested using the following packages/versions.  Other versions may also work.    \n",
    " - python (3.8.8)\n",
    " - numpy (1.20.1)\n",
    " - pandas (1.2.4)\n",
    " - scipy (1.6.2)\n",
    " - matplotlib (3.3.4)\n",
    " - ipython (7.22.0)\n",
    " - ipympl (0.7.0)\n",
    " - ALS (1.2.0)\n",
    " - StaticCell (1.0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import display\n",
    "import ALS\n",
    "import StaticCell as SC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup: User Model\n",
    "\n",
    "#### Step 1: Define the user model. \n",
    "\n",
    "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.\n",
    "\n",
    "**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.\n",
    "\n",
    "*function* **user_model**(*t, model_params*)\n",
    "\n",
    "**Parameters:**\n",
    "\n",
    ">**t : *ndarray***\n",
    ">>Time axis (ms) over which to integrate the model.  You may assume that the points are evenly spaced in ascending order.\n",
    "\n",
    ">**model_params : *dict***\n",
    ">>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.\n",
    "\n",
    ">>Only one parameter is required: 'X0' is the key and the initial radical concentration immediately after photolysis is its value.\n",
    "\n",
    "**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.\n",
    "\n",
    ">**m : *dict***\n",
    ">>Keys (strings) are the names of species returned by the model; values (floats) are their masses (amu).\n",
    "\n",
    ">**c : *dict***\n",
    ">>Keys (strings) are the names of species returned by the model; values (ndarray) are concentrations (molc/cm<sup>3</sup>) corresponding to times in *t*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ex_model_2 import model_H2O2_depletion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Instantiate the StaticCell object.\n",
    "\n",
    "*class* **StaticCell**(*user_model*)\n",
    "\n",
    "**Parameters:**\n",
    "\n",
    ">**user_model : *function***\n",
    ">>The user model function defined in Step 1.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_StaticCell = SC.StaticCell(model_H2O2_depletion)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Setup: Model\n",
    "\n",
    "\n",
    "#### Step 1: Calculations for fixed parameters.\n",
    "\n",
    "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.\n",
    "\n",
    "**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."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fix initial H2O2 concentration and H2O2 scale factor - ignore uncertainties for now\n",
    "c_H2O2_0 = 5.25e14\n",
    "c_H2O2_0_err = 0\n",
    "\n",
    "T = 298 # K\n",
    "P = 30  # Torr\n",
    "M = (P*133.3224)/(1.38e-23*T)/1e6 # molc/cm3\n",
    "\n",
    "# Helper function to compute bimolecular rate constants\n",
    "def calc_k_bi(A, E_R):\n",
    "    k_T = A*np.exp(-E_R/T)\n",
    "    return k_T\n",
    "\n",
    "# Helper function to compute rate constant uncertainty factor (JPL format)\n",
    "def calc_f(f_298, g):\n",
    "    f_T = f_298*np.exp(np.abs(g*(1/T - 1/298.0)))\n",
    "    return f_T\n",
    "\n",
    "# Calculate bimolecular rate constants and their uncertainties (cm3/molc/s)\n",
    "\n",
    "# OH + H2O2 --> HO2 + H2O\n",
    "# JPL 2015, temperature independent over 200 - 300 K\n",
    "k1 = 1.8e-12\n",
    "k1_err = k1*(calc_f(1.15, 45) - 1)\n",
    "\n",
    "# OH + HO2 --> H2O + O2\n",
    "# JPL 2015, T = 252-420 K\n",
    "k2 = calc_k_bi(4.8e-11, -250)\n",
    "k2_err = k2*(calc_f(1.15, 50) - 1)\n",
    "\n",
    "# HO2 + HO2 --> H2O2 + O2\n",
    "# JPL 2015, T = 222-1120 K, M is air\n",
    "k3a = calc_k_bi(3e-13, -460)\n",
    "k3a_err = k3a*(calc_f(1.15, 100) - 1)\n",
    "k3b = calc_k_bi(2.1e-33*M, -920)\n",
    "k3b_err = k3b*(calc_f(1.2, 200) - 1)\n",
    "k3 = k3a + k3b\n",
    "k3_err = (k3a_err**2 + k3b_err**2)**0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Organize model parameters into a DataFrame.\n",
    "\n",
    "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.\n",
    "\n",
    "**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.\n",
    "\n",
    "**Columns:**\n",
    "\n",
    ">**val : *float***\n",
    ">>Value of the parameter.\n",
    "\n",
    ">**err : *float***\n",
    ">>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**\n",
    "\n",
    ">**fit : *bool***\n",
    ">>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.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inputted Model Params DataFrame:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>val</th>\n",
       "      <th>err</th>\n",
       "      <th>fit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>k_OH_wall</th>\n",
       "      <td>1.500000e+01</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k_HO2_wall</th>\n",
       "      <td>4.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k1</th>\n",
       "      <td>1.800000e-12</td>\n",
       "      <td>2.700000e-13</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k2</th>\n",
       "      <td>1.110663e-10</td>\n",
       "      <td>1.665994e-11</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>k3</th>\n",
       "      <td>1.449221e-12</td>\n",
       "      <td>2.108588e-13</td>\n",
       "      <td>False</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     val           err    fit\n",
       "k_OH_wall   1.500000e+01  0.000000e+00   True\n",
       "k_HO2_wall  4.000000e+00  0.000000e+00   True\n",
       "k1          1.800000e-12  2.700000e-13  False\n",
       "k2          1.110663e-10  1.665994e-11  False\n",
       "k3          1.449221e-12  2.108588e-13  False"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df_model_params = {}\n",
    "df_model_params['k_OH_wall'] =  {'val':15,       'err':0,            'fit':True }\n",
    "df_model_params['k_HO2_wall'] = {'val':4,        'err':0,            'fit':True }\n",
    "df_model_params['k1'] =         {'val':k1,       'err':k1_err,       'fit':False}\n",
    "df_model_params['k2'] =         {'val':k2,       'err':k2_err,       'fit':False}\n",
    "df_model_params['k3'] =         {'val':k3,       'err':k3_err,       'fit':False}\n",
    "df_model_params = pd.DataFrame.from_dict(df_model_params, orient='index')\n",
    "\n",
    "print('Inputted Model Params DataFrame:')\n",
    "display(df_model_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3: Set Initial Concentrations\n",
    "\n",
    "This 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)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_concentrations = {\n",
    "    'H2O2': c_H2O2_0,\n",
    "    'OH':   0,\n",
    "    'HO2':  0,\n",
    "    'H2O':  0,\n",
    "    'O2' :  0,\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4: Photolysis Parameters\n",
    "\n",
    "Set up your photolysis parameters in the next block.\n",
    "\n",
    "**df_photolysis_params : *DataFrame***\n",
    ">Each row is a species to be photolyzed, with the following columns:\n",
    "\n",
    ">**xsn : *float***\n",
    ">>The cross section of the species, in cm^2.\n",
    "\n",
    ">**products : *list of strings***\n",
    ">>List of products formed. Must correspond in order with *qyields*.\n",
    "\n",
    ">**qyields : *list of floats***\n",
    ">>List of quantum yields. Must correspond in order with *products*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Inputted Photolysis Params DataFrame:\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>xsn</th>\n",
       "      <th>products</th>\n",
       "      <th>qyields</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>H2O2</th>\n",
       "      <td>8.300000e-20</td>\n",
       "      <td>[OH]</td>\n",
       "      <td>[2.0]</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               xsn products qyields\n",
       "H2O2  8.300000e-20     [OH]   [2.0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fluence = 3.2e16# in photons cm^-2\n",
    "\n",
    "df_photolysis_params = {} ## Cross sections in cm^2, \n",
    "df_photolysis_params['H2O2'] = {'xsn': 8.3e-20, 'products': ['OH'], 'qyields': [2.0]}\n",
    "df_photolysis_params = pd.DataFrame.from_dict(df_photolysis_params, orient='index')\n",
    "\n",
    "print('Inputted Photolysis Params DataFrame:')\n",
    "display(df_photolysis_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Method: Plot the model\n",
    "\n",
    "Integrates 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*.\n",
    "\n",
    "\n",
    "See the file 'ex_init_model_2.csv' for an example of how the output is formatted if *save_fn* is specified.\n",
    "\n",
    "*method* **StaticCell.StaticCell.plot_model**(*t_start, t_end, tbin, df_model_params, , delta_xtick=20.0, save_fn=None*)\n",
    "\n",
    "**Parameters:**\n",
    "\n",
    ">**t_prephoto : *float***\n",
    ">>Amount of time to plot pre-photolysis.  Must be an integer multiple of *tbin* \\* 0.02 ms and cannot be less than -20 ms.\n",
    "\n",
    ">**t_react : *float***\n",
    ">>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*.\n",
    "\n",
    ">**tbin : *int***\n",
    ">>The time axis step size will be *tbin* \\* 0.02 ms.\n",
    "\n",
    ">**df_model_params : *DataFrame***\n",
    ">>Contains the model parameters.  See setup above for formatting.\n",
    "\n",
    ">**initial_concentrations : *dictionary***\n",
    ">>Contains the pre-photolysis concentrations of all species, including products and intermediates.\n",
    "\n",
    ">**df_photolysis_params : *DataFrame***\n",
    ">>Contains the photolysis parameters.  See setup above for formatting.\n",
    "\n",
    ">**fluence : *float***\n",
    ">>Laser fluence in photons cm^-2\n",
    "\n",
    ">**photolysis_cycles : *int***\n",
    ">>Number of photolysis pulses.\n",
    "\n",
    ">**delta_xtick : *float, optional***\n",
    ">>Tick marks and labels for the time axis include zero and are spaced by *delta_xtick* (ms).\n",
    "\n",
    ">**save_fn : *str, optional***\n",
    ">>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/cm<sup>3</sup>) for each species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "75823984ab444c9e9a0aa4b503d58379",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "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')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}