examples/process_simulation_output.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "7968866b-24b8-4ed6-bc13-6035310cb4f5",
   "metadata": {},
   "source": [
    "[![nbviewer](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg)](https://nbviewer.jupyter.org/github/open-atmos/PyPartMC/blob/main/examples/process_simulation_output.ipynb)   \n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/open-atmos/PyPartMC/blob/main/examples/process_simulation_output.ipynb)    \n",
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/open-atmos/PyPartMC.git/main?urlpath=lab/tree/examples/process_simulation_output.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "159edeb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This file is a part of PyPartMC licensed under the GNU General Public License v3\n",
    "# Copyright (C) 2023 University of Illinois Urbana-Champaign\n",
    "# Authors:\n",
    "#  - https://github.com/compdyn/partmc/graphs/contributors\n",
    "#  - https://github.com/open-atmos/PyPartMC/graphs/contributors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4f8359c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "if \"google.colab\" in sys.modules:\n",
    "    !pip --quiet install open-atmos-jupyter-utils\n",
    "    from open_atmos_jupyter_utils import pip_install_on_colab\n",
    "    pip_install_on_colab('PyPartMC')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b494ea6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from open_atmos_jupyter_utils import show_plot\n",
    "import PyPartMC as ppmc\n",
    "from PyPartMC import si"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "df8d1a84",
   "metadata": {},
   "outputs": [],
   "source": [
    "gas_data = ppmc.GasData((\"H2SO4\",\"HNO3\",\"HCl\",\"NH3\",\"NO\",\"NO2\", \"NO3\",\n",
    "                        \"N2O5\", \"HONO\", \"HNO4\", \"O3\", \"O1D\", \"O3P\", \"OH\",\n",
    "                        \"HO2\", \"H2O2\", \"CO\", \"SO2\", \"CH4\", \"C2H6\", \"CH3O2\", \n",
    "                        \"ETHP\", \"HCHO\", \"CH3OH\", \"ANOL\", \"CH3OOH\", \"ETHOOH\",\n",
    "                        \"ALD2\", \"HCOOH\", \"RCOOH\", \"C2O3\", \"PAN\", \"ARO1\", \"ARO2\",\n",
    "                        \"ALK1\", \"OLE1\", \"API1\", \"API2\", \"LIM1\", \"LIM2\", \"PAR\", \"AONE\",\n",
    "                        \"MGLY\", \"ETH\", \"OLET\", \"OLEI\", \"TOL\", \"XYL\", \"CRES\", \"TO2\",\n",
    "                         \"CRO\", \"OPEN\", \"ONIT\", \"ROOH\", \"RO2\", \"ANO2\", \"NAP\", \"XO2\",\n",
    "                         \"XPAR\", \"ISOP\", \"ISOPRD\", \"ISOPP\", \"ISOPN\", \"ISOPO2\", \"API\",\n",
    "                         \"LIM\", \"DMS\", \"MSA\", \"DMSO\", \"DMSO2\", \"CH3SO2H\", \"CH3SCH2OO\", \n",
    "                         \"CH3SO2\", \"CH3SO3\", \"CH3SO2OO\", \"CH3SO2CH2OO\", \"SULFHOX\"\n",
    "                        ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "a2d7bad8",
   "metadata": {},
   "outputs": [],
   "source": [
    "env_state = ppmc.EnvState(\n",
    "    {\n",
    "        \"rel_humidity\": 0.95,\n",
    "        \"latitude\": 0,\n",
    "        \"longitude\": 0,\n",
    "        \"altitude\": 0 * si.m,\n",
    "        \"start_time\": 21600 * si.s,\n",
    "        \"start_day\": 200,\n",
    "    }\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "815283a0",
   "metadata": {},
   "outputs": [],
   "source": [
    "aero_data = ppmc.AeroData(\n",
    "    (\n",
    "        #         density  ions in soln (1) molecular weight    kappa (1)\n",
    "        #         |                     |   |                   |\n",
    "        {\"SO4\": [1800 * si.kg / si.m**3, 1, 96.0 * si.g / si.mol, 0.00]},\n",
    "        {\"NO3\": [1800 * si.kg / si.m**3, 1, 62.0 * si.g / si.mol, 0.00]},\n",
    "        {\"Cl\": [2200 * si.kg / si.m**3, 1, 35.5 * si.g / si.mol, 0.00]},\n",
    "        {\"NH4\": [1800 * si.kg / si.m**3, 1, 18.0 * si.g / si.mol, 0.00]},\n",
    "        {\"MSA\": [1800 * si.kg / si.m**3, 0, 95.0 * si.g / si.mol, 0.53]},\n",
    "        {\"ARO1\": [1400 * si.kg / si.m**3, 0, 150.0 * si.g / si.mol, 0.10]},\n",
    "        {\"ARO2\": [1400 * si.kg / si.m**3, 0, 150.0 * si.g / si.mol, 0.10]},\n",
    "        {\"ALK1\": [1400 * si.kg / si.m**3, 0, 140.0 * si.g / si.mol, 0.10]},\n",
    "        {\"OLE1\": [1400 * si.kg / si.m**3, 0, 140.0 * si.g / si.mol, 0.10]},\n",
    "        {\"API1\": [1400 * si.kg / si.m**3, 0, 184.0 * si.g / si.mol, 0.10]},\n",
    "        {\"API2\": [1400 * si.kg / si.m**3, 0, 184.0 * si.g / si.mol, 0.10]},\n",
    "        {\"LIM1\": [1400 * si.kg / si.m**3, 0, 200.0 * si.g / si.mol, 0.10]},\n",
    "        {\"LIM2\": [1400 * si.kg / si.m**3, 0, 200.0 * si.g / si.mol, 0.10]},\n",
    "        {\"CO3\": [2600 * si.kg / si.m**3, 1, 60.0 * si.g / si.mol, 0.00]},\n",
    "        {\"Na\": [2200 * si.kg / si.m**3, 1, 23.0 * si.g / si.mol, 0.00]},\n",
    "        {\"Ca\": [2600 * si.kg / si.m**3, 1, 40.0 * si.g / si.mol, 0.00]},\n",
    "        {\"OIN\": [2600 * si.kg / si.m**3, 0, 1.0 * si.g / si.mol, 0.10]},\n",
    "        {\"OC\": [1400 * si.kg / si.m**3, 0, 1.0 * si.g / si.mol, 0.10]},\n",
    "        {\"BC\": [1800 * si.kg / si.m**3, 0, 1.0 * si.g / si.mol, 0.00]},\n",
    "        {\"H2O\": [1000 * si.kg / si.m**3, 0, 18.0 * si.g / si.mol, 0.00]},\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c822823d",
   "metadata": {},
   "outputs": [],
   "source": [
    "gas_state = ppmc.GasState(gas_data)\n",
    "\n",
    "input_gas_state = (\n",
    "    {\"NO\": [0.1E+00]},\n",
    "    {\"NO2\": [1.0E+00]},\n",
    "    {\"HNO3\": [1.0E+00]},\n",
    "    {\"O3\": [5.0E+01]},\n",
    "    {\"H2O2\": [1.1E+00]},\n",
    "    {\"CO\": [2.1E+02]},\n",
    "    {\"SO2\": [0.8E+00]},\n",
    "    {\"NH3\": [0.5E+00]},\n",
    "    {\"HCl\": [0.7E+00]},\n",
    "    {\"CH4\": [2.2E+03]},\n",
    "    {\"C2H6\": [1.0E+00]},\n",
    "    {\"HCHO\": [1.2E+00]},\n",
    "    {\"CH3OH\": [1.2E-01]},\n",
    "    {\"CH3OOH\": [0.5E+00]},\n",
    "    {\"ALD2\": [1.0E+00]},\n",
    "    {\"PAR\": [2.0E+00]},\n",
    "    {\"AONE\": [1.0E+00]},\n",
    "    {\"ETH\": [0.2E+00]},\n",
    "    {\"OLET\": [2.3E-02]},\n",
    "    {\"OLEI\": [3.1E-04]},\n",
    "    {\"TOL\": [0.1E+00]},\n",
    "    {\"XYL\": [0.1E+00]},\n",
    "    {\"ONIT\": [0.1E+00]},\n",
    "    {\"PAN\": [0.8E+00]},\n",
    "    {\"RCOOH\": [0.2E+00]},\n",
    "    {\"ROOH\": [2.5E-02]},\n",
    "    {\"ISOP\": [0.5E+00]}\n",
    ")\n",
    "\n",
    "gas_state.mix_rats = input_gas_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "58706963",
   "metadata": {},
   "outputs": [],
   "source": [
    "times = [0 * si.s]\n",
    "back_gas = [{\"time\": times},\n",
    "            {\"rate\": [1.5e-5 / si.s]},\n",
    "            {\"NO\": [0.1E+00]},\n",
    "            {\"NO2\": [1.0E+00]},\n",
    "            {\"HNO3\": [1.0E+00]},\n",
    "            {\"O3\": [5.0E+01]},\n",
    "            {\"H2O2\": [1.1E+00]},\n",
    "            {\"CO\": [2.1E+02]},\n",
    "            {\"SO2\": [0.8E+00]},\n",
    "            {\"NH3\": [0.5E+00]},\n",
    "            {\"HCl\": [0.7E+00]},\n",
    "            {\"CH4\": [2.2E+03]},\n",
    "            {\"C2H6\": [1.0E+00]},\n",
    "            {\"HCHO\": [1.2E+00]},\n",
    "            {\"CH3OH\": [1.2E-01]},\n",
    "            {\"CH3OOH\": [0.5E+00]},\n",
    "            {\"ALD2\": [1.0E+00]},\n",
    "            {\"PAR\": [2.0E+00]},\n",
    "            {\"AONE\": [1.0E+00]},\n",
    "            {\"ETH\": [0.2E+00]},\n",
    "            {\"OLET\": [2.3E-02]},\n",
    "            {\"OLEI\": [3.1E-04]},\n",
    "            {\"TOL\": [0.1E+00]},\n",
    "            {\"XYL\": [0.1E+00]},\n",
    "            {\"ONIT\": [0.1E+00]},\n",
    "            {\"PAN\": [0.8E+00]},\n",
    "            {\"RCOOH\": [0.2E+00]},\n",
    "            {\"ROOH\": [2.5E-02]},\n",
    "            {\"ISOP\": [0.5E+00]}\n",
    "           ]\n",
    "\n",
    "gas_emit_times = [0, 3600, 7200, 10800, 14400, 18000, 21600, 25200, 28800, 32400, 36000,\n",
    "                  39600, 43200, 46800, 50400, 54000, 57600, 61200, 64800, 68400, 72000,\n",
    "                  75600, 79200, 82800, 90000, 93600, 97200, 100800, 104400, 108000]\n",
    "\n",
    "gas_emit_rates = np.zeros(len(gas_emit_times))\n",
    "gas_emit_rates[0:12] = .5\n",
    "\n",
    "SO2 = [4.234E-09, 5.481E-09, 5.089E-09, 5.199E-09, 5.221E-09, 5.284E-09, 5.244E-09,\n",
    "       5.280E-09, 5.560E-09, 5.343E-09, 4.480E-09, 3.858E-09, 3.823E-09, 3.607E-09,\n",
    "       3.533E-09, 3.438E-09, 2.866E-09, 2.667E-09, 2.636E-09, 2.573E-09, 2.558E-09,\n",
    "       2.573E-09, 2.715E-09, 3.170E-09, 4.2344E-09, 5.481E-09, 5.089E-09, 5.199E-09,\n",
    "       5.221E-09, 5.284E-09]\n",
    "\n",
    "NO2 = [3.024E-09, 3.334E-09, 3.063E-09, 3.281E-09, 3.372E-09, 3.523E-09, 3.402E-09,\n",
    "       3.551E-09, 3.413E-09, 3.985E-09, 3.308E-09, 2.933E-09, 2.380E-09, 1.935E-09,\n",
    "       1.798E-09, 1.537E-09 , 9.633E-10, 8.873E-10, 7.968E-10, 6.156E-10, 5.920E-10,\n",
    "       6.320E-10, 9.871E-10, 1.901E-09, 3.024E-09, 3.334E-09, 3.063E-09, 3.281E-09,\n",
    "       3.372E-09, 3.523E-09]\n",
    "\n",
    "NO = [5.749E-08, 6.338E-08, 5.825E-08, 6.237E-08, 6.411E-08, 6.699E-08, 6.468E-08,\n",
    "      6.753E-08, 6.488E-08, 7.575E-08, 6.291E-08, 5.576E-08, 4.524E-08, 3.679E-08,\n",
    "      3.419E-08, 2.924E-08, 1.832E-08, 1.687E-08, 1.515E-08, 1.171E-08, 1.125E-08,\n",
    "      1.202E-08, 1.877E-08, 3.615E-08, 5.749E-08, 6.338E-08, 5.825E-08, 6.237E-08,\n",
    "      6.411E-08, 6.699E-08]\n",
    "\n",
    "CO = [7.839E-07, 5.837E-07, 4.154E-07, 4.458E-07, 4.657E-07, 4.912E-07, 4.651E-07,\n",
    "      4.907E-07, 6.938E-07, 8.850E-07, 8.135E-07, 4.573E-07, 3.349E-07, 2.437E-07,\n",
    "      2.148E-07, 1.662E-07, 8.037E-08, 7.841E-08, 6.411E-08, 2.551E-08, 2.056E-08,\n",
    "      3.058E-08, 1.083E-07, 3.938E-07, 7.839E-07, 5.837E-07, 4.154E-07, 4.458E-07,\n",
    "      4.657E-07, 4.912E-07]\n",
    "\n",
    "NH3 = [8.93E-09, 8.705E-09, 1.639E-08, 1.466E-08, 1.6405E-08, 1.8805E-08, 1.65E-08,\n",
    "       1.8045E-08, 1.347E-08, 6.745E-09, 5.415E-09, 2.553E-09, 2.087E-09, 2.2885E-09,\n",
    "       2.7265E-09, 2.7338E-09, 9.96E-10, 2.707E-09, 9.84E-10, 9.675E-10, 9.905E-10,\n",
    "       1.0345E-09, 1.0825E-09, 2.7465E-09, 8.93E-09, 8.705E-09, 1.639E-08, 1.466E-08,\n",
    "       1.6405E-08, 1.8805E-08]\n",
    "\n",
    "HCHO = [4.061E-09, 3.225E-09, 2.440E-09, 2.639E-09, 2.754E-09, 2.888E-09, 2.741E-09,\n",
    "        2.885E-09, 4.088E-09, 5.186E-09, 4.702E-09, 2.601E-09, 1.923E-09, 1.412E-09,\n",
    "        1.252E-09, 9.776E-10, 4.687E-10, 4.657E-10, 3.836E-10, 1.717E-10, 1.448E-10,\n",
    "        1.976E-10, 6.193E-10, 2.090E-09, 4.061E-09, 3.225E-09, 2.440E-09, 2.639E-09,\n",
    "        2.754E-09, 2.888E-09]\n",
    "\n",
    "ALD2 = [1.702E-09, 1.283E-09, 9.397E-10, 1.024E-09, 1.076E-09, 1.132E-09, 1.068E-09,\n",
    "        1.130E-09, 1.651E-09, 2.132E-09, 1.985E-09, 1.081E-09, 7.847E-10, 5.676E-10,\n",
    "        5.003E-10, 3.838E-10, 1.784E-10, 1.766E-10, 1.430E-10, 5.173E-11, 4.028E-11,\n",
    "        6.349E-11, 2.428E-10, 8.716E-10, 1.7022E-09, 1.283E-09, 9.397E-10, 1.024E-09,\n",
    "        1.076E-09, 1.132E-09]\n",
    "\n",
    "ETH = [1.849E-08, 1.391E-08, 1.010E-08, 1.095E-08, 1.148E-08, 1.209E-08, 1.142E-08,\n",
    "       1.205E-08, 1.806E-08, 2.320E-08, 2.149E-08, 1.146E-08, 8.384E-09, 6.124E-09,\n",
    "       5.414E-09, 4.119E-09, 1.953E-09, 1.927E-09, 1.575E-09, 6.164E-10, 4.973E-10,\n",
    "       7.420E-10, 2.653E-09, 9.477E-09, 1.849E-08, 1.391E-08, 1.010E-08, 1.095E-08,\n",
    "       1.148E-08, 1.209E-08]\n",
    "\n",
    "OLEI = [5.948E-09, 4.573E-09, 3.374E-09, 3.668E-09, 3.851E-09, 4.050E-09, 3.841E-09,\n",
    "        4.052E-09, 6.094E-09, 7.795E-09, 7.215E-09, 3.738E-09, 2.718E-09, 1.973E-09,\n",
    "        1.729E-09, 1.338E-09, 6.333E-10, 6.394E-10, 5.126E-10, 2.089E-10, 1.708E-10,\n",
    "        2.480E-10, 8.947E-10, 3.057E-09, 5.948E-09, 4.573E-09, 3.374E-09, 3.668E-09,\n",
    "        3.851E-09, 4.050E-09]\n",
    "\n",
    "OLET = [5.948E-09, 4.573E-09, 3.374E-09, 3.668E-09, 3.851E-09, 4.050E-09, 3.841E-09,\n",
    "        4.052E-09, 6.094E-09, 7.795E-09, 7.215E-09, 3.738E-09, 2.718E-09, 1.973E-09,\n",
    "        1.729E-09, 1.338E-09, 6.333E-10, 6.394E-10, 5.126E-10, 2.089E-10, 1.708E-10,\n",
    "        2.480E-10, 8.947E-10, 3.057E-09, 5.948E-09, 4.573E-09, 3.374E-09, 3.668E-09,\n",
    "        3.851E-09, 4.050E-09]\n",
    "\n",
    "TOL = [6.101E-09, 8.706E-09, 7.755E-09, 8.024E-09, 8.202E-09, 8.410E-09, 8.218E-09,\n",
    "       8.407E-09, 1.020E-08, 1.139E-08, 7.338E-09, 4.184E-09, 3.078E-09, 2.283E-09,\n",
    "       2.010E-09, 1.575E-09, 8.966E-10, 6.705E-10, 5.395E-10, 2.462E-10, 2.106E-10,\n",
    "       2.852E-10, 9.300E-10, 3.144E-09, 6.101E-09, 8.706E-09, 7.755E-09, 8.024E-09,\n",
    "       8.202E-09, 8.410E-09]\n",
    "\n",
    "XYL = [5.599E-09, 4.774E-09, 3.660E-09, 3.909E-09, 4.060E-09, 4.239E-09, 4.060E-09,\n",
    "       4.257E-09, 6.036E-09, 7.448E-09, 6.452E-09, 3.435E-09, 2.525E-09, 1.859E-09,\n",
    "       1.650E-09, 1.302E-09, 6.852E-10, 6.773E-10, 5.437E-10, 2.697E-10, 2.358E-10,\n",
    "       3.059E-10, 8.552E-10, 2.861E-10, 5.599E-09, 4.774E-09, 3.660E-09, 3.909E-09,\n",
    "       4.060E-09, 4.239E-09]\n",
    "\n",
    "AONE = [7.825E-10, 2.858E-09, 2.938E-09, 2.947E-09, 2.948E-09, 2.951E-09, 2.947E-09,\n",
    "        2.954E-09, 3.032E-09, 2.766E-09, 1.313E-09, 1.015E-09, 8.363E-10, 7.040E-10,\n",
    "        6.404E-10, 6.264E-10, 5.661E-10, 1.538E-10, 1.500E-10, 1.395E-10, 1.476E-10,\n",
    "        1.503E-10, 2.256E-10, 4.244E-10, 7.825E-10, 2.858E-09, 2.938E-09, 2.947E-09,\n",
    "        2.948E-09, 2.951E-09]\n",
    "\n",
    "PAR = [1.709E-07, 1.953E-07, 1.698E-07, 1.761E-07, 1.808E-07, 1.865E-07, 1.822E-07,\n",
    "       1.8599E-07, 2.412E-07, 2.728E-07, 2.174E-07, 1.243E-07, 9.741E-08, 7.744E-08,\n",
    "       6.931E-08, 5.805E-08, 3.900E-08, 3.317E-08, 2.956E-08, 2.306E-08, 2.231E-08,\n",
    "       2.395E-08, 4.284E-08, 9.655E-08, 1.709E-07, 1.953E-07, 1.698E-07, 1.761E-07,\n",
    "       1.808E-07, 1.865E-07]\n",
    "\n",
    "ISOP = [2.412E-10, 2.814E-10, 3.147E-10, 4.358E-10, 5.907E-10, 6.766E-10, 6.594E-10,\n",
    "        5.879E-10, 5.435E-10, 6.402E-10, 5.097E-10, 9.990E-11, 7.691E-11, 5.939E-11,\n",
    "        5.198E-11, 4.498E-11, 3.358E-11, 2.946E-11, 2.728E-11, 2.183E-11, 1.953E-11,\n",
    "        1.890E-11, 2.948E-11, 1.635E-10, 2.412E-10, 2.814E-10, 3.147E-10, 4.358E-10,\n",
    "        5.907E-10, 6.766E-10]\n",
    "\n",
    "CH3OH = [2.368E-10, 6.107E-10, 6.890E-10, 6.890E-10, 6.890E-10, 6.889E-10, 6.886E-10,\n",
    "         6.890E-10, 6.890E-10, 5.414E-10, 3.701E-10, 2.554E-10, 1.423E-10, 6.699E-11,\n",
    "         2.912E-11, 2.877E-11, 2.825E-11, 2.056E-12, 2.056E-12, 2.056E-12, 2.435E-12,\n",
    "         2.435E-12, 4.030E-11, 1.168E-10, 2.368E-10, 6.107E-10, 6.890E-10, 6.890E-10,\n",
    "         6.890E-10, 6.889E-10]\n",
    "\n",
    "ANOL = [5.304E-09, 7.960E-09, 7.649E-09, 7.649E-09, 7.432E-09, 7.428E-09,\n",
    "        7.431E-09, 7.434E-09, 7.434E-09, 6.979E-09,5.666E-09, 4.361E-09, 4.148E-09,\n",
    "        3.289E-09, 2.858E-09, 2.856E-09, 1.127E-09, 9.615E-10, 9.616E-10, 9.616E-10,\n",
    "        9.654E-10, 9.654E-10, 1.397E-09, 2.264E-09, 5.304E-09, 7.960E-09, 7.649E-09,\n",
    "        7.649E-09, 7.432E-09, 7.428E-09]\n",
    "\n",
    "emit_gas = [\n",
    "    {\"time\": gas_emit_times},\n",
    "    {\"rate\": list(gas_emit_rates)},\n",
    "    {\"SO2\": SO2},\n",
    "    {\"NO\": NO},\n",
    "    {\"NO2\": NO2},\n",
    "    {\"CO\": CO},\n",
    "    {\"NH3\": NH3},\n",
    "    {\"HCHO\": HCHO},\n",
    "    {\"ALD2\": ALD2},\n",
    "    {\"ETH\": ETH},\n",
    "    {\"OLEI\": OLEI},\n",
    "    {\"OLET\": OLET},\n",
    "    {\"TOL\": TOL},\n",
    "    {\"XYL\": XYL},\n",
    "    {\"AONE\": AONE},\n",
    "    {\"PAR\": PAR},\n",
    "    {\"ISOP\": ISOP},\n",
    "    {\"CH3OH\": CH3OH},\n",
    "    {\"ANOL\": ANOL},\n",
    "]\n",
    "\n",
    "AERO_DIST_BACKGROUND = {\n",
    "    \"back_small\": {\n",
    "        \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
    "        \"diam_type\": \"geometric\",\n",
    "        \"mode_type\": \"log_normal\",\n",
    "        \"num_conc\": 3.2e9 / si.m**3,\n",
    "        \"geom_mean_diam\": 0.02 * si.um,\n",
    "        \"log10_geom_std_dev\": 0.161,\n",
    "    },\n",
    "    \"back_large\": {\n",
    "        \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
    "        \"diam_type\": \"geometric\",\n",
    "        \"mode_type\": \"log_normal\",\n",
    "        \"num_conc\": 2.9e9 / si.m**3,\n",
    "        \"geom_mean_diam\": 0.16 * si.um,\n",
    "        \"log10_geom_std_dev\": 0.217,\n",
    "    },\n",
    "}\n",
    "\n",
    "AERO_DIST_EMIT = {\n",
    "    \"gasoline\": {\n",
    "        \"mass_frac\": [{\"OC\": [0.8]}, {\"BC\": [0.2]}],\n",
    "        \"diam_type\": \"geometric\",\n",
    "        \"mode_type\": \"log_normal\",\n",
    "        \"num_conc\": 5e7 / si.m**3,\n",
    "        \"geom_mean_diam\": 5e-8 * si.m,\n",
    "        \"log10_geom_std_dev\": 0.24,\n",
    "    },\n",
    "    \"diesel\": {\n",
    "        \"mass_frac\": [{\"OC\": [0.3]}, {\"BC\": [0.7]}],\n",
    "        \"diam_type\": \"geometric\",\n",
    "        \"mode_type\": \"log_normal\",\n",
    "        \"num_conc\": 1.6e8 / si.m**3,\n",
    "        \"geom_mean_diam\": 5e-8 * si.m,\n",
    "        \"log10_geom_std_dev\": 0.24,\n",
    "    },\n",
    "    \"cooking\": {\n",
    "        \"mass_frac\": [{\"OC\": [1]}],\n",
    "        \"diam_type\": \"geometric\",\n",
    "        \"mode_type\": \"log_normal\",\n",
    "        \"num_conc\": 9e6 / si.m**3,\n",
    "        \"geom_mean_diam\": 8.64e-8 * si.m,\n",
    "        \"log10_geom_std_dev\": 0.28,\n",
    "    },\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c6a96b7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_timeseries = list(np.linspace(0,24*3600,25))\n",
    "pressure_timeseries = list(np.ones(25) * 1e5)\n",
    "temp_timeseries = [290.016,292.5, 294.5, 296.112, 297.649, 299.049, 299.684, 299.509,299.002,\n",
    "                   298.432, 296.943, 295.153, 293.475, 292.466, 291.972, 291.96, 291.512,\n",
    "                   291.481, 290.5,   290.313, 290.317, 290.362, 290.245, 290.228, 291.466]\n",
    "height_timeseries = [171.045, 228.210, 296.987, 366.002, 410.868, 414.272, 417.807,414.133,\n",
    "                     397.465, 376.864, 364.257, 352.119, 338.660, 322.028, 305.246, 258.497,  \n",
    "                     240.478, 187.229, 145.851, 128.072, 110.679, 97.628, 93.034, 93.034, 93.034]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "920e41e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "scenario = ppmc.Scenario(\n",
    "    gas_data,\n",
    "    aero_data,\n",
    "    {\n",
    "        \"temp_profile\": [{\"time\": time_timeseries}, {\"temp\": temp_timeseries}],\n",
    "        \"pressure_profile\": [\n",
    "            {\"time\": time_timeseries},\n",
    "            {\"pressure\": pressure_timeseries},\n",
    "        ],\n",
    "        \"height_profile\": [{\"time\": time_timeseries}, {\"height\": height_timeseries}],\n",
    "        \"gas_emissions\": emit_gas,\n",
    "        \"gas_background\": back_gas,\n",
    "        \"aero_emissions\": [\n",
    "            {\"time\": [0 * si.s, 12 * 3600 * si.s]},\n",
    "            {\"rate\": [1 / si.s, 0 / si.s]},\n",
    "            {\"dist\": [[AERO_DIST_EMIT],[AERO_DIST_EMIT]]},\n",
    "        ],\n",
    "        \"aero_background\": [\n",
    "            {\"time\": [0 * si.s]},\n",
    "            {\"rate\": [1.5e-5 / si.s]},\n",
    "            {\"dist\": [[AERO_DIST_BACKGROUND]]},\n",
    "        ],\n",
    "        \"loss_function\": \"none\",\n",
    "    },\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6722ba83",
   "metadata": {},
   "outputs": [],
   "source": [
    "T_INITIAL = 0.0\n",
    "scenario.init_env_state(env_state, T_INITIAL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "9781ca2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "AERO_DIST_INIT = [\n",
    "    {\n",
    "        \"init_small\": {\n",
    "            \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
    "            \"diam_type\": \"geometric\",\n",
    "            \"mode_type\": \"log_normal\",\n",
    "            \"num_conc\": 3.2e9 / si.m**3,\n",
    "            \"geom_mean_diam\": 0.02 * si.um,\n",
    "            \"log10_geom_std_dev\": 0.161,\n",
    "        },\n",
    "        \"init_large\": {\n",
    "            \"mass_frac\": [{\"SO4\": [1]}, {\"OC\": [1.375]}, {\"NH4\": [0.375]}],\n",
    "            \"diam_type\": \"geometric\",\n",
    "            \"mode_type\": \"log_normal\",\n",
    "            \"num_conc\": 2.9e9 / si.m**3,\n",
    "            \"geom_mean_diam\": 0.16 * si.um,\n",
    "            \"log10_geom_std_dev\": 0.217,\n",
    "        },\n",
    "    }\n",
    "]\n",
    "\n",
    "aero_dist_init = ppmc.AeroDist(aero_data, AERO_DIST_INIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d8d9c8fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "284"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "run_part_opt = ppmc.RunPartOpt(\n",
    "    {\n",
    "        \"output_prefix\": \"urban_plume\",\n",
    "        \"do_coagulation\": True,\n",
    "        \"coag_kernel\": \"brown\",\n",
    "        \"t_max\": 86400 * si.s,\n",
    "        \"del_t\": 60 * si.s,\n",
    "        \"t_output\": 3600.0,\n",
    "    }\n",
    ")\n",
    "\n",
    "N_PART = 1000\n",
    "aero_state = ppmc.AeroState(aero_data, N_PART, 'nummass_source')\n",
    "aero_state.dist_sample(\n",
    "    aero_dist_init,\n",
    "    sample_prop=1.0,\n",
    "    create_time=0.0,\n",
    "    allow_doubling=True,\n",
    "    allow_halving=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "5c8aef2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "camp_core = ppmc.CampCore()\n",
    "photolysis = ppmc.Photolysis()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f7c29fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "ppmc.run_part(\n",
    "        scenario,\n",
    "        env_state,\n",
    "        aero_data,\n",
    "        aero_state,\n",
    "        gas_data,\n",
    "        gas_state,\n",
    "        run_part_opt,\n",
    "        camp_core,\n",
    "        photolysis,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0130f477",
   "metadata": {},
   "source": [
    "## Example of post processing from netCDF files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47474cd3",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_TIMES = 25\n",
    "rh = np.zeros(N_TIMES)\n",
    "temperature = np.zeros(N_TIMES)\n",
    "height = np.zeros(N_TIMES)\n",
    "time = np.zeros(N_TIMES)\n",
    "\n",
    "for i_time in range(N_TIMES):\n",
    "    ppmc.input_state(f\"urban_plume_0001_{i_time+1:08}.nc\",\n",
    "                               aero_data, aero_state, gas_data, gas_state,\n",
    "                               env_state)\n",
    "    rh[i_time] = env_state.rh\n",
    "    temperature[i_time] = env_state.temp\n",
    "    height[i_time] = env_state.height\n",
    "    time[i_time] = env_state.elapsed_time\n",
    "plt.plot(time,temperature,'r')\n",
    "plt.ylabel('Temperature (K)', color='r')\n",
    "plt.ylim([275,300])\n",
    "plt.xticks(np.linspace(0, time[-1], 5))\n",
    "plt.xlim([0,time[-1]])\n",
    "plt.xlabel('Time (s)')\n",
    "plt.twinx()\n",
    "plt.plot(time,rh*100,'g')\n",
    "plt.ylabel('Relative humidity (%)', color='g')\n",
    "plt.ylim([50,100])\n",
    "show_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e848f60b",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time,height)\n",
    "plt.xticks(np.linspace(0, time[-1], 5))\n",
    "plt.xlim([0,time[-1]])\n",
    "plt.xlabel('Time (s)')\n",
    "plt.ylim([0,500])\n",
    "plt.grid()\n",
    "plt.ylabel('Boundary layer height (m)')\n",
    "show_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e1b89e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_conc = np.zeros(N_TIMES)\n",
    "mass_conc = np.zeros(N_TIMES)\n",
    "for i_time in range(N_TIMES):\n",
    "    filename = ppmc.input_state(f\"urban_plume_0001_{i_time+1:08}.nc\",\n",
    "                               aero_data, aero_state, gas_data, gas_state,\n",
    "                               env_state)\n",
    "    num_conc[i_time] = aero_state.total_num_conc\n",
    "    mass_conc[i_time] = aero_state.total_mass_conc\n",
    "\n",
    "plt.plot(time, mass_conc, \"b\", label=\"mass conc\")\n",
    "plt.ylabel(\"Mass concentration\", color='b')\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.twinx()\n",
    "plt.plot(time, num_conc, \"g\", label=\"num conc\")\n",
    "plt.xticks(np.linspace(0, time[-1], 5))\n",
    "plt.ylabel(\"Number concentration\", color='g')\n",
    "show_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "602a2a01",
   "metadata": {},
   "outputs": [],
   "source": [
    "chi = np.zeros(N_TIMES)\n",
    "for i_time in range(N_TIMES):\n",
    "    filename = ppmc.input_state(f\"urban_plume_0001_{i_time+1:08}.nc\",\n",
    "                               aero_data, aero_state, gas_data, gas_state,\n",
    "                               env_state)\n",
    "    chi[i_time] = aero_state.mixing_state[2]\n",
    "\n",
    "plt.plot(time, chi, \"g\") #, label=\"num conc\")\n",
    "plt.xticks(np.linspace(0, time[-1], 5))\n",
    "plt.ylabel(\"Mixing state\")\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.xticks(np.linspace(0, time[-1], 5))\n",
    "plt.ylim([0,1])\n",
    "plt.xlim([time[0],time[-1]])\n",
    "show_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "386aa7c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "gases = [\"NO\", \"NO2\", \"O3\", \"HNO3\", \"SO2\", \"NH3\"]\n",
    "gas_mix_rat = np.zeros((N_TIMES, gas_data.n_spec))\n",
    "for i_time in range(N_TIMES):\n",
    "    filename = ppmc.input_state(f\"urban_plume_0001_{i_time+1:08}.nc\",\n",
    "                               aero_data, aero_state, gas_data, gas_state,\n",
    "                               env_state)\n",
    "    gas_mix_rat[i_time,:] = gas_state.mix_rats\n",
    "\n",
    "for i_spec, spec in enumerate(gases):\n",
    "    i_spec = gas_data.spec_by_name(spec)\n",
    "    l, = plt.plot(time, gas_mix_rat[:, i_spec], label=spec)\n",
    "plt.xlabel(\"Time (s)\")\n",
    "plt.ylabel(\"Mixing ratio (ppb)\")\n",
    "plt.xticks(np.linspace(0, time[-1], 5))\n",
    "plt.legend()\n",
    "show_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faa1de28",
   "metadata": {},
   "outputs": [],
   "source": [
    "diam_grid = ppmc.BinGrid(60, \"log\", 1e-9, 1e-6)\n",
    "outputs_to_plot = [0,1,6,12,24]\n",
    "for i, i_time in enumerate(outputs_to_plot):\n",
    "    filename = ppmc.input_state(f\"urban_plume_0001_{i_time+1:08}.nc\",\n",
    "                               aero_data, aero_state, gas_data, gas_state,\n",
    "                               env_state)\n",
    "    dry_diameters = aero_state.dry_diameters\n",
    "    num_concs = aero_state.num_concs\n",
    "    dist = ppmc.histogram_1d(diam_grid, dry_diameters, num_concs)\n",
    "    plt.plot(diam_grid.centers, dist, label=f'$t =$ {time[i_time]} s')\n",
    "plt.xscale(\"log\")\n",
    "plt.xlabel(\"Dry diameter (m)\")\n",
    "plt.ylim(bottom=0)\n",
    "plt.legend()\n",
    "show_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7015860b",
   "metadata": {},
   "outputs": [],
   "source": [
    "I_OUTPUT = 24\n",
    "ppmc.input_state(f\"urban_plume_0001_{I_OUTPUT+1:08}.nc\", aero_data, aero_state, \n",
    "                 gas_data, gas_state, env_state)\n",
    "\n",
    "mass_frac_grid = ppmc.BinGrid(100, \"linear\", 0, 1)\n",
    "mass_frac_histogram = []\n",
    "num_concs = aero_state.num_concs\n",
    "dry_diameters = aero_state.dry_diameters\n",
    "species_mass_frac = np.zeros(len(aero_state))\n",
    "species_index = aero_data.spec_by_name(\"BC\")\n",
    "for i_part in range(len(aero_state)):\n",
    "    particle = aero_state.particle(i_part)\n",
    "    masses = particle.species_masses\n",
    "    species_mass_frac[i_part] = masses[species_index] / np.sum(masses[0:19])\n",
    "vals = ppmc.histogram_2d(\n",
    "    diam_grid, dry_diameters, mass_frac_grid, species_mass_frac, num_concs\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07fabfba",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.pcolormesh(\n",
    "    diam_grid.edges,\n",
    "    mass_frac_grid.edges,\n",
    "    np.array(vals).T / 1e6,\n",
    "    norm=matplotlib.colors.LogNorm(),\n",
    ")\n",
    "plt.xscale(\"log\")\n",
    "plt.xlim([1e-8, 1e-6])\n",
    "plt.colorbar(label=\"Number concentration (# cm$^{-3}$)\")\n",
    "plt.xlabel(\"Dry diameter (m)\")\n",
    "plt.ylabel(\"BC mass fraction\")\n",
    "plt.grid()\n",
    "plt.title(f'$t =$ {time[I_OUTPUT]} s')\n",
    "show_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3ef6f38",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}