{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Integration of the FTE\n",
    "\n",
    "Formal transfer equation:\n",
    "$$ \\frac{dI(s)}{ds} = \\alpha(s) (S(s) - I(s)) $$\n",
    "\n",
    "Ignoring scattering, we have $S(s) = B(s)$. Assuming a Gaussian form for $B(s)$:\n",
    "$$ B(s) = B_0 + (B_1 - B_0) \\exp\\left(-\\frac{s^2}{2w^2}\\right) $$\n",
    "\n",
    "and a simple box function for the extinction coefficient:\n",
    "` alpha = A if abs(s) <= L else 0`\n",
    "\n",
    "* At $s < -L$, $I = I_{in}$.\n",
    "* Want to compute $I_{out}$ at $s > L$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set default parameters\n",
    "L = 3\n",
    "B0 = 1\n",
    "B1 = 3\n",
    "w = 1\n",
    "I_init = 0\n",
    "A = 1\n",
    "\n",
    "# setup grid within -L to L with N = 100.\n",
    "N = 100\n",
    "h = 2 * L/N\n",
    "\n",
    "s = np.linspace(-L, L, N)\n",
    "I = np.zeros(N)\n",
    "I[0] = I_init"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 1\n",
    "Total optical depth is [given by](https://en.wikipedia.org/wiki/Optical_depth_(astrophysics))\n",
    "\n",
    "$$ \\tau = \\int_{-L}^{L} \\alpha dz $$\n",
    "\n",
    "which gives $\\tau = 2A L = 6$.\n",
    "\n",
    "Note that we have the following given **parameters**:\n",
    "* L = 3\n",
    "* A = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 2\n",
    "\n",
    "Use simple forward Euler to get $I(s)$ at all grid points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define some helper functions\n",
    "\n",
    "def f(alpha, source, intensity):\n",
    "    \"\"\"\n",
    "    RHS of the FTE\n",
    "    \"\"\"\n",
    "    return alpha * (source  - intensity)\n",
    "\n",
    "def B(B0, B1, s, w):\n",
    "    \"\"\"\n",
    "    Gaussian Planck function\n",
    "    \"\"\"\n",
    "    return B0 + (B1 - B0) * np.exp(-1/2 * s**2/w**2)\n",
    "\n",
    "def alpha(s, A, L):\n",
    "    \"\"\"\n",
    "    Extinction coefficient\n",
    "    \"\"\"\n",
    "    return A if abs(s) <= L else 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's iterate using explicit Euler\n",
    "# I[0] = I_init\n",
    "\n",
    "for i in range(1, N):\n",
    "    I[i] = h * f(alpha(s[i-1], A, L), B(B0, B1, s[i-1], w), I[i-1]) + I[i-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Numerical solution with the forward Euler method')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xu4XGV59/HvzxggHKMlrWQbDCpSqdIEt4jigVoVBIWIWrBUpYqUV1H0EhSrryKKYOmloGgtVgooAlaRomBRCySA5bBz4ih9IwdzQEjAcIxCyP3+sZ7ZWZnMYc3sWXtOv8917WvPrFmz1rNm1sw9676fZy1FBGZmZgDP6HYDzMysdzgomJnZOAcFMzMb56BgZmbjHBTMzGycg4KZmY1zUOhRkl4j6c4OLOceSW/oRJtyywxJL2zzuR3ZrnY1W7+k2Wn7njmBdXTsNZe0j6T/J+kxSfM6scyylLGvpeWeKOl7nV5uWSby+aixrFJe00aGNiikF/t+Sdvkph0p6eouNmtcRFwTEbt1ux0TVf0B6fZ2Va9/oh86SedI+mJnWlfTScCZEbFtRFxS4npKlV6nJ1Nwq/wt7Xa7JkrS1ZKO7HY7Omlog0LyTODYbjei2kR+pdrAeR5wWztPLGs/msBy/ykFt8rfX3a0YVWUGfbvuJYN+wt2GnCcpOnVD9RKI+R/FUg6QtJ1kr4qaa2kuyS9Kk1fLukBSe/NPXdLSf8s6bfpCOVbkqalx/aVtELSJyX9Dvj3yrTc82dJuljSakkPSjozTX+BpCvTtDWSzq+1PbVIOkDS7ZIelbRS0nG5xz4gaZmkhyRdKmlmnWVs8kspbf+16faCNHlp+mV4aI3tenFaxlpJt0k6KPfYOZK+Iemy1MYbJL2gTjvOlfTxdHskvXcfTPdfmLZD+fVL+i6wM/CT1L5P5BZ5eHqv1kj6dJ11HgUcDnwiPf8nuYfnSLpZ0sOSLpK0Ve55b5G0JG3zryTtUWf5vwGen2vflpJmpvfjofT+fCA3/4mSfijpe5IeAY6WtE7Sjunxz0haL2n7dP+Lkk5Ptw+UtFjSI2n/PTG33Mpn4f2Sfgtcmaa/W9K9ad+r+RoVUb1PpGl1j+Ak7Z1et7WSlkraN/fY1ZJOlnQd8ER6/aqff4+k49P787ik70j6M0k/S/vZLyU9q9n6JJ0MvAY4M70/Z+ZW8wZlab/fp31Y6TnPSO/Dvcq+I86TtENuXR15TSckIobyD7gHeANwMfDFNO1I4Op0ezYQwDNzz7kaODLdPgJYD/w9MAX4IvBb4BvAlsCbgEeBbdP8pwOXAs8GtgN+ApySHts3LevL6bnT0rQV6fEpwFLgq8A2wFbAq9NjLwTemJ43A1gAnF69nXVeg/uA16TbzwL2TLdfD6wB9kzL/TqwIPe8AF5Y/ZrkXpdra82b29bKdk0FlgH/CGyR1vsosFt6/BzgIWAvsqO684EL62zL+4CfpNt/C/wGuCj32H9Wr7/W65N737+d3oe/BP4IvLjOes8h7T9Vy7wRmJne7zuAo9NjewIPAK9I7+t70/xbNtpPc/fnA99M+8AcYDXw1+mxE4GngHlkP/impf3h7enxn6fX5c3p/gLgbbnX5aXpeXsA9wPzql6T88j2v2nA7sBjwGvJ9pGvkO3D9fa1zV6nWvtEre1O2/W9dHsEeBA4ILX1jen+jNz++FvgL8j2mal1XtPrgT9Ly3sAWATMTdtyJfC5FtZ3ZNXyA/gpMJ3sR8dqYP/cvriMLFhtS/b98930WEuvaVl/w36kAPBZ4MOSZrTx3Lsj4t8j4mngImAWcFJE/DEifg48Cbww/Ur4APCxiHgoIh4FvgQcllvWBrId8Y8Rsa5qPXuRfcEcHxGPR8QfIuJagIhYFhG/SM9bTbYjva5g+58Cdpe0fUT8PiIWpemHA2dHxKKI+CPwKeCVkma38uIUsDfZB+PUiHgyIq4k+zC9KzfPxRFxY0SsJwsKc+osaz7wGmXpgtcC/wTskx57XXq8FZ+PiHURsZQsILea6vhaRKyKiIfIfgBU2v0B4F8j4oaIeDoiziULOns3W6CkWcCrgU+mfWAJ8G/Au3Oz/U9EXBIRG9J+NB94nbIj3j2Ar6X7WwEvB64BiIirI+KW9LybgQvYfD86Me1/64B3AD+NiAVpH/m/ZPtwI8elX9uVv3ObbXMNfwdcHhGXp7b+Ahgj+9KuOCcibouI9RHxVJ3lfD0i7o+IlWSvwQ0RsThty4/JAkTR9dVyakSsjYjfAlex8f0/HPhKRNwVEY+RfbYOS+9PO69pxw19UIiIW8m+iE5o4+n3526vS8urnrYt2S/4rYGFlQ8E8F9pesXqiPhDnfXMAu5NX4ybkPSnki5Ulv55BPgesGPB9r+dbOe+V9J8Sa9M02cC91ZmSjvvg2S/mjppJrA8IvI7/r1V6/ld7vYTZK/nZiLiN2S/suaQHdL/FFglaTfaCwqF1tvG858HfDz/5Uj2/tZMz1WZCVR+VFRUv17Lq54zn+yX+J7ALcAvyF6PvYFlEbEGQNIrJF2lLD35MHA0m+9H+WXPzN+PiMfJ9pFG/jkipuf+3ttk/lqeB7yz6vV7NbBTnXbWU/05rfW5Lbq+Wuq9/5t8ttLtZ5IdtbTzmnbc0AeF5HNkv+DyH67H0/+tc9Oe0+by15DtaH+R+0DsEBH5L5pGp6tdDuys2gW+U9Jz94iI7cl+2ahIoyLipog4GPhT4BLgB+mhVWQfBgCU9dD6E2BljcU8Tvuv0SpgljYtBu5cZz1FzCf7tbVF+gU4H3gPWWpsSZ3nTPQ0wa0+fzlwctWX49YRcUGB564Cni1pu9y06teruj2/AnYD3gbMj4jb03MOZNNA+X2y9OasiNgB+Bab70f5Zd9HFswAkLQ12T7Sjk32IUlT2PQHU95ysnRL/vXbJiJOrdPOiWq2vlbXtclni+y9WE8WlDr5mrbNQYEsBUOW/vlIbtpqsg/b30maIul9QM0iZ4HlbyDLUX9V0p/CeDF0v4KLuJFshzlV0jaStpJUSY1sR/YLea2kEeD4IguUtIWkwyXtkA6xHwGeTg9/H/h7SXMkbUmW6rohIu6psaglwCGStlbW9fT9VY/fT41iX3ID2RfCJyRNTQW8twIXFtmGGuYDx5DlyiHL936YrMbxdJ3nNGpfEa0+/9tkBeBXKLONsiLvds2eGBHLyb7kT0n7wB5kr/f5DZ7zBLAQ+BAbg8CvgH9g06CwHdlRyB8k7UVWl2nkh8BbJL1a0hZkXWfb/T75X2Cr9DpMBT5DllOv5XvAWyXtlz6XWykrVD+3zXU302x9rb7/FwAfk7SLpG3JPlsXpSxAJ1/TtjkobHQSWREt7wNkX7IPkhWufjWB5X+SrMB0fUrz/JLsF1xT6QvtrWRF5d8CK4BD08OfJ0sNPAxcRla4KurdwD2pPUeTHWUQEf9Nls/8EVkwegGb1j/yvkpWO7kfOJfNv6BOBM5Nh95/U7VdTwIHAW8mO5r6JvCeiPh1C9uQN5/sy60SFK4l+wW6oO4zsiOtz6T2Hddgvnq+Q1aXWSup6TiCiBgj26/OBH5Ptk8c0cL63kVW+F1Flvv+XMpzNzKfrKh/Y+5+/nUC+CBwkqRHyepsP6CBiLiNLNB8n2wf+T3ZftlIpZdW5W9NWtbDaf3/RvZD7PF6y0qB8WCyzgmryX7JH09J32UF1ncG8I7Uy+hrBRZ5NvBdstf+buAPZD9c2n1NO04RvsiOmZllfKRgZmbjHBTMzGycg4KZmY1zUDAzs3F9d+K1HXfcMWbPnt3tZpiZ9ZWFCxeuiYimZ27ou6Awe/ZsxsbGut0MM7O+Iune5nM5fWRmZjkOCmZmNs5BwczMxjkomJnZOAcFMzMb13e9j8yssUsWr+S0K+5k1dp17DBtKhKsfeKpTW7PnD6N4/fbjXlzO32JDOt3fXdCvNHR0XCXVLPaX/6/f+IpRLGT/FfmG3GAGAqSFkbEaNP5HBTM+kclEKxcu67wl38RDhCDz0HBbMBcsngln7r4FtY9Ve96QZ3hADGYuh4U0kXGzyO7POMG4KyIOKNqnn2B/yS72ARkF2k/qdFyHRRs2OSPDiabA8TgKBoUyiw0rwc+HhGL0qUGF0r6RbpGbN41EfGWEtth1nfKShO1qrLelWvX8amLbwFwYBhwpXVJjYj7ImJRuv0ocAfgvcmsiUqaqHJk0E5AUPo/fdpUnrX1VJS7nX+8FeueepqPXrSEfU69kksWr2xjCdYPJqVLqqTZwFyyC7VXe6WkpWTXnD0uXafUbOi0kyaqHEVMb7G76USORHzUMNhKLzRL2pbsQuEnR8TFVY9tD2yIiMckHQCcERG71ljGUcBRADvvvPPL7r230Mn+zPpGO0XkTuX52w0QUyQ2RHjMQ5/oeqE5NWIq8FPgioj4SoH57wFGI2JNvXlcaLZBtM+pVxY+Qpg2dQqnHPLSUr6E2w0QZbbJOqPrQUGSgHOBhyLio3XmeQ5wf0SEpL2AHwLPiwaNclCwQZEffNbsU9iNXkDtpLPcS6l39ULvo32AdwO3SFqSpv0jsDNARHwLeAfwfyStB9YBhzUKCGaDopV0Ube+aOfNHWHe3JGW2up6Q/8rLShExLU06eQQEWcCZ5bVBrNe08qv715JyVTWXzmqeYbE0w1+u6176mlOu+LOrrfb2uMRzWaTpOgvbkFPF2+LbodTSb2lF9JHZpZz2hV3Fvoive6E109Si9qTP3JodMTjVFJ/8vUUzEp2yeKVhXoXTZs6heP3222SWjUx8+aOcN0Jr+f0Q+cwbeqUuvN5wFv/8ZGCWYkGPdXio4bB4yMFsxI1SxlNmzqF0w+dw3UnvL5vvywrRw0j06c1nK9SgLbe5qBgVoIiKaOR6dN6ondRpxy/324NU0mQHTE4ldTbnD4y67AiKaN+KCi3yqmkweAjBbMOK5Iy6peCcqtaKUA7ldSbHBTMOmzVEKWM6pk3d4RTDnlpwzqDU0m9yekjsw6pjFauNxx0EFNGjVROk9GotuJUUu/xkYJZB1RfGKfaIKeMmmlWgPZYht7iIwWzDmhUR+jXMQid4gJ0f/GRgtkENOt6KujrMQid4rEM/cNBwaxNzVJGkJ3YzjbyWIbe5/SRWZuGuetpu5xK6n0+UjBrk7uetsdjGXqbjxTMWuSup51R5KihUeC1cvhIwawF7nraWc0K0AGuL0wyBwWzFjTreuqUUXsaFaAr9QUHhsnh9JFZC+qlMypdT609zVJJvu7z5PGRglkBlfEI9eoI7no6cZVUkuo87q6qk8NBwawJ1xEmV6MA61RS+RwUzJpwHWFyFTlXkruqlsdBwayJZnUEB4TO8mm3u8tBwawO1xG6p8i5kpxKKoeDglkNriP0BqeSJp+DglkNriP0BqeSJp+DglkNriP0DqeSJpeDglmO6wi9y6mkyeGgYJa4jtDbnEqaHKUFBUmzJF0l6Q5Jt0k6tsY8kvQ1Scsk3Sxpz7LaY9aM6wi9z6mk8pV5pLAe+HhEvBjYG/iQpN2r5nkzsGv6Owr4lxLbY9aQ6wj9w6mk8pQWFCLivohYlG4/CtwBVH+qDgbOi8z1wHRJO5XVJrNaXEfoP0VTSbuccJnTSS2alJqCpNnAXOCGqodGgOW5+yvYPHAg6ShJY5LGVq9eXVYzbQi5jtC/iqSSAqeTWlV6UJC0LfAj4KMR8Uj1wzWestkPtog4KyJGI2J0xowZZTTThpTrCP2vWSoJnE5qRanXU5A0lSwgnB8RF9eYZQUwK3f/ucCqMttklufrI/S//LUYVq1dVzcN6Et7FlNm7yMB3wHuiIiv1JntUuA9qRfS3sDDEXFfWW0yq1avXuA6Qn+ppJLuPvVAX9pzgspMH+0DvBt4vaQl6e8ASUdLOjrNczlwF7AM+DbwwRLbYzauUlxeuXbdZjlM1xH6my/tOTGlpY8i4lpq1wzy8wTwobLaYFZLpbhcqSUE2Y4aZHWE4/fbzXWEPuZLe06MRzTb0KlVXK4EBI9HGAy+tGf7HBRs6NQrOLoQOXh8ac/WOSjY0PAgteHjkc+tc1CwoeBBasPJJ9FrnYOCDQUPUhtePoleaxwUbCj4ZHfmVFIxDgo2FDxIzZxKKqbU01yYddMli1eOn/pgh2lTmTpFPPX0xjKz6wjDZ97cEebNHRkfuFhLJZVUmX/Y+EjBBlK+sBzA2nVPQcCztp6KcB1h2DmVVJ+PFGwg1SosP7Uh2HqLZ7L4s2/qUqusVzQb9QwbU0nDNsLdRwo2kDxAzZpxr6TaHBRsoHiAmrXKqaRNOSjYwPAANWuHeyVtykHBBoYHqFm7nErayEHBBoYHqNlEOZXkoGADwHUE6xSnkhwUrM+5jmCdNuypJAcF62uuI1hZiqSSPnrRkoE7avDgNetrzeoIZu0qMsANBu+0GD5SsL7mE91ZmYqkkmCwCtAOCtaXKsXllWvXbXYdXtcRrNOapZJgcArQTh9Z36kUlyu1hCBLFwVZHWHYzlVj5RumVJIi6nXk602jo6MxNjbW7WZYF9U77fHI9GmuI1jpqn+U1NNrP1AkLYyI0Wbz+UjB+o5PdmfdNOhHDa4pWN/wIDXrFYNcgHZQsL7gQWrWiwaxAO30kfWFZoPUeil3a8NjEFNJLjRbX9jlhMtqpo0E3H3qgZPdHLPN9HoB2oVmGwiXLF7JaVfc6TqC9bxBOWpwTcF6lusI1m9aKUD36nmTSgsKks6W9ICkW+s8vq+khyUtSX+fLast1p98sjvrV0UK0NCbZ1stlD6S9OyIeKjFZZ8DnAmc12CeayLiLS0u14aET3Zn/apoKgk2dlvtlR84RY8UbpD0H5IOkFR9qpmaImIB0GogMfN4BBsIlVTS6YfOKdRtdZcTLuuJdFLRoPAi4Czg3cAySV+S9KIOrP+VkpZK+pmkv6g3k6SjJI1JGlu9enUHVmu9ynUEGzRFruYG2bm7eiGd1HKXVEl/BXwP2AZYCpwQEf9TZ97ZwE8j4iU1Htse2BARj0k6ADgjInZttn53SR1s9c5rBB6PYP2vaLdV6Pz+3tEuqZL+BPg7siOF+4EPA5cCc4D/AHZptYER8Uju9uWSvilpx4hY0+qybHC4jmCDLF9rWLV2Xd0UKXSv62rR9NH/ANsD8yLiwIi4OCLWR8QY8K12VizpOZX6hKS9UlsebGdZ1v9cR7BhUak13H3qgT3ZdbVoUPhMRHwhIlZUJkh6J0BEfLnWEyRdQBZMdpO0QtL7JR0t6eg0yzuAWyUtBb4GHBb9NrzaOsJ1BBtWrXRd/dhFS5g9CcXoQjUFSYsiYs9m0yaDawqDx3UEG2aVUfvNuq7mTZs6peVxOh2pKUh6M3AAMCLpa7mHtgfWF26NWQOuI9gwmzd3hHlzR1oqQpc5tqFZoXkVMAYcBCzMTX8U+FjHW2NDxec1MtuolQFvUN5FpRoGhYhYCiyVdH5E+MjAOqbZryLXEWwYtXLUUNaPpmbpox9ExN8AiyXlf9AJiIjYo5RW2cDz9RHM6qs+ahBsckRd5o+mZumjY9N/n5/IOsp1BLPGKkcNsDHVumrtOmaW/KOpWfrovnRzDbAuIjak01v8OfCzUlpkA811BLPW5QNE2YqOU1gAbCVpBPhv4O/JzoJqVpjHI5j1vqJBQRHxBHAI8PWIeBuwe3nNskHk6yOY9b6il+OUpFcChwPvb/G5ZoDrCGb9oOgX+7HAp4AfR8Rtkp4PXFVes2yQuI5g1j8KBYV0wZwFuft3AR8pq1E2ODweway/FD119ouA44DZ+edEhI/5rSGPRzDrL0XTR/9BdorsfwOan5jDhl6zk3y5jmDWm4oGhfUR8S+ltsQGRjeH6JvZxBTtkvoTSR+UtJOkZ1f+Sm2Z9a1GKSNwHcGslxU9Unhv+n98bloAz+9sc2wQNDp7o+sIZr2taO+jlq/BbMOnWdfTkenTXEcw63GF0keStpb0GUlnpfu7SvJJ8mycT2FhNhiK1hT+HXgSeFW6vwL4Yiktsr7kU1iYDYaiNYUXRMShkt4FEBHrJKnEdlmfcNdTs8FSNCg8KWka6ToPkl4A/LG0VllfcNdTs8FTNCicCPwXMEvS+cA+ZKfPtiHmrqdmg6do76OfS1oI7E2WETg2ItaU2jLrSfkrQNXrZQTuemrWr4qe++i/I+KvgctqTLMhUSRdBO56atbPGgYFSVsBWwM7SnoW2VECwPbAzJLbZj2mWboInDIy63fNjhT+AfgoWQBYyMag8AjwjRLbZT2o0UhlQekXFDez8jUMChFxBnCGpA9HxNcnqU3WYzxS2Wx4FC00f13Sq9j8egrnldQu6xG+SI7ZcClaaP4u8AJgCRuvpxCAg8KA80VyzIZL0XEKo8DuEdGoF+ImJJ0NvAV4ICJeUuNxAWcABwBPAEdExKKiy7dyeaSy2XAqeu6jW4HntLjsc4D9Gzz+ZmDX9HcU4Iv49IhmJ7cDj1Q2G1RFjxR2BG6XdCO501tExEH1nhARCyTNbrDMg4Hz0tHH9ZKmS9opIu4r2CYriUcqmw2vVk5z0WkjwPLc/RVpmoNClzRLGYHrCGaDrmjvo/klrLvWWVZr1iwkHUWWYmLnnXcuoSlWZLSyu56aDb5mI5ofpfYXtYCIiO0nsO4VwKzc/ecCq2rNGBFnAWcBjI6OFi52W3NFjg7AKSOzYdFs8Np2Ja77UuAYSRcCrwAedj1hcrVyLiOnjMyGQ9GaQsskXQDsS3bepBXA54CpABHxLeBysu6oy8i6pPpU3JOsyLmMnDIyGy6lBYWIeFeTxwP4UFnrt/qcMjKzekoLCtabnDIys0YcFIZMkTEIpxzyUgcDsyHloDAkPAbBzIpwUBgCHoNgZkU5KAwwF5TNrFUOCgPKBWUza4eDwoDyGAQza4eDwoBxysjMJsJBYYA4ZWRmE+WgMEA8BsHMJspBYQB4DIKZdYqDQp/zGAQz6yQHhT7lgrKZlcFBoQ+5oGxmZXFQ6EMeg2BmZXFQ6CNOGZlZ2RwU+oRTRmY2GRwUelwrRwceg2BmE+Wg0MN8dGBmk81BoYe5oGxmk81BocdU0kWr1q4jmszrgrKZdZqDQg8pmi4Cp4zMrBwOCj2kSLrIBWUzK5ODQg8o0sNIwEwfHZhZyRwUuswntDOzXuKg0CUenWxmvchBoQs8/sDMepWDQhd4/IGZ9SoHhUnklJGZ9ToHhUnilJGZ9YNnlLlwSftLulPSMkkn1Hj8CEmrJS1Jf0eW2Z5uapYymjZ1CqcfOofrTni9A4KZdU1pRwqSpgDfAN4IrABuknRpRNxeNetFEXFMWe3otiIpIx8dmFmvKDN9tBewLCLuApB0IXAwUB0UBpbHIJhZvykzfTQCLM/dX5GmVXu7pJsl/VDSrFoLknSUpDFJY6tXry6jraUokjJyQdnMekmZQUE1plWf+PMnwOyI2AP4JXBurQVFxFkRMRoRozNmzOhwMzvvksUr2efUK5umjHwOIzPrNWWmj1YA+V/+zwVW5WeIiAdzd78NfLnE9kwKp4zMrJ+VeaRwE7CrpF0kbQEcBlyan0HSTrm7BwF3lNieSeGUkZn1s9KOFCJivaRjgCuAKcDZEXGbpJOAsYi4FPiIpIOA9cBDwBFltads7mVkZoNAEc2u79VbRkdHY2xsrNvN2IRTRmbW6yQtjIjRZvOVOnhtWDhlZGaDwqe56IBVThmZ2YBwUJiASh2hXgLOKSMz6zcOCm1qVkdwysjM+pGDQpsa1RGcMjKzfuWg0KJmXU8FThmZWd9yUGhBka6nM6dPm8QWmZl1lruktsBdT81s0PlIoQXuempmg85BoQB3PTWzYeGg0IS7nprZMHFQaMJdT81smDgo1OGup2Y2jBwUanDXUzMbVu6SWoO7nprZsPKRQg3uempmw8pBIcddT81s2DkoJO56ambmoDDOXU/NzBwUxtWrI7jrqZkNE/c+Sup1MXXXUzMbJkN/pJAfpCbYpMjsOoKZDZuhDgrVxeWA8cDgOoKZDaOhDgq1isuVgOA6gpkNo6GuKdQrLjcavGZmNsiG7kihUkNYtXYdz5B4OjYfqubispkNq6EKCtU1hFoBwcVlMxtmQxUU6g1QmyKxIYKZLi6b2ZAbqqBQr1awIYK7Tz1wkltjZtZ7hqrQ7AFqZmaNlRoUJO0v6U5JyySdUOPxLSVdlB6/QdLsMtpxyeKV7HPqleMD1PJcQzAz26i0oCBpCvAN4M3A7sC7JO1eNdv7gd9HxAuBrwJf7nQ7KsXlymU1KwPUIBuPcMohL3UNwcwsKbOmsBewLCLuApB0IXAwcHtunoOBE9PtHwJnSlJEjW5BbfIANTOz4spMH40Ay3P3V6RpNeeJiPXAw8CfVC9I0lGSxiSNrV69uqVGeICamVlxZQaF6vQ9sNlFzYrMQ0ScFRGjETE6Y8aMlhrh4rKZWXFlBoUVwKzc/ecCq+rNI+mZwA7AQ51sxPH77ca0qVM2mebisplZbWUGhZuAXSXtImkL4DDg0qp5LgXem26/A7iyk/UEgHlzRzjlkJcyMn0awsVlM7NGSis0R8R6SccAVwBTgLMj4jZJJwFjEXEp8B3gu5KWkR0hHFZGW+bNHXEQMDMroNQRzRFxOXB51bTP5m7/AXhnmW0wM7PihmpEs5mZNeagYGZm4xwUzMxsnIOCmZmNU4d7gJZO0mrg3jafviOwpoPN6SZvS28alG0ZlO0Ab0vF8yKi6ejfvgsKEyFpLCJGu92OTvC29KZB2ZZB2Q7wtrTK6SMzMxvnoGBmZuOGLSic1e0GdJC3pTcNyrYMynaAt6UlQ1VTMDOzxobtSMHMzBpwUDAzs3FDFxQkfUHSzZKWSPq5pJndblO7JJ0m6ddpe34saXq329QuSe+UdJukDZL6rvugpP0l3SlpmaQTut2edkk6W9IDkm7tdlsmStIsSVdJuiPtW8d2u03tkLSVpBslLU3b8flS1zdsNQVJ20fEI+n2R4DdI+LoLjerLZLeRHYNivWSvgwQEZ/scrPaIunFwAbgX4HjImKsy00qTNIU4H+BN5JdOOom4F0RcXvDJ/YgSa8FHgPOi4iXdLs9EyFpJ2CniFgkaTtgITCv394XSQK2iYjHJE0FrgWOjYjry1jf0B0pVAKSxeTRAAAEyUlEQVRCsg01Lv/ZLyLi5+na1gDXk13dri9FxB0RcWe329GmvYBlEXFXRDwJXAgc3OU2tSUiFtDhqx92S0TcFxGL0u1HgTvY/DrxPS8yj6W7U9Nfad9bQxcUACSdLGk5cDjw2Wbz94n3AT/rdiOG1AiwPHd/BX345TPIJM0G5gI3dLcl7ZE0RdIS4AHgFxFR2nYMZFCQ9EtJt9b4OxggIj4dEbOA84FjutvaxpptS5rn08B6su3pWUW2pU+pxrS+PQIdNJK2BX4EfLQqU9A3IuLpiJhDlg3YS1Jpqb1Sr7zWLRHxhoKzfh+4DPhcic2ZkGbbIum9wFuAv+709a07rYX3pd+sAGbl7j8XWNWltlhOysH/CDg/Ii7udnsmKiLWSroa2B8opTPAQB4pNCJp19zdg4Bfd6stEyVpf+CTwEER8US32zPEbgJ2lbSLpC3IrjV+aZfbNPRSgfY7wB0R8ZVut6ddkmZUehZKmga8gRK/t4ax99GPgN3IerrcCxwdESu726r2SFoGbAk8mCZd38c9qd4GfB2YAawFlkTEft1tVXGSDgBOB6YAZ0fEyV1uUlskXQDsS3aK5vuBz0XEd7raqDZJejVwDXAL2ecd4B/TteP7hqQ9gHPJ9q1nAD+IiJNKW9+wBQUzM6tv6NJHZmZWn4OCmZmNc1AwM7NxDgpmZjbOQcHMzMY5KNjAkfQcSRdK+o2k2yVdLulFJaznREnHpdsnSWprcJ6kOalLq1nXDeSIZhteacDSj4FzI+KwNG0O8GdkZzJtd7lTIuLpeo9HxETOoTUHGAX6qv+8DSYfKdig+SvgqYj4VmVCRCyJiGuUOS2db+kWSYdCFkjqTN83nY//+2QDoJD06XTdhF+SDYIkTT9H0jvS7XskfV7SorS8P0/T95L0K0mL0//d0gjok4BDlV3j41BJ2yi7rsFNad7Nzg0laSdJC9JzbpX0mtJeURsqPlKwQfMSsvPm13II2a/yvyQbsXuTpAXAq+pMh+y02C+JiLslvYzsFBZzyT47ixqsa01E7Cnpg8BxwJFkpyZ4bbr+xRuAL0XE2yV9FhiNiGMAJH2J7DoZ70unN7hR0i8j4vHc8v8WuCIiTlZ2PYetW3uZzGpzULBh8mrggpQGul/SfODlDaY/AtwYEXen578G+HHlPFOSGp3fqHLytYVkwQhgB+DcdP6tIDsvfi1vAg6q1CuArYCdya4HUHETcHY64dslEbGk+eabNef0kQ2a24CX1Xms1imuG00HeLzqftHzwvwx/X+ajT++vgBcla5o9layL/t67Xl7RMxJfztHRD4gVC6G81pgJfBdSe8p2C6zhhwUbNBcCWwp6QOVCZJeLul1wAKy3P0USTPIvlRvbDC92gLgbZKmKbu841tbbNsOZF/iAEfkpj8KbJe7fwXw4VQ0R9Lc6gVJeh7wQER8m+xMoHu22BazmhwUbKCka0q8DXhj6pJ6G3Ai2fUNfgzcDCwlCx6fiIjfNZhevexFwEXAErJz9F/TYvP+CThF0nVkZ7ysuArYvVJoJjuimArcLOnWdL/avsASSYuBtwNntNgWs5p8llQzMxvnIwUzMxvnoGBmZuMcFMzMbJyDgpmZjXNQMDOzcQ4KZmY2zkHBzMzG/X9i2MOlBYY+uAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s, I, 'o')\n",
    "plt.xlabel(\"Coordinate s\")\n",
    "plt.ylabel(\"Intensity\")\n",
    "plt.title(\"Numerical solution with the forward Euler method\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 3\n",
    "\n",
    "Now set **A = 10**. Note that this affects only the extinction coefficient alpha."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Numerical solution with the forward Euler method')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xu8XGV97/HPlxBgc41KWsk2GKqYU1RMcBcveOEgFkSFiFiwVKGKHA5FsS9BQ/UoUhGUvrxia6FQEBGwijlRsIgCCWC57Fy4l56IYi4UwiVct5CE3/ljPTOsTOa2Z8/ac/u+X6/92jNr1qz1zJo167ee57ee9SgiMDMzA9ii0wUwM7Pu4aBgZmZlDgpmZlbmoGBmZmUOCmZmVuagYGZmZQ4KXUrSWyXd24bl/E7S/u0oU26ZIemVLb63LZ+rVY3WL2lW+nxbTmAdbdvmkvaR9P8kPSVpXjuWWZQi9rW03FMlfb/dyy3KRH4fVZZVyDatZ2CDQtrYD0raLjftGEnXdbBYZRFxfUTM7nQ5JqryB9Lpz1W5/on+6CRdIOlL7SldVacBZ0fE9hGxoMD1FCptp+dScCv93dbpck2UpOskHdPpcrTTwAaFZEvgxE4XotJEzlKt77wcuKuVNxa1H01guV9Nwa3097q2FqyCMoN+jBu3Qd9gZwEnSZpW+UK1ZoT8WYGkoyXdKOnrktZJuk/Sm9P0lZIeknRU7r1bS/oHSb9PNZTvShpKr+0raZWkz0j6b+BfS9Ny758p6XJJayU9IunsNP0Vkq5J0x6WdHG1z1ONpIMk3S3pSUmrJZ2Ue+1jklZIelTSQkkzaixjkzOl9PlvSI8Xp8m3pTPDw6t8rj9Ny1gn6S5JB+deu0DSdyRdkcp4s6RX1CjHhZI+lR4Pp+/u+PT8lelzKL9+SRcBuwI/TeX7dG6RR6bv6mFJn62xzmOBI4FPp/f/NPfyHEm3S3pc0mWStsm97z2SlqfP/GtJe9ZY/m+AP8mVb2tJM9L38Wj6fj6Wm/9UST+S9H1JTwDHSRqTtHN6/XOSNkjaMT3/kqRvpMfvlrRM0hNp/z01t9zSb+Gjkn4PXJOmf0jS/Wnfq7qNmlG5T6RpNWtwkt6Ytts6SbdJ2jf32nWSTpd0I/BM2n6V7/+dpJPT9/O0pPMk/bGkn6f97JeSXtRofZJOB94KnJ2+n7Nzq9lfWbPfY2kfVnrPFul7uF/ZMeJ7knbKrast23RCImIg/4DfAfsDlwNfStOOAa5Lj2cBAWyZe891wDHp8dHABuCvgSnAl4DfA98Btgb+HHgS2D7N/w1gIfBiYAfgp8AZ6bV907K+kt47lKatSq9PAW4Dvg5sB2wDvCW99krgnel904HFwDcqP2eNbfAA8Nb0+EXAXunxfsDDwF5pud8GFufeF8ArK7dJbrvcUG3e3Gctfa6pwArg74Ct0nqfBGan1y8AHgX2JqvVXQxcWuOzfAT4aXr8l8BvgMtyr/3fyvVX2z657/3c9D28DngW+NMa672AtP9ULPMWYEb6vu8Bjkuv7QU8BLwhfa9Hpfm3rref5p4vAv4x7QNzgLXAO9JrpwLrgXlkJ3xDaX94f3r9F2m7vCs9Xwy8L7ddXpvetyfwIDCvYpt8j2z/GwL2AJ4C3ka2j3yNbB+uta9ttp2q7RPVPnf6XN9Pj4eBR4CDUlnfmZ5Pz+2PvwdeTbbPTK2xTW8C/jgt7yFgKTA3fZZrgC+MY33HVCw/gJ8B08hOOtYCB+b2xRVkwWp7suPPRem1cW3Tov4GvaYA8Hng45Kmt/De30bEv0bERuAyYCZwWkQ8GxG/AJ4DXpnOEj4G/G1EPBoRTwJfBo7ILet5sh3x2YgYq1jP3mQHmJMj4umI+ENE3AAQESsi4ur0vrVkO9Lbmyz/emAPSTtGxGMRsTRNPxI4PyKWRsSzwCnAmyTNGs/GacIbyX4YZ0bEcxFxDdmP6YO5eS6PiFsiYgNZUJhTY1mLgLcqay54G/BVYJ/02tvT6+PxxYgYi4jbyALyeJs6vhURayLiUbITgFK5Pwb8c0TcHBEbI+JCsqDzxkYLlDQTeAvwmbQPLAf+BfhQbrb/iIgFEfF82o8WAW9XVuPdE/hWer4N8GfA9QARcV1E3JHedztwCZvvR6em/W8MOAz4WUQsTvvI/yHbh+s5KZ1tl/4ubPSZq/gr4MqIuDKV9WpglOygXXJBRNwVERsiYn2N5Xw7Ih6MiNVk2+DmiFiWPstPyAJEs+ur5syIWBcRvweu5YXv/0jgaxFxX0Q8RfbbOiJ9P61s07Yb+KAQEXeSHYjmt/D2B3OPx9LyKqdtT3YGvy2wpPSDAP49TS9ZGxF/qLGemcD96cC4CUl/JOlSZc0/TwDfB3ZusvzvJ9u575e0SNKb0vQZwP2lmdLO+wjZWVM7zQBWRkR+x7+/Yj3/nXv8DNn23ExE/IbsLGsOWZX+Z8AaSbNpLSg0td4W3v9y4FP5gyPZ91u1ea7CDKB0UlFSub1WVrxnEdmZ+F7AHcDVZNvjjcCKiHgYQNIbJF2rrHnyceA4Nt+P8suekX8eEU+T7SP1/ENETMv9HdVg/mpeDnygYvu9BdilRjlrqfydVvvdNru+amp9/5v8ttLjLclqLa1s07Yb+KCQfIHsDC7/43o6/d82N+2lLS7/YbId7dW5H8ROEZE/0NS7Xe1KYFdVT/Cdkd67Z0TsSHZmo2YKFRG3RsQhwB8BC4AfppfWkP0YAFB2hdZLgNVVFvM0rW+jNcBMbZoM3LXGepqxiOxsa6t0BrgI+DBZ09jyGu+Z6G2Cx/v+lcDpFQfHbSPikibeuwZ4saQdctMqt1dleX4NzAbeByyKiLvTe97NpoHyB2TNmzMjYifgu2y+H+WX/QBZMANA0rZk+0grNtmHJE1h0xOmvJVkzS357bddRJxZo5wT1Wh9413XJr8tsu9iA1lQauc2bZmDAlkTDFnzzydy09aS/dj+StIUSR8BqiY5m1j+82Rt1F+X9EdQToYe0OQibiHbYc6UtJ2kbSSVmkZ2IDtDXidpGDi5mQVK2krSkZJ2SlXsJ4CN6eUfAH8taY6krcmaum6OiN9VWdRy4FBJ2yq79PSjFa8/SJVkX3Iz2QHh05KmpgTee4FLm/kMVSwCTiBrK4esvffjZDmOjTXeU698zRjv+88lSwC/QZntlCV5d2j0xohYSXaQPyPtA3uSbe+L67znGWAJ8De8EAR+DfwvNg0KO5DVQv4gaW+yvEw9PwLeI+ktkrYiu3S21ePJfwHbpO0wFfgcWZt6Nd8H3ivpgPS73EZZovplLa67kUbrG+/3fwnwt5J2k7Q92W/rstQK0M5t2jIHhRecRpZEy/sY2UH2EbLE1a8nsPzPkCWYbkrNPL8kO4NrKB3Q3kuWVP49sAo4PL38RbKmgceBK8gSV836EPC7VJ7jyGoZRMSvyNozf0wWjF7BpvmPvK+T5U4eBC5k8wPUqcCFqer9FxWf6zngYOBdZLWpfwQ+HBH/OY7PkLeI7OBWCgo3kJ2BLq75jqym9blUvpPqzFfLeWR5mXWSGvYjiIhRsv3qbOAxsn3i6HGs74Nkid81ZG3fX0jt3PUsIkvq35J7nt9OAMcDp0l6kizP9kPqiIi7yALND8j2kcfI9st6Sldplf4eTst6PK3/X8hOxJ6utawUGA8huzhhLdmZ/MkUdCxrYn3fBA5LVxl9q4lFng9cRLbtfwv8gezEpdVt2naK8CA7ZmaWcU3BzMzKHBTMzKzMQcHMzMocFMzMrKznbry28847x6xZszpdDDOznrJkyZKHI6LhnRt6LijMmjWL0dHRThfDzKynSLq/8VxuPjIzsxwHBTMzK3NQMDOzMgcFMzMrc1AwM7Oywq4+SoN4LCa72+GWwI8i4gsV82xNNprT68luOnd4jTtxmvWUBctWc9ZV97Jm3Rg7DU1FgnXPrGfGtCFOPmA28+a2e2gKs/Yo7IZ4abSx7SLiqXQ73BuAEyPiptw8x5ONA3CcpCPIhgY8vMYiARgZGQlfkmrdotrB/7Fn1iNq32i/9NqwA4RNIklLImKk0XyFNR9F5qn0dGr6q/ydHEJ2u2XI7iX+jhRMzLregmWrOeXyO1i9bowA1o2t57FnstEf651qlV5bvW6Mv71sObPmX8E+Z17DgmWtji1k1j6Fdl5LIygtIRsH4DsRcXPFLMOk4eciYkMaBvAlZPfWzy/nWOBYgF133bXIIps1VKodrF5XOZT2+OUDxCmX3wHgmoN11KSMpyBpGtmAIB9PYyKXpt8FHBARq9Lz3wB7R0TNcUndfGSdkA8E9ZqGJmqKxPMRzj1Y23W8+SgvItaRDY14YMVLq0hjkqbxh3cCHp2MMpk1K99MBMUFBICNEQQv1BzcpGSTrbCgIGl6qiEgaQjYH6gcZnEhcFR6fBhwTXgoOOsyZ111L2Praw3xXFspOTZtaCov2nbqJtOaMbZ+I5+8bLnzDTapiswp7EI2Nu8UsuDzw4j4maTTgNGIWEg2vu1FklaQ1RBqjQNsNunGmzuY1sSlp600QznfYJOp58Zodk7BJkOpyaiZGsLQ1Cmccehrx33Azl/OuoXExiZ+i76M1VrVbE6h526dbVakZmsH7ehrMG/ucPl9zQYh1xqsaA4KZkmzB+YiztZLy2omII2t38hZV93roGCF8L2PzJJmEsrD04a4cf5+hRyQ580d5sb5+/GNw+cwNHVK3XlXrxtzAtoK4ZqCDbxmm4yGpk7h5ANmF16eZmsNbkqyIrimYAOtsg9CLcPThlpKJreq2VpDqSnJrF1cU7CB1qjJqNUri9qlmVpDqSnJVyVZO7imYANpwbLV7HPmNXVrCJNdO6ilVGsYnjZUcx73gLZ2cVCwgdNMk1GRCeVWnXzAbDclWeHcfGQDp5kmo8lIKI9XM01Ja9pw51YbbK4p2MCpd+DsliajWho1JQX4UlWbEAcFGxilPEKtm0l0Y5NRLfWakpxfsIlwULCB0CiP0K1NRrXMmzvMGYe+tmaNwfkFa5WDgg2EenmEbm8yqqXUlFTrdtzOL1grHBRsINQ6QAp6psmolhnOL1gbOShY3yrlEHabfwVbqPr5dK0Dai9xfsHayUHB+lI+hxBQdayCXssj1OL8grWTg4L1pVo5hCkSonfzCLU4v2Dt4s5r1pdqHQSfj+C3Z757kkszeWZMG6p6hVUpv+D7I1kjrilYX2nUF6Efcgj1OL9gE+WgYH2j3/oitML5BZsoBwXrG/3YF6EVzi/YRDgoWN/o574IrXD/BWuFg4L1vEHPI9Ti/IK1wkHBeprzCLU5v2CtKCwoSJop6VpJ90i6S9KJVebZV9Ljkpanv88XVR7rT84j1Of8go1Xkf0UNgCfioilknYAlki6OiLurpjv+oh4T4HlsD7WKI9gmVr9Fwa1ac1qK6ymEBEPRMTS9PhJ4B5gcE/ZrBC1Dmo+2G2qWn5BZLkFJ50tb1JyCpJmAXOBm6u8/CZJt0n6uaRX13j/sZJGJY2uXbu2wJJarygll1evG9usaWSQ8wi1VOYXBOXEvJPOlqeocqOwtq5A2h5YBJweEZdXvLYj8HxEPCXpIOCbEbF7veWNjIzE6OhocQW2rldKLudzCaWD3PC0Id/KoYFSMK1UGnnO+pOkJREx0mi+Qu99JGkq8GPg4sqAABART+QeXynpHyXtHBEPF1ku623VksulgOCDWmO18jBOOhsUe/WRgPOAeyLiazXmeWmaD0l7p/I8UlSZrD/4oDYx7tRm9RSZU9gH+BCwX+6S04MkHSfpuDTPYcCdkm4DvgUcEUW3Z1nPc3J5YtypzeoprPkoIm6AmpdHl+Y5Gzi7qDJYf1mwbDVnXXVvObmcP3twcrl5pXxLaVtWKnVqc15mMLlHs/WEyp7LwQtnHO6kNn7u1Ga1eJAd6wlOLhfDndqskmsK1hOcXC6GO7VZJdcUrKuV8gi+A2oxKvML1Tq15eez/ueagnUt3wF1cpTyC8PThjYLvr6T6uBxTcG6VqM7oLrncnu5ic7AQcG6mO+AOrmcdDZw85F1MXdSm1xOOhu4pmBdyJ3UOsNJZwPXFKzLuJNaZznpbK4pWFdxJ7Xu4KTz4HJNwbqKD0bdwfmcweWgYF3FB6Pu4KTz4HLzkXUFJ5e7i5POg8s1Bes4J5e7k5POg8k1Bes4J5e7m/M8g8U1Bes4H3S6W618zhYSu82/wjmGPuOgYB3n5HJ3qzV858YIAg/h2W8cFKxjFixbzT5nXlNOZOY5udw95s0d5oxDX8vwtCEETNHm47U5x9A/nFOwjigll0u5hFJyuZRL8B1Qu8u8ucPl72O3+VdUncfNff3BQcE6wsnl3uW7qfY3Nx9ZRzi53Lvcsa2/OShYRzi53LvyOQagasc2B4be5aBgk8rJ5f7gjm39q7CgIGmmpGsl3SPpLkknVplHkr4laYWk2yXtVVR5rPPcc7n/uBmw/xSZaN4AfCoilkraAVgi6eqIuDs3z7uA3dPfG4B/Sv+tDzm53H+cdO4/hdUUIuKBiFiaHj8J3ANUngYeAnwvMjcB0yTtUlSZrLN8Vtl/nHTuP5OSU5A0C5gL3Fzx0jCwMvd8FZsHDiQdK2lU0ujatWuLKqYVpJRHqGx7LvFZZe9y0rn/FB4UJG0P/Bj4ZEQ8UflylbdsduyIiHMiYiQiRqZPn15EMa0glXmESk4u9z4nnftLoZ3XJE0lCwgXR8TlVWZZBczMPX8ZsKbIMtnkqpZHKHHP5f7i5sH+UOTVRwLOA+6JiK/VmG0h8OF0FdIbgccj4oGiymSTr9YBQcCN8/dzQOgjtZoBA5xf6CFFNh/tA3wI2E/S8vR3kKTjJB2X5rkSuA9YAZwLHF9geawD3EltcNS6myo4v9BLCms+iogbqJ4zyM8TwN8UVQbrHA+vOXgqh/CsVMovuHbY3dyj2drOndQGVynpXOts0PmF7ue7pFrbuZOauVNb73JNwdrOV6GYO7X1LtcUrG1KeQR3UrPK/EK1Tm35+ax7uKZgbeFOalbJndp6k2sK1hbupGa1uDmxt7imYG3hTmpWizu19RYHBWsLd1KzWtyprbc0FRQkvbjoglhv8khq1kjlnVQrOb/QXZqtKdws6d/SbSrq9lK2weFOatYsd2rrHc0GhVcB55Ddy2iFpC9LelVxxbJe0KiTmgOCVXJ+ofs1FRTSyGhXR8QHgWOAo4BbJC2S9KZCS2hdy1eV2Hg5v9D9mrokVdJLgL8iqyk8CHyc7LbXc4B/A3YrqoDWXUod1NasG2MLiY2xeVc1J5etFt80r/s123z0H8COwLyIeHdEXB4RGyJiFPhuccWzbpLPIQRUDQhOLlsjzi90t2aDwuci4u8jYlVpgqQPAETEVwopmXWdWh3UpkgIJ5dtfJxf6E7NBoX5Vaad0s6CWPerdQb3fAS/PfPdTi7buDi/0J3q5hQkvQs4CBiW9K3cSzsCG4osmHUP3+jOiuD8QndqVFNYA4wCfwCW5P4WAgcUWzTrBr7RnRXJ+YXuU7emEBG3AbdJujgiXDMYQL7RnU2GWoPylPIL3s8mT6Pmox9GxF8AyyTlWw9E1n1hz0JLZx2TH2O5mtKN7sza4eQDZnPK5XdUPQHx+AuTq1E/hRPT//cUXRDrHqUmo1o1BHAewdrL+YXuUTenEBEPpIcPAysj4n5ga+B1ZPkG60P1mozAeQQrhvML3aHZS1IXA9tIGgZ+Bfw1cEFRhbLOqvfjc18EK5r7L3RWs0FBEfEMcCjw7Yh4H7BHccWyTijdBrvWpae+0Z1NBvdf6Kymg0K68d2RwBVpWqMk9fmSHpJ0Z43X95X0uKTl6e/zzRfb2s2Xnlq38PgLndVsUDiRrAfzTyLiLkl/Alzb4D0XAAc2mOf6iJiT/k5rsixWgEaXnrrJyCZTo/zC6nVjbkoqSFN3SY2IxWR5hdLz+4BPNHqPpFkTKZwVz5eeWjer1X8BfKlqUZodjvNVks6R9AtJ15T+2rD+N0m6TdLPJb26zvqPlTQqaXTt2rVtWK1B4yYj8KWn1ln18gvgpqQiNFVTIBsz4bvAvwC1r1Ucn6XAyyPiKUkHAQuA3avNGBHnkI38xsjISK08qI2TLz21bteo/wK80JTkXs/t0WxOYUNE/FNE3BIRS0p/E1lxRDwREU+lx1cCUyXtPJFlWnNKVxnVqyE4j2DdopRfqJV4Bl+V1E7NBoWfSjpe0i6SXlz6m8iKJb1UktLjvVNZHpnIMq2xZpqMfOmpdSM3JU2OZpuPjkr/T85NC+BPar1B0iXAvsDOklYBXwCmAkTEd4HDgP8taQMwBhwRUWUoL2srNxlZr3JT0uRQrx2HR0ZGYnR0tNPF6DmNrjIC3/XUekej5s+hqVPc/FlB0pKIGGk0X7NXH20r6XOSzknPd5fkm+T1CDcZWb9xU1Jxms0p/CvwHPDm9HwV8KVCSmRt5yYj6zeNej1D1pS02/wr3MltnJoNCq+IiK8C6wEiYgxqdja0LuGrjKyfNXNVUuArk8ar2aDwnKQhsm2MpFcAzxZWKpswNxnZoGjUlARZc9InL1vuWkMTmr366FTg34GZki4G9iG7fbZ1mWYSyuAmI+sf+auS1qwbq3mXX/CtMZrR9NVHkl4CvJGs2eimiHi4yILV4quPamtmxDTwVUbW3xo1mcILteRB0u6rj34VEY9ExBUR8bOIeFjSryZeTGuHUu7gk5ctbyoguMnI+lkzzUm+y2ptjcZE2AbYlqwD2ot4Ibm8IzCj4LJZE5qtHYCbjGwwNNPJDdyUVEvd5iNJJwKfJAsAq3khKDwBnBsRZxdewgpuPso0mzsocZORDSI3qb6g2eajujWFiPgm8E1JH4+Ib7etdDYh460d+JJTG1SuNYzfeBLNbwZmkQskEfG9YopV2yDXFFw7MGtdMwlo6N/fTbM1haaCgqSLgFcAy3lhPIWIiLqjrxVhUIOCawdmEzOe35DIOmX1U4BoS/NRzgiwh+9iOvlcOzBrj2abkoByX4dBbFZqtqbwb8AnIuKB4otUX7/XFEpBYM26MXYamsrTz21g/cbG35FrB2bNG0+toaTXT7ja3Xx0LTAHuIXc7S0i4uCJFLIV/RwUWtlRofd3VrNOGG8tHHq7WandzUenTqw4Vk8rOye4dmA2EfPmDjNv7vC4TsYGoVnJg+x0SD4QlM4+xqMXz1TMutVEfo+98ltsS/ORpCepvn1EdvXRjq0XsTW9HBQmGgjAtQOzovVrs1JbcwrdpNeCwkQDwdQtxPbbbMm6Z9Yzo4t3OLN+02qOr1sDhINCB7WjRgDdt1OZDZqJ/pZL75k2NBWJjp7cOShMgsrLRyV47Jn1EwoE4CYis27U6gUh1XSiNuGgUJB21QIqdWuV08w21WqzUi35JuKdCqxROCi0kQOBmeUVdUyo1M7mp3b3Uxg4tb70iX75DgRmva/UxwGKDRClZa0bW1+eVnQficKCgqTzgfcAD0XEa6q8LuCbwEHAM8DREbG0qPKMR2X10IHAzGqZrACRN7Z+I2dddW9vBQXgAuBsoNbttd8F7J7+3gD8U/rfcWddde+E2wsdCMwGT7UA0e4LUUrWtCHhXU1hQSEiFkuaVWeWQ4DvpTuv3iRpmqRdOnnTvYl0WumGS87MrHvkA0Reqze9rDRj2lA7irmZTuYUhoGVueer0rTNgoKkY4FjAXbddddCCjPo91o3s8lRGSxaqVEUOd56J4OCqkyrug0i4hzgHMiuPiqiMI2ajBwIzKwIzdYoJqslopNBYRUwM/f8ZcCaDpWlbvucA4GZTbZawaJonQwKC4ETJF1KlmB+vJP5hBnThqrmEoanDXHj/P06UCIzs8m3RVELlnQJ8B/AbEmrJH1U0nGSjkuzXAncB6wAzgWOL6oszTj5gNkMTZ2yybQi2+3MzLpRkVcffbDB6wH8TVHrH6/8+K1r1o35CiIzG0gD36M5n8xxIDCzQTfQQaHyMtR+HmLPzKwZheUUekG1y1BL3cfNzAbRQAeFWpehFtV93Mys2w10UKjVTbyo7uNmZt1uoIOCL0M1M9vUQCeafRmqmdmmBjooQOe6kpuZdaOBbj4yM7NNOSiYmVmZg4KZmZU5KJiZWZmDgpmZlTkomJlZmYOCmZmVOSiYmVnZwHVe8/gJZma1DVRQ8PgJZmb1DVTzkcdPMDOrb6CCgsdPMDOrb6CCgsdPMDOrb6CCgsdPMDOrb6ASzR4/wcysvoEKCuDxE8zM6im0+UjSgZLulbRC0vwqrx8taa2k5envmCLLY2Zm9RVWU5A0BfgO8E5gFXCrpIURcXfFrJdFxAlFlcPMzJpXZE1hb2BFRNwXEc8BlwKHFLg+MzOboCKDwjCwMvd8VZpW6f2Sbpf0I0kzqy1I0rGSRiWNrl27toiympkZxQYFVZkWFc9/CsyKiD2BXwIXVltQRJwTESMRMTJ9+vQ2F9PMzEqKDAqrgPyZ/8uANfkZIuKRiHg2PT0XeH2B5TEzswaKDAq3ArtL2k3SVsARwML8DJJ2yT09GLinwPKYmVkDhV19FBEbJJ0AXAVMAc6PiLsknQaMRsRC4BOSDgY2AI8CRxdVHjMza0wRlc383W1kZCRGR0c7XQwzs54iaUlEjDSab6DufWRmZvU5KJiZWZmDgpmZlTkomJlZmYOCmZmVOSiYmVmZg4KZmZU5KJiZWZmDgpmZlTkomJlZmYOCmZmVOSiYmVmZg4KZmZU5KJiZWZmDgpmZlTkomJlZmYOCmZmVOSiYmVmZg4KZmZU5KJiZWZmDgpmZlTkomJlZmYOCmZmVOSiYmVlZoUFB0oGS7pW0QtL8Kq9vLemy9PrNkmYVUY4Fy1azz5nXsNv8K9jnzGtYsGx1EasxM+t5hQUFSVOA7wDvAvYAPihpj4rZPgo8FhGvBL4OfKXd5ViwbDWnXH4Hq9eNEcDqdWOccvkdDgxmZlUUWVPYG1gREfdFxHPApcAhFfMcAlyYHv8IeIcktbMQZ111L2PrN24ybWz9Rs666t52rsbMrC8UGRSGgZW556vStKrzRMQG4HHgJZULknSspFFJo2vXrh1XIdasGxvXdDOzQVZkUKh2xh8tzENEnBMRIxExMn369HEVYsa0oXFNNzOI7hG6AAAG6ElEQVQbZEUGhVXAzNzzlwFras0jaUtgJ+DRdhbi5ANmMzR1yibThqZO4eQDZrdzNWZmfaHIoHArsLuk3SRtBRwBLKyYZyFwVHp8GHBNRGxWU5iIeXOHOePQ1zI8bQgBw9OGOOPQ1zJvbmVLlpmZbVnUgiNig6QTgKuAKcD5EXGXpNOA0YhYCJwHXCRpBVkN4YgiyjJv7rCDgJlZEwoLCgARcSVwZcW0z+ce/wH4QJFlMDOz5rlHs5mZlTkomJlZmYOCmZmVOSiYmVmZ2nwFaOEkrQXub/HtOwMPt7E4neTP0p365bP0y+cAf5aSl0dEw96/PRcUJkLSaESMdLoc7eDP0p365bP0y+cAf5bxcvORmZmVOSiYmVnZoAWFczpdgDbyZ+lO/fJZ+uVzgD/LuAxUTsHMzOobtJqCmZnV4aBgZmZlAxcUJP29pNslLZf0C0kzOl2mVkk6S9J/ps/zE0nTOl2mVkn6gKS7JD0vqecuH5R0oKR7Ja2QNL/T5WmVpPMlPSTpzk6XZaIkzZR0raR70r51YqfL1ApJ20i6RdJt6XN8sdD1DVpOQdKOEfFEevwJYI+IOK7DxWqJpD8nG4Nig6SvAETEZzpcrJZI+lPgeeCfgZMiYrTDRWqapCnAfwHvJBs46lbggxFxd0cL1gJJbwOeAr4XEa/pdHkmQtIuwC4RsVTSDsASYF6vfS9p3PrtIuIpSVOBG4ATI+KmItY3cDWFUkBItqPK8J+9IiJ+kca2BriJbHS7nhQR90TEvZ0uR4v2BlZExH0R8RxwKXBIh8vUkohYTJtHP+yUiHggIpamx08C97D5OPFdLzJPpadT019hx62BCwoAkk6XtBI4Evh8o/l7xEeAn3e6EANqGFiZe76KHjz49DNJs4C5wM2dLUlrJE2RtBx4CLg6Igr7HH0ZFCT9UtKdVf4OAYiIz0bETOBi4ITOlra+Rp8lzfNZYAPZ5+lazXyWHqUq03q2BtpvJG0P/Bj4ZEVLQc+IiI0RMYesNWBvSYU17RU68lqnRMT+Tc76A+AK4AsFFmdCGn0WSUcB7wHe0e7xrdttHN9Lr1kFzMw9fxmwpkNlsZzUBv9j4OKIuLzT5ZmoiFgn6TrgQKCQiwH6sqZQj6Tdc08PBv6zU2WZKEkHAp8BDo6IZzpdngF2K7C7pN0kbUU21vjCDpdp4KUE7XnAPRHxtU6Xp1WSppeuLJQ0BOxPgcetQbz66MfAbLIrXe4HjouI1Z0tVWskrQC2Bh5Jk27q4Sup3gd8G5gOrAOWR8QBnS1V8yQdBHwDmAKcHxGnd7hILZF0CbAv2S2aHwS+EBHndbRQLZL0FuB64A6y3zvA36Wx43uGpD2BC8n2rS2AH0bEaYWtb9CCgpmZ1TZwzUdmZlabg4KZmZU5KJiZWZmDgpmZlTkomJlZmYOC9R1JL5V0qaTfSLpb0pWSXlXAek6VdFJ6fJqkljrnSZqTLmk167i+7NFsgyt1WPoJcGFEHJGmzQH+mOxOpq0ud0pEbKz1ekRM5B5ac4ARoKeun7f+5JqC9Zv/CayPiO+WJkTE8oi4Xpmz0v2W7pB0OGSBpMb0fdP9+H9A1gEKSZ9N4yb8kqwTJGn6BZIOS49/J+mLkpam5f2PNH1vSb+WtCz9n516QJ8GHK5sjI/DJW2nbFyDW9O8m90bStIukhan99wp6a2FbVEbKK4pWL95Ddl986s5lOys/HVkPXZvlbQYeHON6ZDdFvs1EfFbSa8nu4XFXLLfztI663o4IvaSdDxwEnAM2a0J3pbGv9gf+HJEvF/S54GRiDgBQNKXycbJ+Ei6vcEtkn4ZEU/nlv+XwFURcbqy8Ry2Hd9mMqvOQcEGyVuAS1Iz0IOSFgF/Vmf6E8AtEfHb9P63Aj8p3WdKUr37G5VuvraELBgB7ARcmO6/FWT3xa/mz4GDS/kKYBtgV7LxAEpuBc5PN3xbEBHLG398s8bcfGT95i7g9TVeq3aL63rTAZ6ueN7sfWGeTf838sLJ198D16YRzd5LdrCvVZ73R8Sc9LdrROQDQmkwnLcBq4GLJH24yXKZ1eWgYP3mGmBrSR8rTZD0Z5LeDiwma7ufImk62UH1ljrTKy0G3idpSNnwju8dZ9l2IjuIAxydm/4ksEPu+VXAx1PSHElzKxck6eXAQxFxLtmdQPcaZ1nMqnJQsL6SxpR4H/DOdEnqXcCpZOMb/AS4HbiNLHh8OiL+u870ymUvBS4DlpPdo//6cRbvq8AZkm4ku+NlybXAHqVEM1mNYipwu6Q70/NK+wLLJS0D3g98c5xlMavKd0k1M7My1xTMzKzMQcHMzMocFMzMrMxBwczMyhwUzMyszEHBzMzKHBTMzKzs/wM8JpOBO486KwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = 10\n",
    "\n",
    "# reset the intensity values\n",
    "I = np.zeros(N)\n",
    "I[0] = I_init\n",
    "\n",
    "# explicit euler\n",
    "for i in range(1, N):\n",
    "    I[i] = h * f(alpha(s[i-1], A, L), B(B0, B1, s[i-1], w), I[i-1]) + I[i-1]\n",
    "    \n",
    "plt.plot(s, I, 'o')\n",
    "plt.xlabel(\"Coordinate s\")\n",
    "plt.ylabel(\"Intensity\")\n",
    "plt.title(\"Numerical solution with the forward Euler method\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 4\n",
    "\n",
    "Now we try again with **A = 100**. We get now numerical instability. One needs a disproportionately small step size to tackle this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Numerical solution with the forward Euler method')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEWCAYAAAB1xKBvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3X28HHV59/HPl8MBkoBBJbYmJAYV06Y+BBuDFVDaIgQfAJHeJFArCqTc3qB3WxCs3IL4QJS+igawNEoMKoIUMYLEBjWQ+ICQQHhIoHhHBHMSSxLwyNNRTsLVP2Y2TJbdc/bs2dnZh+/79Tqvs/ubnZlrZmf32plr5jeKCMzMzIazS9EBmJlZe3DCMDOzmjhhmJlZTZwwzMysJk4YZmZWEycMMzOriRNGG5J0iKQHGzCdhyUd1oiYMtMMSa+uc9yGLFe9hpu/pKnp8u06ink0bJ1LOkjS/5f0lKRjGjHNvOSxraXTPV/SNxo93byM5vNRYVq5rNOhOGFUkL4Rj0oal2k7RdKtBYa1Q0T8OCKmFR3HaJV/eIpervL5j/YDKWmxpE83JrqKLgAujYg9I2JJjvPJVbqenk0TX+nvnqLjGi1Jt0o6peg4GskJo7pdgY8UHUS50fy6tY7zCmBdPSPmtR2NYrqfTxNf6e8NDQ2sjBL+/hshr7DqLgLOlLR3+YBKhyayvyYknSTpp5IultQv6SFJb0nbN0jaLOn9mXF3l/Qvkn6d7tlcLmlMOuxQSX2Szpb038BXS22Z8SdLul7SFkmPSbo0bX+VpOVp21ZJV1VankokvUPS/ZKelLRR0pmZYadKWi/pcUk3SJpYZRo7/cJKl/8n6eOVafM96S/K4yss15+m0+iXtE7SUZlhiyVdJummNMbbJb2qShxXSvqn9PGk9L37UPr81elyKDt/SV8HpgA3pvF9NDPJE9P3aqukj1eZ5zzgROCj6fg3ZgbPkHSvpN9J+pakPTLjvUvS3eky/0zS66tM/5fAKzPx7S5pYvp+PJ6+P6dmXn++pOskfUPSE8BpkgYk7ZMOP1fSNkkvSp9/WtIX0sfvlLRG0hPp9nt+Zrqlz8LJkn4NLE/b3yfpkXTbq7iOalG+TaRtVff8JL05XW/9ku6RdGhm2K2SPiPpp8Az6forH/9hSWel78/Tkq6Q9EeSvp9uZz+U9OLh5ifpM8AhwKXp+3NpZjaHKTmU+Nt0G1Y6zi7p+/CIku+Ir0kan5lXQ9bpqESE/8r+gIeBw4DrgU+nbacAt6aPpwIB7JoZ51bglPTxScA24ANAD/Bp4NfAZcDuwOHAk8Ce6eu/ANwAvATYC7gRuDAddmg6rc+l445J2/rS4T3APcDFwDhgD+DgdNirgben400AVgJfKF/OKuvgN8Ah6eMXA29MH/8VsBV4YzrdS4CVmfECeHX5Osmsl59Uem1mWUvL1QusB/4Z2C2d75PAtHT4YuBxYBbJ3uBVwDVVluWDwI3p4xOAXwLfygz7bvn8K62fzPv+5fR9eAPwB+BPq8x3Men2UzbNO4CJ6fv9AHBaOuyNwGbgwPR9fX/6+t2H2k4zz1cAX0q3gRnAFuCv02HnA4PAMSQ/FMek28N70+E3p+vlyPT5SuA9mfXyunS81wOPAseUrZOvkWx/Y4DpwFPAW0m2kX8l2YarbWsvWE+VtolKy50u1zfSx5OAx4B3pLG+PX0+IbM9/hr4M5JtprfKOv058Efp9DYDdwEHpMuyHDhvBPM7pWz6AXwP2JvkB8kWYHZmW1xPksj2JPn++Xo6bETrNK8/72EM7RPAGZIm1DHuryLiqxGxHfgWMBm4ICL+EBE3A88Cr05/XZwK/ENEPB4RTwKfBeZkpvUcyUb6h4gYKJvPLJIvn7Mi4umI+H1E/AQgItZHxA/S8baQbGRvqzH+QWC6pBdFxG8j4q60/URgUUTcFRF/AD4G/IWkqSNZOTV4M8mHZn5EPBsRy0k+aHMzr7k+Iu6IiG0kCWNGlWmtAA5RcgjircDngYPSYW9Lh4/EJyNiICLuIUnWIz18siAiNkXE4yQ/Dkpxnwr8e0TcHhHbI+JKkoT05uEmKGkycDBwdroN3A18BXhf5mW3RcSSiHgu3Y5WAG9Tsqf8emBB+nwP4E3AjwEi4taIuC8d717gal64HZ2fbn8DwHHA9yJiZbqN/D+SbXgoZ6a/0kt/Vw63zBX8LbA0Ipamsf4AWE3yhV6yOCLWRcS2iBisMp1LIuLRiNhIsg5uj4g16bJ8hyR51Dq/SuZHRH9E/Bq4heff/xOBf42IhyLiKZLP1pz0/alnnTacE8YQImItyZfUOXWM/mjm8UA6vfK2PUl++Y8F7ix9WID/TNtLtkTE76vMZzLwSPqluRNJL5N0jZJDSk8A3wD2qTH+95Js+I9IWiHpL9L2icAjpRelG/ZjJL+2GmkisCEish+KR8rm89+Zx8+QrM8XiIhfkvw6m0FymOB7wCZJ06gvYdQ03zrGfwXwT9kvTpL3t+IhvzITgdIPjpLy9bWhbJwVJL/g3wjcB/yAZH28GVgfEVsBJB0o6RYlhzx/B5zGC7ej7LQnZp9HxNMk28hQ/iUi9s78vX+Y11fyCuBvytbfwcDLq8RZTfnntNLnttb5VVLt/d/ps5U+3pVkb6eeddpwThjDO4/kl1/2g/d0+n9spu2P65z+VpKN8M8yH5bxEZH9EhqqS+ENwBRVLjZemI77+oh4EckvItUSVESsioijgZcBS4Br00GbSD4oACg5k+ylwMYKk3ma+tfRJmCydi5MTqkyn1qsIPmVtlv6y3EF8Hckh9vurjLOaLtyHun4G4DPlH1xjo2Iq2sYdxPwEkl7ZdrK11d5PD8DpgHvAVZExP3pOO9k5yT6TZJDppMjYjxwOS/cjrLT/g1JogNA0liSbaQeO21DknrY+cdU1gaSQzjZ9TcuIuZXiXO0hpvfSOe102eL5L3YRpKwGrlO6+aEMYyIWE9ySOnDmbYtJB/Ev5XUI+mDQMWCaw3Tf47kmPjFkl4GOwqzR9Q4iTtINqb5ksZJ2kNS6XDLXiS/rPslTQLOqmWCknaTdKKk8elu+xPA9nTwN4EPSJohaXeSw2e3R8TDFSZ1N3CspLFKTp89uWz4o1QoPKZuJ/my+Kik3rSY+G7gmlqWoYIVwOkkx+YhOb58BklNZXuVcYaKrxYjHf/LJMXoA5UYp6TgvNdwI0bEBpIEcGG6DbyeZH1fNcQ4zwB3Av+H5xPEz4C/Z+eEsRfJ3svvJc0iqQMN5TrgXZIOlrQbyem/9X7X/ALYI10PvcC5JMfwK/kG8G5JR6Sfyz2UFM33rXPewxlufiN9/68G/kHSfpL2JPlsfSs9etDIdVo3J4zaXEBS0Ms6leQL+DGSItrPRjH9s0mKXT9PDx39kOSX37DSL7t3kxS4fw30Acengz9Jcrjhd8BNJEW0Wr0PeDiN5zSSvRMi4kckx0+/TZKoXsXO9Zasi0lqNY8CV/LCL6/zgSvT3fn/VbZczwJHAUeS7IV9Cfi7iPivESxD1gqSL75SwvgJyS/XlVXHSPbQzk3jO3OI11VzBUkdqF/SsNdJRMRqku3qUuC3JNvESSOY31ySIvQmkmPt56XH1YeyguQEgzsyz7PrCeBDwAWSniSp613LECJiHUkS+ibJNvJbku1yKKWzyUp/W9Np/S6d/1dIfqQ9XW1aadI8muREiS0kewBnkdP3XA3z+yJwXHo21IIaJrkI+DrJuv8V8HuSHzX1rtOGU4RvoGRmZsPzHoaZmdXECcPMzGrihGFmZjVxwjAzs5p0VEd2++yzT0ydOrXoMKyN3bfxd1WHvW7S+KrDzNrZnXfeuTUihu3RoqMSxtSpU1m9enXRYVgbO2j+cjb2l/e+ApP2HsNPz/mrAiIyy5+kR4Z/VYclDLN6LVmzkYuWPcjG/gHEzpfojunt4awj2v72I2aj5oRhXW/Jmo187Pr7GBhMLvgO2JE0Ju09hrOOmMYxBzS6qyyz9uOEYV3vomUP7kgWJaVk4cNQZs/zWVLW9TZVqFkM1W7WrZwwrOtN3HvMiNrNupUPSVnXcqHbbGScMKwrudBtNnJOGNaVXOg2GznXMKwrudBtNnJOGNaVXOg2G7mWOSQl6RiSewm/DLgsIm6WNIXk7mNbgV+U3ZvXbMRc6DarX657GJIWSdosaW1Z+2xJD0paL+kcgIhYEhGnktySsnSL0dcAN0XEB4HpecZqna9U6C71FVUqdENSu7jw2Ne50G02hLwPSS0GZmcbJPUAl5Hcq3k6MFdSNhmcmw4HWAPMkbQcuCXnWK3DDVfodrIwG1quCSMiVgKPlzXPAtZHxEMR8SxwDXC0Ep8Dvh8Rd6Wv/QDJjez/iuRw1QtImidptaTVW7ZsyWlJrBO40G02OkUUvScBGzLP+9K2M4DDgOMknZYO+0/gw5IuBx6uNLGIWBgRMyNi5oQJw3bnbl3MhW6z0Smi6K0KbRERC4AFZY1rgeOaEpV1LBe6zRqjiITRB0zOPN8X2FRAHNYFfEW3WeMUkTBWAftL2g/YCMwBTiggDusCvqLbrHHyPq32auA2YJqkPkknR8Q24HRgGfAAcG1ErMszDuteLnSbNU6uexgRMbdK+1JgaZ7zNoOkoF3pHt0udJuNXMtc6W3WSC50mzWeE4Z1HBe6zfLhhGEdx4Vus3y4t1rrOC50m+XDCcM6jq/oNsuHD0lZRygVuTf1DzB+TC+9PWJw+/Olbhe6zUbPCcPaXnmRu39gkN5dxIvH9tL/zCATXeg2awgnDGt7lYrcg88FY3fblTWfOLygqMw6j2sY1vZc5DZrDicMa3sucps1hw9JWdvy1dxmzeWEYW3JV3ObNZ8ThrUlX81t1nyuYVhbcqHbrPmcMKwtudBt1nw+JGVtxYVus+I4YVjbcKHbrFhOGNY2XOg2K5ZrGNY2XOg2K1bLJAxJx0j6sqTvSjo80z5O0p2S3lVkfFY8F7rNipVrwpC0SNJmSWvL2mdLelDSeknnAETEkog4FTgJOD7z8rOBa/OM01rbkjUbOWj+8h2F7iwXus2aJ+89jMXA7GyDpB7gMuBIYDowV9L0zEvOTYcj6TDgfuDRnOO0FlUqdG9MDzuVCt2Q1C4uPPZ1LnSbNUmuRe+IWClpalnzLGB9RDwEIOka4GhJDwDzge9HxF3pa/8SGEeSWAYkLY2I57ITkzQPmAcwZcqUvBbFCuJCt1nrKOIsqUnAhszzPuBA4AzgMGC8pFdHxOUR8XEASScBW8uTBUBELAQWAsycOTPKh1t7c6HbrHUUkTDKD0MDREQsABZUGiEiFucakbWsiXuP2XE4qrzdzJqriLOk+oDJmef7ApsKiMPawFlHTGNMb89ObS50mxWjiD2MVcD+kvYDNgJzgBMKiMNaWKkLkE39A4wf08sevbv4/txmBcs1YUi6GjgU2EdSH3BeRFwh6XRgGdADLIqIdXnGYe2lvAuQ/oFBxvT2cPHxM5wozAqU91lSc6u0LwWW5jlva1+VzowaGNzORcsedMIwK1DLXOltVuIzo8xakxOGtRx3AWLWmtxbrbUM3+vCrLU5YVhL8L0uzFqfE4a1BHcBYtb6XMOwluBCt1nrc8KwluBCt1nr8yEpK5QL3WbtwwnDCuNCt1l7ccKwwrjQbdZeXMOwwrjQbdZenDCsMC50m7UXH5KypnOh26w9OWFYU7nQbda+nDCsqVzoNmtfrmFYU7nQbda+nDCsqVzoNmtfPiRlTeFCt1n7c8Kw3LnQbdYZWiZhSDoGeCfwMuCyiLhZ0jjgS8CzwK0RcVWRMVp9XOg26wy51jAkLZK0WdLasvbZkh6UtF7SOQARsSQiTgVOAo5PX3oscF3aflSesVp+XOg26wx5F70XA7OzDZJ6gMuAI4HpwFxJ0zMvOTcdDrAvsCF9vPNPVGsbLnSbdYZcE0ZErAQeL2ueBayPiIci4lngGuBoJT4HfD8i7kpf20eSNKrGKmmepNWSVm/ZsiWHpbB6LVmzkYPmL99R6M5yodus/RRxWu0knt9rgCQpTALOAA4DjpN0WjrseuC9kv4NuLHSxCJiYUTMjIiZEyZMyDFsG4lSoXtjetipVOiGpHZx4bGvc6HbrM0UUfQu/7EJEBGxAFhQ1vg08IGmRGUN5UK3WecpYg+jD5iceb4vsKmAOCxHLnSbdZ4iEsYqYH9J+0naDZgD3FBAHJYjF7rNOk/ep9VeDdwGTJPUJ+nkiNgGnA4sAx4Aro2IdXnGYc3jQrdZ58q1hhERc6u0LwWW5jlvaz5f0W3W2VrmSm9rfy50m3U291ZrDeNCt1lnc8KwhnGh26yz+ZCUjZq7LjfrDk4YNioudJt1DycMGxUXus26h2sYNioudJt1DycMGxUXus26hw9JWV1c6DbrPk4YNmIudJt1p5oShqSXRET5jZCsS7nQbdadaq1h3C7pPyS9Q1Kl+1lYF3Gh26w71ZowXgMsBN4HrJf0WUmvyS8sa2UudJt1p5oSRiR+kPY+ewrwfuAOSSsk/UWuEVrLcNflZt2t1hrGS4G/JdnDeJTk/ts3ADOA/wD2yytAaw0udJtZrWdJ3QZ8HTgmIvoy7aslXd74sKzVuNBtZrXWMM6NiE9lk4WkvwGIiM/lEpm1FBe6zazWhHFOhbaPNTIQa20udJvZkIekJB0JvAOYJGlBZtCLgG15BmbFK13Nval/gPFjeuntEYPbn7+m24Vus+4yXA1jE7AaOAq4M9P+JPAPjQxE0iuBjwPjI+K4tG0X4FMkCWp1RFzZyHladeVF7v6BQXp3ES8e20v/M4NMdKHbrOsMmTAi4h7gHklXRcSI9ygkLQLeBWyOiNdm2mcDXwR6gK9ExPyIeAg4WdJ1mUkcDUwCHgeyxXbLWaUi9+BzwdjddmXNJw4vKCozK9KQNQxJ16YP10i6N/N3n6R7a5j+YmB22TR7gMuAI4HpwFxJ06uMPw24LSL+EfjfNczPGsRFbjMrN9whqY+k/99Vz8QjYqWkqWXNs4D16R4Fkq4h2ZO4v8Ik+oBn08fbKwxH0jxgHsCUKVPqCdMqmLj3GDZWSA4ucpt1ryH3MCLiN+nDrcCGiHgE2B14A0l9ox6TgA2Z530kRfWXptd0HCCpdAbW9cARki4BVlaJcWFEzIyImRMmTKgzJCvx1dxmVk2tF+6tBA6R9GLgRySF8OOBE+uYZ6XOCyMiHgNOK2t8Bji5jnlYHXw1t5kNpdaEoYh4RtLJwCUR8XlJa+qcZx8wOfN8X+rfW7EG8tXcZjaUWi/cU9rJ4InATWlbvTdfWgXsL2k/SbsBc0j6pbKCudBtZkOpNWF8hOTK7u9ExLr0molbhhtJ0tUk/VBNk9Qn6eT09NzTgWXAA8C1EbGuvvCtkXw1t5kNpaa9hIhYSabonJ7h9OEaxptbpX0psLTGGC1nvj+3mdWi1u7NXwOcCUzNjhMRPrDd5lzoNrNa1VqH+A/gcuArVLkewtqTC91mVqtaE8a2iPi3XCOxQrjQbWa1qrXofaOkD0l6uaSXlP5yjcyawoVuM6tVrXsY70//n5VpC+CVjQ3HmsWFbjMbqVrPkvI9uzuIC91mVo9az5IaC/wjMCUi5knaH5gWEd/LNTrLhQvdZlaPWmsYXyXpNfYt6fM+4NO5RGS5c6HbzOpRa8J4VUR8HhgEiIgBKnciaG3AhW4zq0etRe9nJY0hrY1KehXwh9yisly40G1mo1Frwjgf+E9gsqSrgIOAD+QVlDWeC91mNlq1niV1s6Q7gTeTfM98JCK25hqZNZQL3WY2WjXVMCT9KCIei4ibIuJ7EbFV0o/yDs4ax4VuMxutIfcwJO0BjAX2Se+2Vyp0vwiYmHNs1kC+R7eZjdZwexh/D9wJ/En6v/T3XeCyfEOzRjrriGmM6e3Zqc2FbjMbiSH3MCLii8AXJZ0REZc0KSZroNKZUZv6Bxg/ppc9eneh/5lBJrrQbWYjVGvR+xJJb+GF98P4Wk5xWQOUnxnVPzDImN4eLj5+hhOFmY1YrUXvrwP/AhwMvCn9m5ljXNYAlc6MGhjczkXLHiwoIjNrZ7VehzETmB4RMewrrWX4zCgza6RauwZZC/xxnoEASHqlpCskXZdpO0bSlyV9V9LhecfQSdwFiJk1Uq0JYx/gfknLJN1Q+qtlREmLJG2WtLasfbakByWtl3QOQEQ8FBEnZ18XEUsi4lTgJOD4GuPtakvWbOSg+ct3dAGS5TOjzKxeI+kapF6LgUuBHQVyST0kp+W+naTn21WSboiI+4eYzrn4VN5huQsQM8tLrWdJrah3BhGxUtLUsuZZwPqIeAhA0jXA0cALEoYkAfOB70fEXRWGzwPmAUyZMqXeMDuGuwAxs7wMeUhK0pOSnqjw96SkJ0Yx30nAhszzPmCSpJdKuhw4QNLH0mFnAIcBx0k6rXxCEbEwImZGxMwJEyaMIqTO4EK3meVluAv39sppvpXupRER8RhwWlnjAmBBTnF0HHcBYmZ5qbWG0Wh9wOTM832BTQXF0hF8rwszy1tRCWMVsL+k/YCNwBzghIJiaXsudJtZM+SeMCRdDRxK0uNtH3BeRFwh6XRgGdADLIqIdXnH0qlc6DazZsg9YUTE3CrtS4Glec+/G7jQbWbNUOuFe9bCfEW3mTVDUTUMawAXus2smZww2pQL3WbWbE4YbcqFbjNrNtcw2pQL3WbWbE4YbcqFbjNrNh+SajMudJtZUZww2ogL3WZWJCeMNuJCt5kVyTWMNuJCt5kVyQmjjbjQbWZF8iGpNuBCt5m1AieMFudCt5m1CieMFudCt5m1CtcwWpwL3WbWKpwwWpwL3WbWKnxIqkW50G1mrcYJowW50G1mrcgJowW50G1mrailE4akKcClwFbgFxExv+CQmsKFbjNrRU0vektaJGmzpLVl7bMlPShpvaRz0ubXADdFxAeB6c2OtSgudJtZKyriLKnFwOxsg6Qe4DLgSJLEMFfSdGANMEfScuCWJsfZdEvWbOSg+ct3FLqzXOg2s6I1/ZBURKyUNLWseRawPiIeApB0DXA0MAicl45zHfDV8ulJmgfMA5gyZUqOkefLhW4za3WtUsOYBGzIPO8DDgQuB86XdALwcKURI2IhsBBg5syZUek17cCFbjNrda2SMMqPwABERKwFjmt2MEVwodvMWl2rXOndB0zOPN8X2FRQLIVwodvMWl2r7GGsAvaXtB+wEZgDnFBsSPkrXc29qX+A8WN66e0Rg9ufP6rmQreZtZIiTqu9GrgNmCapT9LJEbENOB1YBjwAXBsR65odWzOVitwb+wcIoH9gEAJePLYXkdQuLjz2dS50m1nLKOIsqblV2pcCS5scTmEqFbkHnwvG7rYraz5xeEFRmZlV1yo1jK7jIreZtRsnjIK4yG1m7aZVit5dw92Wm1m7csJoIl/NbWbtzAmjiXw1t5m1M9cwmsiFbjNrZ04YTeRCt5m1Mx+SagIXus2sEzhh5MyFbjPrFE4YOXOh28w6hWsYOXOh28w6hRNGzlzoNrNO4UNSOXGh28w6jRNGDlzoNrNO5ISRAxe6zawTuYaRAxe6zawTOWHkwIVuM+tEPiTVQC50m1knc8JoEBe6zazTtXzCkDQOWAmcFxHfKzqealzoNrNO1/QahqRFkjZLWlvWPlvSg5LWSzonM+hs4NrmRjlyLnSbWacroui9GJidbZDUA1wGHAlMB+ZKmi7pMOB+4NFmBzlSLnSbWadr+iGpiFgpaWpZ8yxgfUQ8BCDpGuBoYE9gHEkSGZC0NCKey44oaR4wD2DKlCn5Bl+BC91m1i1apYYxCdiQed4HHBgRpwNIOgnYWp4sACJiIbAQYObMmVE+PE8udJtZN2mVhKEKbTu+/CNicfNCqZ0L3WbWTVrlwr0+YHLm+b7ApoJiqZkL3WbWTVolYawC9pe0n6TdgDnADQXHNCwXus2smxRxWu3VwG3ANEl9kk6OiG3A6cAy4AHg2ohY1+zYarVkzUYOmr98R6E7y4VuM+tURZwlNbdK+1JgaZPDGTEXus2sW7VK0bttuNBtZt2qVWoYbcOFbjPrVk4YI+RCt5l1KyeMETrriGmM6e3Zqc2FbjPrBq5h1KjUBcim/gHGj+llj95d6H9mkIkudJtZl3DCqEH5mVH9A4OM6e3h4uNnOFGYWdfwIakaVDozamBwOxcte7CgiMzMms8JowY+M8rMzAmjJj4zyszMNYwh+V4XZmbPc8Kowl2AmJntzAmjCncBYma2M9cwqnCh28xsZ97DKFOqW1S716sL3WbWrZwwMsrrFuVc6DazbuaEkVGpblHiQreZdTsnjIxq9QmBC91m1vVc9M7wBXpmZtU5YWS463Izs+pa+pCUpHHAl4BngVsj4qo851eqT5S6MXfX5WZmz2t6wpC0CHgXsDkiXptpnw18EegBvhIR84Fjgesi4kZJ3wJyTRiQJA0nCDOzFyrikNRiYHa2QVIPcBlwJDAdmCtpOrAvsCF9WeXTl8zMrCmavocRESslTS1rngWsj4iHACRdAxwN9JEkjbupktwkzQPmpU+fkjSam1TsA2wdxfitolOWA7wsrapTlqVTlgNGtyyvqOVFrVLDmMTzexKQJIoDgQXApZLeCdxYacSIWAgsbEQQklZHxMxGTKtInbIc4GVpVZ2yLJ2yHNCcZWmVhKEKbRERTwMfaHYwZmb2Qq1yWm0fMDnzfF9gU0GxmJlZBa2SMFYB+0vaT9JuwBzghgLiaMihrRbQKcsBXpZW1SnL0inLAU1YFkVU65c1pxlKVwOHkhRoHgXOi4grJL0D+ALJabWLIuIzTQ3MzMyG1PSEYWZm7alVDkmZmVmLc8IoI+lTku6VdLekmyVNLDqmeki6SNJ/pcvyHUl7Fx1TvST9jaR1kp6T1HanQEqaLelBSeslnVN0PKMhaZGkzZLWFh3LaEiaLOkWSQ+k29ZHio6pXpL2kHSHpHvSZflkbvPyIamdSXpRRDyRPv4wMD0iTis4rBGTdDiwPCK2SfocQEScXXBYdZH0p8BzwL8DZ0bE6oJDqlnai8EvgLeTnA24CpgbEfcXGlidJL0VeAr4WrZrn3Yj6eXAyyPiLkl7AXcCx7Tj+yJJwLiIeEpSL/AT4CMR8fPw7DXtAAAEpklEQVRGz8t7GGVKySI1DqrerbWlRcTNEbEtffpzklOV21JEPBARo7mCv0g7ejGIiGeBUi8GbSkiVgKPFx3HaEXEbyLirvTxk8ADJBcQt51IPJU+7U3/cvnecsKoQNJnJG0ATgQ+UXQ8DfBB4PtFB9GlKvVi0JZfTJ0q7aroAOD2YiOpn6QeSXcDm4EfREQuy9KVCUPSDyWtrfB3NEBEfDwiJpP0jnt6sdFWN9xypK/5OLCNJvT0Oxq1LEubqtiLQdOjsIok7Ql8G/i/ZUcX2kpEbI+IGSRHEmZJyuVwYat0DdJUEXFYjS/9JnATcF6O4dRtuOWQ9H6SruT/Olq8WDWC96TduBeDFpUe7/82cFVEXF90PI0QEf2SbiXpEbzhJyZ05R7GUCTtn3l6FPBfRcUyGun9Rc4GjoqIZ4qOp4u1Si8GlpEWiq8AHoiIfy06ntGQNKF0FqSkMcBh5PS95bOkykj6NjCN5KycR4DTImJjsVGNnKT1wO7AY2nTz9vxbC8ASe8BLgEmAP3A3RFxRLFR1a6TejGo1lNDoUHVQdLBwI+B+0g+6wD/HBFLi4uqPpJeD1xJsn3tAlwbERfkMi8nDDMzq4UPSZmZWU2cMMzMrCZOGGZmVhMnDDMzq4kThpmZ1cQJw7qKpD+WdI2kX0q6X9JSSa/JYT7nSzozfXyBpLouTJQ0Iz0t16xwXXmlt3Wn9GKt7wBXRsSctG0G8EckPcrWO92eiNhebXhEjKY/shnATKDtrg+wzuM9DOsmfwkMRsTlpYaIuDsifqzERWn/VfdJOh6SJFOl/dD0fgrfJLn4C0kfT+978UOSiz9J2xdLOi59/LCkT0q6K53en6TtsyT9TNKa9P+09MrwC4Djldyf5XhJ45Tck2JV+toX9LUl6eWSVqbjrJV0SG5r1LqK9zCsm7yW5L4HlRxL8mv+DSRXMa+StBJ4S5V2SLouf21E/ErSn5N0+3EAyefqriHmtTUi3ijpQ8CZwCkkXTm8Nb1/yWHAZyPivZI+AcyMiNMBJH2W5D4nH0y7g7hD0g8j4unM9E8AlkXEZ5Tcj2PsyFaTWWVOGGaJg4Gr00NLj0paAbxpiPYngDsi4lfp+IcA3yn12yVpqP6iSh3d3UmSqADGA1emfZkFyT0NKjkcOKpUHwH2AKaQ3M+hZBWwKO1cb0lE3D384psNz4ekrJusA/68yrBK3ZAP1Q7wdNnzWvvZ+UP6fzvP/2j7FHBLehe7d5MkgmrxvDciZqR/UyIimyxKNzl6K7AR+Lqkv6sxLrMhOWFYN1kO7C7p1FKDpDdJehuwkqRW0CNpAskX7h1DtJdbCbxH0hglt/x89whjG0/yBQ9wUqb9SWCvzPNlwBlpAR9JB5RPSNIrgM0R8WWSHlnfOMJYzCpywrCukd4T5D3A29PTatcB55Pcn+I7wL3APSSJ5aMR8d9DtJdP+y7gW8DdJPdY+PEIw/s8cKGkn5L0OlpyCzC9VPQm2RPpBe6VtDZ9Xu5Q4G5Ja4D3Al8cYSxmFbm3WjMzq4n3MMzMrCZOGGZmVhMnDDMzq4kThpmZ1cQJw8zMauKEYWZmNXHCMDOzmvwPE1f3xSyEp4oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = 100\n",
    "\n",
    "# reset the intensity values\n",
    "I = np.zeros(N)\n",
    "I[0] = I_init\n",
    "\n",
    "for i in range(1, N):\n",
    "    I[i] = h * f(alpha(s[i-1], A, L), B(B0, B1, s[i-1], w), I[i-1]) + I[i-1]\n",
    "    \n",
    "plt.plot(s, np.abs(I), 'o')\n",
    "plt.ylim(1., 1e30)\n",
    "plt.yscale(\"Log\")\n",
    "plt.xlabel(\"Coordinate s\")\n",
    "plt.ylabel(\"Intensity\")\n",
    "plt.title(\"Numerical solution with the forward Euler method\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 5, 6, and 7\n",
    "\n",
    "Since we have a uniform grid, the optical depth for each grid point is\n",
    "$$ \\tau = h \\alpha = h A$$\n",
    "because the extinction coefficient is constant within the box. Now, the exact integral is\n",
    "$$I_{next} = I_{prev}  e^{-\\tau} + (1 - e^{-\\tau})  B$$\n",
    "assuming $B$ is constant throughout the cell and we simply use the average between two neighbouring grid points. This method gives reasonable emerging intensity $I_{out}$ and is called *first-order integration*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 1\n",
    "tau = h * A\n",
    "\n",
    "# reset the intensity values\n",
    "I = np.zeros(N)\n",
    "I[0] = I_init\n",
    "\n",
    "for i in range(1, N):\n",
    "    B_avg = B(B0, B1, (s[i] + s[i-1])/2, w) # can change this to use the average between two neighbouring grid points\n",
    "    I[i] = I[i-1] * np.exp(-tau) + (1 - np.exp(-tau)) * B_avg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Piecewise Analytic Solution for A = 1')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xm4HGWZ9/HvjxghyhKVKOSQEESIIijoERcQM4KyDELEBdBRcBlUBgVfYIzLKKIIDtelIuowGUFAWX2FvIyAgCIEUYSThTXGiQJDEpawhEWiknC/f1T1odPppbrTdbqr+ve5rnOlu6q66umTpO9+7vt5nlJEYGZmBrBBrxtgZmb9w0HBzMxGOSiYmdkoBwUzMxvloGBmZqMcFMzMbJSDQgFJulLSYX3QjrdKWtzrdrRD0j2S9urwtWdI+rdut6nJ9aZJCknP6/D1H5R0dbfblZ77U5IelPSUpJfkcQ3rDQeFPpV+eK1K/9M9KOlHkjYGiIh9I+KcXrcxIm6IiOl5nFvSNpKelfSDPM6f4fqHS/pN9baI+GREfK2Dc+0u6beSHpf0qKQbJb2he62tH0Ai4ryIeGc3r5NeazzwLeCdEbFxRDzSxXNfJ+kxSRt265w15z9K0oikv0k6O49rFJ2DQn97V0RsDLwOeAPwpR63Zyx9GHgMOCSvD4ixIGlT4OfA6cCLgSHgq8Dfetmu9fQyYCPgznZfqETdzx1J04C3AgEcsB7ta2Y58HXgrJzOX3gOCgUQEcuAK4EdYfTb1Mcr+yV9VNKi9BvWVZK2rtr3aknXpN9QH5T0hXT7BpJmSfqTpEckXSzpxem+cyQdmz4eSr+BHpk+f0V6LkmaIWlp1bU+J2mZpCclLZa0Z6trNfFhkiD4DPCu6h1pez4p6X/S9/x9SUr3bSvp2vQ6D0s6T9LE2pNL2kLS09WpD0mvl7RC0k7AGcCb057aynT/2ZK+XnX8gZIWSnoifW/71Hkf26d/hxdExJqIWBURV0fEbVW/my9JulfSQ5LOlbRZvV9IbepL0gmSfpI+nZv+uTJt85trezuS3iLplrTHcoukt1Ttu07S19JezJOSrpa0eZ02bA8srrrWtRnPfZKkG4GngZfXe38kf+c3AWcDuaRHI+KSiJgDdK13UzYOCgUgaQqwH7Cgzr6ZwBeAg4BJwA3ABem+TYBfAr8AJgOvAH6VvvQzwEzgbem+x4Dvp/uuB2akj98G/Dn9E2AP4IaoWR9F0nTgKOANEbEJsDdwT4Zr1Xu/bwW2Ai4ELib5sKi1P0nv6bXA+9PrAQg4Ob3Oq4ApwAm1L46IB4Dr0tdW/BNwYUTcDnwS+F2aHqkXVHYFzgWOByaS/F7uqT0O+COwJg20+0p6Uc3+w9OffyD5sNwY+F6d87SyR/rnxLTNv6tp74uBy4HvAi8hSf9crrXrAR8APgK8FHg+cFztRSLij8Crq6719ozn/hBwBLAJcG+D9/Bh4Lz0Z29JL2v0ZiX9QNLKBj+3NXqdZRAR/unDH5IPmKeAlST/iX4ATEj3XQd8PH18JfCxqtdtQPJtbGvgUGBBg/MvAvaser4lybfy5wHbptfdgOQb8yeApelx5wD/J308o2r7K4CHgL2A8Vmv1aBtPwTmpI/fnB770qr9Aexe9fxiYFaDc82s/h2kv9e90scHAzemj8cBDwC7ps8PB35Tc66zga+nj/8T+HbGv8tXpa9dCqwGLgNelu77FXBk1bHTq/4epqXv9Xm1bU+fnwD8JH281rG174HkQ/nmmnb9Dji86t/Ul6r2HQn8osH7qW1XlnOf2OJ3tHv6vjdPn/8B+GyO/7++Dpyd1/mL/OOeQn+bGRETI2LriDgyIlbVOWZr4LTKtyTgUZJvy0Mk35L/1ODcWwOXVr1uEbCG5MPqTyQBaWeSHO/PgeVpb+BtJD2JtUTEEuAYkg+qhyRdKGlyq2vVnkfSBOB9JN8WieQb7/+SfIut9kDV46dJvmEj6aXptZdJegL4CbBOGiT1/4AdJL0ceAfweETc3ODYWs1+t2uJiEURcXhEbEWSApwMfCfdPZm1vznfSxIQGn5L7lDtdSrXGqp6Xvd32qVz39fiHIcBV0fEw+nz88kphWTNOSgU333AJ9LgUfmZEBG/Tfdt2+R1+9a8bqNI6heQfPC/F3h+uu16ku79i4CF9U4YEedHxO4kQSCAb2a8VrV3A5sCP5D0gKQHSD5c6qWQ6jk5vfZrImJTkpSQGrT3ryS9jA+SfNv9cfXuFtdp9rttKCL+QNJr2DHdtJzk91UxlaQ38WCdl/8FeEHV8y2qT93i0rXXqVyr3t9Bu7Kcu2H70i8C7wfeVvV3/lngtZJe2+A1Z6S1k3o/bRfA7TkOCsV3BvB5Sa8GkLSZpPel+34ObCHpGEkbStpE0hurXneS0qK0pEmSDqw67/UkNYJKAfM64NMk6Yg1tY2QNF3S25WMFPorsIqkN5DlWtUOIxkZshNJT2VnYDdg57QA3MompGk3SUMkOf9mziVJsxxA0quoeBDYStLzG7zuTOAjkvZMi8VDkl5Ze5CkV0o6VtJW6fMpJGm9m9JDLgA+q2QI7sbAN4CLImJ1nWsuJBmNNV7SMEnQrlgBPEvjIu4VwPaSPiDpeZIOBnYg+Teyvtb33DNJ/q3swHN/568iqY/V/TIQyfDgjRv8vLreawDS9m1Eki4cJ2kjdTgPpKwcFAouIi4l+UZ+YZouuQPYN933JEla5F0kqYH/ISloApxGktu+WtKTJB9Sb6w69fUkH7CVoPAbkm+pc6lvQ+AU4OH0Wi8lKYBnuRaQjHQC9gS+ExEPVP3MIymWZ0knfJVkCO/jJMXPS5odHBE3knyYzo+Ie6p2XUsy5PIBSQ/Xed3NJEXZb6fXup51vy0DPEnyXn8v6S8k7/0O4Nh0/1kkPZS5wN0kAfXTDZr7byS9k8fS93l+VXueBk4CbkzTdG+qae8jJMX5Y0lG3vwrsH9VuqZjXTj3YcCPIuJ/q//eSQruH+zyh/aXSL6wzCLpRa5isIZ6t6QI32THBls6rPL8iPhhr9ti1msOCjbQlMwsvgaYkvaszAaa00c2sCSdQzKP4xgHBLOEewpmZjbKPQUzMxtVuKFYm2++eUybNq3XzTAzK5R58+Y9HBGTWh1XuKAwbdo0RkZGet0MM7NCkdRozam1OH1kZmajHBTMzGyUg4KZmY1yUDAzs1EOCmZmNqpwo4/MrLk5C5Zx6lWLWb5yFZtNGI8EK59+hskTJ3D83tOZuctQ65PYwCrcjObh4eHwkFSz+h/+jz39DKLxzQsq+yY6WAwcSfMiYrjlcQ4KZsVRCQTLVq5q+uHfrsq5hhwgSqvnQSG9mci5JHeHehaYHRGn1Rwzg+SWiHenmy6JiBObnddBwQZNXoGgEQeIcsoaFPKsKawGjo2I+ZI2AeZJuiYi7qo57oaI2D/HdpgV1pwFy/j8Jbez6pnkJnZj0a+vXGPZylV8/pLbARwYBkhuQSEi7gfuTx8/KWkRyb12a4OCmdWo7h300qpn1nDMRQs59arF7jUMiDGpKUiaRnK7wR0j4omq7TOAnwFLSW7+fVxErHPTbUlHAEcATJ069fX33ptpCQ+zQqrtHXSitqDcqgDdzjmdViqmfkgfVRqyMckH/zHVASE1H9g6Ip6StB8wB9iu9hwRMRuYDUlNIecmm/VEJ72DdkYTdTJaqZrTSoMh156CpPHAz4GrIuJbGY6/BxhudsNvF5qtjNrpHXT7G3unhexxEs9GeEhrQfS8pyBJwJnAokYBQdIWwIMREZJ2JZlh/UhebTLrJ9Xf3DeQWJPhC1oeqZuZuwyNnq+d3kqlve45lEue6aPdgA8Bt0tamG77AjAVICLOAN4LfErSamAVcEgUbeKEWQdqewatAsKE8eM4+aCdcv/QrQSIdusaq55Zw6lXLXZQKAFPXjPrgd1OuTZz7aBXhd1O0kouQvevnqePzGxd7aRnxqp30Ei9tFKrVJdTScXnnoLZGMmSkilC8TZrasm9hv7inoJZn8jaO+h1zyCrSvtavSf3GorJPQWzHJX9W3XW2khR31+ZZO0p+CY7Zjk69arFmQLCjbPeXsgPzOP3ns6E8eNaHlfpNcxZsGwMWmXrw0HBLAdzFizL9C16wvhxHL/39DFqVffN3GWIkw/aiaGJE1oeWxm2av3NQcGsyyopo1YBYWjihELUEFqZucsQN856O985eOeWvYZlK1ex2ynXusfQx1xoNuuyVimjohSU2+UCdDm4p2DWZcubfCCWpXfQSNZeg1NJ/ctBwaxLKnWERuP5ilxQbleWWoNTSf3JQcGsC1rVEYpeUO5EpdfQKjB4VFJ/cVAw64JmdYSyp4xaaTVs1amk/uKgYLYeWg09FQxMyqgRp5KKxUHBrENZhp5OzjB+fxA4lVQcDgpmHcoy9HTQ6gitOJXU/xwUzDo0yENPO+VUUv/z5DWzNlVWPW019NTqq9ynoVktxhPcesc9BbM2eOhp9ziV1J8cFMza4KGn3eNUUn9yUDBrQ6M6goeedsajkvqPg4JZBq2WsPDQ0/XjVFL/cFAwa8F1hPw5ldQ/HBTMWnAdYWw4ldQfHBTMWnAdYWw5ldRbDgpmDbiO0BtOJfWWg4JZHa4j9JZTSb3joGBWh+sI/cGppLHnoGBWh+sI/cGppLHnoGBWxXWE/uNU0thyUDBLuY7Q35xKGhu5BQVJUyT9WtIiSXdKOrrOMZL0XUlLJN0m6XV5tcesFdcR+ptTSWMjz6WzVwPHRsR8SZsA8yRdExF3VR2zL7Bd+vNG4D/SP83GXKs6gvWel93OX249hYi4PyLmp4+fBBYBtX9DBwLnRuImYKKkLfNqk1k9riMUj1NJ+RmTmoKkacAuwO9rdg0B91U9X8q6gQNJR0gakTSyYsWKvJppA8h1hGLKmkraZtblTie1KfegIGlj4GfAMRHxRO3uOi9Z5wtbRMyOiOGIGJ40aVIezbQB5TpCcWUZlRR4ZFK7cg0KksaTBITzIuKSOocsBaZUPd8KWJ5nm8yqeT5C8bVKJYHTSe3Ic/SRgDOBRRHxrQaHXQZ8OB2F9Cbg8Yi4P682mVW4jlAe1amkeqmHikZfAGxteY4+2g34EHC7pIXpti8AUwEi4gzgCmA/YAnwNPCRHNtjBjxXR2iUNnIdoXgqo5KAhiOTIt13/N7T3QNsIregEBG/oX7NoPqYAP4lrzaY1dOqjuAPjWI7fu/pDYO+h6u2lmdPwawveT5CuVU+7E+9anHdHkOlvuCgUJ+XubCB06he4DpCeVRGJjVKVXjmc2MOCjYwKsXlZStXrfNh4TpCOTUL9B6qWp+Dgg2E2klqwXMFL89HKC/PfG6fg4INhHrF5SAJCJ6PUF5eRK99Dgo2EBoVlz12vfx8P4b2OChYaVVqCNvMupwNVL/k6OLy4HAqKRsHBSul6hpCAGti3bnLLi4PFqeSsnFQsFJqNEFtnIRwcXlQOZXUmoOClVKjWsGzEdx9yj+6uDzgnEpqzEHBSskT1KwZp5Iac1CwUvEENcvKqaT6HBSsNDxBzTqRJZV0zEULB6bX4AXxrDRaTVAzq6fVAnoVg7LCqnsKVhqeoGadypJKgsEoQDsoWOH5LmrWLVlu7Vn2ArTTR1ZovouadZNTSe4pWMG1uouai8vWrkoq6TsH7zyQcxkcFKzQWt1FzQHBOjWocxkcFKyQXEewsTCIcxkcFKxwaucj1HIdwbptkOYyuNBshdOqjnD83tOdNrKuGqQCtHsKVjiuI1gvDMpcBgcFKxwvdme9VPa5DE4fWWHMWbBstPsuWKvI7DqCjZWyp5IUde5I1c+Gh4djZGSk182wMVZvklolMLiOYL3SavJkRT/8G5U0LyKGWx3nnoIVghe7s35Uxl6DawpWCF7szvpV2QrQDgrW1zxJzYqiLAVop4+sb3mxOyuSsqSScuspSDpL0kOS7miwf4akxyUtTH++nFdbrJi82J0VTTuL6fXrDOg8ewpnA98Dzm1yzA0RsX+ObbACazVJzaxfFbnXkKmnIOnF7Z44IuYCj7bdIht4riNYGbRTgD724lvZZtblfdFzyJo++r2kn0raT5JaH57ZmyXdKulKSa9udJCkIySNSBpZsWJFFy9v/caL3VnZZClAr4kg6I8VV7MGhe2B2cCHgCWSviFp+/W89nxg64h4LXA6MKfRgRExOyKGI2J40qRJ63lZ62euI1jZZLkvQ7Ve1xsyBYVIXBMRhwIfBw4DbpZ0vaQ3d3LhiHgiIp5KH18BjJe0eSfnsvLwYndWRlkL0NV61WvIWlN4iaSjJY0AxwGfBjYHjgXO7+TCkraopKIk7Zq25ZFOzmXF5zqCDYLqXoOAcS2y8b3oNWQdffQ74MfAzIhYWrV9RNIZ9V4g6QJgBrC5pKXAV4DxABFxBvBe4FOSVgOrgEOiaAsxWVd4PoINkpm7DI32eLOunTSWo5QyLYgn6f0RcXHNtvdFxE9za1kDXhCvfHY75dqGheV+WEjMLE/Vq/9m0en/iawL4mUtNM+qs+3zbbXIrAHXEWyQtVtvyLvW0DR9JGlfYD9gSNJ3q3ZtCqzOpUU2MCrfkFxHMMs+4Q2eW1wvjy9MrWoKy4ER4ABgXtX2J4HPdr01NjBcRzBbV6XekKXWkNcKwU2DQkTcCtwq6byIcM/AuqbVfATXEWyQZek15NWTbpU+ujgi3g8skFTdyxfJ9IXX5NIqKz2va2TWXLNeQ5496Vbpo6PTP71onXWF6whm7anuNSxfuYrJOfekW6WP7k8fPgysiohn0+UtXglcmUuLrLRcRzDrTPXchrxlHZI6F9hI0hDwK+AjJEtjm2XmdY3M+l/WGc2KiKclfQw4PSL+XdKCPBtm5eM6gln/y9pTULrw3QeBy9NtvpWnZeJ1jcyKI+sH+9EkM5gvjYg7Jb0c+HV+zbKycB3BrFgyBYX0Lmpzq57/GfhMXo2y8vB8BLNiyRQU0hFHxwHTql8TEU4EW1OuI5gVS9b00U+BM4AfAs3XeDXD8xHMiiprUFgdEf+Ra0usNFxHMCuurEHhvyUdCVwK/K2yMSIezaVVVmiuI5gVV9agcFj65/FV2wJ4eXebY0XW6mYhriOY9b+so4+2ybshVmxZlvp1HcGs/2WavCbpBZK+JGl2+nw7SV4kz0Y1SxmB6whmRZF1RvOPgL8Db0mfLwW+nkuLrJCa3fDD6xqZFUfWmsK2EXGwpEMBImKVJOXYLiuIVkNPhyZOcB3BrECyBoW/S5pAUlxG0rZUjUKyweShp2blkzUonAD8Apgi6TxgN5Lls22AeeipWflkHX10taR5wJtIRhYeHREP59oy60uVdNHylasapow89NSsuLKuffSriNiT55bNrt5mAyLLsFPw0FOzImsaFCRtBLwA2FzSi0i+BAJsCkzOuW3WZ1oNOwXXEcyKrlVP4RPAMSQBYB7PBYUngO/n2C7rI61mKkPyDyPvG4qbWf6aBoWIOA04TdKnI+L0MWqT9ZEsKSMPOzUrj6yF5tMlvYV176dwbk7tsj7hmcpmgyVrofnHwLbAQp67n0IADgollSVl5GGnZuWTdZ7CMLBDRDQahbgOSWcB+wMPRcSOdfYLOA3YD3gaODwi5mc9v+XHKSOzwZV17aM7gC3aPPfZwD5N9u8LbJf+HAH4Jj59wikjs8GVtaewOXCXpJtZ+yY7BzR6QUTMlTStyTkPBM5Nex83SZooacuIuD9jm6zLnDIys3aWuei2IeC+qudL023rBAVJR5D0Jpg6dWoOTTGnjMwMso8+uj6Ha9dbZbVuzSIiZgOzAYaHhzPXNSw7p4zMDFrPaH6S+h/UAiIiNl2Pay8FplQ93wpYvh7nsw44ZWRm1VpNXtskx2tfBhwl6ULgjcDjrieMLaeMzKxW1ppC2yRdAMwgWTdpKfAVYDxARJwBXEEyHHUJyZBUL8U9xpwyMrNauQWFiDi0xf4A/iWv61tjThmZWSO5BQXrT04ZmVkzWSevWUk4ZWRmzbinMCCcMjKzLBwUBoBTRmaWlYNCiWXpHYBTRmb2HAeFksp6P2WnjMysmoNCSWW5n7JTRmZWy0GhZJwyMrP14aBQIk4Zmdn6clAokSxzEE4+aCcHAzNryEGhBDwHwcy6xUGh4DwHwcy6yUGhoFxQNrM8OCgUkAvKZpYXB4UC8hwEM8uLg0KBOGVkZnlzUCgIp4zMbCw4KBSE5yCY2VhwUOhznoNgZmPJQaGPeQ6CmY01B4U+U+kZLF+5ig0k1kQ0PNYFZTPrNgeFPlLbM2gWEJwyMrM8OCj0kSzzD8ApIzPLj4NCH8g6/wCcMjKzfDko9FiWYvI4iWcjmOyUkZnlzEGhxzz/wMz6iYNCj3j+gZn1IweFHvD8AzPrVxv0ugGDKEvKyMVkM+sF9xTGkFNGZtbvcu0pSNpH0mJJSyTNqrP/cEkrJC1Mfz6eZ3t6qZIyahUQbpz1dgcEM+uZ3HoKksYB3wfeASwFbpF0WUTcVXPoRRFxVF7t6BdOGZlZEeSZPtoVWBIRfwaQdCFwIFAbFErNKSMzK5I8g8IQcF/V86XAG+sc9x5JewB/BD4bEffVHiDpCOAIgKlTp+bQ1Hx4lJGZFU2eNQXV2Va7wtt/A9Mi4jXAL4Fz6p0oImZHxHBEDE+aNKnLzcyPU0ZmVjR5BoWlwJSq51sBy6sPiIhHIuJv6dP/Al6fY3vGzJwFy9jtlGtbpow8U9nM+k2e6aNbgO0kbQMsAw4BPlB9gKQtI+L+9OkBwKIc2zMmnDIysyLLLShExGpJRwFXAeOAsyLiTkknAiMRcRnwGUkHAKuBR4HD82rPWHHKyMyKTNHkRi79aHh4OEZGRnrdjHV4lJGZ9TNJ8yJiuNVxntHcBU4ZmVlZeO2jLnDKyMzKwj2F9eCUkZmVjYNCh5wyMrMycvqoQ04ZmVkZuafQJqeMzKzMHBTa4JSRmZWd00dtcMrIzMrOPYU2LHfKyMxKzkEhg0ododHcb6eMzKwsHBRaaFVHcMrIzMrEQaGFZnUEp4zMrGwcFFpoVEcQOGVkZqXjoNBAqzrC5IkTxrQ9ZmZjwUGhDtcRzGxQOSjU4TqCmQ0qB4U6XEcws0HloFDFdQQzG3QOCinXEczMHBRGuY5gZuagMMp1BDMzBwXXEczMqgx0UHAdwcxsbQMdFFxHMDNb20AHBdcRzMzWNnBBoVJDWL5yFRtIrIl1qwmuI5jZoBqooFBbQ6gXEFxHMLNBNlBBoVENYZzEsxFMdh3BzAbcQAWFRjWEZyO4+5R/HOPWmJn1nw163YCx1KhW4BqCmVki16AgaR9JiyUtkTSrzv4NJV2U7v+9pGl5tGPOgmXsdsq1LFu5CtXscw3BzOw5uQUFSeOA7wP7AjsAh0raoeawjwGPRcQrgG8D3+x2OyrF5WVp6ihgNDAMTZzAyQft5BqCmVkqz5rCrsCSiPgzgKQLgQOBu6qOORA4IX38f4HvSVJEnWFBHapXXA6SgOC5CGZma8szfTQE3Ff1fGm6re4xEbEaeBx4Se2JJB0haUTSyIoVK9pqRKPicqPtZmaDLM+gUJu+B9ZZdy7LMUTE7IgYjojhSZMmtdUIF5fNzLLLMygsBaZUPd8KWN7oGEnPAzYDHu1mI47fezoTxo9ba5uLy2Zm9eUZFG4BtpO0jaTnA4cAl9UccxlwWPr4vcC13awnAMzcZYiTD9qJoYkTEC4um5k1k1uhOSJWSzoKuAoYB5wVEXdKOhEYiYjLgDOBH0taQtJDOCSPtszcZchBwMwsg1xnNEfEFcAVNdu+XPX4r8D78myDmZllN1Azms3MrDkHBTMzG+WgYGZmoxwUzMxslLo8AjR3klYA93b48s2Bh7vYnF7ye+lPZXkvZXkf4PdSsXVEtJz9W7igsD4kjUTEcK/b0Q1+L/2pLO+lLO8D/F7a5fSRmZmNclAwM7NRgxYUZve6AV3k99KfyvJeyvI+wO+lLQNVUzAzs+YGradgZmZNOCiYmdmogQsKkr4m6TZJCyVdLWlyr9vUKUmnSvpD+n4ulTSx123qlKT3SbpT0rOSCjd8UNI+khZLWiJpVq/b0ylJZ0l6SNIdvW7L+pI0RdKvJS1K/20d3es2dULSRpJulnRr+j6+muv1Bq2mIGnTiHgiffwZYIeI+GSPm9URSe8kuQfFaknfBIiIz/W4WR2R9CrgWeA/geMiYqTHTcpM0jjgj8A7SG4cdQtwaETc1fSFfUjSHsBTwLkRsWOv27M+JG0JbBkR8yVtAswDZhbt70WSgBdGxFOSxgO/AY6OiJvyuN7A9RQqASH1Qurc/rMoIuLq9N7WADeR3N2ukCJiUUQs7nU7OrQrsCQi/hwRfwcuBA7scZs6EhFz6fLdD3slIu6PiPnp4yeBRax7n/i+F4mn0qfj05/cPrcGLigASDpJ0n3AB4Evtzq+ID4KXNnrRgyoIeC+qudLKeCHT5lJmgbsAvy+ty3pjKRxkhYCDwHXRERu76OUQUHSLyXdUefnQICI+GJETAHOA47qbWuba/Ve0mO+CKwmeT99K8t7KSjV2VbYHmjZSNoY+BlwTE2moDAiYk1E7EySDdhVUm6pvVzvvNYrEbFXxkPPBy4HvpJjc9ZLq/ci6TBgf2DPbt/futva+HspmqXAlKrnWwHLe9QWq5Lm4H8GnBcRl/S6PesrIlZKug7YB8hlMEApewrNSNqu6ukBwB961Zb1JWkf4HPAARHxdK/bM8BuAbaTtI2k55Pca/yyHrdp4KUF2jOBRRHxrV63p1OSJlVGFkqaAOxFjp9bgzj66GfAdJKRLvcCn4yIZb1tVWckLQE2BB5JN91U4JFU7wZOByYBK4GFEbF3b1uVnaT9gO8A44CzIuKkHjepI5IuAGaQLNH8IPCViDizp43qkKTdgRuA20n+vwN8Ib13fGFIeg1wDsm/rQ2AiyPixNyuN2hBwczMGhu49JGZmTXmoGBmZqMcFMzMbJSDgpmZjXJQMDOzUQ4KVjqStpB0oaQ/SbpL0hWSts/hOidIOi59fKKkjibnSdo5HdJq1nOlnNFsgyudsHQpcE5EHJJu2xl4GclKpp2ed1xErGm0PyLWZw2tnYFhoFDj562c3FOwsvmE2oq3AAACaklEQVQH4JmIOKOyISIWRsQNSpyarrd0u6SDIQkkDbbPSNfjP59kAhSSvpjeN+GXJJMgSbefLem96eN7JH1V0vz0fK9Mt+8q6beSFqR/Tk9nQJ8IHKzkHh8HS3qhkvsa3JIeu87aUJK2lDQ3fc0dkt6a22/UBop7ClY2O5Ksm1/PQSTfyl9LMmP3Fklzgbc02A7Jstg7RsTdkl5PsoTFLiT/d+Y3udbDEfE6SUcCxwEfJ1maYI/0/hd7Ad+IiPdI+jIwHBFHAUj6Bsl9Mj6aLm9ws6RfRsRfqs7/AeCqiDhJyf0cXtDer8msPgcFGyS7AxekaaAHJV0PvKHJ9ieAmyPi7vT1bwUurawzJanZ+kaVxdfmkQQjgM2Ac9L1t4JkXfx63gkcUKlXABsBU0nuB1BxC3BWuuDbnIhY2Prtm7Xm9JGVzZ3A6xvsq7fEdbPtAH+peZ51XZi/pX+u4bkvX18Dfp3e0exdJB/2jdrznojYOf2ZGhHVAaFyM5w9gGXAjyV9OGO7zJpyULCyuRbYUNI/VzZIeoOktwFzSXL34yRNIvlQvbnJ9lpzgXdLmqDk9o7varNtm5F8iAMcXrX9SWCTqudXAZ9Oi+ZI2qX2RJK2Bh6KiP8iWQn0dW22xawuBwUrlfSeEu8G3pEOSb0TOIHk/gaXArcBt5IEj3+NiAeabK8993zgImAhyRr9N7TZvH8HTpZ0I8mKlxW/BnaoFJpJehTjgdsk3ZE+rzUDWChpAfAe4LQ222JWl1dJNTOzUe4pmJnZKAcFMzMb5aBgZmajHBTMzGyUg4KZmY1yUDAzs1EOCmZmNur/A4xE2YQCROs4AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(s, I, 'o')\n",
    "plt.xlabel(\"Coordinate s\")\n",
    "plt.ylabel(\"Intensity\")\n",
    "plt.title(\"Piecewise Analytic Solution for A = 1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Piecewise Analytic Solution for A = 10')"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3XmcXGWd7/HPl9BAI0tYopAmIagQRZBEW1wQQVADiBARBWUEXC4yioIXGIM6iDgIyr2OCyoXLwg4rI6QywhOQEE2ZcnGEkM04kISgYCELQ0m4Xf/OKcqlUotp6vrdG3f9+tVr3SdOnXOU+mu8zvP83sWRQRmZmYAG7S6AGZm1j4cFMzMrMhBwczMihwUzMysyEHBzMyKHBTMzKzIQaEDSPqFpGPaoBx7S1rU6nIMh6Q/S3pXg+89X9K/NrtMNc43SVJI2rDB9x8l6cZmlys99j9LekzSc5K2yeMc1h4cFNpEevEaSr90j0n6saTNACLiwIi4pNVljIjbI2JyHseWtJOklyT9II/jZzj/sZLuKN0WEcdHxNcaONbbJf1G0tOS/i7pTklval5pKweQiLgsIt7TzPOk5+oDvgW8JyI2i4gnm3jsX0t6StLGzTpm2fFPkDRb0ouSLq7w+v6SHpK0UtItknbMoxydxEGhvbwvIjYD3gC8Cfhyi8szmo4GngKOzOsCMRokbQH8HPgesDUwAHwVeLGV5RqhVwCbAAuG+0YlKl5nJE0C9gYCOGQE5atlGfBvwEUVzr8tcA3wryS/q9nAVTmVo2M4KLShiFgK/ALYDYp3U58svC7p45IWpndYs0rvbiS9TtJN6R3qY5K+mG7fQNIMSX+U9KSkqyVtnb52iaST058H0jvQT6fPX50eS5L2lbSk5FxfkLRU0rOSFknav965ajiaJAiuAt5X+kJanuMl/SH9zN+XpPS1V0m6OT3PE5IukzS2/OCStkvvBrcp2fZGScsl7Q6cD7w1ramtSF+/WNK/lex/qKT5kp5JP9sBFT7HLunv8IqIWBMRQxFxY0TcX/J/82VJf5H0uKRLJW1Z6T9EZU1fks6Q9B/p09vSf1ekZX5reW1H0tsk3ZvWWO6V9LaS134t6WtpLeZZSTemF8nyMuwCLCo5180Zj32WpDuBlcArK30+kt/5XcDFQC7NoxFxTUTMBCrVbg4DFkTETyPiBeAMYA9Jr8mjLJ3CQaENSZoAHATMq/DadOCLJH/Q44DbgSvS1zYHfgn8NzAeeDXwq/StnwOmA/ukrz0FfD997VZg3/TnfYCH038B3gHcHmXzoUiaDJwAvCkiNgemAX/OcK5Kn3dvYAfgSuBqkotFuYNJak97AB9Kzwcg4Oz0PK8FJpB8udcREY8Cv07fW/BPwJUR8QBwPPDbtHmkUlDZE7gUOBUYS/L/8ufy/YDfA2vSQHugpK3KXj82fbyT5GK5GXBehePU847037FpmX9bVt6tgeuB7wLbkDT/XK918wEfAT4GvBzYCDil/CQR8XvgdSXn2i/jsT8KHAdsDvylymc4GrgsfUyT9IpqH1bSDyStqPK4v9r76ngdcF/hSUQ8D/yRtZ+3JzkotJeZ6V3qHSQX6q9X2OdTwNkRsTAiVqf7TElrCwcDj0bE/46IFyLi2Yi4u+R9X4qIJRHxIsmF83AlbdK3Ansrqea/A/gmsFf6vn3S18utATYGdpXUFxF/jog/ZjhXJccAv4iIp4DLgQMlvbxsn3MiYkVE/BW4BZgCEBGLI+KmiHgxIpaTXKD2obJLSAIBksYAHwZ+UmXfcp8ALkrP9VJELI2Ih8p3iohngLeTNIn8CFgu6bqSC95RwLci4uGIeA44jaTJrKHkcg3vBf4QET+JiNURcQXwEOvWwn4cEb+PiCGSYDylice+OCIWpK+vKj+ApLcDOwJXR8QckovxR6qdMCI+HRFjqzxen7Hc5TYDni7b9jRJIOtZDgrtZXr6R75j+iUYqrDPjsB3CndJwN9J7pYHSO6S/1jhPYX3XVvyvoUkF/ZXpBfz50guCnuTtIkvS2sDFYNCRCwGTiK54D8u6UpJ4+udq/w4kvqBD5LcLZLe8f6V9S8Qj5b8vJLkC42kl6fnXirpGeA/gPWaQVL/jySIvRJ4N/B0RNxTZd9ytf5v15EG7GMjYgeSJsDxwLfTl8ez7p3zX4ANqfB/M0Ll5ymca6DkecX/0yYd+5E6xzgGuDEinkifX05OTUg1PAdsUbZtC+DZUS5HW3FQ6DyPAJ8qu1Pqj4jfpK+9qsb7Dix73yZp/gKSC//hwEbptltJqvdbAfMrHTAiLo+Iwh1fAN/IeK5S7yf5Iv5A0qOSHiW5uFRqQqrk7PTcr4+ILUhqAqpS3hdI7oiPImneKK0l1JsuuNb/bVVpbeJi0vwQSeKztIfLRGA18FiFtz8PbFryfLvSQ9c5dfl5Cueq9DsYrizHrlq+9EbgQ8A+Jb/zz5O05+9R5T3np7mTSo9hJ8BTC0iaIwvneBnJ77jR43UFB4XOcz5wmqTXAUjaUtIH09d+Dmwn6SRJG0vaXNKbS953VtrMhKRxkg4tOe6tJDmCQgLz18BngTsiYk15ISRNlrSfkp5CLwBDJLWBLOcqdQxJz5DdSWoqU0iarqakCeB6Nie541shaYCkzb+WS0na9A8hqVUUPAbsIGmjKu+7EPiYki6MGyhJyK+XkJT0GkknS9ohfT6BpJnqrnSXK4DPK+mCuxlJ899VaVNgufkkTUt9kgZJgnbBcuAlqidxbwB2kfQRSRtKOgLYleRvZKRGeuzpJH8ru7L2d/5akvxYxZuBtHvwZlUeVXMAafk2AcYAYyRtUtJUdy2wm6QPpPucDtxfqVmwlzgodJiIuJbkjvzKtLnkQeDA9LVnSZpF3kfSNPAHkoQmwHeA64AbJT1LcpF6c8mhbyW5wBaCwh0kd6m3UdnGwDnAE+m5Xk6SAM9yLiDp6QTsD3w7Ih4tecwhSZZnaU74KkkX3qdJkp/X1No5Iu4kuZjOjYg/l7x0M8kd4qOSnqjwvntIkrL/np7rVta/W4ak6eHNwN2Snif57A8CJ6evX0RSQ7kN+BNJQP1sleL+K8md61Pp57y8pDwrgbOAO9NmureUlfdJkhzTySQ9b/4FOLikuaZhTTj2MST5jL+W/t5JEu5HNTm/8mWSG5YZJLXIoXQbaQ7qAyT/j0+R/N6ObOK5O5LCi+xYj1HSrfLyiPi/rS6LWbtxULCeomRk8U3AhLRmZWYl3HxkPUPSJSTjOE5yQDCrzDUFMzMrck3BzMyKmj2KMnfbbrttTJo0qdXFMDPrKHPmzHkiIsbV26/jgsKkSZOYPXt2q4thZtZRJFWbg2odbj4yM7MiBwUzMytyUDAzsyIHBTMzK8otKKQTT90j6T5JCyR9tcI+G0u6StJiSXcrWZ7PzMxaJM/eRy8C+0XEc0oW/r5D0i8i4q6SfT4BPBURr5Z0JMlEb0fkWCazUTFz3lLOnbWIZSuG2LK/DwlWrFzF+LH9nDptMtOnDtQ/iFkLjMqIZkmbksy6+c8lK4EhaRZwRkT8Np0Z8VFgXNQo1ODgYLhLqrWLShf/p1auQlRfUKDw2oADhI0iSXMiYrDefrnmFCSNkTQfeBy4qTQgpAZIV2hK55N/mmTNV7O2N3PeUk675gGWrhgigBVDq3hqZbLyZK1brcJrS1cM8fmr5jNpxvXsdc7NzJzXjPVvzEYm18Fr6eIsUySNJVmecbeIeLBkl0orZK33fZJ0HMki4EycODGXspplVagdLF1RabXU4SkNEKdd8wCAaw7WUqM2IZ6krwDPR8T/Ktnm5iPrCKWBoFbT0EiNkXgpwrkHa7qWNx+lSzCOTX/uB94FlC9zdx1rV9c6HLi5VkAwa4XSZiLILyAArIkgWFtzcJOSjbY8cwrbA7dIuh+4lySn8HNJZ0o6JN3nQmAbSYuB/0myZJ5ZWzl31iKGVq23THVdhbbRsf19bLVp3zrbshhatYaTrprvfIONqtxyChFxPzC1wvbTS35+Afhg+T5m7WC4uYOxGbqeNtIM5XyDjaaOW2THOQUbDYUmoyw1hP6+MZx92O7DvmCXdmfdQGJNhu+iu7Fao7LmFDpu6myzPGWtHTRjrMH0qQPF92UNQq41WN4cFMxSWS/MedytF46VJSANrVrDubMWOShYLjwhnlkqS0J5YGw/d87YL5cL8vSpA9w5Yz++fcQU+vvG1Nx36YohJ6AtF64pWM/L2mTU3zeGU6dNzr08WWsNbkqyPLimYD2tfAxCNQNj+xtKJjcqa62h0JRk1iyuKVhPq9dk1GjPombJUmsoNCW5V5I1g2sK1pNmzlvKXufcXLOGMNq1g2oKtYaBsf1V9/EIaGsWBwXrOVmajPJMKDfq1GmT3ZRkuXPzkfWcLE1Go5FQHq4sTUnLmjBzq/U21xSs59S6cLZLk1E19ZqSAtxV1UbEQcF6RiGPUG0yiXZsMqqmVlOS8ws2Eg4K1hPq5RHatcmomulTBzj7sN2r1hicX7BGOShYT6iVR2j3JqNqCk1J1abj9qhna4QTzdYTquURBNw5Y7/RLUyTjR/bX3MMg0c923C4pmBdq5BD2GnG9WygyvfT42v0/e8U7qpqzeSagnWl8hlPK61V0Gl5hGrcVdWayTUF60rVcghjJETn5hGqcVdVaxbXFKwrVbszfimCP53z3lEuzeg5ddrkqmtCOL9gWbimYF2l3liEbsgh1OKuqjZSDgrWNbptLEKj6nVVdX7BanFQsK7RjWMRRqJarajba0s2Mg4K1jXqjUXopYAAlbuqCg9qs9qcaLaOV1hOs1fzCNWUd1UVFP+PnHS2anKrKUiaIOkWSQslLZB0YoV99pX0tKT56eP0vMpj3cl5hNpKu6qWB00nna2SPGsKq4GTI2KupM2BOZJuiojfle13e0QcnGM5rIvVyyN4icpEtaY1J52tXG41hYj4W0TMTX9+FlgI+NtpTeU8QjbVmtA8qM3KjUqiWdIkYCpwd4WX3yrpPkm/kPS60SiPdQ/3sMnG6y9YVrkHBUmbAT8DToqIZ8pengvsGBF7AN8DZlY5xnGSZkuavXz58nwLbB2hMEitkEAt1et5hEo8qM2yyjUoSOojCQiXRcQ15a9HxDMR8Vz68w1An6RtK+x3QUQMRsTguHHj8iyydYDy5HJAMTD04niErDyozbLILdEsScCFwMKI+FaVfbYDHouIkLQnSZB6Mq8yWXeolFwO1i6nabVVW3/BTW4G+dYU9gI+CuxX0uX0IEnHSzo+3edw4EFJ9wHfBY6MqDDHsVkJ96QZGQ9qs1pyqylExB1QtaZa2Oc84Ly8ymDdyXe6I+NBbVaLp7mwjuHkcvN4UJtV42kurCOUr6RWSC4XcgkepNYYN8VZOdcUrCPUSy47IDTGg9qsnIOCdQTf0ebDg9qsnIOCtbVeX0ktbx7UZuUcFKxteQbU0eFBbVbKQcHalldSG12eR8rAQcHamGdAHV0e1GbgoGBtzHeuo6s8v1BpUJsDQ/dzULC240FqreNBbebBa9ZWPEitPbgLcO9yULC24hlQ24Pnl+pdbj6ytuI71PbgpHPvclCwtuLkcntw0rl3OShYW3Byuf046dybnFOwlnNyub25Sa+3OChYyzm53N6cdO4tbj6ylvOdaHurlHTu20Cs/MdqdppxvRPPXcZBwVrOyeX2Vpp0FjC2vw8ET61cReDEc7dxULCWcXK5cxSSzn865728bOMNWbVm3dSzE8/dwzkFawknlzuXm/u6m4OCtYSTy53Liefu5uYjawnfbXYuj3bubg4K1hJOLncuj3bubg4KNqqcXO4OHu3cvXILCpImSLpF0kJJCySdWGEfSfqupMWS7pf0hrzKY61XvuZyIbkMXl6zU7kZsPvkWVNYDZwcEa8F3gJ8RtKuZfscCOycPo4DfphjeazF6iWXHRA6T7XmvgDnFzpUbkEhIv4WEXPTn58FFgLl3/pDgUsjcRcwVtL2eZXJWst3ld2nUtK5wPmFzjQqOQVJk4CpwN1lLw0Aj5Q8X8L6gQNJx0maLWn28uXL8yqm5aSQRyhvey5wcrlzlSedyzm/0HlyDwqSNgN+BpwUEc+Uv1zhLetdOyLigogYjIjBcePG5VFMy0l5HqGck8udr5B0rvRlBtcEO02uQUFSH0lAuCwirqmwyxJgQsnzHYBleZbJRlelPEKBk8vdxd2Mu0OevY8EXAgsjIhvVdntOuDotBfSW4CnI+JveZXJRl+1u0SBk8tdxoPaukOe01zsBXwUeEDS/HTbF4GJABFxPnADcBCwGFgJfCzH8lgLeEqE3lEI8OfOWlQch1I+qK10P2tPuQWFiLiDyjmD0n0C+ExeZbDWmTlvacWLAziP0M2mTx1g+tSB4gDFUoWks4NCe/OEeNZ0ngHV3P24c3maC2s6D1IzD2rrXA4K1nS+SzQPautcDgrWNB6kZgUe1Na5HBSsKTxIzcp5UFtnclCwpvAgNavGg9o6i4OCNYUHqVk1HtTWWdwl1UakMB7BeQSrxoPaOkummoKkrfMuiHUe5xEsK6/U1jmyNh/dLemnkg5K5zQycx7Bhs3dldtf1qCwC3AByVxGiyV9XdIu+RXLOoHzCDZcHtTW/jIFhXRltJsi4sPAJ4FjgHsk3SrprbmW0NqWe5XYcHlQW/vLlGiWtA3wTyQ1hceAz5JMez0F+CmwU14FtPZSSCwvWzHElv199I0Rq9asbSV2HsFqKU86l/Okea2Xtfnot8AWwPSIeG9EXBMRqyNiNnB+fsWzdlKaWA5gxdAqCNhq0z6E8wiWjQe1tbesXVK/HBFXl26Q9MGI+GlEfCOHclkbqpRYXvVSsOlGGzLv9Pe0qFTWqaqttVHIL3g23dbIWlOYUWHbac0siLWvwpxG1bqe+s7OGuH8QnuqWVOQdCDJymgDkr5b8tIWwOo8C2btoXxthEqcWLZGOL/QnurVFJYBs4EXgDklj+uAafkWzdpBrbEI4MSyjYzzC+2nZk0hIu4D7pN0WUS4ZtCDan0pvYqaNYvzC+2jXvPR1RHxIWCepNLR6SIZvvD6XEtnLVNvTqPCKmpmzXDqtMlVmyk9P9Loqtf76MT034PzLoi1j3p5BDcZWbM5v9A+auYUIuJv6Y9PAI9ExF+AjYE9SPIN1oU8p5G1Qr38gqfaHh1ZxyncBuwtaSvgVyTJ5yOAo/IqmLVOvTmNzPJULb8AbkoaDVnHKSgiVgKHAd+LiPcDu+ZXLGsFr7Fs7aDW+AXwVNt5yxwU0onvjgKuT7fVS1JfJOlxSQ9WeX1fSU9Lmp8+Ts9ebGs2r41g7WL61AHOPmx3BmrchLgpKT9Zg8KJJCOYr42IBZJeCdxS5z0XAwfU2ef2iJiSPs7MWBbLgfMI1k5KF+WpxqOe85F16uzbIuKQwjxHEfFwRHyu3nuAvzehjJajelNYeG0EayU3JY2+rFNn7wKcAkwqfU9EjDTr+FZJ95H0ZDolIhZUOf9xwHEAEydOHOEprcBTWFi7q9dVFTzqudmyNh/9FJgHfBk4teQxEnOBHSNiD+B7wMxqO0bEBRExGBGD48aNG+FprcBTWFgnqNeU5FXbmitrUFgdET+MiHsiYk7hMZITR8QzEfFc+vMNQJ+kbUdyTMumXpMROI9g7cezqo6OrOMU/kvSp4FrgRcLGyOi4ZyBpO2AxyIiJO1JEqCebPR4lk2WJiNPYWHtyKOeR0fWmsIxJM1Fv2HtTKmza71B0hUkK7ZNlrRE0ickHS/p+HSXw4EH05zCd4EjI6JaF3lrEjcZWSfzqOf8ZaopRMSw12COiA/Xef084LzhHtcaU5jgrl6TkWejtE7gUc/5yVRTkLSppC9LuiB9vrMkT5LXIeoNTIO1TUb+ElknyNJV9eSr72OnGde75jBMWXMKPyZpMnpb+nwJSY+kn+dRKGuOLLUDcJORdZ4sXVXXpK3RrjkMT9acwqsi4pvAKoCIGIKqzXrWBrLUDsC9jKxzZRn1XOBBbtllDQr/kNRP0iUYSa+ipBeStZ96CWVwk5F1h3pNSQVOQmeTtfnoDOC/gQmSLgP2Aj6WV6GscW4ysl5T2pS0bMUQG0jFpqNybkqqT1l7gUraBngLSbPRXRHxRJ4Fq2ZwcDBmz67ZG7ZnZRmDAO5lZN3N34PKJM2JiMF6+2Wd++hXEbE/a6fNLt1mLTac2oHzB9btsiShwbWGamrmFCRtImlrYFtJW0naOn1MAsaPRgGtNieUzdaXNQk9tGoNJ10137mGEvVqCp8CTiIJAHNY2+PoGeD7OZbL6shaOwBPW2G969RpkzM1JbnWsFbNmkJEfCcdzXxKRLwyInZKH3ukI5KtBbLWDsAJZettWVZxK3CtITGcRPPbWH89hUvzKVZ1vZxoHk7tAHovkWZWS9YENHRn/q3ZieafAK8C5gOF/9EARj0o9Kpe/4M2G6msCWhYW2s4d9ainruxyjpOYRDY1bOYjj7XDsyaZ/rUAaZPHch8k9WLuYasQeFBYDvgbzmWxVgbBJatGGLL/j6e/8dqVq2pH4tdOzDLzrWG6jLlFCTdAkwB7mHdRXYOya9olXVzTmE4TUSlXDswa9xwvnciaTfvxO9c1pxC1qCwT6XtEXFrA2UbkW4MCsNtIipw7cCsORr5DnZagGhqUGgn3RIUSv8IC39cw9Epf4hmnaTR2nonBIimBAVJz1L5eiUgImKLxovYmE4OCiMNBODagVneGq25F7RrgHBNoU2MNBD0bSA222RDVqxcxfg2+yMz62aN1hpKtVOAcFBooWbUCKA9/pDMelmzvsvQ+gDhoDDKmvnH4yYis/aTV4B452vGcctDy1m2YijX1gAHhVHQTXcRZpZdM7/75QrHG9vfh0TTmo4dFHLiQGBmpfIMEKVGer1o6txHlihPPDXyy3cgMOsuhakzIN8AUThW3lNv5BYUJF0EHAw8HhG7VXhdwHeAg4CVwLERMTev8oxEt3ZRM7PmGq0AMbRqDefOWtRZQQG4GDiP6jOpHgjsnD7eDPww/betdPNgFjPLT94BYlmDN6n15BYUIuK2dNnOag4FLk1nXr1L0lhJ20dEW026d+6sRZkDggOBmVVSKUAUehsVeh8NN2CMz7BwUCNamVMYAB4peb4k3bZeUJB0HHAcwMSJE0elcAX1orEDgZkNR2mAKFc+S7IET61ctV6wyHNFxVYGBVXYVjFIRsQFwAWQ9D7Ks1AFhV9OrZM5EJhZM1ULGOW1izyvO60MCkuACSXPdwCWtags66iXR/DgMjMbTbVqF822waicpbLrgKOVeAvwdLvkE2rlEQbG9jsgmFnXyrNL6hXAvsC2kpYAXwH6ACLifOAGku6oi0m6pH4sr7IMV7U8goA7Z+w3uoUxMxtFefY++nCd1wP4TF7nH4nxY/srjknIK9tvZtYuPKK5RK2+xHlm+83M2oWDQqrSFBbubmpmvcZBIVUpuVwICM4jmFmvaGXvo7ZSLbmc11ByM7N25KCQqpZEdnLZzHqJg0Lq1GmT6e8bs842J5fNrNc4p5AqJJFHayi5mVk7clAoMZpDyc3M2pGbj8zMrKjnawqjOfugmVm76+mgUD5gLe+1T83M2l1PNx9VGrBWWPvUzKwX9XRQ8IA1M7N19XRQ8IA1M7N19XRQ8IA1M7N19XSi2QPWzMzW1dNBATxgzcysVE83H5mZ2bocFMzMrMhBwczMihwUzMysyEHBzMyKHBTMzKzIQcHMzIpyDQqSDpC0SNJiSTMqvH6spOWS5qePT+ZZHjMzqy23wWuSxgDfB94NLAHulXRdRPyubNerIuKEvMphZmbZ5TmieU9gcUQ8DCDpSuBQoDwojCovqmNmVl2ezUcDwCMlz5ek28p9QNL9kv5T0oRKB5J0nKTZkmYvX7684QIVFtVZumKIYO2iOjPnLW34mGZm3STPoKAK26Ls+X8BkyLi9cAvgUsqHSgiLoiIwYgYHDduXMMF8qI6Zma15RkUlgCld/47AMtKd4iIJyPixfTpj4A35lgeL6pjZlZHnkHhXmBnSTtJ2gg4EriudAdJ25c8PQRYmGN5vKiOmVkduQWFiFgNnADMIrnYXx0RCySdKemQdLfPSVog6T7gc8CxeZUHvKiOmVk9iihv5m9vg4ODMXv27Ibf795HZtaLJM2JiMF6+/XcIjteVMfMrDpPc2FmZkUOCmZmVuSgYGZmRQ4KZmZW5KBgZmZFDgpmZlbkoGBmZkUOCmZmVuSgYGZmRQ4KZmZW5KBgZmZFDgpmZlbkoGBmZkUOCmZmVuSgYGZmRQ4KZmZW5KBgZmZFDgpmZlbkoGBmZkUOCmZmVuSgYGZmRQ4KZmZW5KBgZmZFuQYFSQdIWiRpsaQZFV7fWNJV6et3S5qURzlmzlvKXufczE4zrmevc25m5ryleZzGzKzj5RYUJI0Bvg8cCOwKfFjSrmW7fQJ4KiJeDfw78I1ml2PmvKWcds0DLF0xRABLVwxx2jUPODCYmVWQZ01hT2BxRDwcEf8ArgQOLdvnUOCS9Of/BPaXpGYW4txZixhatWadbUOr1nDurEXNPI2ZWVfIMygMAI+UPF+Sbqu4T0SsBp4Gtik/kKTjJM2WNHv58uXDKsSyFUPD2m5m1svyDAqV7vijgX2IiAsiYjAiBseNGzesQowf2z+s7WZmvSzPoLAEmFDyfAdgWbV9JG0IbAn8vZmFOHXaZPr7xqyzrb9vDKdOm9zM05iZdYU8g8K9wM6SdpK0EXAkcF3ZPtcBx6Q/Hw7cHBHr1RRGYvrUAc4+bHcGxvYjYGBsP2cftjvTp5a3ZJmZ2YZ5HTgiVks6AZgFjAEuiogFks4EZkfEdcCFwE8kLSapIRyZR1mmTx1wEDAzyyC3oAAQETcAN5RtO73k5xeAD+ZZBjMzy84jms3MrMhBwczMihwUzMysyEHBzMyK1OQeoLmTtBz4S4Nv3xZ4oonFaSV/lvbULZ+lWz4H+LMU7BgRdUf/dlxQGAlJsyNisNXlaAZ/lvbULZ+lWz4H+LMMl5uPzMysyEHBzMyKei0oXNDqAjSRP0t76pbP0i2fA/xZhqWncgpmZlZbr9UUzMysBgcFMzMr6rmgIOlrku6XNF/SjZLGt7pMjZI6xc9YAAAFlElEQVR0rqSH0s9zraSxrS5ToyR9UNICSS9J6rjug5IOkLRI0mJJM1pdnkZJukjS45IebHVZRkrSBEm3SFqY/m2d2OoyNULSJpLukXRf+jm+muv5ei2nIGmLiHgm/flzwK4RcXyLi9UQSe8hWYNitaRvAETEF1pcrIZIei3wEvB/gFMiYnaLi5SZpDHA74F3kywcdS/w4Yj4XUsL1gBJ7wCeAy6NiN1aXZ6RkLQ9sH1EzJW0OTAHmN5pv5d03fqXRcRzkvqAO4ATI+KuPM7XczWFQkBIvYwKy392ioi4MV3bGuAuktXtOlJELIyIRa0uR4P2BBZHxMMR8Q/gSuDQFpepIRFxG01e/bBVIuJvETE3/flZYCHrrxPf9iLxXPq0L33kdt3quaAAIOksSY8ARwGn19u/Q3wc+EWrC9GjBoBHSp4voQMvPt1M0iRgKnB3a0vSGEljJM0HHgduiojcPkdXBgVJv5T0YIXHoQAR8aWImABcBpzQ2tLWVu+zpPt8CVhN8nnaVpbP0qFUYVvH1kC7jaTNgJ8BJ5W1FHSMiFgTEVNIWgP2lJRb016uK6+1SkS8K+OulwPXA1/JsTgjUu+zSDoGOBjYv9nrWzfbMH4vnWYJMKHk+Q7AshaVxUqkbfA/Ay6LiGtaXZ6RiogVkn4NHADk0hmgK2sKtUjaueTpIcBDrSrLSEk6APgCcEhErGx1eXrYvcDOknaStBHJWuPXtbhMPS9N0F4ILIyIb7W6PI2SNK7Qs1BSP/Aucrxu9WLvo58Bk0l6uvwFOD4ilra2VI2RtBjYGHgy3XRXB/ekej/wPWAcsAKYHxHTWluq7CQdBHwbGANcFBFntbhIDZF0BbAvyRTNjwFfiYgLW1qoBkl6O3A78ADJ9x3gi+na8R1D0uuBS0j+tjYAro6IM3M7X68FBTMzq67nmo/MzKw6BwUzMytyUDAzsyIHBTMzK3JQMDOzIgcF6zqStpN0paQ/SvqdpBsk7ZLDec6QdEr685mSGhqcJ2lK2qXVrOW6ckSz9a50wNK1wCURcWS6bQrwCpKZTBs97piIWFPt9YgYyRxaU4BBoKP6z1t3ck3Bus07gVURcX5hQ0TMj4jblTg3nW/pAUlHQBJIqmzfN52P/3KSAVBI+lK6bsIvSQZBkm6/WNLh6c9/lvRVSXPT470m3b6npN9Impf+OzkdAX0mcISSNT6OkPQyJesa3Jvuu97cUJK2l3Rb+p4HJe2d2/+o9RTXFKzb7EYyb34lh5Hcle9BMmL3Xkm3AW+rsh2SabF3i4g/SXojyRQWU0m+O3NrnOuJiHiDpE8DpwCfJJma4B3p+hfvAr4eER+QdDowGBEnAEj6Osk6GR9Ppze4R9IvI+L5kuN/BJgVEWcpWc9h0+H9N5lV5qBgveTtwBVpM9Bjkm4F3lRj+zPAPRHxp/T9ewPXFuaZklRrfqPC5GtzSIIRwJbAJen8W0EyL34l7wEOKeQrgE2AiSTrARTcC1yUTvg2MyLm1//4ZvW5+ci6zQLgjVVeqzTFda3tAM+XPc86L8yL6b9rWHvz9TXglnRFs/eRXOyrlecDETElfUyMiNKAUFgM5x3AUuAnko7OWC6zmhwUrNvcDGws6X8UNkh6k6R9gNtI2u7HSBpHclG9p8b2crcB75fUr2R5x/cNs2xbklzEAY4t2f4ssHnJ81nAZ9OkOZKmlh9I0o7A4xHxI5KZQN8wzLKYVeSgYF0lXVPi/cC70y6pC4AzSNY3uBa4H7iPJHj8S0Q8WmN7+bHnAlcB80nm6L99mMX7JnC2pDtJZrwsuAXYtZBoJqlR9AH3S3owfV5uX2C+pHnAB4DvDLMsZhV5llQzMytyTcHMzIocFMzMrMhBwczMihwUzMysyEHBzMyKHBTMzKzIQcHMzIr+P+53bWT7NDuFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = 10\n",
    "tau = h * A\n",
    "\n",
    "# reset the intensity values\n",
    "I = np.zeros(N)\n",
    "I[0] = I_init\n",
    "\n",
    "for i in range(1, N):\n",
    "    B_avg = B(B0, B1, (s[i] + s[i-1])/2, w) # can change this to use the average between two neighbouring grid points\n",
    "    I[i] = I[i-1] * np.exp(-tau) + (1 - np.exp(-tau)) * B_avg\n",
    "    \n",
    "plt.plot(s, I, 'o')\n",
    "plt.xlabel(\"Coordinate s\")\n",
    "plt.ylabel(\"Intensity\")\n",
    "plt.title(\"Piecewise Analytic Solution for A = \" + str(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Piecewise Analytic Solution for A = 100')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xm4HFW57/HvjxBhyxSVqGQTCSpGGYToFgdEEVEGESKioF4Fh8NBDwpeRIMDzoKH53qcjxePCHiYHICbI2hAg0zKsBPCZIxGREkCskEChERIwnv/qOpO0emhdu+uHn+f5+lnd3dVV60edr211rvWKkUEZmZmAJt0ugBmZtY9HBTMzKzMQcHMzMocFMzMrMxBwczMyhwUzMyszEGhC0n6haSjuqAce0ta0ulyjIekuyTt1+RrvyfpM60uU539zZAUkjZt8vXvknR5q8uVbvuDkv4uaZWkZxSxD+tODgodkh681qT/dH+X9ENJWwJExIERcXanyxgR10TEzCK2LWlHSU9I+m4R28+x/6MlXZt9LiKOjYgvNrGtV0v6raSHJP1D0nWSXta60lYPIBFxbkS8sZX7Sfc1Gfga8MaI2DIiHmjhtn8j6UFJm7VqmxXbP07SqKTHJJ1VZfnrJf1B0mpJV0raIbNsM0lnSnpY0r2S/ncRZex2Dgqd9eaI2BJ4CfAy4NMdLk87vQd4EDiyqANEO0jaGvg58C3g6cAw8HngsU6Wa4KeBWwO3DHeFypR9bgiaQawNxDAIRMoXz0rgC8BZ1bZ/7bARcBnSL6rUeDCzCqfA3YCdgBeB3xc0gEFlbN7RYRvHbgBdwH7ZR6fDvw8vf8b4AOZZe8DFpMcROcBO2SW7QJcAfwD+DvwyfT5TYA5wJ+BB4AfA09Pl50NnJjeHyb5J/1Q+vj56bYE7AMsy+zrE8By4BFgCfD6Rvuq8/7/DHwwLfPhFcsCOBb4U/qevwMoXfY8YH66n/uBc4EplZ8r8GxgNfCMzLKXAmPAbsA/gfXAKmBluvws4EuZ9Q8FFgEPp+U9oMr7GCm9vsb73IQk2P8VuA84B9gmXTYjfa+b1vhNfA747/T+39J1V6W3VwJHA9dm1n8VcBPwUPr3VZllvwG+CFyXfn+XA9tWKe8LgEcz+5qfc9tfTre9Bnh+jc/ilHSdr5H+1gv8//oScFbFc8cAv8083iIt7wvTx8tJakel5V8ELuj0saLdN9cUuoCk6cBBwM1Vls0GPgkcBkwFrgHOT5dtBfwK+CUwjeSA/uv0pR8BZgOvTZeVDq4AV5Ec8EmX35n+BXgNcE2k/xWZcswEjgNeFhFbAfuTHMQa7ava+90b2B64gCSAvKfKageT1J52B96e7g+SYHVqup8XAdNJDp5PEhH3khys3p55+n+R/JPfRhJ0fhdJ88iUKmXck+QAfhIwheRzuatyPeCPwHpJZ0s6UNLTKpYfnd5eBzwX2BL4dpXtNPKa9O+UtMy/qyjv04FLgW8CzyA58F5akQ94J/Be4JnAU4CPVe4kIv5IcqJR2te+Obf9bpKD7lYkAbCa95AE8XOB/SU9q9ablfRdSStr3G6t9boGdgFuKT2IiEdJgv0u6fc2Lbs8vb8LA8ZBobMukbQSuJbkQP2VKuv8K3BqRCyOiHXpOnukbaEHA/dGxP+JiH9GxCMRcUPmdZ+KiGUR8RjJgfPwtE36KmDvtJr/GuDfgb3S1702XV5pPbAZsLOkyRFxV0T8Oce+qjkK+EVEPAicBxwo6ZkV65wWESsj4m/AlcAeABGxNCKuiIjHImKM5AD1Wqo7myQQIGkS8A7gRzXWrfR+4Mx0X09ExPKI+EPlShHxMPBqkjPr7wNjkuZmDnjvAr4WEXdGxCrgZJIms6aSy3W8CfhTRPwoItZFxPnAH4A3Z9b5YUT8MSLWkATjPVq47bMi4o50+drKDUh6NUmzzI8jYgHJwfidtXYYER+KiCk1bi/OWe5KW5LUdLIeIglkW2YeVy4bKA4KnTU7/ZHvkP4TrKmyzg7AN0pnSWxo2hkmOUv+c5XXlF53ceZ1i0kO7M9KD+arSA4Ke5O0ia9IawNVg0JELAVOIDng3yfpAknTGu2rcjuShoC3kZwtkp7x/o2NDxD3Zu6vJv2nlfTMdN/LJT0M/DewbY3P4P+RBLHnAm8AHoqIG2usW6neZ/skacA+OiK2B3YlOeP8erp4Gk8+c/4rsClVPpsJqtxPaV/DmcdVP9MWbfvuBts4Crg8Iu5PH5+XPtdOq4CtK57bmqQ5bVXmceWygeKg0P3uBv614kxpKCJ+my57Xp3XHVjxus0jYnm6/CrgcOAp6XNXkVTvn0bSjr6RiDgvIkpnfAF8Nee+st5C8s/23bSHx70kB5dqTUjVnJru+8URsTVJTUA1yvtPkjPid5E0b2RrCY2mB6732daU1ibOIgkOkCQ+d8is8hxgHUkupdKjwFMzj5+d3XSDXVfup7Svat/BeOXZds3ypScCbwdem/nOPwrsLmn3Gq/5Xtozr9pt3Anw1B0kzZGlfWxB8h3fkdZa78kuT+83u6+e5aDQ/b4HnCxpFwBJ20h6W7rs58CzJZ2QdqfbStLLM6/7cqnLnaSpkg7NbPcqkhzB1enj3wAfJklcrq8shKSZkvZNewr9kyRBV1qv0b6yjiLpGbIbSU1lD5Kmqz0k7Zbj89iKNDksaZikzb+ec0ja9A8hqVWU/B3YXtJTarzuB8B70y6Mm0galvTCypUkvVDSiZK2Tx9PJ2mmuj5d5Xzgo0q64G5J0vx3YdoUWGkRSdPSZEkjJEG7ZAx4giQvUc1lwAskvVPSppKOAHYm+Y1M1ES3PZvkt7IzG77zF5Hkx6qeDETSPXjLGrea7fxp+TYHJgGTJG2eaaq7GNhV0lvTdU4Bbs00C54DfFrS09Lv+l9IAvxAcVDochFxMckZ+QVpc8ntwIHpskdImkXeTNI08CeShCbAN4C5wOWSHiE5SL08s+mrSA6wpaBwLclZ6tVUtxlwGkmPn3tJkpWfzLkvANKD+OuBr0fEvZnbApJkeZ7mhM+TdOF9iCT5eVG9lSPiOpKD6cKIuCuzaD7JWeC9ku6v8robSZKy/5Hu6yo2PluGpHnh5cANkh4lee+3Ayemy88kqaFcDfyFJKB+uEZxP0Ny5vpg+j7Py5RnNWkPn7SZ7hUV5X2AJMd0IknPrI8DB2eaa5rWgm0fRZLP+Fv2eydJuL+rxfmVT5OcsMwhqUWuSZ8jzUG9leRzfJDkezsy89rPkjQZ/pXk+z49In7ZwrL1hFI3P7O+JWk+cF5E/Feny2LW7RwUrK8pGVl8BTA9rVmZWR1uPrK+JelsknEcJzggmOXjmoKZmZW5pmBmZmWtHlVZuG233TZmzJjR6WKYmfWUBQsW3B8RUxut13NBYcaMGYyOjna6GGZmPUVSrTmpnsTNR2ZmVuagYGZmZQ4KZmZW5qBgZmZlDgpmZlZWWO+jdBbCq0kmUtsU+GlEfLZinc1IZiZ8KclEW0dUTFpm1tUuuXk5p89bwoqVa9hmaDISrFy99kn3p00Z4qT9ZzJ71nDjDZp1WGEjmiUJ2CIiVkmaTDIL5/ERcX1mnQ+RzIt/rKQjgbdExBH1tjsyMhLukmqdVAoEy1euQTS+0AFQXm/YAcI6RNKCiBhpuF47prmQ9FSSoPDBzOUikTQP+FxE/C6dPvdeYGrUKZSDgnVCM4GgFgcI64S8QaHQnIKkSZIWAfcBV2QDQmqY9DJ+6UVHHiK5MHjldo6RNCppdGxsrMgim23kkpuXc/JFt7F8ZXK11ImeRpVev3zlGk6+6DYuubkVF0cza41Cg0JErI+IPYDtgT0l7VqxSrXLKG70PxcRZ0TESESMTJ3acJS2WUtccvNy9jptPidcuIg1aze6GF1LrFm7nhMuXMRep813cLCu0JbeRxGxkuRyjwdULFpGcoF00uajbUguTG/WUZW1g6K51mDdosjeR1OBtRGxMr1w935suNB7yVySS/X9juR6tPPr5RPMipbNHeRVyhFMqdL76MHVa3PnINasXc+JP76Fj164yD2WrGOKnBBvO+BsSZNIaiQ/joifS/oCMBoRc0kujv4jSUtJaghH1t6cWbFKtYM8TUXjSRaPJ0m9Pj0nKtUcAAcGa6ueu8iOex9ZUfY6bX6uGsJEeg2NtyYyPGWI6+bsO+79mFXqqi6preSgYK2W90A9NHkSpx62W0vO3MdTK3HXVWuFvEGh566nYNZKeQ/OrT4wl7ZTGg29iVRuOqrkpiRrJwcFG2inz1tSNyC0snZQafas4fJ2GwWnNWvXc/q8JQ4KVjhPiGcDqTQGoV6T0fCUocICQqXZs4Y59bDdGJ4yVHOd5SvXeDyDFc41BRs4eZqMOpHgLdUc6gUrNyVZ0VxTsIGTp8nopP1ntrFET3bS/jMZmjyp5vJSU5JZERwUbGB0W5NRLW5Ksk5y85ENhG5tMqrFTUnWKa4p2EDo9iajWtyUZO3mmoL1rexV0eoN0ezmwWHZ8Qz1agx7nTa/a9+D9RYHBetL4xmU1i1NRrW4Kcnayc1H1pcaNRdB9zYZ1eKmJGsH1xSsL62o08NI0JNTU+dpSqr3vs3y8IR41lcaTW7XC81FedRrSurmHIl1Tldco9msnRpdLa3XmovqqdeU5Ku42UQ4KFjfqJdH6IZBaa3UaICb8wvWLOcUrG/Uak8X9EWTUaVSr6Qd51xatcut8wvWDNcUrG9Mq3HWXOv5flHr/QV4KgwbNwcF63nZOY1Usayf8gi1OL9greSgYD2tMrkcUA4M/ZZHqMX5BWsl5xSsp1VLLgf90/U0L+cXrFVcU7CeVutgN6gHwUHNq1jrOChYTyrlEWoNvRzUg2C1/ILw9RcsPzcfWc9pNNndICSXa6mcCkNQDpyeNM/yKKymIGm6pCslLZZ0h6Tjq6yzj6SHJC1Kb6cUVR7rH4M0SK0Zs2cNc92cfRmeMrRRTcpJZ2ukyJrCOuDEiFgoaStggaQrIuL3FetdExEHF1gO6zODNkitWc63WDMKqylExD0RsTC9/wiwGBjc0zebMOcRxseD2qwZbUk0S5oBzAJuqLL4lZJukfQLSbvUeP0xkkYljY6NjRVYUutWgzTZXat4UJs1o/CgIGlL4GfACRHxcMXihcAOEbE78C3gkmrbiIgzImIkIkamTp1abIGtKzmPMH4e1GbNKDQoSJpMEhDOjYiLKpdHxMMRsSq9fxkwWdK2RZbJelOjPIIDQnWlpHPl9B8lzi9YpSJ7Hwn4AbA4Ir5WY51np+shac+0PA8UVSbrXR6UNTH+/CyvImsKewHvBvbNdDk9SNKxko5N1zkcuF3SLcA3gSOj1y4FZ4Ua9MnuWsWD2iyvwrqkRsS1ULPWWlrn28C3iyqD9bbKQWqlye5Kcxv5kpP5eVCb5eVpLqxrNZrszgew8fGgNsvDQcG6lgdfFcOfq9XjoGBdy8nRYvhztXocFKzrOLlcLCedrR7Pkmpdxcnl4jnpbPW4pmBdxcnl9nDS2WpxULCu4iRoe/nztkoOCtZVnARtL3/eVslBwbqCk8ud4aSzVXKi2TrOyeXOcdLZKrmmYB3n5HJnOelsWQ4K1nFOdnYHfw8GDgrWBZzs7A7+HgwcFKyDnFzuLk46GzjRbB3i5HL3cdLZwDUF6xAnl7uTk87moGAd4aRmd/P3M7gcFKytSnmEWtdcdVKzO9T6HgKcX+hzDgrWNqU8wvIaZ5tOLnePaknnklJ+wYGhPzkoWNtUyyOUDE8Z4tTDdnMuoUvMnjXMqYftxnCNGoPzC/3LvY+sbWq1Rwu4bs6+7S2MNTR71jCzZw2z45xLqzb3Ob/Qn1xTsLbx4Kje5O9tsDgoWOE8SK23VcsvTN5ErH58HTvOudSJ5z7j5iMrlAep9b7soLYVK9ewzdBkHn18HQ+uXgt4YFu/KaymIGm6pCslLZZ0h6Tjq6wjSd+UtFTSrZJeUlR5rDM8SK0/lAa1/eW0N7HFZpuydv2TswxOPPePImsK64ATI2KhpK2ABZKuiIjfZ9Y5ENgpvb0c+M/0r/UJD4LqP/5O+1thNYWIuCciFqb3HwEWA5WnhYcC50TiemCKpO2KKpO1n5OU/cffaX9rS6JZ0gxgFnBDxaJh4O7M42VsHDiQdIykUUmjY2NjRRXTWsjJ5f7l2VT7W+FBQdKWwM+AEyLi4crFVV6yUZfoiDgjIkYiYmTq1KlFFNNaqHLkcim5DB6k1g8qB7ZVm03VgaF3FRoUJE0mCQjnRsRFVVZZBkzPPN4eWFFkmax4Ti73P8+m2r+K7H0k4AfA4oj4Wo3V5gLvSXshvQJ4KCLuKapM1h5ORA4Of9f9p8iawl7Au4F9JS1KbwdJOlbSsek6lwF3AkuB7wMfKrA8VjDPgDp4PJtq/ymsS2pEXEv1nEF2nQD+ragyWPtUDlKr5ORyfzpp/5k1v3cPautNnubCWsIzoA4mz6bafzzNhbWEZ0AdXJ5Ntb+4pmAt4QFN5t9Af3BQsAnxIDUr8aC2/uDmI2uaZ0C1rOxsqqWThMpBbdn1rDu5pmBN8yA1q+RBbb3PQcGa5oFLVot/G73LQcHGpZRD2HHOpWyi6sNQnFg0D2rrXbmCgqSnF10Q637Zie4CWB8bd0B0ctmgetK5xJPmdbe8NYUbJP0knaai7ihl61+1BqhNkhAepGYbeFBb78rb++gFwH7A+4BvSboQOCsi/lhYyazr1GoPfiKCv5z2pjaXxrqdB7X1plw1hfTKaFdExDuADwBHATdKukrSKwstoXWcJ7qziXB+obfkzSk8Q9LxkkaBjwEfBrYFTgTOK7B81mGVF8yp5ByCNeL8Qm/Jm1P4HbA1MDsi3hQRF0XEuogYBb5XXPGs0zzRnU2U8wu9JW9Q+HREfDEilpWekPQ2gIj4aiEls67QaKI7BwTLozSorVYvFecXukfeoDCnynMnt7Ig1l2cR7AiOL/Q/er2PpJ0IHAQMCzpm5lFWwPriiyYdY4vmGNF8UV5ul+jmsIKYBT4J7Agc5sL7F9s0axTnEewoji/0P3q1hQi4hbgFknnRoRrBgPCF8yxInn8QnerW1OQ9OP07s2Sbs3cbpN0axvKZ23kPIK1k/ML3anRiObj078HF10Q6yznEazdnF/oTnVrChFxT3r3fuDuiPgrsBmwO0m+wfqE8wjWbs4vdKe8XVKvBjaXNAz8GngvcFZRhbL2yV5OsxqPR7AiNRq/4Et5tl/eoKCIWA0cBnwrIt4C7FxcsawdGk1hAc4jWHvU+515Koz2yh0U0onv3gVcmj7XaIzDmZLuk3R7jeX7SHpI0qL0dkr+Ylsr1GsyAucRrH3qzY8Ebkpqp7xTZx9PMoL54oi4Q9JzgSsbvOYs4NvAOXXWuSYinMTukHpd/4anDHHS/jPdbGRtUfqdnT5vSc2aq7uqtkfeqbOvjohDSvMcRcSdEfGRRq8B/tGCMlqLNep6OjxlyHkEa7tSfqFW4tldVdsj79TZL5B0hqTLJc0v3Vqw/1dKukXSLyTtUmf/x0galTQ6NjbWgt0OLk+Fbd3OU213Vt7mo5+QTJH9X0DtRujxWQjsEBGrJB0EXALsVG3FiDgDOANgZGSk1gmu5dCo66mbjKzTGjUllfIL/p0WQ1Hl4usbrSQtiIiXjnvj0gzg5xGxa4517wJGIuL+euuNjIzE6OjoeIsy8C65eXnd9lqBL6lpXafWVBjgk5jxSo/jI43Wy9v76H8kfUjSdpKeXrpNsIDPlqT0/p5pWR6YyDatOnc9tV7lrqrtlzcoHAWcBPyWDTOl1j1dl3Q+yRXbZkpaJun9ko6VdGy6yuHA7ZJuAb4JHBl5qi02bu56ar3KXVXbL1fzUTdx81F+jZqMwFVw637+HbdGS5uPJD1V0qclnZE+3kmSxxd0sTxNRu56ar2gUVdVcFNSK+VtPvoh8DjwqvTxMuBLhZTIJqQ0BuGECxe5ycj6Sp6mpBMuXOSxDBOUt0vq8yLiCEnvAIiINaUksXWPRtNfl7iqbb0oz6hn8LTbE5W3pvC4pCGSQYVIeh7wWGGlsqY0SiiDm4yst+VpSgInoCcib1D4HPBLYLqkc0mmz/5EUYWy8Wk0/XWJm4ysXzRqSoKkxrDjnEvdnDROuZqPIuJySQuAV5CMczq+0SAzaw83GdkgytuUFLg5abzyjmj+dUS8vtFz7eAuqYk83fQgqR34qmnWz/KeGMFgnxzl7ZLa6JoImwNPBbaV9DQoXyBpa2DahEtpTXHtwGyDbK1hxco1NafFANca8qhbU5B0PHACSQBYzoag8DDw/Yj4duElrDDINYW8tQPYkFA2GzR58msweCdNLRm8FhHfiIgdgY9FxHMjYsf0tnsnAsIgyzMYrcQJZRtkeZLQ4AFvteSe5kLSq4AZZJqcIqLeVdUKMWg1hfHUDmDwzn7MqvH/zcZaklPIbOxHwPOARWy4nkJQ/1KbNkHjSaA5oWy2wexZw8yeNZz7f8i5hg3y9j5aDOzcDbOY9ntNoXSGs2LlGjaRWJ/jIx+EsxyzZrnWkMhbU8gbFH4CfCQi7mlF4SaiH4NC9kcrqNt7Isu1A7P8xlPzLv0f9lOAaGnzEbAt8HtJN5KZ3iIiDmmyfJaq/KHmDQj99GM1a4e8A95gw//hIDYr5a0pvLba8xFxVctL1EC/1BTGW6Utce3AbOLGU2so6fUTsZY2H3WTXg4KzTYTTZJ4IoJpPf6jNOsmzZyY9XKzUkuCgqRHqH7sEhARsXXzRWxOrwWFZgNBiWsGZsVqptYAvRcgXFPooIkGgl77sZn1ukH4n3VQaLOJ/qhKuvlHZTYIms33lXRrgHBQKEh2HME2Q5OR4MHVaycUCMDNRGbdptlmpazScWFKeqxYuXptx3KDDgoTVNTBP6tbzyjMLNGqFoBK2f/9171wKlf+YYwVK9cUGjAcFHJqx8E/y4HArDcVFSAqVatdbNOCmoaDQg2VQeDRx9exdn2xn4EDgVl/aVeAqKWZ5uaOBwVJZwIHA/dFxK5Vlgv4BnAQsBo4OiIWNtruRIJCK9oI83IgMBsMnQoQ471mSqunuWjGWcC3qT2T6oHATunt5cB/pn8Lc/q8JYUGBAcCs8FTmpEV2hsgVjTZO6qRwoJCRFwtaUadVQ4FzklnXr1e0hRJ2xU56V6rPsRu6lFgZt2jWoAoKl85bcpQC7aysSJrCo0MA3dnHi9Ln9soKEg6BjgG4DnPeU7TO5w2ZaipIe0++JvZeGUDRFY2WEzL9D4aT+2iyKsrdjIoqMpzVT+PiDgDOAOSnEKzOzxp/5kb5RQmbyK23HzTlmX4zczqqRUsoHrtot3Hpk4GhWXA9Mzj7YEVRe4wO3Vu0X2CzczGq17AaJdOBoW5wHGSLiBJMD/Ujov4dMOHbmbWrQoLCpLOB/YBtpW0DPgsMBkgIr4HXEbSHXUpSZfU9xZVFjMzy6fI3kfvaLA8gH8rav9mZjZ+m3S6AGZm1j0cFMzMrMxBwczMyhwUzMyszEHBzMzKHBTMzKzMQcHMzMocFMzMrMxBwczMyhwUzMyszEHBzMzKHBTMzKzMQcHMzMocFMzMrMxBwczMyhwUzMyszEHBzMzKHBTMzKzMQcHMzMocFMzMrMxBwczMyhwUzMyszEHBzMzKHBTMzKys0KAg6QBJSyQtlTSnyvKjJY1JWpTePlBkeczMrL5Ni9qwpEnAd4A3AMuAmyTNjYjfV6x6YUQcV1Q5zMwsvyJrCnsCSyPizoh4HLgAOLTA/ZmZ2QQVGRSGgbszj5elz1V6q6RbJf1U0vRqG5J0jKRRSaNjY2NFlNXMzCg2KKjKc1Hx+H+AGRHxYuBXwNnVNhQRZ0TESESMTJ06tcXFNDOzkiKDwjIge+a/PbAiu0JEPBARj6UPvw+8tMDymJlZA0UGhZuAnSTtKOkpwJHA3OwKkrbLPDwEWFxgeczMrIHCeh9FxDpJxwHzgEnAmRFxh6QvAKMRMRf4iKRDgHXAP4CjiyqPmZk1pojKZv7uNjIyEqOjo50uhplZT5G0ICJGGq3nEc1mZlbmoGBmZmUOCmZmVuagYGZmZQ4KZmZW5qBgZmZlDgpmZlbmoGBmZmUOCmZmVuagYGZmZQ4KZmZW5qBgZmZlDgpmZlbmoGBmZmUOCmZmVuagYGZmZQ4KZmZW5qBgZmZlDgpmZlbmoGBmZmUOCmZmVuagYGZmZQ4KZmZW5qBgZmZlhQYFSQdIWiJpqaQ5VZZvJunCdPkNkmYUUY5Lbl7OXqfNZ8c5l7LXafO55OblRezGzKznFRYUJE0CvgMcCOwMvEPSzhWrvR94MCKeD/wH8NVWl+OSm5dz8kW3sXzlGgJYvnINJ190mwODmVkVRdYU9gSWRsSdEfE4cAFwaMU6hwJnp/d/CrxeklpZiNPnLWHN2vVPem7N2vWcPm9JK3djZtYXigwKw8DdmcfL0ueqrhMR64CHgGdUbkjSMZJGJY2OjY2NqxArVq4Z1/NmZoOsyKBQ7Yw/mliHiDgjIkYiYmTq1KnjKsS0KUPjet7MbJAVGRSWAdMzj7cHVtRaR9KmwDbAP1pZiJP2n8nQ5ElPem5o8iRO2n9mK3djZtYXigwKNwE7SdpR0lOAI4G5FevMBY5K7x8OzI+IjWoKEzF71jCnHrYbw1OGEDA8ZYhTD9uN2bMqW7LMzGzTojYcEeskHQfMAyYBZ0bEHZK+AIxGxFzgB8CPJC0lqSEcWURZZs8adhAwM8uhsKAAEBGXAZdVPHdK5v4/gbcVWQYzM8vPI5rNzKzMQcHMzMocFMzMrMxBwczMytTiHqCFkzQG/LXJl28L3N/C4nSS30t36pf30i/vA/xeSnaIiIajf3suKEyEpNGIGOl0OVrB76U79ct76Zf3AX4v4+XmIzMzK3NQMDOzskELCmd0ugAt5PfSnfrlvfTL+wC/l3EZqJyCmZnVN2g1BTMzq8NBwczMygYuKEj6oqRbJS2SdLmkaZ0uU7MknS7pD+n7uVjSlE6XqVmS3ibpDklPSOq57oOSDpC0RNJSSXM6XZ5mSTpT0n2Sbu90WSZK0nRJV0panP62ju90mZohaXNJN0q6JX0fny90f4OWU5C0dUQ8nN7/CLBzRBzb4WI6CT4bAAAFOElEQVQ1RdIbSa5BsU7SVwEi4hMdLlZTJL0IeAL4v8DHImK0w0XKTdIk4I/AG0guHHUT8I6I+H1HC9YESa8BVgHnRMSunS7PREjaDtguIhZK2gpYAMzute8lvW79FhGxStJk4Frg+Ii4voj9DVxNoRQQUltQ5fKfvSIiLk+vbQ1wPcnV7XpSRCyOiCWdLkeT9gSWRsSdEfE4cAFwaIfL1JSIuJoWX/2wUyLinohYmN5/BFjMxteJ73qRWJU+nJzeCjtuDVxQAJD0ZUl3A+8CTmm0fo94H/CLThdiQA0Dd2ceL6MHDz79TNIMYBZwQ2dL0hxJkyQtAu4DroiIwt5HXwYFSb+SdHuV26EAEfGpiJgOnAsc19nS1tfovaTrfApYR/J+ulae99KjVOW5nq2B9htJWwI/A06oaCnoGRGxPiL2IGkN2FNSYU17hV55rVMiYr+cq54HXAp8tsDiTEij9yLpKOBg4PWtvr51q43je+k1y4DpmcfbAys6VBbLSNvgfwacGxEXdbo8ExURKyX9BjgAKKQzQF/WFOqRtFPm4SHAHzpVlomSdADwCeCQiFjd6fIMsJuAnSTtKOkpJNcan9vhMg28NEH7A2BxRHyt0+VplqSppZ6FkoaA/SjwuDWIvY9+Bswk6enyV+DYiFje2VI1R9JSYDPggfSp63u4J9VbgG8BU4GVwKKI2L+zpcpP0kHA14FJwJkR8eUOF6kpks4H9iGZovnvwGcj4gcdLVSTJL0auAa4jeT/HeCT6bXje4akFwNnk/y2NgF+HBFfKGx/gxYUzMystoFrPjIzs9ocFMzMrMxBwczMyhwUzMyszEHBzMzKHBSs70h6tqQLJP1Z0u8lXSbpBQXs53OSPpbe/4KkpgbnSdoj7dJq1nF9OaLZBlc6YOli4OyIODJ9bg/gWSQzmTa73UkRsb7W8oiYyBxaewAjQE/1n7f+5JqC9ZvXAWsj4nulJyJiUURco8Tp6XxLt0k6ApJAUuP5fdL5+M8jGQCFpE+l1034FckgSNLnz5J0eHr/Lkmfl7Qw3d4L0+f3lPRbSTenf2emI6C/AByh5BofR0jaQsl1DW5K191obihJ20m6On3N7ZL2LuwTtYHimoL1m11J5s2v5jCSs/LdSUbs3iTpauBVNZ6HZFrsXSPiL5JeSjKFxSyS/52FdfZ1f0S8RNKHgI8BHyCZmuA16fUv9gO+EhFvlXQKMBIRxwFI+grJdTLel05vcKOkX0XEo5ntvxOYFxFfVnI9h6eO72Myq85BwQbJq4Hz02agv0u6CnhZnecfBm6MiL+kr98buLg0z5SkevMblSZfW0ASjAC2Ac5O598Kknnxq3kjcEgpXwFsDjyH5HoAJTcBZ6YTvl0SEYsav32zxtx8ZP3mDuClNZZVm+K63vMAj1Y8zjsvzGPp3/VsOPn6InBlekWzN5Mc7GuV560RsUd6e05EZANC6WI4rwGWAz+S9J6c5TKry0HB+s18YDNJ/1J6QtLLJL0WuJqk7X6SpKkkB9Ub6zxf6WrgLZKGlFze8c3jLNs2JAdxgKMzzz8CbJV5PA/4cJo0R9Ksyg1J2gG4LyK+TzIT6EvGWRazqhwUrK+k15R4C/CGtEvqHcDnSK5vcDFwK3ALSfD4eETcW+f5ym0vBC4EFpHM0X/NOIv378Cpkq4jmfGy5Epg51KimaRGMRm4VdLt6eNK+wCLJN0MvBX4xjjLYlaVZ0k1M7My1xTMzKzMQcHMzMocFMzMrMxBwczMyhwUzMyszEHBzMzKHBTMzKzs/wN6xcFfQ2a56wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = 100\n",
    "tau = h * A\n",
    "\n",
    "# reset the intensity values\n",
    "I = np.zeros(N)\n",
    "I[0] = I_init\n",
    "\n",
    "for i in range(1, N):\n",
    "    B_avg = B(B0, B1, (s[i] + s[i-1])/2, w) # can change this to use the average between two neighbouring grid points\n",
    "    I[i] = I[i-1] * np.exp(-tau) + (1 - np.exp(-tau)) * B_avg\n",
    "    \n",
    "plt.plot(s, I, 'o')\n",
    "plt.xlabel(\"Coordinate s\")\n",
    "plt.ylabel(\"Intensity\")\n",
    "plt.title(\"Piecewise Analytic Solution for A = \" + str(A))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question 8, 9, 10\n",
    "\n",
    "However, first order integration does not account for radiative diffusion and can produce wrong results for scattering problems in the regimes of high optical depth. To tackle this problem, let's consider linear interpolation of $B$ within a cell, using optical depth as the coordinate: $ d\\tau = \\alpha(s) ds$ such that the FTE becomes\n",
    "\n",
    "$$ \\frac{dI(\\tau)}{d\\tau} = S(\\tau) - I(\\tau) $$\n",
    "\n",
    "with the interpolated values for the source:\n",
    "\n",
    "$$ S(\\tau) = \\left(1 - \\frac{\\tau}{\\Delta\\tau}\\right) S_i + \\left(\\frac{\\tau}{\\Delta\\tau}\\right) S_{i+1} $$\n",
    "\n",
    "where $\\Delta\\tau$ is the grid spacing. The iterative formula for $I$ won't be typed here but directly implemented below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'OK87 for A = 100')"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi40LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcv7US4rQAAIABJREFUeJzt3Xu4XHV97/H3h7CBDQJBkkeTTUKAIj2ImMguXlBLvYHUQkAUPLaAl6ZeUPAoFbTlZhWUPp5W5EjTQgWLigrkSSsa0VCjHrnkxiWE2HjhkAQ1XMJFdiUJ3/PHWjMMw6yZtWfP2msun9fzzJOZtdbMfGdn7/Wd3+/7+/2WIgIzMzOAHcoOwMzMuoeTgpmZVTkpmJlZlZOCmZlVOSmYmVmVk4KZmVU5KZjVUOJfJT0i6bay4zGbbE4KVjpJp0m6S9KTkn4t6UuSptbsP1/Sv9U8HpF0r6QvpCfx50u6VtKD6e0aSXukx86W9ETdLSR9NCOcVwNvBPaJiMM7+Bl3S9/7xk69ZoP3WChpnaSnJZ3WYP9H0p/vo5KulLRzzb45km5O/w/ulfSGouK07uakYKVKT86fBc4C9gReAewL3CRppwbH7wssAxZHxIcjmX35d8BewP7AAcALgPMBIuL/RcTzKjfgJcDTwHUZIe0L/CoiftfGZ9mxye4Tgd8Db5I0Y7yvndMdwAeAlQ1iOwo4G3g9MIfkZ3VBzSFfA1YBewOfBL4laXpBcVoXc1Kw0qTf5i8APhQR342IrRHxK+DtJCfnP687/gCShPDViPjrml37AYsi4rGIeBS4AXhxxtueAixL36c+nvcA/wK8Mv1Wf0G6/S8lrZf0sKTFkmbWPCckfVDSfwH/1eTjngpcDtwJvLPJcW2LiMsi4gfAf2e8/xURsSYiHgE+BZwGIOlFwMuA8yJiLCKuA+4C3lpEnNbdnBSsTK8CdgGur90YEU8A3yHpxqnYnyQh/FNE/G3d61wGvEXSXpL2IjmZfSfjPU8Brmq0IyKuAN4H/DRtWZwn6XXARSSJagZwH/D1uqfOB14OHNzodSXNBo4Erklvp2TEVjn+TklbMm7/p9lzm3gxSUui4g7gBZL2Tvf9IiIer9uflVitjzVr7poVbRrwYERsa7DvAeCwmseHkHT7XNvg2JXATsBD6eMfAM85eUp6DUnX0rfGEeM7gSsjYmX6GucAj0iaU9PauCgiHm7yGqcAd0bEPZK2AJ+TNC8iVjU6OCIOHUd8eT0PeLTmceX+7g32VfaPFBCHdTm3FKxMDwLTMvriZ6T7KxYDVwJL07pCrW8CPyM5we0B/Bz4N57rVOC6tCWS10yS1gFQbcU8xLNPmPe3eI1TSFoIRMQm4IdpLJPpCZKfTUXl/uMN9lX2P44NHCcFK9NPSYqvJ9RulLQb8GaSb/xVEfG/gP8gSQy1J+WXknQr/S49aV8OHFP3msPA28joOmpiE0l9oza2vYGNtaFlPVnSq4ADgXPSkT+/JulqekdWYVrSmgYjpiq3y8cZf8Uakp9TxUuB30TEQ+m+/SXtXrd/TZvvZT3MScFKkxaFLwAulXS0pCFJc0i++W8AvtLgaacDS4EfSHpBuu124L2ShtOT/wKe3X8OcDywBbh5nGF+FXiXpLnpEM7PALc2KlRnOBW4iaTeMDe9HQLsSpL4niMiXlw7Yqru9r6sN5K0k6RdAAFDknaRVPkbvxp4j6SD07rL3wBfTt/vZ8Bq4Lz0OccDh5I9Qsv6mJOClSoiPgd8Avh74DHgVpLumNdHxO8bHB/AXwG3Ad+XNA14N8kwyw0k3+D3Jx1ZU+NU4OoY5wVE0tE8f0tygnyAZMjryXmem56g3w5cGhG/rrn9kiThdboL6XvAGEkBf2F6/7Xp5/gu8DmSpHhfejuv5rknA6PAI8DFwIkRsbnD8VkPkC+yY2ZmFW4pmJlZlZOCmZlVOSmYmVmVk4KZmVX13IzmadOmxZw5c8oOw8ysp6xYseLBiGi5yGHPJYU5c+awfPnyssMwM+spku5rfZS7j8zMrIaTgpmZVTkpmJlZlZOCmZlVOSmYmVlVYaOP0sXAlgE7p+/zrYg4r+6YnUlWbzyMZI36k8ax+qRZ6Rat2sglS9axacsYew4PIcGWJ7cyc+owZx11EPPn+To11lsKWxBPkoDdIuIJSUPAj4EzIuKWmmM+ABwaEe+TdDJwfESc1Ox1R0dHw0NSrUyVRLBxyxgi+2IKlX0jThDWBSStiIjRlsdNxiqpknYlSQrvj4hba7YvAc6PiJ+mFxz5NTC92fLGTgpWhryJIIsThJWtK5KCpCnACuAPgMsi4uN1++8Gjo6IDenjnwMvj4gH645bQHLhFGbPnn3YffflmoNhNiETTQRZnCCsDHmTQqGF5ojYHhFzgX2AwyUdUneIGj2twessjIjRiBidPr3lLG2zCVu0aiPnXH8XG7eMAZ1LCLWvtXHLGOdcfxeLVm1serzZZJqU0UcRsQX4T+Doul0bgFkAaffRnsDDkxGTWSOLVm3kiIuXcua1qxnbur3w9xvbup0zr13NERcvdXKwrlBYUpA0XdLU9P4w8Abg3rrDFvPMJQlPBJaO93KJZp1S3zrIq9LcnTo8xF67Dj1rW15uNVi3KHJBvBnAVWldYQfgGxHxH5IuBJZHxGLgCuArktaTtBByXfvWrAiXLFmXu3XQqi7QTj1ibOt2PvqNO/jItas9pNVK03PXaPboI+u02hN4M+0WiNstWA8PTeGiE17ixGAdkbfQ3HNLZ5t1UqXLqFULYSIjhebPG6k+L28CgqTlcMmSdU4KNqncUrCBdsTFS5ueoIv6tp43GYGHrlpnuKVg1kSeb+xFnowrr1lZImMHie0ZX9AqReja55kVxUnBBk6eb+kjU4f5ydmvKzSO+m6lZjG5K8kmi1dJtYHTapTR8NAUzjrqoEmMKEkQF53wEkamDmces3HLmOczWOGcFGxgVCamteoyKmvEz/x5I/zk7Ne1TAyez2BFclKwgZBnYlqly6jsLpqzjjqI4aEpmfsrXUlmRXBSsIHQjV1GWdyVZGVyodn6Vu0FcJoNvO7GIZ+VInSz7i6PSrIiuKVgfam2u6hVQuiGLqMs7kqyyeaWgvWlPOsYdVOXUZba+QxZLYZN41zAz6wZtxSsLzU7UYpyRxmNV6tRSQGuL1jHuKVgfaVSR8jqMpqMSWlFOeuogzInuLm+YJ3iloL1jVbDTnuhu6iZVqOSXF+wTnBSsL7RrI7QS91FzVS6krIu4uP6gk2Uk4L1jawToqCrRxi1Y6brC1YQJwXreZXlK7LqCFkn0F7WbKiql8KwiXBSsJ7W73WELK4vWFGcFKynDUIdIYvrC1YEJwXraYNUR8ji+oJ1kpOC9aRBrCNkcX3BOslJwXrOoNYRsri+YJ3kpGA9Z5DrCFlcX7BOcVKwnuM6QjbXF2yiCksKkmZJulnSWklrJJ3R4JgjJT0qaXV6O7eoeKx/ZJ34BqmOkMX1BZuoIhfE2wZ8NCJWStodWCHppoi4p+64H0XEWwqMw/pEZbG7jVvGEDyryDxodYQsrZbartQXBrk1Zc0V1lKIiAciYmV6/3FgLeDfRGtLfXE5oNp/Pqh1hCyuL9hETMrS2ZLmAPOAWxvsfqWkO4BNwMciYk2D5y8AFgDMnj27uECtazUqLge9vRR20WZOHW7YWnA3mzVTeKFZ0vOA64AzI+Kxut0rgX0j4qXApcCiRq8REQsjYjQiRqdPn15swNaVsr7d+ltvtkb1BZHUFlx0tiyFJgVJQyQJ4ZqIuL5+f0Q8FhFPpPdvBIYkTSsyJustnqTWvvr5C7V1GBedLUuRo48EXAGsjYjPZxzzwvQ4JB2exvNQUTFZb/EktYmrvZRnfWL1pDZrpMiawhHAXwB3SVqdbvsEMBsgIi4HTgTeL2kbMAacHBFZXwptwLSapHbWUQe5uJyTu98sr8KSQkT8GDIHQFSO+SLwxaJisN7WapKa5eeis+XlGc3WtTxJrXNcdLa8JmVIqtl4eJJa59VPamtUdK49zgaXWwrWVTxJrTguOlsebilYV/EkteK56GzNuKVgXcUnrOJ5JVVrxknBuoqLy8XzSqrWjJOCdYXKzOVKEbSWi8ud5Su1WTNOClY6F5cnn1dStSwuNFvpXFwujye1WT23FKx0Li6Xx5ParJ6TgpXOxeXyeCVVq+ekYKVxcbk7eFKb1XJNwUpRKS5XagmV4nKlluAVUCefu/EMnBSsJC4udx8XnQ3cfWQl8bfS7uOis4GTgpXExeXu46KzgZOCTTIXl7ubi87mmoJNGheXe4e79waXk4JNGheXe4eLzoPL3Uc2afzts3e46Dy43FKwwlUur1nfR13hb5/dx5fvHFxuKVih6ldAreficvdy0XkwuaVghWpUR6hwcbk3uNtvsDgpWKGyThwCF5d7hIvOg6Ww7iNJsyTdLGmtpDWSzmhwjCR9QdJ6SXdKellR8Vg5PEmt9zUqOg/tIJ58ahv7nf1tF577TJEthW3ARyNipaTdgRWSboqIe2qOeTNwYHp7OfCl9F/rcZXicn2RElxH6DW1RedNW8bYc3iI3z21jUee3Aq48NxvCmspRMQDEbEyvf84sBao/405Drg6ErcAUyXNKCommxy+vGb/qRSdf3nxn7LbzjuydfuzS88uPPePSakpSJoDzANurds1Atxf83hDuu2BuucvABYAzJ49u6gwrUM8Sa2/ufDc3wofkirpecB1wJkR8Vj97gZPec5w9ohYGBGjETE6ffr0IsK0DvJJo79l1YMCXF/oA4UmBUlDJAnhmoi4vsEhG4BZNY/3ATYVGZMVz8Xl/tao8Fzh1VR7X5GjjwRcAayNiM9nHLYYOCUdhfQK4NGIeCDjWOtyXgF1MNQvsV3P9YXeVmRN4QjgL4C7JK1Ot30CmA0QEZcDNwLHAOuBJ4F3FRiPFcgroA6W+fNGmD9vhP3O/nbD5UvcVdi7CksKEfFjGtcMao8J4INFxWCTx8XlweSJbf3Hax9ZR7i4PJi8mmr/8TIXNiFeAXWweTXV/uOWgrXNK6AaeDXVfuOWgrXNK6BaLXch9ge3FKxtrVZAdUIYLJ7U1h+cFKxtnqRmtTyprT+4+8jGpVJYrqyWOTRFz1oczXWEwVVfdK5XqS+4Bdnd3FKw3GoLywFsGdsKAXvtOoTwCqj2TNE5a4KS6wvdL1dLQdLzI+LhooOx7taosLz16WDXnXZk1blvKikq60ZZk9oq9QUPQuheeVsKt0r6pqRj0jWNbAB5dInl5fpC78qbFF4ELCRZy2i9pM9IelFxYVk3qSx05wlqlpcXzetduZJCemW0myLiHcB7gVOB2yT9UNIrC43QSuUJatYu1xd6U96awt7An5O0FH4DfIhk2eu5wDeB/YoK0MrlCWo2UV40r7fk7T76KbAHMD8i/jQiro+IbRGxHLi8uPCsbJ6gZhPlRfN6S955Cn8TEd+o3SDpbRHxzYj4bAFxWcm80J11ihfN6y15WwpnN9h2TicDse7hOoJ1mhfN6x1NWwqS3kxyZbQRSV+o2bUHsK3IwKw8riNYUbK6IytdSf7dKl+r7qNNwHLgWGBFzfbHgY8UFZSVq1UdwaxdWUVncFdSt2jafRQRd0TEVcABEXFVze36iHhkkmK0SeL5CFa0ZpPawF1J3aBV99E3IuLtwCpJtecKkUxfOLTQ6GzSVOoIWd1GriNYJ7RaNA88f6Fsisj6XgiSZkTEA5L2bbQ/Iu4rLLIMo6OjsXz58sl+2753xMVLM/9IXUewIvh3bnJJWhERo62Oa9V99EB690Hg/jQJ7Ay8lKTeYD2u0mWU9cfp+QhWFK+P1J3yDkldBuwiaQT4AfAu4MtFBWWTo9XQU3AdwYrj9ZG6U96koIh4EjgBuDQijgcOLi4smwzNhp6C6whWPK+P1H1yJ4V04bt3At9Ot7UqUl8p6beS7s7Yf6SkRyWtTm/n5g/bOqHZH5wvmGOTydd37h55k8IZJDOYb4iINZL2B25u8ZwvA0e3OOZHETE3vV2YMxaboFZDT0emDruOYJPK9YXukWvto4hYRlJXqDz+BfDhVs+RNGciwVnneeipdSNf37l75GopSHqRpIWSvidpaeXWgfd/paQ7JH1H0oubvP8CScslLd+8eXMH3nZwtVrCwl1GVpZW9QWvqjo58q6S+k2SJbL/BciuTI7PSmDfiHhC0jHAIuDARgdGxEKSK78xOjqaPbHCMlVWPW019NSsbF4Ko1x5awrbIuJLEXFbRKyo3CbyxhHxWEQ8kd6/ERiSNG0ir2mNeeip9RIvhVGuvEnh3yV9QNIMSc+v3CbyxpJeKEnp/cPTWB6ayGtaYx56ar2k1fwFcFdSkfJ2H52a/ntWzbYA9s96gqSvAUcC0yRtAM4DhgAi4nLgROD9krYBY8DJ0WzNDRu3Vl1G4OUErDvNnzfC/HkjTWfbuyupGE3XPupGXvson1ajjOCZoadm3cq/x53TkbWPal5sV0l/I2lh+vhASW+ZaJBWHHcZWT9wV9Lky1tT+FfgKeBV6eMNwN8VEpFNSKsF7sBDT6231F7KM4snuHVO3qRwQER8DtgKEBFjkDmc2EqSZ5SRZytbr8ozKunMa1e71TBBeQvNT0kaJikuI+kA4PeFRWXjkqegDO4yst6W5wI94AL0ROVtKZwPfBeYJekakuWzP15UUJZfntYBuMvI+kOeriTwXIaJyD36SNLewCtIuo1uiYgHiwwsi0cfJfK2DsCjM6z/5BmVBMnJaqaHXQP5Rx/l6j6S9IOIeD3PLJtdu80mWd4/CHCXkfWnvF1JgbuTxqtp95GkXdKZy9Mk7VUzm3kOMHMyArTnajXctMJdRtbPKl1J/3DS3KYFaHARejxatRT+CjiTJAGs4JkRR48BlxUYlzUwnoKyk4ENitpWw6YtY5nXCQG3GvLIVVOQ9KGIuHQS4mlp0GoKtYlA0PQXHrxshVmreToVg/a30tGaQkRcKulVwJza50TE1W1HaC3V1w6aJQS3DswSZx11UK6am1sNjeUtNH8FOABYzTPXUwjASaHDKi2DTVvG2EFie46W3KB94zFrJm8RGp6pNVyyZJ3/hlJ5u4/WAgd3wyqm/dx9NJ5RRRUebmqWbTx/U5Xu2X79ktXR7iPgbuCFwAMTisoaGs+cg1oebmrW3HhaDZVvvIPerZS3pXAzMBe4jZrlLSLi2OJCa6xfWgrjLSBX9Pu3GbOitNsS75e/s7wthbxJ4Y8bbY+IH7YR24T0Q1IY7y/nFImnIzwz02yC2mmV98sXsY4mhW7Sy0mhnV9Ijyoy67x2Wg3Q2wmiI0lB0uM07tkQEBGxR/shtqfXkkK73UTQm794Zr1iIn+b0HsJwi2FEk30l82tA7PJ1e5gj4peSBBOCpOgdk7BnsNDSPDIk1sH4luHWT9qt1upVrf+LTspFGSirYBGuu2Xx2yQdfJvvPL8qemXxi1Pbi1twIiTwgR1shWQxd1EZt2tiC+B8OzWxJ/84XRuvnczm7aMFZownBQy1J7sZ9b9hxR18q/VrU1LM2uuqARRr1HrYs8OtDScFBroRH9hO5wIzPrLZCWILO30MnR6mYtxk3Ql8BbgtxFxSIP9Av4ROAZ4EjgtIlYWFQ/kvzhNJzgRmPWv+fNGqn/TZSSIyjWoizivFJYUgC8DXyR7JdU3Awemt5cDX0r/LcymNoebNdNNhSQzm3xlJYgizmdQYFKIiGXpZTuzHAdcna68eoukqZJmRERhi+7NnDrc9jjkWm4FmFkjjRJEUfXKmVOHO/Aqz1VkS6GVEeD+mscb0m3PSQqSFgALAGbPnt32G+a9+Eb1fXErwMzaU5sgamUNdhlP66LIFZLLTApqsK3hzyMiFgILISk0t/uG9ddyzRp95JO/mRUlK1lA49ZFp0Yf5VVmUtgAzKp5vA+wqeg3bfYfYmZWpm44P+1Q4nsvBk5R4hXAo0XWE8zMrLUih6R+DTgSmCZpA3AeMAQQEZcDN5IMR11PMiT1XUXFYmZm+RQ5+ugdLfYH8MGi3t/MzMavzO4jMzPrMk4KZmZW5aRgZmZVTgpmZlblpGBmZlVOCmZmVuWkYGZmVU4KZmZW5aRgZmZVTgpmZlblpGBmZlVOCmZmVuWkYGZmVU4KZmZW5aRgZmZVTgpmZlblpGBmZlVOCmZmVuWkYGZmVU4KZmZW5aRgZmZVTgpmZlblpGBmZlVOCmZmVlVoUpB0tKR1ktZLOrvB/tMkbZa0Or29t8h4zMysuR2LemFJU4DLgDcCG4DbJS2OiHvqDr02Ik4vKg4zM8uvyJbC4cD6iPhFRDwFfB04rsD3MzOzCSoyKYwA99c83pBuq/dWSXdK+pakWY1eSNICScslLd+8eXMRsZqZGcUmBTXYFnWP/x2YExGHAt8Hrmr0QhGxMCJGI2J0+vTpHQ7TzMwqikwKG4Dab/77AJtqD4iIhyLi9+nDfwYOKzAeMzNrocikcDtwoKT9JO0EnAwsrj1A0oyah8cCawuMx8zMWihs9FFEbJN0OrAEmAJcGRFrJF0ILI+IxcCHJR0LbAMeBk4rKh4zM2tNEfXd/N1tdHQ0li9fXnYYZmY9RdKKiBhtdZxnNJuZWZWTgpmZVTkpmJlZlZOCmZlVOSmYmVmVk4KZmVU5KZiZWZWTgpmZVTkpmJlZlZOCmZlVOSmYmVmVk4KZmVU5KZiZWZWTgpmZVTkpmJlZlZOCmZlVOSmYmVmVk4KZmVU5KZiZWZWTgpmZVTkpmJlZlZOCmZlVOSmYmVmVk4KZmVUVmhQkHS1pnaT1ks5usH9nSdem+2+VNKeIOBat2sgRFy9lv7O/zREXL2XRqo1FvI2ZWc8rLClImgJcBrwZOBh4h6SD6w57D/BIRPwB8L+Bz3Y6jkWrNnLO9XexccsYAWzcMsY519/lxGBm1kCRLYXDgfUR8YuIeAr4OnBc3THHAVel978FvF6SOhnEJUvWMbZ1+7O2jW3dziVL1nXybczM+kKRSWEEuL/m8YZ0W8NjImIb8Ciwd/0LSVogabmk5Zs3bx5XEJu2jI1ru5nZICsyKTT6xh9tHENELIyI0YgYnT59+riCmDl1eFzbzcwGWZFJYQMwq+bxPsCmrGMk7QjsCTzcySDOOuoghoemPGvb8NAUzjrqoE6+jZlZXygyKdwOHChpP0k7AScDi+uOWQycmt4/EVgaEc9pKUzE/HkjXHTCSxiZOoyAkanDXHTCS5g/r74ny8zMdizqhSNim6TTgSXAFODKiFgj6UJgeUQsBq4AviJpPUkL4eQiYpk/b8RJwMwsh8KSAkBE3AjcWLft3Jr7/w28rcgYzMwsP89oNjOzKicFMzOrclIwM7MqJwUzM6tSh0eAFk7SZuC+Np8+DXiwg+GUyZ+lO/XLZ+mXzwH+LBX7RkTL2b89lxQmQtLyiBgtO45O8GfpTv3yWfrlc4A/y3i5+8jMzKqcFMzMrGrQksLCsgPoIH+W7tQvn6VfPgf4s4zLQNUUzMysuUFrKZiZWRNOCmZmVjVwSUHSpyTdKWm1pO9Jmll2TO2SdImke9PPc4OkqWXH1C5Jb5O0RtLTknpu+KCkoyWtk7Re0tllx9MuSVdK+q2ku8uOZaIkzZJ0s6S16e/WGWXH1A5Ju0i6TdId6ee4oND3G7SagqQ9IuKx9P6HgYMj4n0lh9UWSW8iuQbFNkmfBYiIj5ccVlsk/Q/gaeCfgI9FxPKSQ8pN0hTgZ8AbSS4cdTvwjoi4p9TA2iDptcATwNURcUjZ8UyEpBnAjIhYKWl3YAUwv9f+X9Lr1u8WEU9IGgJ+DJwREbcU8X4D11KoJITUbjS4/GeviIjvpde2BriF5Op2PSki1kbEurLjaNPhwPqI+EVEPAV8HTiu5JjaEhHL6PDVD8sSEQ9ExMr0/uPAWp57nfiuF4kn0odD6a2w89bAJQUASZ+WdD/wTuDcVsf3iHcD3yk7iAE1Atxf83gDPXjy6WeS5gDzgFvLjaQ9kqZIWg38FrgpIgr7HH2ZFCR9X9LdDW7HAUTEJyNiFnANcHq50TbX6rOkx3wS2EbyebpWns/So9RgW8+2QPuNpOcB1wFn1vUU9IyI2B4Rc0l6Aw6XVFjXXqFXXitLRLwh56FfBb4NnFdgOBPS6rNIOhV4C/D6Tl/futPG8f/SazYAs2oe7wNsKikWq5H2wV8HXBMR15cdz0RFxBZJ/wkcDRQyGKAvWwrNSDqw5uGxwL1lxTJRko4GPg4cGxFPlh3PALsdOFDSfpJ2IrnW+OKSYxp4aYH2CmBtRHy+7HjaJWl6ZWShpGHgDRR43hrE0UfXAQeRjHS5D3hfRGwsN6r2SFoP7Aw8lG66pYdHUh0PXApMB7YAqyPiqHKjyk/SMcA/AFOAKyPi0yWH1BZJXwOOJFmi+TfAeRFxRalBtUnSq4EfAXeR/L0DfCK9dnzPkHQocBXJ79YOwDci4sLC3m/QkoKZmWUbuO4jMzPL5qRgZmZVTgpmZlblpGBmZlVOCmZmVuWkYH1H0gslfV3SzyXdI+lGSS8q4H3Ol/Sx9P6FktqanCdpbjqk1ax0fTmj2QZXOmHpBuCqiDg53TYXeAHJSqbtvu6UiNietT8iJrKG1lxgFOip8fPWn9xSsH7zJ8DWiLi8siEiVkfEj5S4JF1v6S5JJ0GSSDK2H5mux/9VkglQSPpket2E75NMgiTd/mVJJ6b3fyXpAkkr09f7w3T74ZL+r6RV6b8HpTOgLwROUnKNj5Mk7abkuga3p8c+Z20oSTMkLUufc7ek1xT2E7WB4paC9ZtDSNbNb+QEkm/lLyWZsXu7pGXAqzK2Q7Is9iER8UtJh5EsYTGP5G9nZZP3ejAiXibpA8DHgPeSLE3w2vT6F28APhMRb5V0LjAaEacDSPoMyXUy3p0ub3CbpO9HxO9qXv9/Aksi4tNKruew6/h+TGaNOSnYIHk18LW0G+g3kn4I/FGT7Y8Bt0XEL9Pnvwa4obLOlKRm6xtVFl9bQZKMAPYErkrX3wqSdfEbeRNwbKVeAewCzCa5HkDF7cCV6YJviyJideuPb9aau4+s36wBDsvY12iJ62bbAX5X9zjvujC/T//dzjNfvj4F3Jxe0ezPSE7hFOHWAAABA0lEQVT2WfG8NSLmprfZEVGbECoXw3ktsBH4iqRTcsZl1pSTgvWbpcDOkv6yskHSH0n6Y2AZSd/9FEnTSU6qtzXZXm8ZcLykYSWXd/yzcca2J8lJHOC0mu2PA7vXPF4CfCgtmiNpXv0LSdoX+G1E/DPJSqAvG2csZg05KVhfSa8pcTzwxnRI6hrgfJLrG9wA3AncQZI8/joift1ke/1rrwSuBVaTrNH/o3GG9zngIkk/IVnxsuJm4OBKoZmkRTEE3Cnp7vRxvSOB1ZJWAW8F/nGcsZg15FVSzcysyi0FMzOrclIwM7MqJwUzM6tyUjAzsyonBTMzq3JSMDOzKicFMzOr+v+nC/DtSr+prgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = 100\n",
    "tau = h * A\n",
    "\n",
    "# reset the intensity values\n",
    "I = np.zeros(N)\n",
    "I[0] = I_init\n",
    "\n",
    "for i in range(1, N):\n",
    "    coef1 = (1 - (1 + tau) * np.exp(-tau))/tau\n",
    "    coef2 = (tau - 1 + np.exp(-tau))/tau\n",
    "    I[i] = I[i-1] * np.exp(-tau) + coef1 * B(B0, B1, s[i-1], w) + coef2 * B(B0, B1, s[i], w)\n",
    "    \n",
    "plt.plot(s, I, 'o')\n",
    "plt.xlabel(\"Coordinate s\")\n",
    "plt.ylabel(\"Intensity\")\n",
    "plt.title(\"OK87 for A = \" + str(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}