{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0f55d04e8790403b9c149d5f8cbb8e77",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(VBox(children=(HBox(children=(Dropdown(description='Solver for Z', options=('Newton', 'Bisectio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import main\n",
    "main.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The program is written so that an arbitrary root-solver and ODE-step function can only needs to be implemented in the module numerikk.py and added to the respectve dictionary Z_solver_dict or step_func_dict to automatically become availible to the interface and function with the rest of the program. In the solvers used to find $Z$, the tolarance decreases linearly with $\\frac{K}{H}$, because $lim_{Z\\rightarrow 1}Z^2H^2 = H^2$ and $Z$ deviates more from one with decreasing $\\frac{K}{H}$.The differential function \n",
    "\\begin{equation}\n",
    "\\frac{d[\\%H]}{dt} = f([\\%H]) = C_1Z^2[\\%H]^2 + C_2[\\%H] + C_3\n",
    "\\end{equation}\n",
    "is implemented in differensial.py as a callable object, so that $C_3$ (which depends on the input parameter $H_v$) as well as $Z([\\%H])$ for any given $[\\%H]$ and corresponding plots for the function used to find $Z$ (displayed to the right) easily can be retrieved by any part of the program. This also reduces run-time because $C_3$ only needs to be re-evaluated if $H_v$ is changed. All constants are defined in constants.py. $C_1$, $C_2$ and $C_3$ are given by:\n",
    "\\begin{split}\n",
    "    C_0 &= -\\frac{100M_H}{M}\\\\\n",
    "    C_1 &= C_0\\frac{2f_H^2K_H^2G_m}{p_{in}}\\\\\n",
    "    C_2 &= C_0\\frac{k_{ts}A_s\\rho}{100M_H}\\\\\n",
    "    C_3 &= -C_2H_v\n",
    "\\end{split}\n",
    "\n",
    "Z is found by finding the root of the function \n",
    "\\begin{equation}\n",
    "    Z + \\ln{(1-Z)} + \\frac{K}{H} = 0\n",
    "\\end{equation}\n",
    "where \n",
    "\\begin{equation}\n",
    "    K = \\frac{\\rho k_t A p_{in}}{400M_H G_mm f_H^2 K_H^2}\n",
    "\\end{equation}\n",
    "The solvers that have been used to find $Z$ are Newtons method, Bisection method and scipy.optimize.fsolve. ODE-solvers are Runge-Kutta 4, Heuns method (improved explicit Euler) and Implicit euler, where Implicit Euler uses the either the bisection method or Newtons method to solve the resulting root-finding problem. If Newtons method is chosen, $Z$ must be calculated for each iteration, here, the same solver that is chosen for calculating $Z$ when evaluating $[\\%H]$ is used (either, Newtons, Bisection or fsolve). The Runge-Kutta 4 method is an explicit method, implemented with stepsize $h$ as:\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "    k_1 &= f(y_n)\\\\\n",
    "    k_2 &= f(y_n + \\frac{h}{2}k_1)\\\\\n",
    "    k_3 &= f(y_n + \\frac{h}{2}k_2)\\\\\n",
    "    k_4 &= f(y_n + hk_3)\\\\\n",
    "    y_{n+1} &= y_n + \\frac{h}{6}(k_1 +2k_2 +2k_3 +k_4)\\\\\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "Heuns method is also an explicit method, it is implemented with stepsize $h$ as:\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "    y* &= y_n + f(y_n)\\\\\n",
    "    y_{n+1} &= y_n + \\frac{h}{2}\\big[f(y*) + f(y_n)\\big].\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "Implicit Euler is given by:\n",
    "\\begin{equation}\n",
    "    y_{n+1} = y_n + \\frac{h}{2}\\big[f(y_n) + f(y_{n+1})\\big],\n",
    "\\end{equation}\n",
    "Where $y_{n+1}$ must be found numerically by solving the corresponding root-equation\n",
    "\\begin{equation}\n",
    "     g(y_{n+1}) = y_n + \\frac{h}{2}\\big[f(y_n) + f(y_{n+1})\\big] - y_{n+1} = 0\n",
    "\\end{equation}\n",
    "Doing this by Bisection is trivial, however solving by newtons method requires differenciating $g$ with respect to $[\\%H]$. From here on out, $H:=[\\%H]$, for ease of notation. As we already know, $H = H(Z(H))$ where $Z$ is found numericaly. We can still find the derivative of $Z$ as an explicit function of $Z$ and $H$:\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "    Z + \\ln{(1-Z)} + \\frac{K}{H} &= 0\\\\\n",
    "    \\frac{dZ}{dH} - \\frac{\\frac{dZ}{dH}}{1-Z} - \\frac{K}{H^2} &= 0\\\\\n",
    "    \\frac{dZ}{dH} &= \\frac{K(Z-1)}{H^2Z}\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "We can use this to aquire $\\frac{dg}{dH}$ as an explicit function of $H$ and $Z$:\n",
    "\\begin{equation}\n",
    "\\begin{split}\n",
    "    g(H) &= C_1Z^2H^2 + C_2H + C_3\\\\\n",
    "    \\frac{dg}{dH} &= 4C_1ZH\\frac{dZ}{dH} + C_2\\\\\n",
    "    \\frac{dg}{dH} &= 4C_1\\frac{K(1-Z)}{H} + C_2\\\\\n",
    "\\end{split}\n",
    "\\end{equation}\n",
    "When implementing Newtons method:\n",
    "\\begin{equation}\n",
    "     y_{n+1} = y_n - \\frac{g(y_n)}{g'(y_n)}\n",
    "\\end{equation}\n",
    "to solve the root equation $g(y_{n+1})=0$ we now have to find $Z$ numerically in each iteration. Starting value is chosen to be the previous time step, and Newtons method can be written as:\n",
    "\\begin{split}\n",
    "    H_{n+1, k+1} &= H_{n+1, k} - \\frac{g(H_{n+1, k})}{g'(H_{n+1, k})}\\\\\n",
    "    H_{n+1, 0} &= H_n\n",
    "\\end{split}\n",
    "where $n$ are time-steps and $k$ are iterations in Newtons method for finding the next time step."
   ]
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}