EX5_FORCE_BALANCE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dynamics of tectonic plates\n",
    "## Exercises on plate dynamics\n",
    "\n",
    "### Author(s), contact(s), and dates:\n",
    "- Author: Valentina Espinoza Fernández (University of Copenhagen)\n",
    "- Email:  vf@ign.ku.dk\n",
    "- Date:   22nd of January 2024"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5: Forces acting on a simplified subducting plate model\n",
    "\n",
    "As we have seen, there are different main forces acting on tectonic plates [Forsyth & Uyeda, 1975]. The relevant force equations for this exercise can be found in the appendix of the exercise document. In this case, we will explore the forces acting on a highly simplified 2-D plate model - A plate with constant thickness and density ($H_s = 80$ km, $\\rho_{plate} = 2900~\\mathrm{kg/m^3}$) that subducts beneath a mantle with equally uniform density throughout ($\\rho_{mantle} = 3100~\\mathrm{kg/m^3}$). The force enabling subduction is the net weight of the subducting plate ($NW$), but this is partially counteracted by the viscous resistance of the mantle ($R$). To calculate this resistance, a mantle viscosity of $1.5~\\cdot 10^{21}~\\mathrm{Pa \\cdot s}$ is assumed, along with an instantaneous plate velocity of $50~mm/year$ and a dissipation distance equivalent to $200~km$. Together, the force produced by the weight of the subducting plate minus the viscous resistance is known as the net slab pull ($NSP$).\n",
    "\n",
    "Perform the following exercises:\n",
    "\n",
    "1. Considering a subduction with an intermediate angle, i.e., $\\alpha = 45^{\\circ}$, calculate the plate's net weight force ($NW$), viscous resistance force ($R$), and net slab pull force ($NSP$) acting on the plate. Remember that the forces are expressed in Newtons ($N$), for which all physical magnitudes included in the equations must be expressed in units of the International System. The subduction zone is 1 kilometer long, and the deepest portion of the slab is at a depth ($z_{tip}$) of $200~km$. What is the order of magnitude of the calculated forces?\n",
    "\n",
    "2. With the same parameters as in item 1, but in this case with a vertical subduction angle ($\\alpha = 90^{\\circ}$), calculate $NW$, $R$, and $NSP$. How do these results compare to the magnitudes from the previous item?\n",
    "\n",
    "3. Now, consider different values for the subduction angle. For instance, consider values between $10^{\\circ}$ and $90^{\\circ}$. Before plotting, can you make a prediction of how the trend will look like? Graph and compare the force trend with your prediction. How do you expect other parameters to influence the magnitude of these forces?\n",
    "\n",
    "4. Include now the force arising from the friction along the subduction interface ($FS$). Assume a value of $0.3$ for the friction coefficient ($\\mu_f = 0.3$) and a depth of $50~km$ for the brittle-ductile boundary. Are the orders of magnitude as expected? From what we discussed in class, if you were to try out a more realistic friction coefficient value, would that value be lower or higher that $0.3$? Do you expect to force to compare better to the NSP force, or drift apart even more? \n",
    "\n",
    "5. With the forces at hand, do a force balance for the plate. If the residual is not close to zero (bare in mind the orders of magnitude we are dealing with), what other force could be missing to achieve a state of equilibrium? Considering the order of magnitude, of the missing forces, which do you think is the most relevant one (i.e. the one that can account for most of the residual)? \n",
    "\n",
    "\n",
    "\n",
    "<br>\n",
    "\n",
    "***\n",
    "References:  \n",
    "- D. Forsyth and S Uyeda. “On the relative importance of the driving forces on plate motion”. In: Geophysical Journal of the Royal Astronomical Society 43 (1975), pp. 163–200.\n",
    "\n",
    "\n",
    "<br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Public dependencies\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# General parameters\n",
    "g = 9.8         # Gravity [m/s2]\n",
    "rho_m = 3100.0  # Density of the mantle [kg/m3]\n",
    "rho_p = 2900.0  # Density of the plate [kg/m3]\n",
    "\n",
    "\n",
    "# Specific parameters\n",
    "d_rho = rho_m - rho_p            # Density contrast between a uniform mantle and the subducting plate [kg/m3]\n",
    "v = 50e-3 / (365.25 * 24* 3600)  # Velocity of the plate [m/s]\n",
    "\n",
    "Hs = 80e3         # Plate thickness [m]\n",
    "mu_m = 1.5e21     # Mantle uniform viscosity [Pa s]\n",
    "Dx = 200e3        # Distance over which stresses dissipate on the mantle [m]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simplified model start with the assumption that all the parameters are uniform with depth, i.e. they do not depend on the depth $z$. We will solve the force integrals (make sure you understand how) and write them into Python functions that we later use for the many calculations ahead.\n",
    "\n",
    "$$ NW = \\frac{1}{sin(\\theta)} \\cdot g \\cdot \\Delta \\rho \\cdot H_s \\cdot z_{tip} $$\n",
    "\n",
    "$$ R = \\frac{2}{sin(\\theta)} \\cdot \\mu \\cdot \\frac{v}{Dx} \\cdot z_{tip} $$\n",
    "\n",
    "Note the expected units for the paramaters of the equations below, specially the units for the angle `theta_slab`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Slab net weight (NW)\n",
    "def nw_force(theta_slab, g, d_rho, Hs, z_tip):\n",
    "    return 1/np.sin(np.radians(theta_slab)) * g * d_rho * Hs * z_tip\n",
    "\n",
    "# Subduction resistance (R)\n",
    "def r_force(theta_slab, mu_m, v, Dx, z_tip):\n",
    "    return 2/np.sin(np.radians(theta_slab)) * mu_m * v/Dx * z_tip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subduction with a 45° angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>NW</th>\n",
       "      <th>R</th>\n",
       "      <th>NSP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>45 deg</th>\n",
       "      <td>4.434974e+13</td>\n",
       "      <td>6.722059e+12</td>\n",
       "      <td>3.762768e+13</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  NW             R           NSP\n",
       "45 deg  4.434974e+13  6.722059e+12  3.762768e+13"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_slab_45 = 45     # Angle of the subducting slab [degrees]\n",
    "z_tip = 200e3          # Depth of the tip of the slab [m]\n",
    "\n",
    "forces_df = pd.DataFrame(columns=['NW', 'R', 'NSP'])    # We will create a dataframe to store the forces\n",
    "\n",
    "\n",
    "NW = nw_force(theta_slab_45, g, d_rho, Hs, z_tip)\n",
    "R = r_force(theta_slab_45, mu_m, v, Dx, z_tip)\n",
    "NSP = NW - R\n",
    "\n",
    "forces_df.loc[\"45 deg\"] = [NW, R, NSP]  # Add a row to the dataframe\n",
    "forces_df                               # Show the dataframe output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "37627678783179.86"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "forces_df[\"NSP\"][\"45 deg\"]  # You can access the value of a specific cell in the dataframe like so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These forces actually expressed in N/m (Newton per meter), and represent the force for each linear meter of convergent margin. The force acting on a plate is in reality the integral over the length of the boundary. For our simplified plate, one can simply multiply by the length of the margin in question."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Net slab pull force: 3.76277e+16 N\n"
     ]
    }
   ],
   "source": [
    "L_conv = 1e3    # Length of the convergent margin [m]\n",
    "NSP_f = forces_df[\"NSP\"][\"45 deg\"] * L_conv\n",
    "\n",
    "print(f\"Net slab pull force: {NSP_f:.5e} N\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subduction with a 90° angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try now to calculate the forces for a subduction with a 90° angle. Before you do, can you predict what will happen? Will the forces be larger or smaller than for the 45° case?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subduction for a range of angles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What about values between $10^{\\circ}$ and $90^{\\circ}$. Before plotting, try and make a prediction of how the trend will look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_angles = 30   # Number of angles to test\n",
    "theta_slab_array = np.linspace(10, 179, N_angles)    # Array of angles to test\n",
    "NW_array = [nw_force(theta_slab, g, d_rho, Hs, z_tip) for theta_slab in theta_slab_array]\n",
    "\n",
    "R_array = [r_force(theta_slab, mu_m, v, Dx, z_tip) for theta_slab in theta_slab_array]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABNoAAAK7CAYAAAAk8+KPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAByFElEQVR4nO3de5xUdf0/8PdwW1BhEQ1Y5CLeUEERQQUUlUwUL0lqonkBL6W/vEaWkpppGWpmaKZmX5XMRL/fQMWvWGKCaGCJgqYpXyoSxCXywi5i3M/vj30wMe6FhT27M8s+n4/HPGDOfM6Z95mZXd685pzPySRJkgQAAAAAUCfN8l0AAAAAAGwLBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbUCDmDBhQmQymWjdunW8++67lR4/8sgjo0+fPtn7ffr0iX322afSuMcffzwymUwMGjSo0mO/+tWvIpPJxJQpU9ItvhHJZDLxve99L99lAADk2NgLbry1aNEiSkpK4vTTT48FCxbku7y82/j6/OMf/8h3KUAdCdqABrV69eq49tprNztu6NCh8c4778TSpUtzls+YMSO23377mDNnTqxYsaLSY82aNYvDDz881ZoBAEjHgw8+GLNnz47nnnsuLrnkkpgyZUocdthh8fHHH+e7NIBUCNqABnXsscfGI488Eq+//nqN44YOHRoRFeHZpmbMmBEXXHBBZDKZeOmllyo91q9fv2jfvn2aJQMAkJI+ffrEwIED48gjj4xrrrkmrr766li2bFk88cQT+S4NIBWCNqBBffvb346ddtoprrrqqhrHHXnkkZHJZHKCtg8//DD+/Oc/x/HHHx/9+/eP6dOnZx9bvHhx/P3vf88GdFX517/+Fa1atYrrrruu0mPvvPNOZDKZuPPOOyMi4tNPP40rr7wyevbsGa1bt44OHTrEgAEDYuLEiTXW/a9//Su+/vWvx7777hs77LBDdOzYMT7/+c/Hiy++mDPuH//4R2Qymbjtttvi9ttvj549e8YOO+wQgwYNipdffrnSdn/xi1/EXnvtFUVFRbHvvvvGI488EqNHj45dd921xnoiIpYuXRoXXnhhdO3aNVq1ahU9e/aMG264IdatW7fZdQEA6tOAAQMiIuKf//xnjePGjx8fmUwm/vrXv1Z67KqrropWrVrFBx98EBERc+fOjRNOOCE6duwYRUVF0aVLlzj++OPjvffeq/E5pk2bFieddFJ07do1WrduHXvssUdceOGF2e1u9L3vfS8ymUy89dZbccYZZ0RxcXF06tQpzjvvvCgrK8sZu3z58jj//POjQ4cOscMOO8Txxx8ff//732s93cdzzz0XRx11VLRr1y622267OPTQQ+P3v//9ZtcD8kfQBjSotm3bxrXXXhu/+93v4vnnn692XIcOHWL//ffPCdNeeOGFaN68eQwePDiOOOKInBBu47iagrbPfe5zccIJJ8Qvf/nL2LBhQ85jDz74YLRq1SrOPPPMiIgYM2ZM3HPPPXHZZZfFb3/72/jVr34VX/7yl+PDDz+scf8++uijiIi4/vrr4+mnn44HH3wwdttttzjyyCMrHZ0XEfGzn/0spk2bFuPHj49f//rXsXLlyjjuuONymrT77rsvvva1r8X+++8fkydPjmuvvTZuuOGGKrf3WUuXLo2DDz44fve738V3v/vdeOaZZ+L888+PcePGxVe/+tXNrg8AUJ8WLlwYERF77bVXjePOOuusaNWqVUyYMCFn+fr16+Phhx+OE088MXbeeedYuXJlHH300fHPf/4zp8/q3r17pWlHPutvf/tbDBo0KO6555549tln47vf/W788Y9/jMMOOyzWrl1bafwpp5wSe+21V0yaNCmuvvrqeOSRR+Ib3/hG9vENGzbEiSeeGI888khcddVV8fjjj8chhxwSxx57bK1em4cffjiGDRsW7dq1i1/+8pfx3//939GhQ4c45phjhG1QyBKABvDggw8mEZG88soryerVq5PddtstGTBgQLJhw4YkSZLkiCOOSHr37p2zzhVXXJFERPL+++8nSZIkl156aTJw4MAkSZJk6tSpSfPmzZOysrIkSZLk3HPPTZo3b56Ul5fXWMeUKVOSiEieffbZ7LJ169YlXbp0SU455ZTssj59+iQjRoyo836vW7cuWbt2bXLUUUclX/rSl7LLFy5cmEREst9++yXr1q3LLv/Tn/6UREQyceLEJEmSZP369Unnzp2TQw45JGe77777btKyZcukR48eOcsjIrn++uuz9y+88MJkhx12SN59992ccbfddlsSEclbb71V530EANicjb3gyy+/nKxduzZZsWJF8tvf/jbp3Llzcvjhhydr167d7DZOPvnkpGvXrsn69euzy6ZOnZpERPLUU08lSZIkc+bMSSIieeKJJ+pU74YNG5K1a9cm7777bhIRyZNPPpl97Prrr08iIrn11ltz1vn617+etG7dOtvfPv3000lEJPfcc0/OuHHjxlXq2Ta+PgsXLkySJElWrlyZdOjQITnxxBNz1l2/fn3St2/f5OCDD67T/gH1xxFtmzFz5sw48cQTo0uXLpHJZLZ47oBVq1bF6NGjY7/99osWLVrEiBEjKo2ZMWNGzhV4Nt7eeeeddHYCCkyrVq3iBz/4QcyZMyf++7//u9pxn52nbcaMGXHkkUdGRMRhhx0WERU/oxsfGzBgQLRt27bG5x4+fHh07tw5Hnzwweyy3/3ud/H+++/Heeedl1128MEHxzPPPBNXX311zJgxI/7973/Xev/uvffeOPDAA6N169bRokWLaNmyZfz+97+Pt99+u9LY448/Ppo3b569v//++0dEZK/MOn/+/Fi6dGmcdtppOet17949Dj300M3W8r//+78xdOjQ6NKlS6xbty57Gz58eERUHCUI0FQ1RJ+3qT/84Q/RokWLOOCAA7a6ZmjsBg4cGC1btoy2bdvGscceGzvuuGM8+eST0aJFi82ue+6558Z7770Xzz33XHbZgw8+GJ07d872NnvssUfsuOOOcdVVV8W9994bf/nLX2pd27Jly+Kiiy6Kbt26ZXu4Hj16RERU2cd98YtfzLm///77x6pVq2LZsmUR8Z8+67N93BlnnLHZWmbNmhUfffRRjBo1KqeH27BhQxx77LHxyiuvxMqVK2u9b0DDEbRtxsqVK6Nv375x1113bdX669evjzZt2sRll10WX/jCF2ocO3/+/CgtLc3e9txzz616TmgMTj/99DjwwAPjmmuuqfJQ/IiII444Ipo1axbTp0+PDz/8MN5888044ogjIqLiFNR+/frFjBkzYtGiRbFw4cIaTxvdqEWLFnH22WfH448/HsuXL4+Iisupl5SUxDHHHJMdd+edd8ZVV10VTzzxRAwdOjQ6dOgQI0aM2Ozl52+//fb4f//v/8UhhxwSkyZNipdffjleeeWVOPbYY6sM63baaaec+0VFRRER2bEbT1Xt1KlTpXWrWvZZ//znP+Opp56Kli1b5tx69+4dEVFpzhGApqQh+7yysrI455xz4qijjtqq54JtxUMPPRSvvPJKPP/883HhhRfG22+/XavgKaLiC9OSkpLsF6Yff/xxTJkyJc4555zsF5fFxcXxwgsvxAEHHBDf+c53onfv3tGlS5e4/vrrq+05IypO8xw2bFhMnjw5vv3tb8fvf//7+NOf/pSdO3dr+7gWLVpEhw4dcsbVtoeLiDj11FMr9XG33HJLJEmSnbIEKCyb/9qgiRs+fHj225GqrFmzJq699tr49a9/HcuXL48+ffrELbfckj3qZvvtt4977rknIiq+xdz4H/uqdOzY0dUSaTIymUzccsstcfTRR8d9991X5Zji4uJsmDZjxoxo1qxZzlFcRxxxREyfPj3222+/iKh5frZNnXvuufGjH/0oHn300Rg5cmRMmTIlrrjiipwjy7bffvu44YYb4oYbboh//vOf2aPbTjzxxBqPNn344YfjyCOPzP7cb7S5OUGqs7GBq2qC4KVLl252/Z133jn233//uOmmm6p8vEuXLltVF8C2oCH7vAsvvDC+8pWvRPPmzV1dkSZtn332yV4AYejQobF+/fr4r//6r/jNb34Tp556ao3rNm/ePM4+++y48847Y/ny5fHII4/E6tWr49xzz80Zt99++8Wjjz4aSZLEG2+8ERMmTIgbb7wx2rRpE1dffXWV237zzTfj9ddfjwkTJsSoUaOyy6u6+EJt7bTTTrFu3br46KOPcsK22vZwERE//elPY+DAgVWOqU1gBzQ8R7TV0bnnnht/+MMf4tFHH4033ngjvvzlL8exxx672aNeqtKvX78oKSmJo446KmcCeNhWfeELX4ijjz46brzxxvjkk0+qHDN06NBYsGBBPPLII9G/f/+cU0OPOOKImDdvXjzxxBPRsmXLWp1KGVHR4B1yyCHx4IMPVtugbapTp04xevToOOOMM2L+/Pnx6aefVjs2k8lkv83c6I033ojZs2fXqrbP6tWrV3Tu3LnSKbaLFi2KWbNmbXb9E044Id58883YfffdY8CAAZVugjaA6qXV5z344IPxt7/9La6//vp6qhQar1tvvTV23HHH+O53v1vpYlVVOffcc2PVqlUxceLEmDBhQgwaNCj23nvvKsdmMpno27dv/OQnP4n27dvHa6+9Vu12M5lMRESlPu7nP//5FuxNro1nYjz22GM5yx999NHNrnvooYdG+/bt4y9/+UuVPdyAAQOiVatWW10bUH8c0VYHf/vb32LixInx3nvvZf+zeuWVV8Zvf/vbePDBB+OHP/xhrbZTUlIS9913X/Tv3z9Wr14dv/rVr+Koo46KGTNmxOGHH16fuwB5d8stt0T//v1j2bJl2dMZNzV06NC47bbb4vHHH48rr7wy57EhQ4ZERMSTTz4ZgwcPju23377Wz3veeefFhRdeGO+//34MHjw4evXqlfP4IYccEieccELsv//+seOOO8bbb78dv/rVr2LQoEGx3XbbVbvdE044Ib7//e/H9ddfH0cccUTMnz8/brzxxujZs2esW7eu1vVt1KxZs7jhhhviwgsvjFNPPTXOO++8WL58edxwww1RUlISzZrV/H3JjTfeGNOmTYvBgwfHZZddFr169YpVq1bFP/7xj5g6dWrce++90bVr1y2uC2Bbl1aft2DBgrj66qvjxRdfrNUcVNDU7LjjjjF27Nj49re/HY888kicddZZNY7fe++9Y9CgQTFu3LhYvHhxpTMj/vd//zfuvvvuGDFiROy2226RJElMnjw5li9fHkcffXSN2919993j6quvjiRJokOHDvHUU0/FtGnTtnrfjj322Dj00EPjm9/8ZpSXl0f//v1j9uzZ8dBDD0VE1NjH7bDDDvHTn/40Ro0aFR999FGceuqp0bFjx/jXv/4Vr7/+evzrX/+qdAYFUBj8a18Hr732WiRJUulS1KtXr650vn5NevXqlfOf/EGDBsXixYvjtttuE7SxzevXr1+cccYZ8cgjj1T5+JAhQ6JFixaxbt267LeCG7Vv3z7233//mDdvXvY0nto6/fTT44orroj33nuvyiMMPv/5z8eUKVPiJz/5SXz66aexyy67xDnnnBPXXHNNjdu95ppr4tNPP437778/br311th3333j3nvvjccffzx7UYct9bWvfS0ymUzceuut8aUvfSl23XXXuPrqq+PJJ5+MRYsW1bhuSUlJzJkzJ77//e/Hj370o3jvvfeibdu20bNnz+wExABUlkaft379+vjKV74SN9xwQ6XtAP9x6aWXxl133RU33nhjnHHGGTnTeVTl3HPPja997WvRpk2bGDlyZM5je+65Z7Rv3z5uvfXWeP/996NVq1bRq1evSqeEflbLli3jqaeeissvvzwuvPDCaNGiRXzhC1+I5557Lrp3775V+9WsWbN46qmn4pvf/GbcfPPNsWbNmjj00EPj4YcfjoEDB2522qCzzjorunfvHrfeemtceOGFsWLFiujYsWMccMABMXr06K2qCah/mSRJknwX0VhkMpl4/PHHs1eUeuyxx+LMM8+Mt956q9I/BjvssEN07tw5Z9no0aNj+fLltZqX46abboqHH364yqvbACxfvjz22muvGDFiRLVz3AFQe/XR5y1fvjx23HHHnPU3bNgQSZJE8+bN49lnn43Pf/7z9bZPQGF65JFH4swzz4w//OEPMXjw4HyXA6TMEW110K9fv1i/fn0sW7YsewpbWubOnRslJSWpbhNonJYuXRo33XRTDB06NHbaaad499134yc/+UmsWLEiLr/88nyXB7BNSqPPa9euXfz5z3/OWXb33XfH888/H7/5zW+iZ8+eaZQKFLCJEyfGkiVLYr/99otmzZrFyy+/HD/60Y/i8MMPF7LBNkrQthmffPJJzpVmFi5cGPPmzYsOHTrEXnvtFWeeeWacc8458eMf/zj69esXH3zwQTz//POx3377xXHHHRcREX/5y19izZo18dFHH8WKFSti3rx5ERFxwAEHRETE+PHjY9ddd43evXvHmjVr4uGHH45JkybFpEmTGnp3gQJUVFQU//jHP+LrX/96fPTRR7HddtvFwIED4957761yXjsAaqe++7xmzZpFnz59cp6zY8eO0bp160rLgW1T27Zt49FHH40f/OAHsXLlyigpKYnRo0fHD37wg3yXBtQTp45uxowZM2Lo0KGVlo8aNSomTJgQa9eujR/84Afx0EMPxZIlS2KnnXaKQYMGxQ033BD77bdfRETsuuuu8e6771baxsaX/tZbb4377rsvlixZEm3atInevXvH2LFjsw0cAADpa4g+77O+973vxRNPPJEN5ACAbYugDQAAAABSUP31hAEAAACAWhO0AQAAAEAKXAyhChs2bIj3338/2rZtG5lMJt/lAACNRJIksWLFiujSpUs0a+b7zEKkzwMAtkZt+zxBWxXef//96NatW77LAAAaqcWLF0fXrl3zXQZV0OcBAHWxuT5P0FaFtm3bRkTFi9euXbs8VwMANBbl5eXRrVu3bC9B4dHnAQBbo7Z9nqCtChtPI2jXrp0GDADYYk5JLFz6PACgLjbX55k8BAAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIQYt8FwAA0CDWr4948cWI0tKIkpKIIUMimjfPd1UAANRRIbV5gjYAYNs3eXLE5ZdHvPfef5Z17Rpxxx0RJ5+cv7oAAKiTQmvznDoKAGzbJk+OOPXU3O4rImLJkorlkyfnpy4AAOqkENs8QRsAsO1av77iK84kqfzYxmVXXFExDgCARqNQ2zxBGwCw7XrxxcpfcW4qSSIWL64YBwBAo1GobZ6gDQDYdpWWpjsOAICCUKhtnqANANh2lZSkOw4AgIJQqG2eoA0A2HYNGVJx2alMpurHM5mIbt0qxgEA0GgUapsnaAMAtl3Nm1dc2z2iche28f748RXjAABoNAq1zRO0AQDbtpNPjvjNbyJ22SV3edeuFctPPjk/dQEAUCeF2Oa1aPinBABoYCefHHHSSRWXnSotrZisY8gQR7IBADRyhdbmCdoAgKahefOII4/MdxUAAKSskNo8p44CAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKchr0DZz5sw48cQTo0uXLpHJZOKJJ56ocfzo0aMjk8lUuvXu3Ts7ZsKECVWOWbVqVT3vDQAAAABNWV6DtpUrV0bfvn3jrrvuqtX4O+64I0pLS7O3xYsXR4cOHeLLX/5yzrh27drljCstLY3WrVvXxy4AAAAAQEREtMjnkw8fPjyGDx9e6/HFxcVRXFycvf/EE0/Exx9/HOeee27OuEwmE507d06tTgAAAADYnEY9R9v9998fX/jCF6JHjx45yz/55JPo0aNHdO3aNU444YSYO3dujdtZvXp1lJeX59wAAAAAYEs02qCttLQ0nnnmmbjgggtylu+9994xYcKEmDJlSkycODFat24dhx56aCxYsKDabY0bNy57tFxxcXF069atvssHAAAAYBvTaIO2CRMmRPv27WPEiBE5ywcOHBhnnXVW9O3bN4YMGRL//d//HXvttVf89Kc/rXZbY8eOjbKysuxt8eLF9Vw9AAAAANuavM7RtrWSJIkHHnggzj777GjVqlWNY5s1axYHHXRQjUe0FRUVRVFRUdplAgAAANCENMoj2l544YX461//Gueff/5mxyZJEvPmzYuSkpIGqAwAAACApiqvR7R98skn8de//jV7f+HChTFv3rzo0KFDdO/ePcaOHRtLliyJhx56KGe9+++/Pw455JDo06dPpW3ecMMNMXDgwNhzzz2jvLw87rzzzpg3b1787Gc/q/f9AQAAAKDpymvQNmfOnBg6dGj2/pgxYyIiYtSoUTFhwoQoLS2NRYsW5axTVlYWkyZNijvuuKPKbS5fvjy+9rWvxdKlS6O4uDj69esXM2fOjIMPPrj+dgQAAACAJi+TJEmS7yIKTXl5eRQXF0dZWVm0a9cu3+UAAI2EHqLweY8AgK1R2x6iUc7RBgAAAACFRtAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAQF6MGzcuDjrooGjbtm107NgxRowYEfPnz69xnRkzZkQmk6l0e+eddxqoagCA6gnaAADIixdeeCEuvvjiePnll2PatGmxbt26GDZsWKxcuXKz686fPz9KS0uztz333LMBKgYAqFmLfBcAAEDT9Nvf/jbn/oMPPhgdO3aMV199NQ4//PAa1+3YsWO0b99+s8+xevXqWL16dfZ+eXn5VtUKAFAbjmgDAKAglJWVRUREhw4dNju2X79+UVJSEkcddVRMnz692nHjxo2L4uLi7K1bt26p1QsA8FmCNgAA8i5JkhgzZkwcdthh0adPn2rHlZSUxH333ReTJk2KyZMnR69eveKoo46KmTNnVjl+7NixUVZWlr0tXry4vnYBAMCpowAA5N8ll1wSb7zxRrz00ks1juvVq1f06tUre3/QoEGxePHiuO2226o83bSoqCiKiopSrxcAoCqOaAMAIK8uvfTSmDJlSkyfPj26du26xesPHDgwFixYUA+VAQBsGUe0AQCQF0mSxKWXXhqPP/54zJgxI3r27LlV25k7d26UlJSkXB0AwJYTtAEAkBcXX3xxPPLII/Hkk09G27ZtY+nSpRERUVxcHG3atImIijnWlixZEg899FBERIwfPz523XXX6N27d6xZsyYefvjhmDRpUkyaNClv+wEAsJGgDQCAvLjnnnsiIuLII4/MWf7ggw/G6NGjIyKitLQ0Fi1alH1szZo1ceWVV8aSJUuiTZs20bt373j66afjuOOOa6iyAQCqlUmSJMl3EYWmvLw8iouLo6ysLNq1a5fvcgCARkIPUfi8RwDA1qhtD+FiCAAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAryGrTNnDkzTjzxxOjSpUtkMpl44oknahw/Y8aMyGQylW7vvPNOzrhJkybFvvvuG0VFRbHvvvvG448/Xo97AQAAAAB5DtpWrlwZffv2jbvuumuL1ps/f36UlpZmb3vuuWf2sdmzZ8fIkSPj7LPPjtdffz3OPvvsOO200+KPf/xj2uUDAAAAQFaLfD758OHDY/jw4Vu8XseOHaN9+/ZVPjZ+/Pg4+uijY+zYsRERMXbs2HjhhRdi/PjxMXHixLqUCwAAAADVapRztPXr1y9KSkriqKOOiunTp+c8Nnv27Bg2bFjOsmOOOSZmzZpV7fZWr14d5eXlOTcAAAAA2BKNKmgrKSmJ++67LyZNmhSTJ0+OXr16xVFHHRUzZ87Mjlm6dGl06tQpZ71OnTrF0qVLq93uuHHjori4OHvr1q1bve0DAAAAANumvJ46uqV69eoVvXr1yt4fNGhQLF68OG677bY4/PDDs8szmUzOekmSVFq2qbFjx8aYMWOy98vLy4VtAAAAAGyRRnVEW1UGDhwYCxYsyN7v3LlzpaPXli1bVukot00VFRVFu3btcm4AAAAAsCUafdA2d+7cKCkpyd4fNGhQTJs2LWfMs88+G4MHD27o0gAAAABoQvJ66ugnn3wSf/3rX7P3Fy5cGPPmzYsOHTpE9+7dY+zYsbFkyZJ46KGHIqLiiqK77rpr9O7dO9asWRMPP/xwTJo0KSZNmpTdxuWXXx6HH3543HLLLXHSSSfFk08+Gc8991y89NJLDb5/AAAAADQdeQ3a5syZE0OHDs3e3zhP2qhRo2LChAlRWloaixYtyj6+Zs2auPLKK2PJkiXRpk2b6N27dzz99NNx3HHHZccMHjw4Hn300bj22mvjuuuui9133z0ee+yxOOSQQxpuxwAAAABocjJJkiT5LqLQlJeXR3FxcZSVlZmvDQCoNT1E4fMeAQBbo7Y9RKOfow0AAAAACoGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAIC/GjRsXBx10ULRt2zY6duwYI0aMiPnz5292vRdeeCH69+8frVu3jt122y3uvffeBqgWAGDzBG0AAOTFCy+8EBdffHG8/PLLMW3atFi3bl0MGzYsVq5cWe06CxcujOOOOy6GDBkSc+fOje985ztx2WWXxaRJkxqwcgCAqrXIdwEAADRNv/3tb3PuP/jgg9GxY8d49dVX4/DDD69ynXvvvTe6d+8e48ePj4iIffbZJ+bMmRO33XZbnHLKKfVdMgBAjRzRBgBAQSgrK4uIiA4dOlQ7Zvbs2TFs2LCcZcccc0zMmTMn1q5dW2n86tWro7y8POcGAFBfBG0AAORdkiQxZsyYOOyww6JPnz7Vjlu6dGl06tQpZ1mnTp1i3bp18cEHH1QaP27cuCguLs7eunXrlnrtAAAbCdoAAMi7Sy65JN54442YOHHiZsdmMpmc+0mSVLk8ImLs2LFRVlaWvS1evDidggEAqmCONgAA8urSSy+NKVOmxMyZM6Nr1641ju3cuXMsXbo0Z9myZcuiRYsWsdNOO1UaX1RUFEVFRanWCwBQHUe0AQCQF0mSxCWXXBKTJ0+O559/Pnr27LnZdQYNGhTTpk3LWfbss8/GgAEDomXLlvVVKgBArQjaAADIi4svvjgefvjheOSRR6Jt27axdOnSWLp0afz73//Ojhk7dmycc8452fsXXXRRvPvuuzFmzJh4++2344EHHoj7778/rrzyynzsAgBADkEbAAB5cc8990RZWVkceeSRUVJSkr099thj2TGlpaWxaNGi7P2ePXvG1KlTY8aMGXHAAQfE97///bjzzjvjlFNOyccuAADkMEcbAAB5sfEiBjWZMGFCpWVHHHFEvPbaa/VQEQBA3TiiDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFeQ3aZs6cGSeeeGJ06dIlMplMPPHEEzWOnzx5chx99NHxuc99Ltq1axeDBg2K3/3udzljJkyYEJlMptJt1apV9bgnAAAAADR1eQ3aVq5cGX379o277rqrVuNnzpwZRx99dEydOjVeffXVGDp0aJx44okxd+7cnHHt2rWL0tLSnFvr1q3rYxcAAAAAICIiWuTzyYcPHx7Dhw+v9fjx48fn3P/hD38YTz75ZDz11FPRr1+/7PJMJhOdO3dOq0wAAAAA2KxGPUfbhg0bYsWKFdGhQ4ec5Z988kn06NEjunbtGieccEKlI94+a/Xq1VFeXp5zAwAAAIAt0aiDth//+MexcuXKOO2007LL9t5775gwYUJMmTIlJk6cGK1bt45DDz00FixYUO12xo0bF8XFxdlbt27dGqJ8AAAAALYhmSRJknwXEVFxuufjjz8eI0aMqNX4iRMnxgUXXBBPPvlkfOELX6h23IYNG+LAAw+Mww8/PO68884qx6xevTpWr16dvV9eXh7dunWLsrKyaNeu3RbtBwDQdJWXl0dxcbEeooB5jwCArVHbHiKvc7RtrcceeyzOP//8+J//+Z8aQ7aIiGbNmsVBBx1U4xFtRUVFUVRUlHaZAAAAADQhje7U0YkTJ8bo0aPjkUceieOPP36z45MkiXnz5kVJSUkDVAcAAABAU5XXI9o++eST+Otf/5q9v3Dhwpg3b1506NAhunfvHmPHjo0lS5bEQw89FBEVIds555wTd9xxRwwcODCWLl0aERFt2rSJ4uLiiIi44YYbYuDAgbHnnntGeXl53HnnnTFv3rz42c9+1vA7CAAAAECTkdcj2ubMmRP9+vWLfv36RUTEmDFjol+/fvHd7343IiJKS0tj0aJF2fE///nPY926dXHxxRdHSUlJ9nb55Zdnxyxfvjy+9rWvxT777BPDhg2LJUuWxMyZM+Pggw9u2J0DAAAAoEkpmIshFBKT5AIAW0MPUfi8RwDA1qhtD9Ho5mgDAAAAgEIkaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFLQojaDpkyZssUbPvroo6NNmzZbvB4AAIVFLwgAUDu1CtpGjBixRRvNZDKxYMGC2G233bamJgAACoheEACgdmp96ujSpUtjw4YNtbptt9129VkzAAANTC8IALB5tQraRo0atUWH/p911lnRrl27rS4KAIDCoRcEAKidTJIkSb6LKDTl5eVRXFwcZWVlmkQAoNb0EIXPewQAbI3a9hCuOgoAAAAAKajVxRA2tWrVqvjpT38a06dPj2XLlsWGDRtyHn/ttddSKw4AgMKiFwQAqN4WB23nnXdeTJs2LU499dQ4+OCDI5PJ1EddAAAUIL0gAED1tjhoe/rpp2Pq1Klx6KGH1vnJZ86cGT/60Y/i1VdfjdLS0nj88cc3e/n4F154IcaMGRNvvfVWdOnSJb797W/HRRddlDNm0qRJcd1118Xf/va32H333eOmm26KL33pS3WuFwCgqUuzFwQA2NZs8Rxtu+yyS7Rt2zaVJ1+5cmX07ds37rrrrlqNX7hwYRx33HExZMiQmDt3bnznO9+Jyy67LCZNmpQdM3v27Bg5cmScffbZ8frrr8fZZ58dp512Wvzxj39MpWYAgKYszV4QAGBbs8VXHX3mmWfizjvvjHvvvTd69OiRXiGZzGaPaLvqqqtiypQp8fbbb2eXXXTRRfH666/H7NmzIyJi5MiRUV5eHs8880x2zLHHHhs77rhjTJw4scrtrl69OlavXp29X15eHt26dXM1KgBgizSFK1rWVy/YUJrCewQApK/erjo6YMCAWLVqVey2227Rtm3b6NChQ86tPs2ePTuGDRuWs+yYY46JOXPmxNq1a2scM2vWrGq3O27cuCguLs7eunXrln7xAADbgDR7wZkzZ8aJJ54YXbp0iUwmE0888USN42fMmBGZTKbS7Z133qnDHgEApGeL52g744wzYsmSJfHDH/4wOnXq1KAT4C5dujQ6deqUs6xTp06xbt26+OCDD6KkpKTaMUuXLq12u2PHjo0xY8Zk7288og0AgFxp9oIbpxE599xz45RTTqn1evPnz8/5Jvlzn/vcVtcAAJCmLQ7aZs2aFbNnz46+ffvWRz2b9dlmbuOZr5sur2pMTU1gUVFRFBUVpVglAMC2Kc1ecPjw4TF8+PAtXq9jx47Rvn37Oj8/AEDatvjU0b333jv+/e9/10ctm9W5c+dKR6YtW7YsWrRoETvttFONYz57lBsAAFsun73gRv369YuSkpI46qijYvr06TWOXb16dZSXl+fcAADqyxYHbTfffHN885vfjBkzZsSHH37YoI3LoEGDYtq0aTnLnn322RgwYEC0bNmyxjGDBw+u19oAAJqCfPaCJSUlcd9998WkSZNi8uTJ0atXrzjqqKNi5syZ1a5jLl4AoCFt8VVHmzWryOaqOz1z/fr1td7WJ598En/9618jouKbydtvvz2GDh0aHTp0iO7du8fYsWNjyZIl8dBDD0VExMKFC6NPnz5x4YUXxle/+tWYPXt2XHTRRTFx4sTsvB6zZs2Kww8/PG666aY46aST4sknn4xrr702XnrppTjkkENqVZerUQEAW6Mp9BBp9oKbqs0V6Kty4oknRiaTiSlTplT5uKvLAwBpqG2ft8VztD3//POpXQBhzpw5MXTo0Oz9jRckGDVqVEyYMCFKS0tj0aJF2cd79uwZU6dOjW984xvxs5/9LLp06RJ33nlnzuS5gwcPjkcffTSuvfbauO6662L33XePxx57rNYhGwAA1UuzF0zDwIED4+GHH672cXPxAgANqdZHtN13333xxS9+MTp37lzfNeVdU/g2GgBI37bcQ9R3L7i1R7Sdeuqp8dFHH8Xzzz9fq/Hb8nsEANSf2vYQtZ6jbeLEibHrrrvGIYccEj/84Q/jrbfeSqVQAAAKX330gp988knMmzcv5s2bFxEV04TMmzcve0bD2LFj45xzzsmOHz9+fDzxxBOxYMGCeOutt2Ls2LExadKkuOSSS+pcCwBAGmodtE2fPj1KS0vj0ksvjXnz5sXgwYNj9913jzFjxsSMGTNiw4YN9VknAAB5VB+94Jw5c6Jfv37Rr1+/iKiYRqRfv37x3e9+NyKi0jQia9asiSuvvDL233//GDJkSLz00kvx9NNPx8knn5zOTgIA1NEWXwxhozVr1sTzzz8fU6ZMiaeeeio+/fTTOP744+OLX/xiDB8+PLbffvu0a20wTikAALZGU+ohGmsv2JTeIwAgPbXtIbY6aPusOXPmxJQpU+LJJ5+MU089Na677ro0NpsXGjAAYGs05R6isfSCTfk9AgC2XoMHbZtau3ZttGzZMu3NNhgNGACwNfQQFQq5F/QeAQBbo7Y9RIvabvDGG2/c7JhMJhPXXXddwTZWAABsHb0gAMDm1fqIto2T1Fa5kUwm5s+fH6tWrYr169enVly++KYTANga23IPsa30gtvyewQA1J/Uj2ibO3dulcvnzZsXV199dbz55pvx1a9+dcsrBQCg4OkFAQA2r9nWrrhw4cI466yz4qCDDori4uJ466234t57702zNgAACpReEACgsi0O2j744IO49NJLY++9947S0tKYNWtWPPbYY7HnnnvWR30AABQQvSAAQPVqferoypUr47bbbovbb7899thjj3jqqadi2LBh9VkbAAAFQi8IALB5tQ7adt9991ixYkVceumlccYZZ0Qmk4k33nij0rj9998/1QIBAMg/vSAAwObV+qqjzZr95yzTTCYTm6628X4mkyn4K03VhqtRAQBbY1vuIbaVXnBbfo8AgPqT+lVHFy5cmEphAAA0PnpBAIDNq3XQ1qNHj/qsAwCAAqYXBADYvFpddfSNN96IDRs21Hqjb731Vqxbt26riwIAoHDoBQEAaqdWQVu/fv3iww8/rPVGBw0aFIsWLdrqogAAKBx6QQCA2qnVqaNJksR1110X2223Xa02umbNmjoVBQBA4dALAgDUTq2CtsMPPzzmz59f640OGjQo2rRps9VFAQBQOPSCAAC1U6ugbcaMGfVcBgAAhUovCABQO7Waow0AAAAAqJmgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEjBVgVtv/rVr+LQQw+NLl26xLvvvhsREePHj48nn3wy1eIAACg8ekEAgKptcdB2zz33xJgxY+K4446L5cuXx/r16yMion379jF+/Pi06wMAoIDoBQEAqrfFQdtPf/rT+MUvfhHXXHNNNG/ePLt8wIAB8ec//znV4gAAKCx6QQCA6m1x0LZw4cLo169fpeVFRUWxcuXKVIoCAKAw6QUBAKq3xUFbz549Y968eZWWP/PMM7HvvvumURMAAAVKLwgAUL0WW7rCt771rbj44otj1apVkSRJ/OlPf4qJEyfGuHHj4r/+67/qo0YAAAqEXhAAoHpbHLSde+65sW7duvj2t78dn376aXzlK1+JXXbZJe644444/fTT66NGAAAKhF4QAKB6mSRJkq1d+YMPPogNGzZEx44d06wp78rLy6O4uDjKysqiXbt2+S4HAGgkmloP0Rh7wab2HgEA6ahtD7HFR7QtXLgw1q1bF3vuuWfsvPPO2eULFiyIli1bxq677rpVBQMAUPj0ggAA1dviiyGMHj06Zs2aVWn5H//4xxg9enQaNQEAUKD0ggAA1dvioG3u3Llx6KGHVlo+cODAKq9ABQDAtkMvCABQvS0O2jKZTKxYsaLS8rKysli/fn0qRQEAUJj0ggAA1dvioG3IkCExbty4nEZq/fr1MW7cuDjssMNSLQ4AgMKiFwQAqN4WXwzhlltuiSOOOCJ69eoVQ4YMiYiIF198McrLy+P5559PvUAAAAqHXhAAoHpbfERb796944033ojTTjstli1bFitWrIhzzjkn3nnnnejTp0991AgAQIHQCwIAVG+Ljmhbu3ZtDBs2LH7+85/HD3/4w/qqCQCAAqQXBACo2RYd0dayZct48803I5PJ1Fc9AAAUKL0gAEDNtvjU0XPOOSfuv//++qgFAIACpxcEAKjeFl8MYc2aNfFf//VfMW3atBgwYEBsv/32OY/ffvvtqRUHAEBh0QsCAFRvi4O2N998Mw488MCIiPi///u/nMecRgAAsG3TCwIAVG+Lg7bp06fXRx0AADQCekEAgOpt8Rxtm3rvvfdiyZIladUCAEAjohcEAMi1xUHbhg0b4sYbb4zi4uLo0aNHdO/ePdq3bx/f//73Y8OGDfVRIwAABUIvCABQvS0+dfSaa66J+++/P26++eY49NBDI0mS+MMf/hDf+973YtWqVXHTTTfVR50AABQAvSAAQPUySZIkW7JCly5d4t57740vfvGLOcuffPLJ+PrXv75NnD5QXl4excXFUVZWFu3atct3OQBAI9EUeojG3gs2hfcIAEhfbXuILT519KOPPoq999670vK99947Pvrooy3dHAAAjYheEACgelsctPXt2zfuuuuuSsvvuuuu6Nu3bypFAQBQmPSCAADV2+I52m699dY4/vjj47nnnotBgwZFJpOJWbNmxeLFi2Pq1Kn1USMAAAVCLwgAUL0tPqLtiCOOiP/7v/+LL33pS7F8+fL46KOP4uSTT4758+fHkCFD6qNGAAAKhF4QAKB6tb4Ywt///vfo2bNnZDKZ+q4p70ySCwBsjW25h9hWesFt+T0CAOpP6hdD2HPPPeNf//pX9v7IkSPjn//8Z92qBACgUdALAgBsXq2Dts8e+DZ16tRYuXJl6gUBAFB49IIAAJu3xXO0AQAAAACV1Tpoy2QylebkaOxzdAAAUDt6QQCAzWtR24FJksTo0aOjqKgoIiJWrVoVF110UWy//fY54yZPnpxuhQAA5J1eEABg82odtI0aNSrn/llnnZV6MQAAFCa9IADA5tU6aHvwwQfrsw4AAAqYXhAAYPNcDAEAAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFOQ9aLv77rujZ8+e0bp16+jfv3+8+OKL1Y4dPXp0ZDKZSrfevXtnx0yYMKHKMatWrWqI3QEAAACgicpr0PbYY4/FFVdcEddcc03MnTs3hgwZEsOHD49FixZVOf6OO+6I0tLS7G3x4sXRoUOH+PKXv5wzrl27djnjSktLo3Xr1g2xSwAAAAA0UXkN2m6//fY4//zz44ILLoh99tknxo8fH926dYt77rmnyvHFxcXRuXPn7G3OnDnx8ccfx7nnnpszLpPJ5Izr3LlzQ+wOAAAAAE1Y3oK2NWvWxKuvvhrDhg3LWT5s2LCYNWtWrbZx//33xxe+8IXo0aNHzvJPPvkkevToEV27do0TTjgh5s6dW+N2Vq9eHeXl5Tk3AAAAANgSeQvaPvjgg1i/fn106tQpZ3mnTp1i6dKlm12/tLQ0nnnmmbjgggtylu+9994xYcKEmDJlSkycODFat24dhx56aCxYsKDabY0bNy6Ki4uzt27dum3dTgEAAADQZOX9YgiZTCbnfpIklZZVZcKECdG+ffsYMWJEzvKBAwfGWWedFX379o0hQ4bEf//3f8dee+0VP/3pT6vd1tixY6OsrCx7W7x48VbtCwAAAABNV4t8PfHOO+8czZs3r3T02rJlyyod5fZZSZLEAw88EGeffXa0atWqxrHNmjWLgw46qMYj2oqKiqKoqKj2xQMAAADAZ+TtiLZWrVpF//79Y9q0aTnLp02bFoMHD65x3RdeeCH++te/xvnnn7/Z50mSJObNmxclJSV1qhcAAAAAapK3I9oiIsaMGRNnn312DBgwIAYNGhT33XdfLFq0KC666KKIqDilc8mSJfHQQw/lrHf//ffHIYccEn369Km0zRtuuCEGDhwYe+65Z5SXl8edd94Z8+bNi5/97GcNsk8AAAAANE15DdpGjhwZH374Ydx4441RWloaffr0ialTp2avIlpaWhqLFi3KWaesrCwmTZoUd9xxR5XbXL58eXzta1+LpUuXRnFxcfTr1y9mzpwZBx98cL3vDwAAAABNVyZJkiTfRRSa8vLyKC4ujrKysmjXrl2+ywEAGgk9ROHzHgEAW6O2PUTerzoKAAAAANsCQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAQF7MnDkzTjzxxOjSpUtkMpl44oknNrvOCy+8EP3794/WrVvHbrvtFvfee2/9FwoAUEuCNgAA8mLlypXRt2/fuOuuu2o1fuHChXHcccfFkCFDYu7cufGd73wnLrvsspg0aVI9VwoAUDst8l0AAABN0/Dhw2P48OG1Hn/vvfdG9+7dY/z48RERsc8++8ScOXPitttui1NOOaXKdVavXh2rV6/O3i8vL69TzQAANXFEGwAAjcLs2bNj2LBhOcuOOeaYmDNnTqxdu7bKdcaNGxfFxcXZW7du3RqiVACgiRK0AQDQKCxdujQ6deqUs6xTp06xbt26+OCDD6pcZ+zYsVFWVpa9LV68uCFKBQCaKKeOAgDQaGQymZz7SZJUuXyjoqKiKCoqqve6AAAiHNEGAEAj0blz51i6dGnOsmXLlkWLFi1ip512ylNVAAD/IWgDAKBRGDRoUEybNi1n2bPPPhsDBgyIli1b5qkqAID/ELQBAJAXn3zyScybNy/mzZsXERELFy6MefPmxaJFiyKiYn61c845Jzv+oosuinfffTfGjBkTb7/9djzwwANx//33x5VXXpmP8gEAKjFHGwAAeTFnzpwYOnRo9v6YMWMiImLUqFExYcKEKC0tzYZuERE9e/aMqVOnxje+8Y342c9+Fl26dIk777wzTjnllAavHQCgKplk4wyyZJWXl0dxcXGUlZVFu3bt8l0OANBI6CEKn/cIANgate0hnDoKAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbQAAAACQAkEbAAAAAKRA0AYAAAAAKRC0AQAAAEAKBG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJCCvAdtd999d/Ts2TNat24d/fv3jxdffLHasTNmzIhMJlPp9s477+SMmzRpUuy7775RVFQU++67bzz++OP1vRsAAAAANHF5Ddoee+yxuOKKK+Kaa66JuXPnxpAhQ2L48OGxaNGiGtebP39+lJaWZm977rln9rHZs2fHyJEj4+yzz47XX389zj777DjttNPij3/8Y33vDgAAAABNWCZJkiRfT37IIYfEgQceGPfcc0922T777BMjRoyIcePGVRo/Y8aMGDp0aHz88cfRvn37Krc5cuTIKC8vj2eeeSa77Nhjj40dd9wxJk6cWKu6ysvLo7i4OMrKyqJdu3ZbtlMAQJOlhyh83iMAYGvUtofI2xFta9asiVdffTWGDRuWs3zYsGExa9asGtft169flJSUxFFHHRXTp0/PeWz27NmVtnnMMcfUuM3Vq1dHeXl5zg0AAAAAtkTegrYPPvgg1q9fH506dcpZ3qlTp1i6dGmV65SUlMR9990XkyZNismTJ0evXr3iqKOOipkzZ2bHLF26dIu2GRExbty4KC4uzt66detWhz0DAAAAoClqke8CMplMzv0kSSot26hXr17Rq1ev7P1BgwbF4sWL47bbbovDDz98q7YZETF27NgYM2ZM9n55ebmwDQAAAIAtkrcj2nbeeedo3rx5pSPNli1bVumItJoMHDgwFixYkL3fuXPnLd5mUVFRtGvXLucGAAAAAFsib0Fbq1aton///jFt2rSc5dOmTYvBgwfXejtz586NkpKS7P1BgwZV2uazzz67RdsEAAAAgC2V11NHx4wZE2effXYMGDAgBg0aFPfdd18sWrQoLrroooioOKVzyZIl8dBDD0VExPjx42PXXXeN3r17x5o1a+Lhhx+OSZMmxaRJk7LbvPzyy+Pwww+PW265JU466aR48skn47nnnouXXnopL/sIAAAAQNOQ16Bt5MiR8eGHH8aNN94YpaWl0adPn5g6dWr06NEjIiJKS0tj0aJF2fFr1qyJK6+8MpYsWRJt2rSJ3r17x9NPPx3HHXdcdszgwYPj0UcfjWuvvTauu+662H333eOxxx6LQw45pMH3DwAAAICmI5MkSZLvIgpNeXl5FBcXR1lZmfnaAIBa00MUPu8RALA1attD5G2ONgAAAADYlgjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AgLy5++67o2fPntG6devo379/vPjii9WOnTFjRmQymUq3d955pwErBgConqANAIC8eOyxx+KKK66Ia665JubOnRtDhgyJ4cOHx6JFi2pcb/78+VFaWpq97bnnng1UMQBAzQRtAADkxe233x7nn39+XHDBBbHPPvvE+PHjo1u3bnHPPffUuF7Hjh2jc+fO2Vvz5s0bqGIAgJoJ2gAAaHBr1qyJV199NYYNG5azfNiwYTFr1qwa1+3Xr1+UlJTEUUcdFdOnT69x7OrVq6O8vDznBgBQXwRtAAA0uA8++CDWr18fnTp1ylneqVOnWLp0aZXrlJSUxH333ReTJk2KyZMnR69eveKoo46KmTNnVvs848aNi+Li4uytW7duqe4HAMCmWuS7AAAAmq5MJpNzP0mSSss26tWrV/Tq1St7f9CgQbF48eK47bbb4vDDD69ynbFjx8aYMWOy98vLy4VtAEC9cUQbAAANbuedd47mzZtXOnpt2bJllY5yq8nAgQNjwYIF1T5eVFQU7dq1y7kBANQXQRsAAA2uVatW0b9//5g2bVrO8mnTpsXgwYNrvZ25c+dGSUlJ2uUBAGwVp44CAJAXY8aMibPPPjsGDBgQgwYNivvuuy8WLVoUF110UURUnPa5ZMmSeOihhyIiYvz48bHrrrtG7969Y82aNfHwww/HpEmTYtKkSfncDQCALEEbAAB5MXLkyPjwww/jxhtvjNLS0ujTp09MnTo1evToERERpaWlsWjRouz4NWvWxJVXXhlLliyJNm3aRO/evePpp5+O4447Ll+7AACQI5MkSZLvIgpNeXl5FBcXR1lZmXk8AIBa00MUPu8RALA1attDmKMNAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUtAi3wUAAGzW+vURL74YUVoaUVISMWRIRPPm+a4KAIA62tbaPEEbAFDYJk+OuPzyiPfe+8+yrl0j7rgj4uST81cXAAB1si22eU4dBQAK1+TJEaeemtt9RUQsWVKxfPLk/NQFAECdbKttnqANAChM69dXfMWZJJUf27jsiisqxgEA0Ghsy22eoA0AKEwvvlj5K85NJUnE4sUV4wAAaDS25TYv70Hb3XffHT179ozWrVtH//7948UaXsXJkyfH0UcfHZ/73OeiXbt2MWjQoPjd736XM2bChAmRyWQq3VatWlXfuwIApKm0NN1xAAAUhG25zctr0PbYY4/FFVdcEddcc03MnTs3hgwZEsOHD49FixZVOX7mzJlx9NFHx9SpU+PVV1+NoUOHxoknnhhz587NGdeuXbsoLS3NubVu3bohdgkASEtJSbrjAAAoCNtym5dJkqrOiG0YhxxySBx44IFxzz33ZJfts88+MWLEiBg3blytttG7d+8YOXJkfPe7342IiiParrjiili+fPlW11VeXh7FxcVRVlYW7dq12+rtAAB1sH59xK67VsyIW1W7kslUXJZq4cKCuQa8HqLweY8AIP8aYZtX6x4ib0e0rVmzJl599dUYNmxYzvJhw4bFrFmzarWNDRs2xIoVK6JDhw45yz/55JPo0aNHdO3aNU444YRKR7x91urVq6O8vDznBgDkWfPmFdd2j6jotja18f748YXTfQEAUCvbcpuXt6Dtgw8+iPXr10enTp1ylnfq1CmWLl1aq238+Mc/jpUrV8Zpp52WXbb33nvHhAkTYsqUKTFx4sRo3bp1HHroobFgwYJqtzNu3LgoLi7O3rp167Z1OwUApOvkkyN+85uIXXbJXd61a8Xyk0/OT10AANTJttrm5e3U0ffffz922WWXmDVrVgwaNCi7/Kabbopf/epX8c4779S4/sSJE+OCCy6IJ598Mr7whS9UO27Dhg1x4IEHxuGHHx533nlnlWNWr14dq1evzt4vLy+Pbt26OaUAAArF+vUVl50qLa2YrGPIkIL8itNpiYXPewQAhaWRtHm17iFaNGBNOXbeeedo3rx5paPXli1bVukot8967LHH4vzzz4//+Z//qTFki4ho1qxZHHTQQTUe0VZUVBRFRUW1Lx4AaFjNm0cceWS+qwAAIGXbWpuXt1NHW7VqFf37949p06blLJ82bVoMHjy42vUmTpwYo0ePjkceeSSOP/74zT5PkiQxb968KGmMl6oAAAAAoNHI2xFtERFjxoyJs88+OwYMGBCDBg2K++67LxYtWhQXXXRRRESMHTs2lixZEg899FBEVIRs55xzTtxxxx0xcODA7NFwbdq0ieLi4oiIuOGGG2LgwIGx5557Rnl5edx5550xb968+NnPfpafnQQAAACgSchr0DZy5Mj48MMP48Ybb4zS0tLo06dPTJ06NXr06BEREaWlpbFo0aLs+J///Oexbt26uPjii+Piiy/OLh81alRMmDAhIiKWL18eX/va12Lp0qVRXFwc/fr1i5kzZ8bBBx/coPsGAAAAQNOSt4shFDKT5AIAW0MPUfi8RwDA1qhtD5G3OdoAAAAAYFsiaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEiBoA0AAAAAUiBoAwAAAIAUtMh3AU3O+vURL74YUVoaUVISMWRIRPPm+a4KAIA60uYBAIK2hjR5csTll0e8995/lnXtGnHHHREnn5y/ugAAqBNtHgAQ4dTRhjN5csSpp+Z2XxERS5ZULJ88OT91AQBQJ9o8AGAjQVtDWL++4ivOJKn82MZlV1xRMQ4AgEZDmwcAbErQ1hBefLHyV5ybSpKIxYsrxgEA0Gho8wCATQnaGkJpabrjAAAoCNo8AGBTgraGUFKS7jgAAAqCNg8A2JSgrSEMGVJx2alMpurHM5mIbt0qxgEA0Gho8wCATQnaGkLz5hXXdo+o3IVtvD9+fMU4AAAaDW0eALApQVtDOfnkiN/8JmKXXXKXd+1asfzkk/NTFwAAdaLNAwA2apHvApqUk0+OOOmkistOlZZWTNYxZIivOAEAGjltHgAQIWhreM2bRxx5ZL6rAAAgZdo8AEDQBgDUv/XrHeoDALAN0ublErQBAPVr8uSIyy+PeO+9/yzr2rViBnmTVwEANFravMpcDAEAqD+TJ0ecempu9xURsWRJxfLJk/NTFwAAdaLNq5qgDQCoH+vXV3zFmSSVH9u47IorKsYBANBoaPOqJ2gDAOrHiy9W/opzU0kSsXhxxTgAABoNbV71BG0AQP0oLU13HAAABUGbVz1BGwBQP0pK0h0HAEBB0OZVT9AGANSPIUMqLjuVyVT9eCYT0a1bxTgAABoNbV71BG0AQP1o3rzi2u4RlbuwjffHj68YBwBAo6HNq56grTFavz5ixoyIiRMr/myKl/EAoHE4+eSI3/wmYpddcpd37Vqx/OST81MXFChtHgCNhTavai3yXQBbaPLkimvobnp5j65dK6LkpvopBqCwnXxyxEknVVx2qrS0YrKOIUOa5lecUANtHgCNjTavskySJEm+iyg05eXlUVxcHGVlZdGuXbt8l/MfkydHnHpqxXVyN7XxuMymHBkDQAEo2B6CrEJ9j7R5AFDYattDOHW0sVi/vuIrzqpy0Y3LrrjC+QUAAI2MNg8Ath2CtsbixRdzzyP4rCSJWLy4YhwAAI2GNg8Ath3maGssSkvTHQcAtbV+vYk3oB5p8wDIF21e+gRtjUVJSbrjAKA2zM4O9U6bB0A+aPPqh1NHG4shQyo+8RtnxP2sTCaiW7eKcQCQho2zs3/2nLYlSyqWT56cn7pgG6PNA6ChafPqj6CtsWjevCJWjqjchW28P368YzwBSIfZ2aHBaPMAaEjavPolaGtMTj654truu+ySu7xrV9d8ByBdZmeHBqXNA6ChaPPqlznaGpuTT4446aS6zVZotkMANsfs7NDg0mjzIrR6ANRMm1e/BG2NUfPmEUceuXXrmu0QgNowOzvkRV3avAitHgCbp82rX04dbUrMdghAbZmdHRodrR4AtaHNq1+CtqbCbIcATdP69REzZkRMnFjxZ21/z5udHRoVrR5A07Q1rZ42r34J2poKsx0CND2TJ0fsumvE0KERX/lKxZ+77lr7w1rMzg6NhlYPoOmpS6unzas/5mhrKtKe7dAsuwCFbeM5ZJ89vGXjOWS17aDSmp0dqFdptnraPIDCl0arp82rH4K2piLN2Q7NsgtQ2DZ3DlkmU3EO2Ukn1a6Tquvs7EC9S6vV0+YBFL40Wz1tXvqcOtpUpDXboVl2ARrW1ky84RwyaHLSaPW0eQANa2un0tXqFTZBW1ORxmyHZtkFaFhbO/FG2tMFAAWvrq2eNg+gYdVlfjWtXmETtDUldZ3tMO3YfGvje4CmoC6HlqQ5XQDQaNSl1dPmATScuh5BrNUrbOZoa2rqMtthmrG5CUCApmBrZxSv68QbG88hW7Kk6m1kMhWPb266AKDR2dpWT5sHsGXy1eZFaPUKnaCtKdra2Q7TnGU3jSvhuSQWUMjq8j/NLTm0pKrf5xvPITv11IpOa9Pft7WdLgBotLam1dPmAdRePtu8CK1eoXPqKLWXxiy7aU0AUpcT2gFqa2vPfarr+QBpHFpS1+kCgCZFmwc0NY25zYvQ6hUyQRu1l8YFFdKYACTNS2KZQAS2XXX9+d7a/+ml8T/NtA4tOfnkiH/8I2L69IhHHqn4c+FCnRdQiTYPaEy0eRW0egUqoZKysrIkIpKysrJ8l1KYJk1Kkq5dk6Tid0nFrVu3iuWb88gjuetVd3vkkarXX7eu8nNvestkKmpZt27r9qNr19rtx6b1TJ9eUe/06bV7XqB26vLzVdef70mTKn6fVPU7JpOpeTvTp9fu99z06TXve9euVdewpb/raFB6iMLnPaqZNi+3Hm0e1A9tnjavMaptD+GINrZcXWLzusb3aV0SK42vS9M6r6GuX8f4upZClK+vGTeuW5ef77p+VZnG+QBpHFoCsBW0ef/ZhjYPqqbN2zxtXhPXQMFfo+KbznpU1/i+rl+VblpDXb4urcvXIJ/dTl2/jqnr17VJks5Xtr72TUchvBd1XT+fXzOm8fNd168q0/iqc9PXYmsPLSEv9BCFz3tUf7R5lbejzWNThfJeNNajybR55Ftte4hooHoaFQ1YPdv42/Wzv2Eb6ljdum4jrfMa6trFFUoXmNY28t01FEINhfBe5LN7SpK6/3yl8TsirXOf0jofwP9uGhU9ROHzHtUvbV7u66DNq5DvFqsQaiiU96Iu29DmafOaOkFbHWjAGsDWxvdp/Gar62/XNE/M39p/ZQqlC0xzG/nuPPJdQyG8F/nunpKk7j9faRwOkcbPeF3+p0mjpocofN6j+qfN0+Z9dhvavMJ5Lxrz0WTaPPJN0FYHGrAGsrXxfV1/s20Lv+ELoQtMaxv57hoKoYZCeC8KoXtKksL6H1Jdv6p0PkCTpIcofN6jhqHN0+YlSf5brEKooVDei23haDJtHvkmaKsDDVgjUJffbHX97VoI/8oUQheYxjYKoWsohBoK4b0ohM91GnWk2T2l8VWl8wGaHD1E4fMeFT5tXv7/OU1jG4XQYhVCDYXwXqSxjUL4XGrzyDdXHWXbVpdLYtX1Mi9DhkR07Vp53U230a1bxbjq1PWyXHVdPyKdS+bUdRtpXF6srtsohBoK4b1Io4Y0Ppd1/flK6zJOJ58c8ZvfROyyS+7yrl0rltfmd83Geo48MuKMMyr+dPkogM3S5tVt/YjCaC0KocUqhBoK4b1IYxvavKrr0eZRFUEbjVddfrPV5bdrGr/h6/qvTCF0gWlsoxC6hkKooRDei0LoniLS+flKq3uqy//0AKgTbZ42L41tFEINhfBepLENbR5sgQY6wq5RcUpBE5L2pYO25MT8uh6zXNf10zj2uhDOzyiEUyYL4Tj4um6j0I7FT2PiC8fzkwd6iMLnPWo6tHnavEKooRDei7S2oc2jqTNHWx1owKi1uv6GT6OLy2cXWNdtFELXUAg11PV1TGsbhdQ9JYkOikZJD1H4vEfUljZPm5fWNvL9XqS9DW0eTZWgrQ40YDSouv4rk+8usK7bKISuoRBq2LiNfL4XadWQJLonmiw9ROHzHtGQtHn5b7EKoYaN28h3m5fWNrR5NFW17SEySZIk+TlptXCVl5dHcXFxlJWVRbt27fJdDtS/9esrZnAtLa2YmGHIkC2fzbMu25g8OeLyy3Nnmu3WrWKShtpOlFDXbRRCDRH5fy/SqgGaKD1E4fMe0dTku7UohBarEGqIyP97keY2oCmqbQ8haKuCBgzyoBC6hkKoAWjU9BCFz3sEDa8QWqxCqAFo3GrbQ+T9qqN333139OzZM1q3bh39+/ePF2u6tnJEvPDCC9G/f/9o3bp17LbbbnHvvfdWGjNp0qTYd999o6ioKPbdd994/PHH66t8IC1pXB+7rtsohBoAmpj66AWBwlIILVYh1AA0DXkN2h577LG44oor4pprrom5c+fGkCFDYvjw4bFo0aIqxy9cuDCOO+64GDJkSMydOze+853vxGWXXRaTJk3Kjpk9e3aMHDkyzj777Hj99dfj7LPPjtNOOy3++Mc/NtRuAQBQC/XRCwIA5FNeTx095JBD4sADD4x77rknu2yfffaJESNGxLhx4yqNv+qqq2LKlCnx9ttvZ5dddNFF8frrr8fs2bMjImLkyJFRXl4ezzzzTHbMscceGzvuuGNMnDixVnU5pQAA2Bp6iC1TH73g5niPAICtUfCnjq5ZsyZeffXVGDZsWM7yYcOGxaxZs6pcZ/bs2ZXGH3PMMTFnzpxYu3ZtjWOq22ZExOrVq6O8vDznBgBA/amvXvCz9HkAQEPKW9D2wQcfxPr166NTp045yzt16hRLly6tcp2lS5dWOX7dunXxwQcf1Dimum1GRIwbNy6Ki4uzt27dum3NLgEAUEv11Qt+lj4PAGhIeb8YQiaTybmfJEmlZZsb/9nlW7rNsWPHRllZWfa2ePHiWtcPAMDWq49ecFP6PACgIbXI1xPvvPPO0bx580rfWC5btqzSN5Ubde7cucrxLVq0iJ122qnGMdVtMyKiqKgoioqKtmY3AADYCvXVC36WPg8AaEh5O6KtVatW0b9//5g2bVrO8mnTpsXgwYOrXGfQoEGVxj/77LMxYMCAaNmyZY1jqtsmAAANr756QQCAfMrrqaNjxoyJ//qv/4oHHngg3n777fjGN74RixYtiosuuigiKg71P+ecc7LjL7roonj33XdjzJgx8fbbb8cDDzwQ999/f1x55ZXZMZdffnk8++yzccstt8Q777wTt9xySzz33HNxxRVXNPTuAQBQg/roBQEA8ilvp45GRIwcOTI+/PDDuPHGG6O0tDT69OkTU6dOjR49ekRERGlpaSxatCg7vmfPnjF16tT4xje+ET/72c+iS5cuceedd8Ypp5ySHTN48OB49NFH49prr43rrrsudt9993jsscfikEMOafD9AwCgevXRCwIA5FMm2TiDLFnl5eVRXFwcZWVl0a5du3yXAwA0EnqIwuc9AgC2Rm17iLxfdRQAAAAAtgWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2AAAAAEhBi3wXUIiSJImIiPLy8jxXAgA0Jht7h429BIVHnwcAbI3a9nmCtiqsWLEiIiK6deuW50oAgMZoxYoVUVxcnO8yqII+DwCoi831eZnEV66VbNiwId5///1o27ZtZDKZfJeTV+Xl5dGtW7dYvHhxtGvXLt/lFBSvTfW8NtXz2lTPa1M9r031Cu21SZIkVqxYEV26dIlmzczQUYj0ef9RaD8/hcRrUz2vTfW8NtXz2lTPa1O9QnttatvnOaKtCs2aNYuuXbvmu4yC0q5du4L4YBcir031vDbV89pUz2tTPa9N9QrptXEkW2HT51VWSD8/hcZrUz2vTfW8NtXz2lTPa1O9QnptatPn+aoVAAAAAFIgaAMAAACAFAjaqFFRUVFcf/31UVRUlO9SCo7Xpnpem+p5barntame16Z6XhvYen5+que1qZ7Xpnpem+p5barntaleY31tXAwBAAAAAFLgiDYAAAAASIGgDQAAAABSIGgDAAAAgBQI2gAAAAAgBYI2Yty4cXHQQQdF27Zto2PHjjFixIiYP39+zpjRo0dHJpPJuQ0cODBPFTec733ve5X2u3PnztnHkySJ733ve9GlS5do06ZNHHnkkfHWW2/lseKGs+uuu1Z6bTKZTFx88cUR0bQ+MzNnzowTTzwxunTpEplMJp544omcx2vzOVm9enVceumlsfPOO8f2228fX/ziF+O9995rwL2oHzW9NmvXro2rrroq9ttvv9h+++2jS5cucc4558T777+fs40jjzyy0mfp9NNPb+A9Sd/mPje1+RnaVj83EZt/far6/ZPJZOJHP/pRdsy2+tmBLaHPq54+r3r6vP/Q59VMr1c9vV71tvU+T9BGvPDCC3HxxRfHyy+/HNOmTYt169bFsGHDYuXKlTnjjj322CgtLc3epk6dmqeKG1bv3r1z9vvPf/5z9rFbb701br/99rjrrrvilVdeic6dO8fRRx8dK1asyGPFDeOVV17JeV2mTZsWERFf/vKXs2Oaymdm5cqV0bdv37jrrruqfLw2n5MrrrgiHn/88Xj00UfjpZdeik8++SROOOGEWL9+fUPtRr2o6bX59NNP47XXXovrrrsuXnvttZg8eXL83//9X3zxi1+sNParX/1qzmfp5z//eUOUX68297mJ2PzP0Lb6uYnY/Ouz6etSWloaDzzwQGQymTjllFNyxm2Lnx3YEvq8munzqqbP+w99Xs30etXT61Vvm+/zEviMZcuWJRGRvPDCC9llo0aNSk466aT8FZUn119/fdK3b98qH9uwYUPSuXPn5Oabb84uW7VqVVJcXJzce++9DVRh4bj88suT3XffPdmwYUOSJE33MxMRyeOPP569X5vPyfLly5OWLVsmjz76aHbMkiVLkmbNmiW//e1vG6z2+vbZ16Yqf/rTn5KISN59993ssiOOOCK5/PLL67e4PKvqtdncz1BT+dwkSe0+OyeddFLy+c9/PmdZU/jswJbS5/2HPq/29HkV9Hk10+tVT69XvW2xz3NEG5WUlZVFRESHDh1yls+YMSM6duwYe+21V3z1q1+NZcuW5aO8BrdgwYLo0qVL9OzZM04//fT4+9//HhERCxcujKVLl8awYcOyY4uKiuKII46IWbNm5avcvFizZk08/PDDcd5550Umk8kub6qfmU3V5nPy6quvxtq1a3PGdOnSJfr06dPkPktlZWWRyWSiffv2Oct//etfx8477xy9e/eOK6+8skkcTRBR88+Qz81//POf/4ynn346zj///EqPNdXPDlRHn5dLn7d5+rzq6fO2nF4vl15v8xpjn9ci3wVQWJIkiTFjxsRhhx0Wffr0yS4fPnx4fPnLX44ePXrEwoUL47rrrovPf/7z8eqrr0ZRUVEeK65fhxxySDz00EOx1157xT//+c/4wQ9+EIMHD4633norli5dGhERnTp1ylmnU6dO8e677+aj3Lx54oknYvny5TF69Ojssqb6mfms2nxOli5dGq1atYodd9yx0piN6zcFq1atiquvvjq+8pWvRLt27bLLzzzzzOjZs2d07tw53nzzzRg7dmy8/vrr2dNYtlWb+xnyufmPX/7yl9G2bds4+eSTc5Y31c8OVEefl0ufVzv6vOrp87aMXi+XXq92GmOfJ2gjxyWXXBJvvPFGvPTSSznLR44cmf17nz59YsCAAdGjR494+umnK33gtyXDhw/P/n2//faLQYMGxe677x6//OUvsxNVbvrNXkRFE/vZZdu6+++/P4YPHx5dunTJLmuqn5nqbM3npCl9ltauXRunn356bNiwIe6+++6cx7761a9m/96nT5/Yc889Y8CAAfHaa6/FgQce2NClNpit/RlqSp+bjR544IE488wzo3Xr1jnLm+pnB6qjz8ulz6sdfd7m6fM2T69XmV6vdhpjn+fUUbIuvfTSmDJlSkyfPj26du1a49iSkpLo0aNHLFiwoIGqKwzbb7997LfffrFgwYLsVak++23CsmXLKn2rtS17991347nnnosLLrigxnFN9TNTm89J586dY82aNfHxxx9XO2Zbtnbt2jjttNNi4cKFMW3atJxvOKty4IEHRsuWLZvcZ+mzP0NN/XOz0Ysvvhjz58/f7O+giKb72YEIfV5t6PMq0+fVTJ9XO3q92tHrVdZY+zxBG5EkSVxyySUxefLkeP7556Nnz56bXefDDz+MxYsXR0lJSQNUWDhWr14db7/9dpSUlGQPU9300NQ1a9bECy+8EIMHD85jlQ3rwQcfjI4dO8bxxx9f47im+pmpzeekf//+0bJly5wxpaWl8eabb27zn6WNjdeCBQviueeei5122mmz67z11luxdu3aJvdZ+uzPUFP+3Gzq/vvvj/79+0ffvn03O7apfnZo2vR5tafPq0yfVzN93ubp9WpPr1dZo+3z8nMNBgrJ//t//y8pLi5OZsyYkZSWlmZvn376aZIkSbJixYrkm9/8ZjJr1qxk4cKFyfTp05NBgwYlu+yyS1JeXp7n6uvXN7/5zWTGjBnJ3//+9+Tll19OTjjhhKRt27bJP/7xjyRJkuTmm29OiouLk8mTJyd//vOfkzPOOCMpKSnZ5l+XjdavX5907949ueqqq3KWN7XPzIoVK5K5c+cmc+fOTSIiuf3225O5c+dmr6ZUm8/JRRddlHTt2jV57rnnktdeey35/Oc/n/Tt2zdZt25dvnYrFTW9NmvXrk2++MUvJl27dk3mzZuX8/tn9erVSZIkyV//+tfkhhtuSF555ZVk4cKFydNPP53svffeSb9+/bbp16a2P0Pb6ucmSTb/c5UkSVJWVpZst912yT333FNp/W35swNbQp9XPX1ezfR5FfR5NdPrVU+vV71tvc8TtJFERJW3Bx98MEmSJPn000+TYcOGJZ/73OeSli1bJt27d09GjRqVLFq0KL+FN4CRI0cmJSUlScuWLZMuXbokJ598cvLWW29lH9+wYUNy/fXXJ507d06KioqSww8/PPnzn/+cx4ob1u9+97skIpL58+fnLG9qn5np06dX+TM0atSoJElq9zn597//nVxyySVJhw4dkjZt2iQnnHDCNvF61fTaLFy4sNrfP9OnT0+SJEkWLVqUHH744UmHDh2SVq1aJbvvvnty2WWXJR9++GF+dywFNb02tf0Z2lY/N0my+Z+rJEmSn//850mbNm2S5cuXV1p/W/7swJbQ51VPn1czfV4FfV7N9HrV0+tVb1vv8zJJkiRbezQcAAAAAFDBHG0AAAAAkAJBGwAAAACkQNAGAAAAACkQtAEAAABACgRtAAAAAJACQRsAAAAApEDQBgAAAAApELQBAAAAQAoEbcBWmTFjRmQymVi+fPkWrTdhwoRo3759vdS0qdGjR8eIESPq/Xkayq677hrjx4/f6vX/8Y9/RCaTiUwmEwcccECNY+vjtdv4eclkMtvU+wIA2yJ9XsPS58G2RdAGTdCyZcviwgsvjO7du0dRUVF07tw5jjnmmJg9e3a+S9tiGxuLefPm5Sy/4447YsKECXmpqZA999xz8fvf/77Bn3fw4MFRWloap512WoM/NwA0Jfq8pkufB4WhRb4LABreKaecEmvXro1f/vKXsdtuu8U///nP+P3vfx8fffRRvktLTXFxcb5LKEg77bRT7LTTTg3+vK1atYrOnTtHmzZtYvXq1Q3+/ADQVOjzmi59HhQGR7RBE7N8+fJ46aWX4pZbbomhQ4dGjx494uCDD46xY8fG8ccfHxFVf3u4fPnyyGQyMWPGjJzt/eEPf4i+fftG69at45BDDok///nPOY9PmDAhunfvHtttt1186Utfig8//DDn8aoOX7/iiiviyCOPzN7fsGFD3HLLLbHHHntEUVFRdO/ePW666aaIiOjZs2dERPTr1y8ymUx2vc9ud/Xq1XHZZZdFx44do3Xr1nHYYYfFK6+8kn184yHvv//972PAgAGx3XbbxeDBg2P+/Pk1vp5XXXVV7LXXXrHddtvFbrvtFtddd12sXbs2+/j3vve9OOCAA+JXv/pV7LrrrlFcXBynn356rFixIjtmxYoVceaZZ8b2228fJSUl8ZOf/CSOPPLIuOKKK6p93rKysvja174WHTt2jHbt2sXnP//5eP3112ustSrr16+PMWPGRPv27WOnnXaKb3/725EkSc6YJEni1ltvjd122y3atGkTffv2jd/85jc5Y6ZMmRJ77rlntGnTJoYOHRq//OUvt+qUEwBg6+nz9Hmb0udBfgjaoInZYYcdYocddognnngilW+cvvWtb8Vtt90Wr7zySnTs2DG++MUvZhuQP/7xj3HeeefF17/+9Zg3b14MHTo0fvCDH2zxc4wdOzZuueWWuO666+Ivf/lLPPLII9GpU6eIiPjTn/4UERWHypeWlsbkyZOr3Ma3v/3tmDRpUvzyl7+M1157LfbYY4845phjKn27e80118SPf/zjmDNnTrRo0SLOO++8Gmtr27ZtTJgwIf7yl7/EHXfcEb/4xS/iJz/5Sc6Yv/3tb/HEE0/E//7v/8b//u//xgsvvBA333xz9vExY8bEH/7wh5gyZUpMmzYtXnzxxXjttdeqfc4kSeL444+PpUuXxtSpU+PVV1+NAw88MI466qgt/rb6xz/+cTzwwANx//33x0svvRQfffRRPP744zljrr322njwwQfjnnvuibfeeiu+8Y1vxFlnnRUvvPBCRFQ07KeeemqMGDEi5s2bFxdeeGFcc801W1QHAFB3+jx93qb0eZAnCdDk/OY3v0l23HHHpHXr1sngwYOTsWPHJq+//nr28YULFyYRkcydOze77OOPP04iIpk+fXqSJEkyffr0JCKSRx99NDvmww8/TNq0aZM89thjSZIkyRlnnJEce+yxOc89cuTIpLi4OHt/1KhRyUknnZQz5vLLL0+OOOKIJEmSpLy8PCkqKkp+8YtfVLkvVdX62e1+8sknScuWLZNf//rX2cfXrFmTdOnSJbn11ltz9ue5557Ljnn66aeTiEj+/e9/V/ncVbn11luT/v37Z+9ff/31yXbbbZeUl5dnl33rW99KDjnkkOz+tWzZMvmf//mf7OPLly9Ptttuu+Tyyy/PLuvRo0fyk5/8JEmSJPn973+ftGvXLlm1alXOc+++++7Jz3/+8yrrqu51KikpSW6++ebs/bVr1yZdu3bNee1at26dzJo1K2e9888/PznjjDOSJEmSq666KunTp0/O49dcc00SEcnHH3+cs7yq9xsASI8+T5+3kT4P8sMRbdAEnXLKKfH+++/HlClT4phjjokZM2bEgQceuFWTyg4aNCj79w4dOkSvXr3i7bffjoiIt99+O+fxz46vjbfffjtWr14dRx111BbXttHf/va3WLt2bRx66KHZZS1btoyDDz44W+tG+++/f/bvJSUlEVExqXB1fvOb38Rhhx0WnTt3jh122CGuu+66WLRoUc6YXXfdNdq2bZuz3Y3b/Pvf/x5r166Ngw8+OPt4cXFx9OrVq9rnfPXVV+OTTz6JnXbaKfvN9Q477BALFy6Mv/3tbzW9FDnKysqitLQ05z1p0aJFDBgwIHv/L3/5S6xatSqOPvronOd66KGHss81f/78OOigg3K2ven+AAANR5+nz4vQ50E+uRgCNFGtW7eOo48+Oo4++uj47ne/GxdccEFcf/31MXr06GjWrCKDTzaZw2HT+Sg2J5PJVFq/Os2aNas0btPnatOmTa2ftzobt7+xrk2Xf3ZZy5Yts3/f+NiGDRuq3O7LL78cp59+etxwww1xzDHHRHFxcTz66KPx4x//uNptbtzuxm3WVFt1NmzYECUlJZXmUYmIaN++fbXrbY2NdT799NOxyy675DxWVFQUEVW/jrV57wGA+qHP0+fVhj4P6ocj2oCIiNh3331j5cqVERHxuc99LiIiSktLs49/9rLqG7388svZv3/88cfxf//3f7H33ntnt7np458dv/G5Nn2ezz7XxolXq7tUeatWrSKiYrLX6uyxxx7RqlWreOmll7LL1q5dG3PmzIl99tmn2vU25w9/+EP06NEjrrnmmhgwYEDsueee8e67727RNnbfffdo2bJldg6SiIjy8vJYsGBBtesceOCBsXTp0mjRokXsscceObedd9651s9dXFwcJSUlOe/JunXr4tVXX83e33fffaOoqCgWLVpU6bm6desWERF77713zoTDERFz5sypdR0AQP3S5205fZ4+D7aWI9qgifnwww/jy1/+cpx33nmx//77R9u2bWPOnDlx6623xkknnRQRFd8uDhw4MG6++ebYdddd44MPPohrr722yu3deOONsdNOO0WnTp3immuuiZ133jl7FajLLrssBg8eHLfeemuMGDEinn322fjtb3+bs/7nP//5+NGPfhQPPfRQDBo0KB5++OF48803o1+/fhFR8Y3sVVddFd/+9rejVatWceihh8a//vWveOutt+L888+Pjh07Rps2beK3v/1tdO3aNVq3bl3pku/bb799/L//9//iW9/6VnTo0CG6d+8et956a3z66adx/vnnb/Vruccee8SiRYvi0UcfjYMOOiiefvrpShPMbk7btm1j1KhR2do6duwY119/fTRr1qzSt4cbfeELX4hBgwbFiBEj4pZbbolevXrF+++/H1OnTo0RI0bknBKwOZdffnncfPPNseeee8Y+++wTt99+e84VpNq2bRtXXnllfOMb34gNGzbEYYcdFuXl5TFr1qzYYYcdYtSoUXHhhRfG7bffHldddVWcf/75MW/evOzpKdXtAwCQPn2ePm9T+jzIk4aeFA7Ir1WrViVXX311cuCBBybFxcXJdtttl/Tq1Su59tprk08//TQ77i9/+UsycODApE2bNskBBxyQPPvss1VOkvvUU08lvXv3Tlq1apUcdNBBybx583Ke7/7770+6du2atGnTJjnxxBOT2267LWeS3CRJku9+97tJp06dkuLi4uQb3/hGcskll2QnyU2SJFm/fn3ygx/8IOnRo0fSsmXLpHv37skPf/jD7OO/+MUvkm7duiXNmjXLrvfZyVj//e9/J5deemmy8847J0VFRcmhhx6a/OlPf8o+vnF/Np3Ude7cuUlEJAsXLqz29fzWt76V7LTTTskOO+yQjBw5MvnJT36Ss3/XX3990rdv35x1fvKTnyQ9evTI3i8vL0++8pWvJNttt13SuXPn5Pbbb08OPvjg5Oqrr86O2XSS3I3rXHrppUmXLl2Sli1bJt26dUvOPPPMZNGiRVXWWd0kuWvXrk0uv/zypF27dkn79u2TMWPGJOecc07Oa7dhw4bkjjvuSHr16pW0bNky+dznPpccc8wxyQsvvJAd8+STTyZ77LFHUlRUlBx55JHJPffcU+UEwybJBYD6o8/T521Knwf5kUkSJ1gDFJKVK1fGLrvsEj/+8Y/r9E3spv7xj39Ez549Y+7cuXHAAQekss2a3HTTTXHvvffG4sWLc5aPHj06li9fHk888US91wAAUGj0ebDtM0cbQJ7NnTs3Jk6cGH/729/itddeizPPPDMiInuKR5oGDx4cgwcPTn27d999d7zyyivx97//PX71q1/Fj370oxg1alT28RdffDF22GGH+PWvf536cwMAFCp9HjQ95mgDKAC33XZbzJ8/P1q1ahX9+/ePF198cYsmvN2crl27Zife3XgVqTQtWLAgfvCDH8RHH30U3bt3j29+85sxduzY7OMDBgzITn68ww47pP78AACFSp8HTYtTRwEAAAAgBU4dBQAAAIAUCNoAAAAAIAWCNgAAAABIgaANAAAAAFIgaAMAAACAFAjaAAAAACAFgjYAAAAASIGgDQAAAABS8P8BhH7C9V4DQ2IAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1500x800 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize =(15,8))\n",
    "ax1.set(title = \"NW vs angle\",\n",
    "        xlabel = \"Subduction angle [deg]\",\n",
    "        ylabel = \"Force [N/m]\",\n",
    "        );\n",
    "\n",
    "ax1.plot(theta_slab_array, NW_array, 'or');\n",
    "\n",
    "\n",
    "ax2.set(title = \"R vs angle\",\n",
    "        xlabel = \"Subduction angle [deg]\",\n",
    "        ylabel = \"Force [N/m]\",\n",
    "        );\n",
    "ax2.plot(theta_slab_array, R_array, 'ob');\n",
    "\n",
    "# Fill the other plot (ax2) with the resistance force results!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Friction along the subducting interface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets consider now the input of the friction along the subduction interface ($FS$) onto the force balance of our imaginary plate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Friction along the subduction interface (FS)\n",
    "def fs_force(f, rho_up, g, Hbdt, theta_interface):\n",
    "    return f * rho_up * g * (Hbdt**2) / (2 * np.tan(np.radians(theta_interface)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\nbt571\\AppData\\Local\\Temp\\ipykernel_7132\\3634833124.py:8: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  forces_df[\"FS\"][\"45 deg\"] = fs_force(f, rho_p, g, Hbdt, theta_interface)    # Fill the new column with our results\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>NW</th>\n",
       "      <th>R</th>\n",
       "      <th>NSP</th>\n",
       "      <th>FS</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>45 deg</th>\n",
       "      <td>4.434974e+13</td>\n",
       "      <td>6.722059e+12</td>\n",
       "      <td>3.762768e+13</td>\n",
       "      <td>7105000000000.000977</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>90 deg</th>\n",
       "      <td>3.136000e+13</td>\n",
       "      <td>4.753213e+12</td>\n",
       "      <td>2.660679e+13</td>\n",
       "      <td></td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  NW             R           NSP                    FS\n",
       "45 deg  4.434974e+13  6.722059e+12  3.762768e+13  7105000000000.000977\n",
       "90 deg  3.136000e+13  4.753213e+12  2.660679e+13                      "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# More parameters\n",
    "theta_interface = 45    # Angle of the subduction interface [degrees]\n",
    "Hbdt = 50e3             # Depth of the brittle-ductile transition [m]\n",
    "f = 0.2                 # Friction coefficient\n",
    "\n",
    "\n",
    "forces_df[\"FS\"] = \"\"                                                        # Add a new column to the dataframe\n",
    "forces_df[\"FS\"][\"45 deg\"] = fs_force(f, rho_p, g, Hbdt, theta_interface)    # Fill the new column with our results\n",
    "\n",
    "forces_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From what we've discussed in class, if you were to try out a other friction coefficient values, what range of values would you use? Do you expect to force to compare better to the NSP force, or drift apart even more?\n",
    "\n",
    "As a last task, do a force balance for the plate (the summation of all forces must add to zero). If the residual is not close to zero (bare in mind the orders of magnitude we are dealing with), what other force could be missing to achieve a state of equilibrium? Considering the order of magnitude, of the missing forces, which do you think is the most relevant one (i.e. the one that can account for most of the residual)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'1.95018e+13 N/m'"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L_conv = 1e3\n",
    "balance = NSP - forces_df[\"FS\"][\"45 deg\"]\n",
    "f\"{balance:.5e} N/m\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py38",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}