{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Code calibration Example\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simply supported beam is to be designed with a given target reliability. A self weight $G$ is uniformly distributed over the span $l$ of the beam. The time-variable load is acting at mid-span and represented by the yearly extreme value distribution $Q$.  The basic random variables that represent the situation are given in the table below. The reliability requirement is introduced as a target reliability index with a one year reference period $\\beta_t=4.7$. The resulting beam span and cross-section modulus are dependent on the design philosophy that is followed. The following formats are to be considered: reliability based design, load and resistance factor design and global partial safety factor.\n",
    "\n",
    "![Image of the beam](img/simplebeam1.pdf)\n",
    "\n",
    "|                                    | $G$        | $Q$  | $F_m$  |\n",
    "|------------------------------------|------------|---------|---------|\n",
    "| Distribution                       | Normal     | Normal  | Normal  |\n",
    "| Mean ($\\mu_X$)                 | 2 N/mm     | 10000 N | 200 MPa |\n",
    "| Standard dev. ($\\sigma_X$)    | 0.30 N/mm | 4000 N  | 20 MPa  |\n",
    "| Coeff. of var. ($COV_X$\\)      | 15\\%      | 50\\%   | 10\\\\%   |\n",
    "| Characteristic fractile ($k_X$) | 0.50      | 0.98   | 0.05   |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import scipy.stats\n",
    "## some general plotting settings:\n",
    "fontsizes=18\n",
    "plt.rcParams.update({'font.size': fontsizes})\n",
    "plt.rcParams.update({\"font.family\": \"serif\"})\n",
    "plt.rcParams.update({\"mathtext.fontset\" : \"cm\"})\n",
    "plt.rcParams.update({'font.serif': 'Times New Roman'})\n",
    "plt.close('all')\n",
    "\n",
    "## RANDOM VARIABLES\n",
    "#G: Self-weight\n",
    "mG = 2        #mean\n",
    "sG = 0.3      #standard deviation\n",
    "VG = sG/mG    #coeff of variation\n",
    "kG = 0.50     #characteristic fractile\n",
    "chG = sp.stats.norm.ppf(kG,mG,sG)  # characteristic value\n",
    "\n",
    "#Q: Variable load\n",
    "mQ = 10000    #mean\n",
    "sQ = 5000     #standard deviation\n",
    "VQ = sQ/mQ    #coeff of variation\n",
    "kQ = 0.98     #characteristic fractile\n",
    "chQ = sp.stats.norm.ppf(kQ,mQ,sQ)  # characteristic value\n",
    "\n",
    "#F: Material strength\n",
    "mF = 200      #mean\n",
    "sF = 20       #standard deviation\n",
    "VF = sF/mF    #coeff of variation\n",
    "kF = 0.05     #characteristic fractile\n",
    "chF = sp.stats.norm.ppf(kF,mF,sF)  # characteristic value\n",
    "\n",
    "#OTHER INPUT\n",
    "bT = 4.7          # Target reliability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The limit state function for bending moment at mid-span is (with $W$ the elastic section modulus):\n",
    "\\begin{equation*}\n",
    "g(f_m,g,q)=W\\cdot f_m-g\\dfrac{l^2}{8}-q\\dfrac{l}{4}\\le 0\n",
    "\\label{eq:lsf} \\tag{1}\n",
    "\\end{equation*}\n",
    "\n",
    "The corresponding reliability index is:\n",
    "\n",
    "\\begin{equation*}\n",
    "\\beta = \\dfrac{W\\mu_{F_m}-\\mu_G\\dfrac{l^2}{8}-\\mu_Q\\dfrac{l}{4}}{\\sqrt{\\left(W\\sigma_{F_m}\\right)^2+\\left(\\sigma_G\\dfrac{l^2}{8}\\right)^2+\\left(\\sigma_Q\\dfrac{l}{4}\\right)^2}}\n",
    "\\label{eq:beta} \\tag{2}\n",
    "\\end{equation*}\n",
    "\n",
    "## Reliability Based Design\n",
    "\n",
    "The beam cross-section $W$ can be chosen such that the reliability index of the beam corresponds to the target reliability, i.e. $\\beta=\\beta_t$. This corresponds to the so-called **reliability based design**. E.g. for $l=6$m, the section elastic modulus giving $\\beta=\\beta_t$ is found to be $W=3.01\\cdot 10^5$. The cross-section elastic section modulus is plotted against the beam length:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Define a function that gives W = W (beta_t) (solving eq.2 for W)\n",
    "def reldesW(l):\n",
    "    a = mF**2 - bT**2*sF**2\n",
    "    b = 2*mF*(mG*l**2/8+mQ*l/4)\n",
    "    c = (mG*l**2/8+mQ*l/4)**2 - ((sG*l**2/8)**2+(sQ*l/4)**2)*bT**2 \n",
    "    return (b + np.sqrt(b**2-4*a*c))/(2*a)\n",
    "\n",
    "## Compute W(beta_t) for several lengths l\n",
    "ls = np.linspace(6e3,30e3,50)\n",
    "W_bT = []\n",
    "for l in ls:\n",
    "    W_bT.append(reldesW(l))\n",
    "\n",
    "W_bT = np.array(W_bT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEHCAYAAABiAAtOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd7hU1dnG4d8rgiIgAir2EPNFsYECRmwEE8GGHQuKiEaKNRbsvXfBHtGgBKUoFixIET2CFCk2LFiJChppQUCQ+n5/rCEexgOeM7Nn9uyZ576uc41nz5qZdznow96rbHN3REREorRe3AWIiEjxUbiIiEjkFC4iIhI5hYuIiERO4SIiIpFTuIiISOTWj7sAqZxNN93UGzVqFHcZGfnpp5+oVatW3GXkTSn1t5T6CqXV3ylTpsxx980yfb3CJSEaNWrE5MmT4y4jI2VlZbRu3TruMvKmlPpbSn2F0uqvmX2dzet1WUxERCKncBERkcgpXEREJHIKFxERiZzCRUREIqdwERGRyClcRERkDb17Z/8eChcREfmfu++Gbt2yfx+Fi4iI4A7XXgs9esBxx2X/fgoXEZES5w4XXgg33ACnnw4DBmT/ngoXEZEStnIldOkCvXrB3/8Ojz4K1apl/74KFxGRErVsGZx0Evzzn3DNNdCzJ6wXUSpo40oRkRK0eDG0bw+vvgp33hnGWqKkcBERKTHz50O7djBuXJh23KVL9J+hcBERKSE//AAHHwwffQSDBkUzM6wiChcRkRLx9dfQpg3MmAEvvQQHHZS7z1K4iIiUgGnTQrAsXAgjR8K+++b28xQuIiJFbvJkOOSQMMX4zTehadPcf6amIouIFLFRo+CAA6B2bRgzJj/BAgoXEZGiNXgwHHooNGoEY8fCH/+Yv89WuIiIFKHeveH446FFCxg9GrbaKr+fr3ARESki7nDLLWFn40MOCYP39erlvw6Fi4hIkVi1KmxAeeWV0LEjvPACbLRRPLVotpiISBFYtgxOOw369w8bUN5zT3T7hGVC4SIiknALF4Z9wkaMgFtvhUsvBbN4a1K4iIgk2KxZcNhh8O670KdPOHspBAoXEZGEmj4d2raFmTPD+Eq7dnFX9AuFi4hIAr33XpgNtnQpvPYa7LNP3BWtSbPFREQS5vXXoVUrWH99eOutwgsWULiIiCTKgAFhy/zttgv3Y9l557grqpjCRUQkIe6+O9yWeO+9wxnLttvGXdHaKVxERArc6sWRPXqEm3sNHw6bbBJ3VeumAX0RkQK2dCmcemq4a+R550HPnvEujqysgi3RgqvNbIaZzTezIWa2dVqbncxssJndl3psVsH7JK6NiAiEe90ffHAIljvugF69khEsUNhnLj2AecDhwF+A24EBQCsAM9sWKAOOcvfxZrYjMM7M9nL3L5LaRkQE4JtvwlTjzz+HJ5+Ek0+Ou6KqKcgMNDMDxrj7g+7+rrvfTQiWPco1uw341t3HA7j7p8D7wF0JbyMiJe6996Bly7A4cvjw5AULFGi4eDAh7fBs4BkAM6sJHANMTGszEWhnZvWS2ObX/yZEpNSMGAH77//LGpYDDoi7oswU8mWx/zGz+sBWwBmpQ82BDQmBU973QDWgGbA0gW1GISIl64knoEuXsHZl6FDYeuvffEnBKuhwMbPawAnAlcBGQBvgeWCLVJO5aS9ZmHrcHFiewDZrMLOuQFeAhg0bUlZWlt4kERYtWpTY2jNRSv0tpb5C7vrrDn37NqJv30a0aDGP6677iM8/X8nnn0f+UXlT0OECLAFGAjWA64CnzaxxuecXp7WvlnpcBlgC26zB3XsDvQFatGjhrVu3Tm+SCGVlZSS19kyUUn9Lqa+Qm/4uWwZnnAH9+kHnztC7d32qV98/0s+IQ0GHi7uvBL4BHjazyaTGJ4DV4zHp4xQbpx5n8cv/rJPURkRKyH//C8ceC2+8ATfeGO4gGfd9WKJS0OFSnrtPMrOpwM/AJ4Szmi3Tmm2Ten4KoW9JayMiJWL6dDj0UPjyy3DW0rFj3BVFKzHhklIdeNPdF5jZYFJrXsppBrzg7osBkthGRIrfxIlw+OGwfDmMHAl//nPcFUWvIKcim1ktM7vUzHYvd+xsYLC7T0sdugH4v9VjMGa2C7AzcEW5t0piGxEpYs8/D61bQ61aYVfjYgwWKNwzlzpAB+AGM5tIuAxW5u4Prm7g7l+Y2UHAzWb2JbA1cIC7T09yGxEpTu5w113h/vZ/+hMMGQING8ZdVe4UZLi4+3+A3SvR7m3g2GJrIyLFZflyOPtsePTRsKtx375Qs2bcVeVWQV4WExEpFvPnh4H7Rx+FK66AgQOLP1igQM9cRESKwfTpcNhh8MUX8PjjYR1LqVC4iIjkwLhxcNRRsGJF2C+shNaaArosJiISuSefDBtO1q0L48eXXrCAwkVEJDKrVoVV9qecAvvsAxMmwI47xl1VPHRZTEQkAosXQ6dO8OyzYa+wBx+EGjXirio+ChcRkSzNnAlHHgnvvAN33w0XXFA8e4RlSuEiIpKFyZNDsCxYAC++CO3axV1RYdCYi4hIhgYODHeNrF4dxo5VsJSncBERqaJVq+Dqq6FDB2jRAiZNgiZN4q6qsOiymIhIFfz0Uxi4f+45+Nvf4KGHSnvgfm0ULiIilfTDDxuw774wdSr06gXnnaeB+7VRuIiIVMKYMdC9e3PcYehQOOiguCsqbBpzERH5Db17w1//CnXqrGDCBAVLZShcRETWYvVW+d26hXB56KF3aNw47qqSQeEiIlKB2bPhwAPDgP3FF8PLL0Pt2iviLisxNOYiIpLm/ffDwsgffgibUJ58ctwVJY/OXEREyhk4EPbeO2yVP3q0giVTChcREUKYXHJJWBjZvHnY1mXPPeOuKrl0WUxESt68eXDiiTByJJx5ZljDooWR2alSuJhZnyw+y939b1m8XkQkch98EO4YOXMmPPZYWHUv2avqmUs74KMMP2vnDF8nIpITgwbB6afDJpvAm29Cy5ZxV1Q8qhoub7r7cZl8kJk9k8nrRESitmIFXHZZuPfKPvvA4MGw5ZZxV1VcNOYiIiVl1iw44QQoK4NzzgkBo/GV6FU1XG7N4rOyea2ISNbefhvat4c5c6Bv37C7seRGlaYiu/s7mX5QNq8VEcnWo49Cq1aw/vowbpyCJdd+M1zMbD0z65T6qZ06doiZvWRmN5iZ1sqISMFasiTMAOvaFVq3DutX9tgj7qqKX2WC4R7gQqAHMNHM/gzcBSwE/grcmbvyREQy99VXsO++0KcPXHll2Cq/QYO4qyoNlRlz+d7ddwcws0bALUBTd1+ROnZzzqoTEcnQK69Ax47hn196Sfe3z7fKnLnMMLPfA7j7v4H+q4Ml5T+5KExEJBMrV8I114Qw+d3vYMoUBUscKhMuC4DPzKwpgLu/DGBmn5jZw6nnRURiN2cOHHoo3HgjdO4M48fD9tvHXVVp+s3LYu7+kpk1cPf0ELkAWOXuI3JTmohI5Y0fD8cfH9axPPIIdOmi+9vHqVLrXCoIFtx9WPTliIhUjTvcdx/06AHbbhtCplmzuKuSnE0jNrNTcvXeIiIAP/4Ixx0H558fLodNmaJgKRRZbf9iZg2By4HGQPkNFKoBTYF+2by/iMjavP9+WG0/fTrccUc4c9FlsMKR7d5iLwBNgPeBpeWOe5bvKyJSIfewNf5550G9evDGG7D//nFXJemyDZemwH4Vbe1iZl2yfG8RkTUsXAjdu0P//tCmDfTrBw0bxl2VVCTbMZfRwOy1PDcoy/cWEfmfDz6AFi3CPe5vugmGDVOwFLJsw6U7cNJantOAvohkbfVlsL32ggULYNSosJXLetrVsKBle1msHtDBzLqlHV8P2Ap4MMv3F5ESpstgyZVtuPQHGgCvA0vKHa9G2NRSRCQj77wTbur11Vdhxf0VV+hsJUmyDZeGwD7uPi39CTM7Jsv3FpES5A4PPBCmFm+2WbhjpGaDJU+2fw94mrVPOx6d5XuLSImZNw+OOSZMM27bNqxlUbAkU7bhchFwcPpBMzPgzCzfW0RKyNix4SZer7wCPXvCiy/q3itJlm24fA3cY2Yry/8AK4Drsq5ORIreypVhTKX8LYjPP1+r7ZMu2zGXp4DfE1bol7/HS3XgqEzf1MxqAjcCHYANCZfYLnL3r8q12SnV5jvCzLRb0hdzJrGNSCn59ttwQ6/Ro+Hkk+Ghh2DjjeOuSqKQbbg8Bixx9y/SnzCzbLbi70UYy7mQsL1MD6CZmTV19/lmti1QBhzl7uPNbEdgnJnttbqWJLYRKSXPPx/ubb98OfzrX3CKVsYVlawui7n71HX8j/F3mbynmdUHprt7d3cf5O5XAucA2/HL2dBtwLfuPj5Vx6eEs6e7yr1VEtuIFL0lS+Css8LA/fbbhynHCpbik1W4mFlDM+tlZsPM7PVyP28C92fx1vem/f5M6rFB6pLZMcDEtDYTgXZmVi+JbdL/BYgUo/feg+bN4eGHw1TjcePgj3+MuyrJhYLbFdnd51VwuHrqcTTQnDAOk76n2feExZvNUrUkrc2o9E6LFItVq6BXL7j8cqhfH4YPD1ONpXglZVfkw4BX3X2SmbVPHZub1mZh6nFzYHkC2/yKmXUFugI0bNiQsrKyipoVvEWLFiW29kyUUn8r09c5c2pw++2NmTy5PvvsM4dLLvmUGjWWk8R/RaX03WYr23DJ+a7IZrYBcBZh5lh5i9N+r5Z6XAZYAtv8irv3BnoDtGjRwlu3bl1Rs4JXVlZGUmvPRCn197f6OmRI2Bts8WL4xz+ga9dNMds0fwVGrJS+22wlYVfk24Br3f3L1O/fph7TxylWT2CcldA2IkVj0SLo2hWOOgq22y4M2nfrprUrpaSgd0VOve9H7v5qucOfEDbJ3DKt+TbAz8AUQr+S1kakKEyYEGZ/ffklXHopXH89bLBB3FVJvhXsrshmdhJQ0917lTtWB1gJDAZapb2kGfCCuy9OtU1cG5EkW74cbr453Mhr663DhpOt0v/ES8koyF2RzawTcCJwn5mt3rtsE6AjYUrvDcBkM2vs7tPMbBdgZ6BzubdJYhuRRPr887DSfuLEcNZy//1Qt27cVUmcsg2XyHdFNrPOQB/CIPghaU8/4O7LgC/M7CDgZjP7EtgaOMDdp69u6O6JayOSNO5hoP6ii8Klr6efhuOOi7sqKQTZhstFwBnAp+UPltsV+caqvqG7PwE8UYl2bwPHFlsbkaSYORMuvbQJkyaFNSt9+oTLYSKgXZFFJAMDB8Juu8HUqXV56CEYNkzBImsqyF2RRaQwzZ0LZ58NgwZBy5Zw9tmT6dhxr7jLkgJUqLsii0iBefll6NIlBMwtt8DFF8Nbby357RdKSarSZTEzW2M3oHXtiuzuawzop79WRJLhxx/htNPg8MPDPe0nTgx7hK2f7V9NpahVdczljCw+K5vXikgMRoyAXXeFfv3gyith8mTYffe4q5IkyHZAX0SK0MKFYbuWgw6COnVg/PiwOLJGjbgrk6So6ontoWb21W83q9BmGb5ORPLotdfgjDPgm2/CuMoNN8CGG8ZdlSRNVcPlzpxUISKxW7AghEnv3rDDDvDWW7DPPnFXJUlVpXBx9+tzVYiIxGfEiHC2MnNmCJjrr4eaNeOuSpJMYy4iJezHH0OoHHQQ1KoFY8fCHXcoWCR7CheREvXSS7DzzvD442Fr/HffDQsjRaJQ1XUut+aqEBHJj9mz4aST4IgjoEGDcP+V227ToL1Eq6pnLqem7qkiIgnjDgMGhLOVwYPDuMrkybDnnnFXJsWoquHSAHjQzP6Ui2JEJDdmzIAjjwxnLNtvHy6BXXON1q1I7lQ1XF52907AVmZ2pZml3xNeRArIqlXw4IPhbOW11+Duu2HcONhll7grk2JX1anIx6YeX0htTHmOmX3n7k/mpDoRydjHH4eNJseNgzZtwk29tt8+7qqkVGQ8W8zdF7v7HcC7ZnarmTWJsC4RydDSpXDddWEPsGnToG9fGD5cwSL5VdXZYr8aa3H3j9z9cqCJmV1uZrUjq05EqmTMGNhjjzBYf/zxIVw6dQKzuCuTUlPVM5fLKzqYCpQxwGRgiJmdkG1hIlJ58+aFxZCtWsHixTB0KDz5ZNgiXyQOVd1brK2ZDQbqE2aObZp6rF6ujQEHmFlndz8kmjJFpCLu0L8/XHBBCJhLLgmzwGrVirsyKXVVDZcNgOaEQPkcGA/MBeakHueW+/2H6MoUkXRffAFnnhlmgf3pTzByJDRtGndVIkFVw2Woux9hZpsCfwc+dPdBOahLRNZi6VK4/fZwq+EaNeCBB6B7d6hWLe7KRH5R1TGXHgDuPsfdrwb+bWZ3mlnj6EsTkXSvvQa77QbXXhsWRU6bBmefrWCRwlOlcHH3z9J+fxu4BNjPzC41s42iLE5Egu+/hw4dwnqVVatg2DAYNAi22iruykQqVtWpyBelH/PgMaAPcKmZHR1VcSKlbsUKuP9+aNwYnnsunLF8+GHYIl+kkFV1zKW9mfUmDOhX9FMfuNXMugDnuvuXURYrUkrGjYOzzoL334cDDwzbuOywQ9xViVROVcNlL2B+2rHyy7M89bwBvYDDMy9NpDTNmhXur/LEE7DNNvDMM3DssVoIKclS1XCZDTzBmlOO55T753nuvirKAkVKxcqVYf+vq66CRYtCwFx1FdTWnheSQFUNlxfc/dKcVCJSwt56C849F957D/761zDOstNOcVclkrmqzhbrlqtCRErRd99Bx46w//4wZ06YATZypIJFki/jXZFFJHPLlsEdd8COO4YxlSuvDGtWjj9eYytSHKp6WUxEsjR0aNgL7LPP4PDDoWdP+MMf4q5KJFo6cxHJk2nT4NBD4bDDwoaTQ4fCiy8qWKQ4KVxEcmz+/HCmsttuMHZsuNXwhx/CIdozXIqYLouJ5MjKlfDYY2E68dy54X4rN90Em28ed2UiuaczF5EcGDky3BGye3fYeWeYMgV691awSOnQmYtIhKZNg8sv340JE+D3v4enn4b27TUDTEqPzlxEIjB3blgEueuuMHVqXe64Az7+GI47TsEipUlnLiJZ+PnncLOum2+GBQuga1do2/Ztjj5637hLE4mVzlxEMuAOAweGlfQXXwwtW4bdix9+GOrVWx53eSKxU7iIVNGYMSFMOnSAunVhxAh49dVwSUxEAoWLSCVNmwZHHw2tWsHMmWFL/ClTwt0hRWRNCheR3/Ddd2EsZZddYNSosFbls8/g1FN173qRtdGAvsha/Phj2FyyZ89wu+FzzgkLIjfbLO7KRAqfwkUkzc8/h5t23XRTmGLcoUP45+23j7sykeQo6HAxs9bA5UB/d++b9tyWwD2EO2BuCTzi7iOT3kbis2IF9OsH110H33wTbtp1++3QvHnclYkkT8GOuZjZIcBpQFvA0p6rA4wGhrj7uUA34Ckz2zfJbSQe7vDcc9CkCZx+OjRsCK+9Fn4ULCKZKdhwcfdXgdvX8vRlwCbAoFTbucDLwIMJbyN55B4CpGVLOPbY8Puzz8Lbb4ezFhHJXMGGS8qStRw/GZjk7l7u2ESgqZk1SXAbyZOxY+EvfwnTiL//Hvr0galT4ZhjtF2LSBQKPVw8/YCZbQP8Dpid9tT3qcc9k9gGyYspU8INu/bbDz75BO69N0wrPu00WL+gRyBFkiWJ/zltkXqcm3Z8Yepx84S2+RUz6wp0BWjYsCFlZWUVNSt4ixYtir326dNr8fjjjRgzZjPq1FlO167fcNRRM6lZcxUTJkT7WYXQ33wppb5C6fU3G0kMl9UWp/2+ejnbsoS3+R937w30BmjRooW3bt26omYFr6ysjLhq/+QTuP76sPV97dpw7bVwwQXVqVv3D0Bu7i8cZ3/zrZT6CqXX32wkMVy+TT3WSzu+cepxVkLbSIQ+/RRuuAEGDICNNoLLL4cLL4QGDeKuTKQ0JDFcZgEzCetEytsm9fhWQttIBD77LCx4fOop2HDDsGNxjx5aVS+Sb4U+oP8rqdlWfYD9zNaY19MMmODu05PYJvN/IwJhU8mOHcMW+IMHwwUXwPTpYRGkgkUk/wo9XKqnHtPPsHoRxinaAJhZQ+Bg4PyEt5Eq+vhjOOmkcJ/655//JVTuukv3qxeJU8FeFjOzvYHTU792MrNZ7v4igLvPM7O/ALeY2YGEy0vt3f3t1a9PYhupvA8+CHd/fOaZMKZy8cVw0UUKFJFCUbDh4u7jgfFAl7U8/xnQ/jfeI3FtZN0mTQqhMmRImP112WVhoH7TTeOuTETKK9hwESlv9OgQKiNGQL16YXPJc8+F+vXjrkxEKqJwkYLlDsOGwa23hlsLb755GKA/80yoUyfu6kRkXRQuUnBWrgxjKbfdBu+/D9tuG7ZpOeOMML4iIoVP4SIFY+lS6Ns33P3xyy+hcWN4/PEwG6xGjbirE5GqULhI7ObPD3d+vPde+M9/YM894c474cgjYb1CnywvIhVSuEhsZsyAXr3gkUdg0SJo2xaefDJsha9t70WSTeEieffhh2GR41NPhUH7E04I61R23z3uykQkKgoXyQt3GDkS7r47TCfeaCM466ywor5Ro7irE5GoKVwkp5YuhWHDtuC888KdHrfYIqxX6d5da1REipnCRXJizpwwlvLgg/D9943Zddcw86tDB9hgg7irE5FcU7hIpD76KMz66tcPfv45DNJfcMH79OjRVIP0IiVEEz0la6tWwauvhiDZddcQLJ06haAZPhz23PO/ChaREqMzF8nYggXwxBPwwAPw+eew1VZwyy3QpYs2khQpdQoXqbJp00Kg9O0b1qe0bBk2kmzfXivpRSRQuEilrFwJr7wSBuhHjAghcuKJYWfiFi3irk5ECo3CRdZp9mx47LGwPcs338DWW8ONN0LXrroxl4isncJFfsUdJkyAhx6Cp5+GZcvCliw9e8IRR8D6+lMjIr9B/5uQ/1mwIOzt9Y9/hAWPG28M3bqF+6fstFPc1YlIkihchClTwoLH/v3hp59gjz2gd++w4LF27birE5EkUriUqIULYcCAECJTpkDNmiFMuncPA/RalyIi2VC4lBB3mDgRHn0UBg4MZym77hpW1HfqBJtsEneFIlIsFC4lYO7csL39P/8JH3wQdiQ+8cSw2HGvvXSWIiLRU7gUqZUrYdSoECgvvBBmfDVrFgbrO3QIg/UiIrmicCkyX30VVs4//jh8+23Y1r57dzj9dGjaNO7qRKRUKFyKwMKFMHhw2Odr9Ohwmatt23BjriOO0Bb3IpJ/CpeEWrUK3ngjnKU8+ywsXgw77BA2juzYEbbdNu4KRaSUKVwS5qOPwpb2Tz0FM2ZA3bpwyinQubMG50WkcChcEuKHH8KA/LvvQrVqcPDBcNdd4bJXzZpxVycisiaFS0LMmAENG0KvXmG2lzaNFJFCpnBJiF12gcmT465CRKRydJvjhNhww7grEBGpPIWLiIhETuEiIiKRU7iIiEjkFC4iIhI5hYuIiERO4SIiIpFTuIiISOTM3eOuQSrBzGYDX8ddR4Y2BebEXUQelVJ/S6mvUFr93dHd62T6Yq3QTwh33yzuGjJlZpPdvUXcdeRLKfW3lPoKpdVfM8tqTxBdFhMRkcgpXEREJHIKF8mH3nEXkGel1N9S6iuUVn+z6qsG9EVEJHI6cxERkcgpXCQyZtbazIab2akVPLelmQ0ws/vNbLCZtYmjxiitq7+p53uZmZf7eTnfNUbBzGqa2V1mNtPM5prZ82a2fVqbnVLf632px2Zx1ZuNyvQ11e6FtO/2gTjqzZYFV5vZDDObb2ZDzGzrtDYZfbeaiiyRMLNDgBOBtsCAtOfqAKOBq919oJk1AD4xs6PdfWz+q83euvqben4LoDFwcbnDQ/NTXeR6AQ5cCDQBegDNzKypu883s22BMuAodx9vZjsC48xsL3f/IraqM7POvgKY2e7AKtb8bgfmu9CI9ADmAYcDfwFuJ/x5bgWQ1Xfr7vrRTyQ/wM6E/zA7px2/GZhNaowvdawP8F7cNeeiv6nn7gF2i7vGCPpYH7gs7ViX8v0GngImp7V5HXgh7vqj7mu5/m4ed70R9NeAlmnH+gEL0/qa0Xery2ISpSVrOX4yMMlTfzJTJgJNzaxJ7svKmQr7mzpr6Qb0NrOrzGzT/JYVuXvTfn8m9djAzGoCxxC+z/ImAu3MrF6ui4vYWvsKYGZNgeOAAWZ2gZnVymdxUfJgQtrh2aT6nO13q3CRKP1q6qGZbQP8jvCHtrzvU4975rqoHFrbVMtmwEhgS+BGYKqZtc1bVRFy93nunh6i1VOPo4HmwIZU/P1WI/y7SIRK9BVgd2A4sAvh7PRdM2uepxJzyszqA1sB56UOZfXdKlwk17ZIPc5NO74w9bh5HmvJC3cf6u5HuXsj4FBgBTDEzHaKt7LIHAa86u6TKP7vt3xfcfe+7n444S8OHQn9G25mSd6eqbaZ/Q2YDLQGVk+2yeq7VbhIvixO+71a6nFZvgvJJ3d/FTgAWAmcG3M5WTOzDYCz+HVfiu77XUdfV19SeooQPvWB0/JcXpSWEM607ySMwzxtZn8o93xG363CRXLt29Rj+vXZjVOPs/JYSyw8zKrpA/zht9omwG3Ate7+Zer3Yv5+0/v6Kx5mO75Mgr9bd1/p7t+4+8NAO8Is4nZk+d1qKrLk2ixgJuEyQnnbpB7fym85sfkU2CjuIrJhZt2Aj1JnY6t9Qvibb0Xf78/AlDyVF6m19HVtPgV+zHFJeeHuk8xsKuG7y+q71ZmL5FRqhlgfYD8zs3JPNQMmuPv0eCrLu+bAP+IuIlNmdhJQ090fK3esDuFy32BS6yLKaUaYrpp+SaXgra2v65gZthPwr7wUlx/VgTfdfQFZfLcKF4nS6pk16WfEvQjXZ9sAmFlD4GDg/PyVlhMV9tfMHjezi8xsAzNbz8y6A1PdPav7Y8TFzDoRBq+nmdnBqZ8TCYvtlgM3AP9nZo1T7XchrAG6Iq6aM/Ubfa1mZs+a2WlmVs3MapjZtUA/d/8m1sIzYGa1zOzS1KLQ1cfOBga7+7TUocy/27gX8uinOH6AvYFHCdNzRwNHpD2/A+FvQXcA/YH94645V/0lhOkCwpTNIcABcdebRT87E9XUN3oAAAJ/SURBVFajewU/95drtxfwbOr7fQrYNe7ao+4rYSB7EPAT8G/gaWD3uOvOor9bAO8BS4ExhF2QT6qgXUbfrXZFFhGRyOmymIiIRE7hIiIikVO4iIhI5BQuIiISOYWLiIhETuEiIiKRU7iIiEjkFC4iIhI5hYuIiERO4SJSwMxsEzN7xsxmm9l9a2lznZn9ZGaDUz8HZ/F5zcu9z2wzeyLj4qWkact9kQLm7vPN7ATC3QDHr6PpbHdvH8HnTQHaA5hZWbbvJ6VLZy4iha8ZsAnwetyFiFSWwkWk8LUBPnT3H+IuRKSyFC4iha8N8FpVXmBmTczsXjOba2a1zewRM/uvmU1Ljav83syGm9kSM5tgZul3GxTJisJFpICZ2UbAPsCoKr50CbA9UB/oQrjHzL5AA8KdQU8EugL7A02AyyMqWQTQgL5IoWtFuEnV6Kq8yN0/N7P3gXbu3nP1cTN7HdjL3W9NHfo61a5xVAWLgM5cRApdG2CSh/uZV9WKCo4tqeDYUqBGBu8vslYKF5HC1oaqXxITiZ3CRaRAmVlDYDdglJk1M7OWcdckUlkKF5HCtV3q8W3gUODDKr5+fQAzq1bu2Hr8+r97I4zriERG4SJSuN4FhgC3AtPdfVFlX2hmfwaOTv16lZltbmbHAgcCW5nZuamtZboCewBNzezUiOuXEmbuHncNIpIFM7sO6OzujSJ+3zLg3+7eOcr3ldKgMxcREYmcwkVERCKnRZQixWEzMxuc+ufH3H1YJm9iZs35ZbX+LsC/I6hNSpDGXEREJHK6LCYiIpFTuIiISOQULiIiEjmFi4iIRE7hIiIikVO4iIhI5BQuIiISOYWLiIhE7v8BJRwl+VungG8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(ls/1e3,W_bT,'-b')\n",
    "ax.set_xlabel('$l$ [m]')\n",
    "ax.set_ylabel('$W$ [mm$^3$]')\n",
    "ax.set_xlim(ls[0]/1e3,ls[-1]/1e3)\n",
    "plt.tight_layout()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Design value method and partial factors\n",
    "\n",
    "Although modern structural design codes (like Eurocodes) allow reliability based design, they propose simpler design procedures in the **semi-probabilistic design format**, where design equations are written in the **load and resistance factor design format (LRFD)**. This format is referred to as **semi-probabilistic**, since random variables (represented by distribution functions) are reduced to design values, which are obtained from characteristic values as $r_d=r_k/\\gamma_R$ (for resistance) and $s_d=s_k/\\gamma_S$ (for actions). The format is simpler, since the designer can check the design just by comparing design values, without the necessity of using reliability analysis.\n",
    "\n",
    "The design equation of the beam above is: \n",
    "\n",
    "\\begin{equation*}\n",
    "W\\cdot f_{m,d}-g_d\\dfrac{l^2}{8}-q_d\\dfrac{l}{4}\\le 0$\n",
    "\\label{eq:deseq} \\tag{3}\n",
    "\\end{equation*}\n",
    "\n",
    "where:\n",
    "\n",
    "- $f_{m,d}=\\dfrac{f_{m,k}}{\\gamma_R}$; $g_{d}=g_k\\gamma_G$ and $q_{d}=q_k\\gamma_Q$ are the **design values**;\n",
    "\n",
    "- $f_{m,k}=F_{F_m}^{-1}(k_{F_m})$; $g_{k}=F_{G}^{-1}(k_{G})$ and $q_{k}=F_{Q}^{-1}(k_{Q})$ are the **characteristic values** corresponding to the characteristic fractiles $k_{F_m}$, $k_{G}$ and $k_{Q}$, which are defined arbitrarily and\n",
    "\n",
    "- $\\gamma_G$; $\\gamma_R$ and $\\gamma_Q$ are the **partial safety factors (p.s.f.)**.\n",
    "\n",
    "The minimum cross-section plastic modulus satisfying the design criterion above is:\n",
    "\\begin{equation*}\n",
    "W_{min}=\\dfrac{g_k\\gamma_G \\dfrac{l^2}{8}+q_k\\gamma_Q \\dfrac{l}{4}}{\\dfrac{f_{m,k}}{\\gamma_R}}\n",
    "\\label{eq:deseq2} \\tag{4}\n",
    "\\end{equation*}\n",
    "\n",
    "## Calibration of partial factors\n",
    "\n",
    "The partial safety factors for a given span, say, $l=6$m, are computed based on the $\\alpha-$values (FORM sensitivity factors), which are:\n",
    "\n",
    "$\n",
    "\\alpha_{F_m}=\\dfrac{-W\\sigma_{F_m}}{k} \\, \\, ; \\, \\,\\alpha_{G}=\\dfrac{\\left(\\sigma_G\\dfrac{l^2}{8}\\right)}{k} \\, \\, ; \\, \\, \\alpha_{Q}=\\dfrac{\\left(\\sigma_Q\\dfrac{l}{4}\\right)}{k}  \n",
    "$\n",
    "\n",
    "with:\n",
    "$\n",
    "k=\\sqrt{\\left(W\\sigma_{F_m}\\right)^2+\\left(\\sigma_G\\dfrac{l^2}{8}\\right)^2+\\left(\\sigma_Q\\dfrac{l}{4}\\right)^2}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The FORM sensitivity factors for l = 6 m are:\n",
      " alpha_F = -0.70\n",
      " alpha_Q = 0.71\n",
      " alpha_G = 0.13\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def alphafun(w,l):\n",
    "    aF = -w*sF/np.sqrt((w*sF)**2+(sG*l**2/8)**2+(sQ*l/4)**2)\n",
    "    aQ = sQ*l/4/np.sqrt((w*sF)**2+(sG*l**2/8)**2+(sQ*l/4)**2)\n",
    "    aG = sG*l**2/8/ np.sqrt((w*sF)**2+(sG*l**2/8)**2+(sQ*l/4)**2)\n",
    "    return np.array([aF,aQ,aG])\n",
    "\n",
    "alpha_display = alphafun(W_bT[0],ls[0])\n",
    "print(\"The FORM sensitivity factors for l = 6 m are:\")\n",
    "print(\" alpha_F = {af:.2f}\\n\".format(af=alpha_display[0]),\"alpha_Q = {aq:.2f}\\n\".format(aq=alpha_display[1]),\"alpha_G = {ag:.2f}\\n\".format(ag=alpha_display[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The partial factors can now be computed based on the $\\alpha-$values, the reliability target, and the definition of the fractile value.\n",
    "\\begin{equation*}\n",
    "\\gamma_m=\\frac{\\left(1+ \\Phi^{-1}(k_F) V_F \\right)}{\\left(1+ \\alpha_F\\beta_t V_F \\right)} \n",
    "\\label{eq:gF} \\tag{5.1}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "\\gamma_G=\\frac{\\left(1+ \\alpha_G\\beta_t V_G \\right)}{\\left(1+ \\Phi^{-1}(k_G) V_G \\right)} \n",
    "\\label{eq:gG} \\tag{5.2}\n",
    "\\end{equation*}\n",
    "\n",
    "\\begin{equation*}\n",
    "\\gamma_Q=\\frac{\\left(1+ \\alpha_Q\\beta_t V_Q \\right)}{\\left(1+ \\Phi^{-1}(k_Q) V_Q \\right)}\n",
    "\\label{eq:gQ} \\tag{5.3}\n",
    "\\end{equation*}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The partial safety factors for l = 6 m are:\n",
      " gamma_F = 1.24\n",
      " gamma_Q = 1.31\n",
      " gamma_G = 1.09\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def psf_alpha(w,l):\n",
    "    alpha = alphafun(w,l)\n",
    "    gR = (1 + sp.stats.norm.ppf(kF)*VF)/(1+alpha[0]*bT*VF)\n",
    "    gQ = (1+alpha[1]*bT*VQ)/(1+sp.stats.norm.ppf(kQ)*VQ)\n",
    "    gG = (1+alpha[2]*bT*VG)/(1+sp.stats.norm.ppf(kG)*VG)\n",
    "    return [gR,gQ,gG]\n",
    "\n",
    "psf_display = psf_alpha(W_bT[0],ls[0])\n",
    "print(\"The partial safety factors for l = 6 m are:\")\n",
    "print(\" gamma_F = {af:.2f}\\n\".format(af=psf_display[0]),\"gamma_Q = {aq:.2f}\\n\".format(aq=psf_display[1]),\"gamma_G = {ag:.2f}\\n\".format(ag=psf_display[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the partial factors and the characteristic values from above, **a one to one correspondence between reliability based design and partial factor design can be established**. However, this is only possible for one **specific span of the beam**. For different spans, the $\\alpha-$values would change and with them the partial factors, as demonstrated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEHCAYAAABiAAtOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUVfrA8e9JCAQCCYSEAAECoSRK7036giBgQeQnqJRVsGBvYHd1VewF7GtfhEV01UUREGnSQXpPgEDoPQk9yfn98aYTQkgmc2eS9/M889zMnTszby4T3jn3nPMeY61FKaWUciUfpwNQSilV/GhyUUop5XKaXJRSSrmcJhellFIup8lFKaWUy2lyUUop5XKlnA7AHUJCQmzt2rWdDqPInDx5koCAAKfD8Ch6TnKn5yV3el4utHLlysPW2tCCPr9EJJfatWuzYsUKp8MoMnPnzqVr165Oh+FR9JzkTs9L7vS8XMgYE1eY5+tlMaWUUi6nyUUppZTLaXJRSinlcppclFJKuZwmF6WUUi6nyUUppZTLlYihyEopVdKdSznH0dNHOXLqCIdPHebI6SMcOXVE9qX/fOZoxr7C0uSilFJeJjk1maOnj3Lo5CEOnzrMoVOyTU8ch08flm2WfYnnEi/6emV8y1C5XGUql61McNlgokKi2MCGQsWoyUUppRx2PuU8h08d5uDJgxw8eZBDpw5x6OQhDp06dMH9w6cOc+z0MSy5L/RYoXQFQsqFULlcZULKhRBVOUrul62csS/95/Rt2VJlMcZkex3zfybX188vTS5KKeVi1loSzyVyIOkAB08e5MDJtG2O++nJ42KXoXyMD6HlQgkNCCW0XChNw5pm3A8pF0JoOdmGlAshNCCUymUrU6ZUGTf/trnT5KKUUvl06vwp9iftz3Y7kHRAtiezb88kn8n1NYLLBlMloAphAWE0rtKYKgFVst2yJpNKZSvhY7xz3JUmF6VUiWat5cT5E6w/uJ59ifvYl7Qvc5v2c3oiya3fwmAIDQglLCCMsPJh1K9cX35Oux8WECbJpHwYIeVCKO1b2oHf0v00uSiliq2EswnsTdzLnoQ97Encw97EvexN3Mu+pH2ZPyfu42zKWViU/bkBfgFUq1CNauWr0bxac6oGVKVq+apUq1CNquXl57CAMEIDQinlo/+V5qRnRCnldVJSUzh48iB7EvcQnxDPnoS0bdr9vYl72ZO4h6RzSRc8N7BMINUrVKd6hep0rNmR6hWqc/LASTo160T1CtWpVl6SR4UyFRz4zYoPTS5KKY+SalM5kHSA3Qm7iU+IZ/eJtG1C5nZv4l6SU5OzPa+UTymqV6hOeIVwmoQ1oXe93oRXCCc8MDxjf/UK1QkofeG6LXPnzqVro65u+g1LBk0uSim3SjybyK4Tu4g7EceuE7vYdWIXuxN2yzYtkZxPPZ/tOWV8y1AjsAY1g2rSOaIzNSrIz+EVwqkRWIPwwHCqBFTx2s7v4kiTi1LKZay1HD51mJ3HdxJ3Io6443EZP6cnlONnjmd7jq/xzUgc7Wq0o1ZQLWoG1qRmUM2MbeWylS+Yh6E8myYXpVS+ZU0eO4/vZMfxHRk/pyeRU+dPZXtOYJlAIoIiiKgYQceaHYmoGEGtoFrUCqpFRFAEVctXxdfH16HfSBUVTS5KqWxOnjvJjuM72HFsB9uPbWfH8cztjmM7OHn+ZLbjg8sGU7tiba4IvYI+9foQUTGCiKAIalesTUTFCCr6V3ToN1FO0uSiVAljreXAyQPEHo0l9lgs249tz9wejeXAyQPZjg/wCyCyUiSRlSLpUacHdSrWoXbF2tSpVIeIoAiC/IMc+k2UJ9PkolQxlGpTiU+IJ+ZoDLFHY4k5GkPMsRhW71rNgUUHsrU+DIYagTWoG1yXvvX7ElkpkrrBdalTsQ6RlSIJKRei/R3qsmlyUcpLWWvZl7SPrUe2su3INrYdTbsd2Ubssdhs5UdK+5YmslIkVcpUoX+j/tStVJe6wXWpW6kutSvW9ph6VKr40OSilIc7ceYEW45sYeuRrRfcsrZASvuWpm6lutSvXJ/e9XpTP7g+9YLrUS+4HjUCa+Dr4yvzObp2de6XUSWGJhelPEBKagpxJ+LYfHgzWw5vke0R2WbtA/ExPtSuWJuoylF0juhMg8oNqB9cn/qV61MzsKaOulIeQ5OLUm50+vxpth7ZyqbDm9h8eDObDm9i06FNbD2yVepbpQkuG0x0SDTX1L+GqMpRRIVEEVU5ishKkXoJS3kFTS5KFYGkc0lsOrSJjYc2svHQRjYdlp+3H9uesciTj/GhTsU6RIdEc3Xdq4kOiSY6JJqoEFncSSlvpslFqUI4ee4kmw5vYsPBDaw/uJ4Nhzaw4dAGdp3YlXGMn48fUSFRtKzeklub3MqVoVcSHRJNg8oN8C/l72D0ShWdkpFc4uJg/36oWtXpSJSXOp9yni1HtrDuwDrWH1zPuoOy3Xl8Z0ZLpLRvaaJDoulYsyOjWoziytAruTL0SuoG19WS7KrEKRmf+MOHoX59GDMGHn4YypVzOiLloay17Encw9oDazNu6w6uY8vhLRnFFH2NL1EhUbQOb83wZsNpGNqQRlUaaRJRKgvH/hKMMV2BJ4BvrbVfXeLYmsAEoBtwFPgUeMVam5qvN2vUSJLLM8/Axx/DSy/BrbeCj1ZQLclOnz/N+oPrWXNgTbZkcuzMsYxjagXVonGVxvSr349GVRrROKwxUZWjtFNdqUtwJLkYY/oANwO9gEmXONYfeAX4ABgHPA38E9gHfJ6vNyxTBn74AebPh0cegWHD4N134c03Qcf8lwj7k/azev9q1uxfw+oDst1yZAupad9PypcuT+MqjRnUcBBNwprQJKwJjao00rpYShWQI8nFWjvdGBMHDM3H4ZHAaGvtCQBjzCDgCND8st+4c2dYuhQmTYInnoBu3eC66+C116BBg8t+OeV5Um0q249tZ+6hucycPZNV+1exat+qbHNFagXVolnVZgy8ciDNqjajaVhT6lSqo2uBKOVCTl4gPp2fg6y1G3PcP2mMSQKmFuhdfXzglltgwAB4+2145RVo2BBGj4Znn4Xg4AK9rHK/lNQUthzZwsq9K/lr31+s3LeS1ftXk3guEZCVCa8MvZKr611N86rNMxJJpbKVHI5cqeLPyeRiC/KktJbLS9baeYV697Jl4ckn4fbbJamMHw9ffy0/33MPlC5dqJdXrpWSmsKmw5tYuXclK/au4K/9f7F6/+qMtUPKlipL06pNua3JbTSv1pyUPSkM6zNMh/oq5RBjbYH+jy/8GxtTG9gBjLDWfpmP46OAkcADwCLg/6y1+/M4fhQwCiAsLKzl5MmT83z9gO3bqfvhhwSvWMGp8HC233UXhzt2BC+oBpuUlET58uWdDsNlUm0q8afj2ZK4JeMWkxTDmVQpxOjv40/98vVpUKFBxrZWuVr4mszSJ8XtnLiKnpfc6Xm5ULdu3VZaa1sV9PnelFwCgAbAnUjSmGet7Zaf92rVqpVdsWLFpQ+0Fn77TTr9N22CLl3grbegRYv8vI1jvLkYobWW+IR4lu9dzrI9y1i+dzkr9q4g4WwCIC2S5tWa06paK1pVl1uDyg0uWUPLm89JUdLzkjs9LxcyxhQquXjNoHxr7UlgFXCXMSYZGG2MCUrv6HcJY6BPH+jZEz79VC6RtWolo8teegmqV3fZW5VUJ86cYPne5SyNX8qyvctYtmcZ+5OkAern40fTqk25pfEttK7emtbhrYkOida5I0p5IW/9q/0cuBs4XySvXqoU3H03DBkCL78M77wDU6bIJMxHH9VJmPmUnJrMhoMbWBK/hKV7lrJ0z1I2HdqUMaO9QeUG9IzsSZvwNrSu3pqmVZtqH4lSxYS3Jhd/YLG19lSRvktQELz6Ktx5J4wdC889B598IiPMbrlFJ2HmcOjkIZbEL2Fx/GKWxC9h2Z5lGeuNhJQLoW14W25ueDNta7SldfXWOmpLqWLMyeTil1sMxpjBwBikw36LMaYV0AH42lp73BhTPu3x+9wWaWSktFwWLoSHHoKhQ+G992Qo81VXuS0MT5KSmsLGQxtZuHshi3YvYnH8YmKOxgAyBLhZ1WaMaDaCdjXa0a5GOyIrRepSuUqVIE7N0G8P/D3t7lBjzEFr7c9p90OAOkCFtPtRwFPA88aYOcBuYIy1drM7YwagY0dYsgS+/VYmYXbqBAMHSusmMtLt4bhT4tlElu1Zli2ZpHe6hwWE0b5me0a2GEn7Gu1pWb0l5fz00qFSJZlTM/QXA4uRocU5HxsPjM9yfyIw0X3RXYKPj9QlGzBAyseMGwc//wwPPABPPSWX0oqB/Un7+XPXnxm31ftXk2JTMBgaVWnEkEZD6FCzAx1rdaROxTraKlFKZeOtfS7OK1dOCmHefrsklTfegC+/hBdegDvukEEBXsJaS8zRGObHzefP3X+yIG4BscdiARkK3K5GO57s9CQda3akXY12BPkXjwSqlCo63vM/oKeqXh2++ALuu0/6Y+6+GyZMkPkxvXo5HV2uUm0q6w+uZ37cfBbsWsD8uPkZw4FDyoVwVa2ruKf1PVxV6yqaV22On6/fJV5RKaWy0+TiKi1awNy58N//wmOPwdVXwzXXSIvmiiscDS0lNYU1B9Ywb+c85sbNZUHcgoyy8jUDa9KjTg86R3Smc0RnoipH6SUupVShaXJxJWOkL6ZvX6lV9uKL0Lgx3HUXPP88hLhnXfSU1BRW71/N3J1zmRc3j/lx8zlxVuaa1guux4ArBtAlogudIzoTUTHCLTEppUoWTS5FoUwZmWw5bJjMjfnwQ5g4Ufpo7r3X5UUxU20qa/avYc7OOczZOYd5O+dlJJMGlRswqOEgutbuSpeILoQHhrv0vZVSKjeaXIpSaCh88IGU83/0UalZ9uGH8Prrso5MAS8/WWvZcmQLf+z4gz92/MGsbbNImC/DgutWqsughoPoVrsbXWp3oXoFLVmjlLq448chJkZusbGZPxeWJhd3aNgQpk+XopgPPww33CALlb31FjRrlq+XiE+IZ/b22czeIbe9iXsB6TNpX7k9Q9oPoVvtbtQMqlmUv4lSygsdPw7btmXeYmIyt0eOZD82PBzq1i38e2pycafeveFvf5MSMs8+K4MA/v53+Oc/oWrVbIeeOHOCOTvnMCt2Fr/v+J2tR7YCMpqre53u9KjTg+51ulO3Ul3mzZtH16ZdHfiFlFKe4uRJSRZbt0ri2Lo18+fDh7MfW7Mm1K8PN94o23r15BYZmVk6sbDjejS5uFupUrIY2ZAh0uE/fjz85z+cH/s4S25qz6w985m1fRbL9iwj1aYS4BdAl9pduLPlnfSo04PGYY11OV6lSqiUFIiLgy1b5LZ1a+Y2Pj77sdWry+rtN9wg23r1JJFERspaiUVNk4tDbFAQMU/exczu5Zkx+xPmnHyWpEngg6F1eBuevOpJetbtSbsa7Sjtq6tiKlWSnDghSWPz5szt5s3SMjl3LvO4ihUhKkquskdFSRJJb4k4vfaZJhc3OnHmBH/s+IMZsTOYGTuTHcd3ABBZO5Jb/fvT88d1dJuzk0rNfeDtfhDR1uGIlVJFxVrYs0fWJdy0SZJH+nZ/ljV2S5WSPpDoaOjXLzOJREXJ7AZPnZamyaUIpdpUVu9fzfRt0/kt9jcW715Mik2hQukKdK/TnUc7PMrVda+mbnBa79ntKfDVV1JOpl07uXT2yitQq5azv4hSqsBSUmDHDti4UZLHxo1y27wZkpIyjwsKkvnWvXtLIkm/RUaCnxcWydDk4mJHTh1hZuxMfov9jRkxMzhw8gAALau1ZEzHMVxd72ra12ife0kVX1/p4L/pJqm0/Oab8MMPMoR57Fjn27lKOSAhIYGDBw9y/nzRrA0IEBQUxKZNmwr1GtZCcjKcP3/hzVqZ3ta0qYzj8fO78Oaby8rdqamuGRaczs/PjypVqhAYGOi6F70ITS6FlN46+XXbr/yy7ReWxi/FYqlctjK96vaiT70+9Krbi7DyYfl/0QoVZATZyJFS2v+ll+Czz2Q7bFjun0KliqGEhAQOHDhAeHg4ZcuWLbLSRImJiVSoUOHSByKJ4uxZOHMGTp/OvJ05I4+BXMoqV046zv39M7f+/s7VtLXWcvr0afbs2QNQ5AlGk0sBnDhzglnbZ/Hrtl+ZHjOd/Un7MRhah7fm2S7P0qdeH1pVb4WvTyGTQESErB1z//1SFPP222V02VtvSQ+eUsXcwYMHCQ8Pp5wDS4tbK62OrAkkPYmkpmYeV7q0JI/AQNmmJxJP+w5ojKFcuXKEh4ezd+9eTS6eYuuRrUzbOo1pW6exYNcCklOTqehfkavrXs019a+hd73eVAmoUjRv3q4dLFoEkyfL5bHu3WWG/+uvy9AQpYqp8+fPU9YN42ZTU6X/I2ciSU7OPMbPTxJHaGhma6RsWc9LIpdStmzZIr3EmE6Ty0WcTznPgl0LMhLKtqPbAGhUpRGPtH+Efg360a5GO0r5uOkUGgODB8P118vyyq+8IjP/R4+mVPfu7olBKQe48lJY+iWt9ORx6pRsz57NvCTm4yNJo1KlzARStqxXLdGUJ3dVPS8mp8s1jp0+xvSY6fy85Wemx0wn4WwCpX1L071Odx5o+wB9G/SldsXazgZZtiw8+aR0/D/9NLz7Lm0//1z6aO66yzuHlShVBFJSsieR9ESS9ZKWv7/0jVSocJaKFctQtqxc5vLU4b3epMQnl9ijsfxv6//4ecvPzI+bT4pNoUpAFQZeMZBro66lR2QPypf2wFFaVavCv/4F991H0ogRVLr/fnj/fVk/pm9f/etQJcr585kJJD2JnDmT+bivrySRkJDsHe3pl7QSE89RoUIZZ4Ivpkpcckm1qazYu4KfNv/ET1t+YsOhDYBc7hrTcQz9o/rTJryN95RYadqUNW++SdekJKm83L8/9Oghw5ibNnU6OqVcylqZoX7yZPYWSdYuhNKlJYEEB2cmEm2NuF+JSC4Wy4yYGfy4+Ud+3vozexP34mt86RzRmZEtRtI/qj+RlSKdDrPgjJGk0ru3lPR//vnMopgvvnhBUUylvEFKiiSNI0eyt0pSUjKPSR+lVa5cZiIpLn0j3q5E/DOs3r+a3hN7E+AXQO96vbku6jr6NuhLcNlgp0NzLT8/GbZ8663SBzN+fOYIs4cfdk+1OqUKIDVVii+uXAkrVsj2r79g6lRpqRiTvTWSnkh8vOQCQ0lUIpJLsH8wXwz+gh51elDWrwT8BxscLHNh7r4bxoyRjv+PP5YRZoMH61+kcpS1sijV8uWSSFaskESSXgrF3x+aN4cRI6ByZbjyStmnH1vvUiKSS0TFCPo16Od0GO5Xv76Uj5k7V1out94K770niadjR6ejUyXE3r2wbJkkk/SEcuyYPFamjKyXN2wYtGwJrVpJfa30S1ubNmWuL1JSJCUlERAQ4LYhw0WlRCSXEq9rV/mL/uYbGcZ81VUwcCCMG+eaJeeUSpOQIJe0li6VhLJsmVT+BRmZ1bixfPRatYLWraFRo4KNnn/wQVi92nVxp6SULfRkyGbN4J13Cv782NhYhg8fzp9//km1atWYOnUqHTp0AKQ8zVdffcW9995buCDdSJNLSeHjI18PBw6UkWSvvgo//SR9NE8/LQtDKHUZUlJgwwZYskSSydKlUu03vb5WvXrQpQu0aSO3Zs202+9ijh07Rr9+/RgwYAAPPPAA69evZ+jQoaxdu5Zy5crx2WefMWDAgEu+TmpqKj4ecv1Qk0tJExAgSyzfcQc884xcIvvyS3juOZ2EqfJ04IAkkvTbihWZ/SSVK0PbtjBokCSS1q1lX1EpTAshN4mJp/NduLIoTJgwgUmTJtGsWTMABg4cSHR0NNOmTeO6664jPj6eWlmW3tiyZQvjx4+nfPnyVKpUCX9/f9q0acOGDRu44447nPo1stHkUlJVry6Vlu+7T+bH3H8/TJgAr70G116rkwJKuORkWLdOStotXiy37dvlsVKlpBUyfLiUvWvXTtYc0Y9MwfXv3z8jsaQbNGgQY8eO5ezZswwZMiRj/8SJE/n000+ZPHkyVdOmGezfv5+GDRsyZ84ct8adF00uJV2zZjBrFvz6qySZ66+XPpo33pAeVlUiHDsmrZGFCyWhLFsmExUBqlWD9u1l8GH79jKFSi9vuVbOxALg4+NDamoqa9as4bbbbgNg9uzZPPjgg2zevJnKWZqGVatW5W9/+xuNGzd2W8yXoslFyVfOvn2hVy/49FO5RNaqFdx2m6whU7Om0xEqF7JWWiF//imJZOFC6TsB6XRv1kzm33boIMmkVi1tlTjl6NGjtEz7kmetZfTo0dx///3ZEku6Bx980KNGmGlyUZn8/OCee+CWW2Qk2dtvw3ffyTDmsWNlETPldZKTYe1aSSY//HAlW7ZkrtEeFCRJ5OabZXR6mzbSLac8Q3JyMoMGDQJgxYoVbNmyhZEjR+Z6bPv27d0Z2iVpclEXCgqSCZd33QVPPQUvvyxFMl94QRYs0/oaHu3MGbmsNX8+LFggrZP0jvewsED+9jcZjX7VVTJB0UMGF6lcBAYGEhoaCshQ5dDQ0Ix+FoDly5fz448/snTpUpo1a8Zrr72mo8WUF4iIgH//Gx54AB55RJLNe+/JImV9+ui1Eg+RmCgd7vPny23pUimZAjKP5LbboFMnSSaxsUvo2rWro/Gq/ImPj+fs2bMZ98PDw0nOunoZ0Lp1a8qXL8/HH3/MjBkzPCaxgCYXlR+tW8O8eTIv5vHHpX+mRw9JMs2bOx1diZOQIJe45s2T4gsrV8qcE19fGYNx//3QubNc5grOUT4vNtaRkFUBzJw5kzZt2mTc79ixI5GRkUydOpWBAwdm7P/mm2/o0qULvh62JKYmF5U/xshIsr59pU7Z88/L/2S33SZFMrXTv8ikJ5O5czOTSWqqdJG1bSvl47p0kc537RYrPubPn8/YsWMz7vv4+DB9+nTGjh3LggULqFmzJv7+/vTu3ZurrrrKwUhzp8lFXR4/P7j3XqlTNm6czGabMkU6/ceMkfrnqlBOnpQRXHPmwB9/ZLZMSpeWZPLUU5nJpKTV3SpJWrduTXR0dLZ9oaGhfPbZZw5FdHk0uaiCqVhRksvdd0u9spdfzhzGPGqUzvS/DOfOyRyT2bMloSxZIuuYlColyeSJJ6BbN0kmOr+k5Bg9erTTIRSKJhdVOBERMHEiPPQQPPaYtGree09ql113nXb65yI1Fdasgd9/l4SyYIEsguXjI1caH35YkknHjlDeA1fYVio/HBtaYIzpaoyZYYwZlo9jGxpjZhljEowxscaYZ40x+tXYk7RqJddwfv5Z/pe84QbpVV661OnIPMKOHfDJJ1J7q0oVmeX++OOwe7eM7v7xR1lxcdkyaRBefbUmFuXdHGm5GGP6ADcDvYBJlzg2BPga+Bb4ErgV+AcQAtxfpIGqy5O+3HKfPlK37LnnpPDUTTfJvJkSVN7/+HHJtbNmyS19lFZ4uJyiHj2ge3cp8aZUceRIcrHWTjfGxAFD83H4AGCQtTYWwBjzLbAAuNMY86i19lwRhqoKolQpuPNOGDJEapS98YZ8Nb/7bqnEHBLidIQul5wsrY4ZM2DmTPk5NVVaH127yvDgXr0gKkqvFKqSwckZN6fzedwv6YkFwFprgalAaUAHXnqyChXgH/+AmBhZs3bCBGm9jBsHp/P7z++54uJkDMPAgRAaKn0k//yn1O566imZ0Hj0KPzvf5JcoqM1saiSw8kOfZuvg6zdk8tuP2CjtfaIa0NSRaJaNZkb88ADUqPsiSfg/ffhxRdlnoyHTf66mLNnJWFMny63zZtlf3g4DBgg/SR/+9uFExeVKomMtfn6P971b2xMbWAHMMJa++VlPncu8IG1dkoex4wCRgGEhYW1nDx5ckFD9XhJSUmU96Le36DVq6n78ccEbt5MUmQk20eN4mibNi79Wu+qc7Jvnz9LlwazbFkwq1ZV4swZX/z8Umna9Dht2hyldeujRESc8poWidd9VoKCqFevXpG/T0pKisfNcC9KMTExnDhxIs9junXrttJa26rAb2KtzfUGXJt263ixYwpzA2ojrZfhl/m87sD/Luc5LVu2tMXZnDlznA7h8qWmWvuf/1hbt661YG23btYuX+6yly/oOTl71trZs619+GFro6MlNLA2MtLa0aOtnTbN2qQkl4Xpdt72Wdm4caNb3ichIcEt7+Mp8nNegRW2EP/H59Xn8iPQElhV4MzlYsaYYOAx4BanY1GFZIyMy924UebFrFsnNcwGD85c8tBN9u+HL76QvpOQEBnJNWGCVLR5+23YskW6jSZMkOo3WpJeqUvLq89lkbX2ufQ7xphXgZuARcBUa+2PRR1cVsaYMsB7wF3W2gR3vrcqQqVLy1LLw4ZJIcy33oLvv5eRZU8/LT3lLmYtrF4N06ZJZ/vy5bI/PFzWNUmvy+lFV4+U8jh5tVyy1U+11o4BTgJDHUgspYH3gZestXFZ9ldzZxyqCAUGSgf/tm2yOHv6yLJ//jNzvd1COH1aVnK++25ZWbFFC5mG4+Mjb7FqlUxo/OQTKSygiUWpwskruaTmsi/GWnvBfmNMVAHeO32GfbbWkzFmsDFmdfprGmP8gSnAViDCGNPbGHONMeZBwLuL76gLVa8u/8OvXy9Dr555RpLMhx9Kwa3LcPgwfPmlFAsICZEWyTffyNW3zz+HffukjtdTT8nSvt7SKa+UN8jrslhuw8guNrRsIPBSft/UGNMe+Hva3aHGmIPW2p/T7ocAdcicwzId6Apcl8tLFXwkg/JsV1wBP/wgq2A9/rgsv/z22/DSS9I5cpFMsG2bLDvzzTfNWL9eJjLWqCGNof79ZUKjv79bfxOlSqS8kss1xpj5OfZdkcs+f6Apl5FcrLWLgcXABYtBW2vHA+Oz3O+W39dVxVD79jK5ZNo0mSMzaJA0PcaNg+7dsVZK0v/3v1IEYONGeVpkZCmeekoucbVooa0S5T2SkpIICAjAePmHNq/LYlWA+kDNLLekHPdrAhFodWVVlNJrlq1dC59/jt23D3r0YFPE1VxTbVVGrqlSRZaX2b4dPvtsBS+8IFWGvfxvVJiAPcoAACAASURBVJUQsbGxdOrUiQoVKhAeHs6iRYsyHktMTGTChAkORnf58koKr1lrx+bxeAZjzD9cFI9SuTp7Fn7/3ZfvF4xg5qmb+T8+4MldLzOdFmxvO5iK418kuHVmYcy4uDxeTHm/Bx+UIX8uUjZ9nejCaNZMvt0UwLFjx+jXrx8DBgzggQceYP369QwdOpS1a9dSrlw5PvvsMwYMGFC4+Nwsr+Ty7WW8ztTCBqJUTqdOSSHIqVPlqlhCggwq69evLB0GPEKZtrfDh68T+fbb0OE7KZb5zDMQFuZ06EpdlgkTJjBp0iSaNWsGwMCBA4mOjmbatGlcd911xMfHU6tWrYzj161bx6RJk6hfvz4JCQns27ePZ555hgAPmoR10eRirV2b1xONMaHASWvtKWvtOpdHpkqkpCT45RdJKL/+KgkmOFj68G+8UeaflCmTfnRF6eAfPRpeeAE++kiGhz34IL7t2jn4W6giV8AWwsWcTkykQgXn6uD2798/I7GkGzRoEGPHjuXs2bMMGTIkY//EiRP57rvv+PrrrwlMW1Z82rRpjBo1iokTJ7o17rxcNLkYY67NcjfVWjstbX9VYDLQCThvjHnbWvtE0YapirP0hPLdd5JQTp+WxsfQoZJUunSRKv4XVb26JJZHHpGWy0sv0S4wUH4ePVrXBlYeL2diAfDx8SE1NZU1a9Zw2223AbBw4UKefvpp1q1bl61G3DXXXMPgwYM5duwYlSpVclvceblU+ZdvgShgA4CR4Qs/Ap2BzcCHwE3GmBFFHKcqZk6ehClTZB2xKlVkZvyff8Lf/w5z58KePTK1pUePSySWrOrXh8mT4a+/SIiOlmWX69eHf/1LFlxRysscPXqUOnXqZNy/9957ue+++y4oPurj44O/vz/btm1zd4gXdan1XAZZa1+31u5Iuz8CaAMsB1pYax8COgCXXKpYqTNnZLjwzTdLQvm//5P147MmlAkTpKVSqL7V5s1Z9+qrMGeOFAgbORIaNZKmUWpuc4OV8kzJyckMGjQIgDVr1rB27VpuueXC0ooHDx7k8OHDBAUFuTvEi8oruWyy1v6afietBMtzyMz9e6y1ZwGstQeB40UapfJa58/L2ifDhsmlrhtugNmz5f6cOS5MKLnp2hUWLZKMVqpU5hyZ336TAmNKebjAwEBC0+rrbdmyhZo1axKWy4CVH3/8kVq1atGgQQN3h3hReSWXnO2r+5F5LZOstStzPKaVmFSG1FSZ93jXXVC1Klxzjcyav/FGGf21bx988IH831/kS2gYIzMp16yBr7+GY8egTx9584ULi/jNlSq4+Ph4zp49m3E/ODiYatUuLKd47tw53nzzTR577DGPmniZV3IpZ4wJAjDG1AOeRQpXZuu8N8ZUBtoXWYTKK1grxR8fewwiIqQl8s03sjrjzz/DgQNSz6tXr8voQ3ElX19Z9XLzZlkFc+tWuOoq6NfPpfMllHKVmTNn0qZNm4z7Xbp0oUzaUMm1a9cyadIkDh06xD333EOtWrW48847nQo1V3kll4+AxcaY14G5SOtkjLU2Pv0AY0wA8DVSAkaVQNu3SzHjK6+UMivvvCNzyb79Fg4elG3//lmHDzusdGmpUxYTA6+8Iq2X5s2lI2jLFqejUyrD/Pnz6dSpU8Z9Pz8/Pv/8c/bs2cO3337LuXPnaN68OUePHuWLL74g2cMGrVw0uVhrfwCeB9oBB4A7rbUfpD9ujBkN/AcIRNZ4USXEkSMy8rdjRylY/Oyz0p/y8cey8Nb//idrfnnQfK4LBQRIrbIdO+DJJ2WW5pVXyugCnd6vPEDr1q2Jjo7Oti8yMpLw8HDGjRvHsGHDuPrqq/nhhx9YunRpRqvGU+Q5WszKGvUPIUUp5+d47H1rbT9rbSdrbadcX0AVG2fOyMTG66+HatVkXZTjx+XLf1ycjPYaNQoqV3Y60stUMW0i5vbtcP/90tSqX18WMNu/3+noVAk2enTeK4pYa/n999958803Wb16NT4+lxr8614XjcYY42eM+Q5YCnwHbDTGfOy2yJTjrJXBVnfeKR3zN90Ey5bJ/8GrVsmSK2PHyuJbXq9KFSnpv20bjBghk2wiI2HMGGmqKeVhjDFER0fzxRdf8OCDDzodzgXySnWPAjcCa5HFulYBdxhjBrsjMOWcuDjpR2nQQC59/fvf0m8yY4as1vjGG8V4ca2aNeX63ubNMGCALL1cpw48/zycOOF0dEplM2PGDNavX09lD7xkkFdy+T9ghLW2ubV2sLW2FXAn0Ns9oSl3SkqCr76Cbt2gdm3pR6lRA774Qq4OffONjPQq8qHDnqJePcmq69bJL/6Pf0iSGTfOJcsuK1Xc5ZVcjllrv8q6w1r7L+B0zgONMTVcHZgqetZmzpCvWlVWa4yPl1bLzp0yyXH4cHCwnp/zGjaUzqa//pJm3BNPyOWyd96RjiilVK7ySi4XW7B8Ty77tDXjRXbvlj7sBg2gc2epipJe22vrVnj6aZmrorJo3lyGwS1eDE2awEMPyVC5Dz6QxWaUUtnkNZ2ttTHmKyDnlfVGxpj6We4HAN2Bf7k6OOU6587JZMZ//QtmzpRWS9euUjj4xhs9fNiwJ2nXDmbNkuFxzz4rVZdffVUy8vDh4OfndIRKeYS8kksQcNtFHstZH1oLNXmojRvhs8+k8snhw9KPkv7/YGSk09F5sa5dYd48+P13ydCjRsm47GefhVtvdagMgVKeI6/LYh8CFay1PnndkJn7n7gnXJUfSUmSUDp0kC6D8eOlHMv06dKX8sILmlhcwhjo2VMulf3yi6xqNmIEXHGFDAZISXE6QqUck1dy+cZae8lhMdbaU8BXlzpOFb1Vq2RyY/XqcMcdMsnxjTekk37qVOjduwSN9nInY6Q65/LlUoE5IEDqmDVsKJMyNcmoEiiv8i9L8vsil3Oscq3ERJg2rRqtW0ttry+/lOkZf/4JGzbI4oxVqjgdZQmRXoH5r78km/v5wS23QOPG8J//6FoyqkTxrHoBKt9WrZKZ89Wrw5tvRnH2rFz+2rtXEkzHjsV0kqM38PGRURJr1shymz4+MhyvSRNdsEyVGJpcvMipU5I42raVVso330hJlvff/4s1a+Dee8FDls9WIEnlpptg7VpZfjk1VRYsa9pUWjaaZFQxpsnFC2zeLNMqwsOlvzgxEd59V1opn38OV16ZoK0UT+bjI2s6r1sHEydCcrIknaZNtSWjii1NLh7q/Hn5ctu9uww+ev99WXhr7lzpS7n/finoq7yIry8MGSIVP9M7+gcN0stlqljS5OJh9u+XocK1a8uX2+3b4eWXZVb95MkypFhbKV7O11cWvFm3DiZNyrxc1qSJdPzr6DJVDGhy8QDWyuiuwYOlKO9zz0GjRrLufGyslLMKC3M6SuVyvr7S0b9unXxzsFbuN26sQ5iV19Pk4qDTp2WyY/Pm0KmTTHK8915ZbXfGDLj2Wp2XUiL4+mb2yUyZIvdvuUVWxvz6a+mjUSVGUlIS1np/0ROtUeGAXbuk3uGnn8LRo/JF9eOP5f8TrfFVgqWPLrvxRmm2vvACDBsm2yeflImZWrsMgAd/e5DV+1e77PVSUlLwLeQ3uWZVm/FO73cK/PzY2FiGDx/On3/+SbVq1Zg6dSodOnQAIDExka+++op77723UDG6k7Zc3MRaKUV1442yLMjrr0t5qrlzZTrEqFGaWFQaHx+44QaZjPnTTzJy4/bbZfnljz7SKszF0LFjx+jXrx+dO3fmu+++Y9SoUQwdOpRTp04B8Nlnn3Httdc6HOXl0ZZLETtzRi6nv/OOJJHgYHjsMSnTomXtVZ6MkWuj/fvLNdMXX5QPzosvwuOPw8iRUK6c01E6ojAthNwkJiZSwcGFiyZMmMCkSZNo1kxqAg8cOJDo6GimTZvGddddR3x8PLW8bD1xbbkUkQMHZGXciAiZm5KcDJ98IqO+xo3TxKIuQ3rtskWLpApz/frw4IPSBH7tNZn4pLxa//79MxJLukGDBrFixQqmTJnCkCFDHIqs4DS5uNiaNZJMatWSlXFbt5blP9atK9FfNJUrGAM9esi11PnzZRLmmDHyTeWFF+DYMacjVAWUM7EA+Pj4kJqaypo1a2jRooUDURWOXhZzgdRU+O03ePNN+OMPSSB33CETHaOinI5OFUudOsmqb0uXyrKizz0nJbDvuQceflirlRYTR48epWXLlgDExMQwbdo04uPj6dChA4mJifz222+MGTOG9evXExsbS0REBMOHD3c26DTacimEM2dkxFejRtC3rwwhfvVVKXH//vuaWJQbtG0rS4yuWSOXzl57TWbgPvCAfBCVV0tOTmbQoEEALFq0iJtvvpmvv/6aTp06MWzYMKy1zJo1i1tvvZUbbriBqVOnOhxxJk0uBXDokFzyqlVLRnn5+8vaUDt2SD+rFo9UbtekiYwc2bRJ5sx88IGsCDdyJMTEOB2dKqDAwEBCQ0MBGDBgAIsXL2bgwIEZ+zZt2pTRUlm/fj1RHvSN1rHkYozpaoyZYYwZls/jyxlj7jPGxBV1bBcTEyODdWrVks76Nm3kMtjKlTJHRacgKMdFRcEXX8iHdeRIKZ0dFSXlH9audTo6dRni4+M5m2XYefny5VmwYAGdOnUC4Pjx45w/fz4j0fz888/07duXEydOOBJvTo4kF2NMH2AE0Au4ZKUsY0wgMBi4F3D7eLylS2HgQGjQQP5ub7tN1qafNg26ddNaX8oDRUTItdmdO2Xs+y+/QNOmNH7iCRl1pjzezJkzadOmTbZ9WZPLwoUL6dixIyCTQJcsWUKnTp2YPHmy22PNjSPJxVo7HXj1Mo5PsNZ+BvxUdFFll5oqf49dukC7djB7ttT42rlThhRfcYW7IlGqEKpWlbHvcXHw4osEbtwoK8l17So1hopBmZHiav78+RmJBCA1NRU/Pz9q1KgBwIYNG+jTpw8gI8tatWrFp59+yo033uhIvDk5OVrstJuec1nOn5dCta+9JqXta9aEt9+WCdIOzrFSqnAqVYKnn2Zxy5Z03rpVSkT07g3NmsHYsdI010J2HqV169ZER0dn3Pfx8WFRllbn448/nvGzMcajOvPB2Q79gnxlKrKvWSdPwnvvQb16Us7Jx0cuV8fGynw1TSyqOEgtW1ZGkm3fLivNnT4tlZijoqTA3ZkzToeo0owePdrpEArFOFV90xhTG9gBjLDWfpnP5zwPPGetzU8/zShgFEBYWFjLi12HTEgoxX//G84PP9QgIcGPJk2OM3jwLtq2Peo1fSlJSUmUL1/e6TA8ip6T3F1wXlJTCVm4kFrffkvg5s2cq1SJ3TfdxN7+/UnxgPMXFBREvXr1ivx9XFG40pvExMRcsuO/W7duK621rQr8JtZaR25AbaQlMvwynvO8hHx579WyZUub05491j78sLUBAdaCtf37W/vnnxcc5hXmzJnjdAgeR89J7i56XlJTrZ0929qePeUPIjDQ2scft3bvXrfGl9PGjRvd8j4JCQlueR9PkZ/zCqywhfg/vsTNc9mxA+66S8oyvfuuFJ9dt07moaUNvFCq5DFG1tSeORNWrJD+mDfekAmZI0fKDGGlLkOJSS6bNsHQoVLz74svpP7X1q3Sr9KokdPRKeVBWraU5Za3bpWRLP/+twyPHDAAlixxOjrlJUpEcomNhYYN4fvvpd7X9u2yLEZkpNORKeXB6taVmf5xcfDUU1Iws3176NwZ/vc/Ga/vBlaHS7uUu86nk8klfT57tuHQxpjBxpjVxpjc6hj4pR1zWUOoExPlbyMuDt56C8LDCxawUiVSlSqyhsyuXTIuPy5O1plp1EjW6S7CxctKlSpFsi7z7FLJycmUKlX0s1CcmqHfHkgfpD3UGJN1ibUQoA5QIcvxPsaYW4Hr03b9wxjTML/v17ix/G2EhBQycKVKsvLlZVx+TIxcKitdWsp/164tEzWPH3f5W/r7+5OUlOTy1y3JEhMT8ff3L/L3cWqG/mJr7UhrrbHWdrbW/pzlsfHW2iBr7Yos+1Kttf+21jZMe85T1toN+X2/EjTCUKmi5+cnxfRWrZIBAI0bS/mKmjXhoYekjIWLhIaGcujQIU6dOqWXxwrJWsupU6c4fPhwRj2yoqTruSilCsYY6NlTbqtXy+iyCRNkNvJNN8Ejj8hqeYXg7+9PWFgY+/fvz1bE0dXOnDnjlm/zTitTpgxhYWFu+V01uSilCq9ZM7lUNm6cJJePP5YRZ507S5Lp10/KXhRAUFAQQUFBLg44u7lz59K8efMifY+SpkSMFlNKuUmNGlKYb/duGT2zcydcd50MZf7oIzh1yukIlZtoclFKuV5goPS/xMTAt99Kcb6775Z+maefhn37nI5QFTFNLkqpouPnJwuVLV8O8+fLZbKXX5b1ZoYPl+WZVbGkyUUpVfSMgU6d4L//lZn/d90FU6dKX02PHrLynpsmZSr30OSilHKvevWk03/3bhkAsGUL9O8P0dGyeqbOaykWNLkopZxRqRKMGSPVZCdNkvv33iv9Mo8/LhUBlNfS5KKUcpafnyxYtnQpLFok82befFOK/w0aBAsX6nLMXkiTi1LKc7RvD1OmSHXZhx6CWbPgqqtkMuY33xRpHTPlWppclFKeJyICXn8d4uOlMvPJk7JmRkQEPP887N/vdITqEjS5KKU8V0CAzI/ZsAFmzJC1Zv7xD6hVC267TYY4K4+kyUUp5fl8fKBXL/jlFxlddued8OOP0KYNtGsHEyfCuXNOR6my0OSilPIuDRrA+PGwZ48MaT52DG69NfOSmc7+9wiaXJRS3ikwEO67T9Ywnz4dWrSQS2YRETBkiIw801FmjtHkopTybj4+0Lu3XDLbtg1Gj4Zff4WOHaWP5vPP4fRpp6MscTS5KKWKj3r1ZCnmPXuk7P/583D77VKt+bHHZIizcgtNLkqp4icgAEaNgrVrYd48qV/29tuSfPr1k5aN1jIrUppclFLFlzFSiXnKFFlb5plnYOVK6NtXEs3rr8Phw05HWSxpclFKlQw1akiHf1ycrJKZXsOsRg2iX3lFys/oAACX0eSilCpZSpeWmmXz5sG6dXD77YQsWCDzZVq2hE8/lYoAqlA0uSilSq5GjeD991k8daqUmUlOlr6a6tWlQvOGDU5H6LU0uSilSryUcuWkzMyaNVKF+dprpQXTqJH02Xz7rRbNvEyaXJRSKp0x0KGDVGDes0c6/PfuhVtuyRzOvG2b01F6BU0uSimVm5AQePRRWZZ55kzo0kWGMzdoIEObv/tO65nlQZOLUkrlxcdHFjCbOlWWZv7nPyE2VgYF1KwJTzwh91U2mlyUUiq/qlWDp56SZPLLLzLC7LXXZM5Mz57amslCk4tSSl0uX1+45hr46SfYtQtefFEunw0aJH0zjz9e4vtmNLkopVRhhIfD009L3bLp02VZ5rfekr6Zbt1krZkzZ5yO0u00uSillCv4+kp15h9+kL6Zl1+WVs2tt8q8mfvvl1pnJYQmF6WUcrVq1aSjf9s2mD1bks7HH0PTprJ65iefQEKC01EWKU0uSilVVHx8oHt3mYS5dy+88w6cOiXLNFerBsOHw4IFxbKmmSYXpZRyh8qV4YEHpJ7ZkiUyMfOHH6QCQHQ0vPoq7N/vdJQuo8lFKaXcyRho21Yuje3bB198AVWqwNixMtLsuutkFNr5805HWiiaXJRSyikBAZmXxrZskYoAy5bB9ddLonn0Ua8tnqnJRSmlPEGDBjBunIw0+9//ZEjzu+9K8cw2beDDD+H4caejzDdNLkop5UlKlZKlmL//XgYBvPWWzJO55x6oWhUGD4YZMyAlxelI86TJRSmlPFVoKDz0kCwFsHw53HGHJJbevaFWLRnuvHmz01HmyrHkYozpaoyZYYwZls/j/26MmWKM+cgY84ExJqCoY1RKKY9gDLRqBRMmyCCA776D5s1lSYArroD27eGjj+DYMacjzeBIcjHG9AFGAL0Ak4/j7wYeAoZYa+8CjgD/KdIglVLKE5UpAwMHwrRp0j/z+uuQmCiLnVWrJvXNpk1zfLSZI8nFWjsdeDU/xxpjKgPjgI+ttclpuz8F+hpjri+iEJVSyvNVqyYjytatg5UrZXLmnDnQv7+MNnv4YVi92pHQnOxzOZ3P464FAoFl6TustbuAA8DQIohLKaW8izHQooWMLtuzR+bJXHWVXEZr3hyaNIE33pABAm7iZHLJb72DjmnbQzn27wNauy4cpZQqBkqXhmuvldFm+/bB++/LfJrHHpPFzXr1kmWck5KKNAxjHappY4ypDewARlhrv8zjuGlAXyDIWpuQZf98oK21tsxFnjcKGAUQFhbWcvLkyS6L3dMkJSVRvnx5p8PwKHpOcqfnJXcl4byU3b2bsFmzCPv9d8ru20eKvz+HOnXiQM+eHG/RAuvrm+34bt26rbTWtiro+5UqdMTucyrHfV/goku+WWs/AT4BaNWqle3atWvRReawuXPnUpx/v4LQc5I7PS+5KzHn5bbbpEjmwoX4fv01VadMoeqsWZnzZ265RS6vmUuOs7okb5jnsjttWynH/kDgoJtjUUop72aM9Md88okUyvz+exnK/P77Mtz5yivhpZcK/TbekFz+SttWy7E/HPjTzbEopVTx4e8PAwZIdeb9+yXhVKkiK2sWkjckl8nASaBz+g5jTB2kJfOtU0EppVSxUqkSjBwJ8+bBzp2Ffjknk4tf2jZbv48xZrAxZrUxJgrAWpsIvAjcbEzGhcC7gB+ttTPcFq1SSpUUERGFfglHOvSNMe2Bv6fdHWqMOWit/TntfghQB6iQfry19tW0vPKlMSYROAPc7MaQlVJKXQZHkou1djGwGBiZy2PjgfG57M/XjH6llFLO84Y+F6WUUl5Gk4tSSimX0+SilFLK5TS5KKWUcjlNLkoppVzOscKV7mSMOQTEOR1HEQoBDjsdhIfRc5I7PS+50/NyoShrbYVLH5Y7bypcWWDW2lCnYyhKxpgVhaleWhzpOcmdnpfc6Xm5kDFmRWGer5fFlFJKuZwmF6WUUi6nyaV4+MTpADyQnpPc6XnJnZ6XCxXqnJSIDn2llFLupS0XpZRSLqfJxYsYY7oaY2YYY4bl8lg1Y8wkY8x4Y8xUY0xPJ2J0t7zOSdrj7xhjbJbbNHfH6E7GmLLGmDeMMXuMMUeMMf81xkTmOOaKtM/Ie2nbFk7F6y75OS9px/2Y4/MywYl43cWIZ4wx8caY48aYn4wx4TmOKdDnpUQMRS4OjDF9kGUGegGTcjxWAZgPPGOtnWyMqQxsMsbcYK1d6P5o3SOvc5L2eFUgGngsy+5f3ROdY94BLPAw0AR4FGhhjGlqrT1ujKkJzAWut9YuTls3aZExpq21NsaxqItenucFwBjTDEgl++dlsrsDdbNHgaNAf6A78Cryt9QZoFCfF2ut3rzkBlyJ/IEMz7H/JeAQaX1oafs+B1Y7HbNT5yTtsbeAxk7H6MZzEQyMzbFvZNbzA0wEVuQ45g9k8T3HfwenzkuWc1PF6XjdeF4M0C7Hvm+AxBznpECfF70s5l1OX2T/LcBym/Yvn2YZ0NQY06Tow3JUruckrdVyJ/CJMeZpY0yIe8NyzLs57n+Xtq1sjCkLDEA+G1ktA/oZYyoVdXAOuuh5ATDGNAVuAiYZYx4yxgS4MzgnWLEkx+5DpJ2bwn5eNLl4lwuG9hljagARyIciq31p29ZFHZTDLjbcsQUwC6iGLJO9zhjTy21ROcBae9RamzPZpi8nPh9oCfiT+2fFFzlnxU4+zgtAM2AG0BBp8a4yxrR0U4gewRgTDFQH7k/bVajPiyYX71c1bXskx/7EtG0VN8biMay1v1prr7fW1gauAZKBn4wxVzgbmdv1BaZba5ejn5Wssp4XrLVfWWv7I19GbkXOxQxjTLEuHQVgjClvjLkdWAF0BdIHAxXq86LJpfg4leO+b9r2nLsD8TTW2ulANyAFuM/hcNzGGFMGuIcLf+cS/VnJ47ykXyqaiCSfYGCEm8Nzwmmklf860g8zxRhTN8vjBfq8aHLxfrvTtjmvfwambQ+6MRaPZWVky+dA3UsdW4yMA56z1sam3dfPish5Xi5gZZTlNErA58Vam2Kt3WWt/RDoh4wi7kchPy86FNn7HQT2IM35rGqkbf90bzgebQtQzukg3MEYcyewIa3Vlm4T8i01t8/KGWClm8JzzEXOy8VsAU4UcUgexVq73BizDvk8FOrzoi0XL5c2Quxz4CpjjMnyUAtgibV2hzOReaSWwEdOB1HUjDFDgLLW2n9l2VcBuSw4lbQ5DFm0QIaW5rz8Uaxc7LzkMTLsCuBrtwTnWfyAedbaBArxedHk4l3SR7jkbHG+g1z/7AlgjAkDegMPui80x+R6TowxXxhjHjHGlDHG+Bhj7gLWWWsLtUaFpzPGDEU6pDcbY3qn3W5GJsadB14A6hljotOOb4jMFXrSqZjd4RLnxdcY870xZoQxxtcYU9oY8xzwjbV2l6OBFyFjTIAxZkza5NH0faOBqdbazWm7Cv55cXoij97yPeGpPfApMvR2PnBtjscbIN8yXgO+BTo5HbOT5wRJuAnIsMmfgG5Ox+uG8zEcmWFuc7mNz3JcW+D7tM/KRKCR07E7eV6QDur/ACeBncAUoJnTcbvhvFQFVgNngQVIFeQhuRxXoM+LVkVWSinlcnpZTCmllMtpclFKKeVymlyUUkq5nCYXpZRSLqfJRSmllMtpclFKKeVymlyUUkq5nCYXpZRSLqfJRSmllMtpclGqEIwxFY0x3xljDhlj3rvIMc8bY04aY6am3XoX4v1aZnmdQ8aYLwscvFJFSEvuK1UI1trjxpj/Q1brW5zHoYestQNd8H4rgYEAxpi5hX09pYqKtlyUKrwWQEXgD6cDUcpTaHJRqvB6AuuttQecDkQpT6HJRanC6wn8fjlPMMY0Mca8a4w5Yowpb4z52BhzzBizL51z3gAAAZxJREFUOa1fpY4xZoYx5rQxZokxJudqgEp5NE0uShWCMaYc0AGYfZlPPQ1EAsHASGT9mY5AZWRl0ZuBUUAnoAnwhItCVsottENfqcLpjCw2Nf9ynmSt3WaMWQP0s9a+nb7fGPMH0NZa+0rarri046JdFbBS7qAtF6UKpyew3Mp645crOZd9p3PZdxYoXYDXV8oxmlyUKpyeXP4lMaWKPU0uShWQMSYMaAzMNsa0MMa0czompTyFJhelCq5W2nYpcA2w/jKfXwrAGOObZZ8PF/5dGqRfRymvoclFqYJbBfwEvALssNYm5feJxpguwA1pd582xlQxxtwI/A2oboy5L620zCigOdDUGDPMxfErVWSMtdbpGJQq1owxzwPDrbW1Xfy6c4Gd1trhrnxdpVxBWy5KKaVcTpOLUkopl9NJlEq5R6gxZmraz/+y1v5WkBcxxrQkc7Z+Q2CnC2JTyuW0z0UppZTL6WUxpZRSLqfJRSmllMtpclFKKeVymlyUUkq5nCYXpZRSLqfJRSmllMtpclFKKeVymlyUUkq53P8DwAh8GgFe9yoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "[gRa,gQa,gGa] = psf_alpha(W_bT,ls)\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(ls/1e3,gGa,'-b',label=r'$\\gamma_{G}$')\n",
    "ax.plot(ls/1e3,gQa,'-r',label=r'$\\gamma_{Q}$')\n",
    "ax.plot(ls/1e3,gRa,'-g',label=r'$\\gamma_{m}$')\n",
    "ax.set_xlabel('$l$ [m]')\n",
    "ax.set_ylabel('PSF')\n",
    "ax.set_xlim(ls[0]/1e3,ls[-1]/1e3)\n",
    "plt.tight_layout()\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4b8eb43fdedf456084f0b55e88d591b8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(FloatSlider(value=6.0, description='$l$ [m] =', max=30.0, min=6.0, readout_format='2.2f'), Butt…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import display, Markdown, clear_output\n",
    "from matplotlib.figure import figaspect\n",
    "import json\n",
    "import ipywidgets as widgets\n",
    "from matplotlib.figure import figaspect\n",
    "\n",
    "slider_l = widgets.FloatSlider(\n",
    "         value=6,\n",
    "         min=6,\n",
    "         max=30,\n",
    "         step=.1,\n",
    "         description='$l$ [m] =',\n",
    "         orientation='horizontal',\n",
    "         readout=True,\n",
    "         readout_format='2.2f',)\n",
    "\n",
    "button = widgets.Button(description='Refresh plot')\n",
    "\n",
    "out = widgets.Output()\n",
    "\n",
    "def plotslider(l,ax):\n",
    "\n",
    "    W_rbd = reldesW(l)\n",
    "    alpha_display = alphafun(W_rbd,l)\n",
    "    \n",
    "    ax[0].pie(alpha_display**2,labels=[r'$\\alpha_R^2$',r'$\\alpha_Q^2$',r'$\\alpha_G^2$'],autopct='%1.1f%%')\n",
    "    ax[0].axis('equal');\n",
    "    \n",
    "    psf_display = psf_alpha(W_rbd,l)\n",
    "    ax[1].bar([r'$\\gamma_R$',r'$\\gamma_Q$',r'$\\gamma_G$'],psf_display)\n",
    "    plt.tight_layout\n",
    "    plt.show()\n",
    "    \n",
    "def on_button_clicked(b):\n",
    "    with out:\n",
    "        clear_output()   \n",
    "        wp, hp = figaspect(1/3)\n",
    "        fig, ax = plt.subplots(1, 2,figsize=(wp,hp))\n",
    "        plotslider(slider_l.value*1e3 ,ax)\n",
    "\n",
    "button.on_click(on_button_clicked)       \n",
    "display(widgets.VBox([slider_l,button,out]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem is solved by finding the set of p.s.f which ensures the reliability to be as homogeneous as possible over a range of possible scenarios. The set can be obtained by solving the minimization problem.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\min_{\\gamma_R^*,\\gamma_G^*,\\gamma_Q^*}\\left\\{ {\\sum}_{i=1}^n\\left(\\beta_t-\\beta_i(\\gamma_R,\\gamma_G,\\gamma_Q,l_i)\\right)^2\\right\\}\n",
    "\\label{eq:gF} \\tag{6}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 0.043783\n",
      "         Iterations: 37\n",
      "         Function evaluations: 72\n",
      "The optimal partial safety factors for l = [6,30] m are:\n",
      " gamma_F = 1.35\n",
      " gamma_Q = 1.13\n",
      " gamma_G = 1.16\n",
      "\n"
     ]
    }
   ],
   "source": [
    "W_psf = lambda psf,l: ((chG * psf[2] * l**2/8 + chQ * psf[1]*l/4)/ (chF/psf[0]))\n",
    "beta_fun = lambda W,l: (W*mF-mG*l**2/8-mQ*l/4)/np.sqrt((W*sF)**2+(sG*l**2/8)**2+(sQ*l/4)**2)\n",
    "\n",
    "def betadiff(psf):\n",
    "    \n",
    "    beta = []\n",
    "    for l in ls:\n",
    "        W_min = W_psf(psf,l)\n",
    "        beta.append(beta_fun(W_min,l))\n",
    "        \n",
    "    return np.sum((np.array(beta) - bT)**2)\n",
    "\n",
    "res = sp.optimize.minimize(betadiff,[1,1,1],method='nelder-mead',options={'xatol': 1e-3, 'disp': True})\n",
    "PSF_opt = res.x\n",
    "\n",
    "print(\"The optimal partial safety factors for l = [6,30] m are:\")\n",
    "print(\" gamma_F = {af:.2f}\\n\".format(af=res.x[0]),\"gamma_Q = {aq:.2f}\\n\".format(aq=res.x[1]),\"gamma_G = {ag:.2f}\\n\".format(ag=res.x[2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this set of partial factors is used for the design of beams with span $l=[6,30]$ the obtained reliability indices are deviating from $\\beta_t$ but the sum of squared differences is minimized. It is interesting to note that the obtained set of partial factors is not very different to the partial factors that are calculated with the design value method above. (see after Equation 5.1-3). The effect on the reliability indices, however, is significant, as seen below where the \"optimal solution\" is compared with the \"design value factors\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEHCAYAAABiAAtOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhV1bn48e+bk3lOGMIckElkECQoCkqYbC2gRav1ZwenipVeW9vqbSu9aq/21vmqrdcWrbW1VqlYQGmrggxBBAVlBtGIzAoBMjJkfH9/rJOczAnJOTkn4f08z3p2zjprn7zZbPJmrb322qKqGGOMMf4UFuwAjDHGdDyWXIwxxvidJRdjjDF+Z8nFGGOM31lyMcYY43eWXIwxxvhdyCUXEVkoIlqt/K6J9r8Qkb+KyB9F5Dci4mmrWI0xxtQvPNgBVCciI4EK4K5q1a800v5BYKSqftX7+kXgKeAHgYzTGGNM4ySUbqIUkZeAH6vq4Wa0PRvYCnxdVRd76y4BVgLnqeqGgAZrjDGmQSEzLCYi5wJXAy+LyI9FJK6JXf4f4AE+qFa3Htfz+U5gojTGGNMcIZNcgJHAW8BQ4HFgg4iMbqT9OECBo5UVqnoCKADGBDBOY4wxTQipYTEAERHgOuBpoAwYoqo59bTbCvRQ1dRa9XuBk6o6uC3iNcYYU1dIXdAHUJftXhKR3cAq4Ebg4Qaan6inzgOU1K4UkVnALICYmJjRvXv39ku8oaaiooKwsFDqkAafHZO67JjUZcekpk8++eSIqnZp6f4hl1wqqepqEVkM9G+gyT6gXz31icDH9XzeXGAuQEZGhq5fv95foYaUFStWkJmZGewwQoodk7rsmNRlx6QmEdnTmv1DPU3vxCWR+nwExIpIYmWF9+t44N02iM0YY0wDQj25DAH+0sB7z+Mu6F9SrW40brbYvADHZYwxphEhkVxEJFpEXhORG0XEIyKRInIv8KKq7vW2uVNE3heRzgCq+hnwB9zF/0q3Ak+p6vY2/yGMMcZUCYnkApTiZob9DvgM+CuwSFWr90DScNdYYqrV3Q4cFJE/iMhcYDPwk7YJ2RhjTENC4oK+qpYD32yizV3UXBYGVS0D7gxgaMYYY1ogVHouxhhjOhBLLsYYY/zOkosxxhi/s+RijDHG7yy5GGOM8TtLLsYYY/zOkosxxhi/s+RijDHG7yy5GGOM8TtLLsYYY/zOkosxxhi/s+RijDHG7yy5GGOM8TtLLsYYY/zOkosxxhi/s+RijDHG7yy5GGOM8TtLLsYYY/zOkosxxhi/s+RijDHG7yy5GGOM8TtLLsYYY/zOkosxxhi/s+RijDHG7yy5GGOM8TtLLsYYY/wuPNgB1EdEngGGqGpmI236AP8F7AcUOAt4QFWz2yRIY4wxDQq55CIiU4FbgaxG2oQD/wZuUdX3vHUXA++IyGBVPdUmwRpjjKlXSA2LiUgScAewtommw4BzgM+q1a0D+njrjTHGBFFIJRfgEWAOUNJEuyO4obBfVaubgBsi2xGY0IwxxjRXyAyLicgMYI+qbhSRRtuq6n4R+T/gByISCzwB3AtMU9WTgY/WGGNMY0IiuYhIJ+AW4MrT2O2HQLl3+w3gfFXdGoDwjDHGnCZR1WDHgIi8ADykqju8r1cANDFbLBr4HZAH3A7k43ou6xpoPwuYBZCWljb6lVde8d8PEEKKioqIj48PdhghxY5JXXZM6rJjUtPEiRM/VNWMlu4f9OQiItcAvVT18Wp1K6DJ5LIQyFLVx0VkPLAYd61mkKrmNfY9MzIydP369X6IPvSsWLGCzMzMYIcRUuyY1GXHpC47JjWJSKuSSyhc0L8NeExEtLLgLs5P8L6+ofYO3mRyBfAigKq+C1wOdAaubrPIjTHG1CsUrrncAtTuiz7n3X4P2FvPPinebdWsMlXNEpENQOOzAYwxxgRc0JNLfXfUi0iR972N3td34nok01T1CLAcl3S+AfzR2yYRiAZeb5vIjTHGNCToyaWZ0oB+QAyAqhZ57+R/SESGAbuBXsDVqvpl0KI0xhgDhGhyqX0hX1XvAu6qVfcJMLMNwzLGGNNMoXBB3xhjTAdjycUYY4zfWXIxxhjjd5ZcjDHG+J0lF2OMMX4XkrPFzJnlVNkp8k/lU1hSSGFxIYUlhRSVFFV9fbzkOCdKT3Cy7KTblp7kRJnbFpcXU1JeQkl5CcVlvq9LyksoqyijXMspOlFExIaIqtflFeUobtmjyuWPKl8DhElYvSU8LLyqRIRFuK3HbaM8UUSFRxHliSLSE1n1dXR4NNHh0cSExxATEVNjGxsRS1xknNtGxBEXGVe1jY+MJ8oTRVMrhBsTqiy5GL8pryjnyIkjHD5+mCMnjtQtJ4+QezKXvFN55J5y27xTeZwqa/6DQ2MjYqt+McdExBAdHk2kJ7KqJEUkEeWJqvqlHx4WTs7hHHp174VHPISHhRMmYYgI4l3MofIXuCAoiqpSoRU1SrmWU67llJaXUlZRRllFGaUV7uvKZFZ0oqgq2RWXFVNcXsypslOcKjvFydKTNRJYc4SHhZMQmUB8ZDwJUQkkRCaQEJVAYlQiSVFJNbfRSSRFJZEcnUxKTArJ0ckkRyeTFJVEhCfitL6vMf5gycU0qUIrOHz8MPsL9rO/YD/78vdxsPAgXx7/ki+LfOXw8cNUaEW9n5EUlUSn2E6kxqSSHJ1Mz8SepETX/CVY/Rdo5TY+Mr7qr/mW/iUfCgsSqiol5SWcLDtZlWxOlJ7geOlxjpccr7EtKimqKpW9t+q9uoOFBykoLqjq7TUlLiKOlJgUUmNSSYl22+K8Yv5Z8k9SY1LpFNuJTjGd6myjwqPa4MiYjsqSi6Gsooz9Bfv5PPdzPs/7nN15u6u2+wv2c6DgAKUVpTX2CQ8Lp1t8N7rHd6d3Ym/G9BhDt/hupMWlkRafRufYzlUlNSaVSE9kkH660CAibqjMz7+wK7SCopIi8k/lk3cqj/zi/KoeYe1e4rGTxzh28hifHvuUL/O+ZNmRZY32GuMj4+kc25kusV3cNq4LnWPctktsF7rGdaVLnNt2jetKXEScDeOZKpZczhBlFWXsydvDp8c+5ZOjn/Dp0U/55Jjb7s3fS7mWV7UNkzB6J/YmPTmd8X3G0yuhF70Se9E7qTe9Et3XnWM7EyY2HyTYwiSMxKhEEqMS6Z3Uu9n7VfbmTpae5OjJoxw9cbTOtnIoM+d4DoePH2ZbzjaOnDjCidIT9X5mTHhMVaJJi08jLS7Nfe39gyMtLo1u8d3oFt+N5OhkS0QdnCWXDqakooQth7awPWc723O2sy1nG9tztpN9LLtG7yMxKpFBnQYxttdYrht+Hf2S+9E3uS/9UvrRO7G3jdOfIWIiYugV4f5gaK4TpSfIOZ5DzgmXdCpLzvEcDh0/VDWE+tEXH3H4+GHKKsrqfEakJ7Iq0XSL70a3OLftntCd7vHdq7bd4rvZudhOnZHJpagI3nsPPB5XwsJ8X4eHQ3S0K1FRvm1YiP2RrqocLDzIxi83sunQJjZ+uZHNhzbz6dFPqVjlrnuESRj9U/oztOtQrhh8BYM6DWJgp4EM6jSILrFd7C9H0yKxEbGkJ6eTnpzeZNsKrSD3ZC6Hjh/iUNEhDh0/VOM63ZdFX7Inbw9r968l53hOvZMeOsd2pkdCD1fie1R93T2hOz0TetIjoQdp8WmEh52Rv85CVtCfRBkMgyVB/8Do09pHBDzVklCNEu6SUmWJqP46AiIjICKi5QlKcdN1C4sLqqbpFpUUUVru+4swOiKa+Ig4PBUeUhM6ERcRS2xErA1dAXl5eSQnJwc7jJASisdE0WrTyiunlHtn33m/Lva+XycHiesNRXmiiPJEElltanhUeGTVVHGPeBr8/qF4TIJJVq5s1ZMoz8hUHxsDIwZ4z091W1X3dYWCVkBF9aJuW15etxSXQPkJKCtzpbFU7fG4JBNRmXAiISoSIiMhMsq7jQTxlHK8rICC4sKqhFKZSESE+Mg4Osd2Ji4invjIeOIj46r+asvLyyM5zv6DmPZHEG9yiCKhkfkfLgmVumRTVpl0fF+fLDtJ3qn8eofjwsM81ZJOzW1xeQnxFWXWA/KTM/Ioes4ZTOr6FX7/XFU35JabC3l5bnvsGOTk1F8OHYJD+5WKinqGp2KOQNJ+Ejrn07NvBYPPiuO8wV0Ze05P+vSPoFcvl6Rq2xgC025DjR2TutrzMREgylsSGml3ovQEBwoOcKDwQNWsx/0F+9lf6L7eV7CPQ0V7aw3F7SQhMoHeSb3pnegtSb5tn6Q+9E7sTUxETEB/xpDQymHzMzK5BIoIJCS40qdP/W3KK8rZ8OUGln2+jFV7V7F6z1pyj3igqDuJpYPpHzGeND2XqBMDKDk2jIMHwtm/Dj5+ExZV+5ywMOjVC/r2rVkKC5MYNAi6d2/1uWFMuxYbEcvATgMZ2Glgg21Kykv4ovAL9hfs5621b5HYK5F9+fvYV+DKxi83cuj4oTr7dY7tTO9El2zqK93iu53xQ9KWXAKsQivYcmgLy3cvZ/nu5azcvZL84nwABncazJXnXMG43uMY12ccA1MHNniR/eRJ2L8f9u2D3btrluXL3Xvu8tkofvQjiImBs86CAQOgf38YOBAGDXKlZ09LPMaAu05TOTmhdFcpmRdl1mlTXFbMgcIDVUlnb/5e9uXvY2/BXnbl7mL57uUUFBfU2CciLIJeib2qkk16UrrbJqdX1cVGxLbRTxkcllwC4IvCL3jrs7d467O3WLprKUdOHAFgQOoArhl6DZP6TSKzbybd4rs1+zNjYlyCGNjAH2ElJbB3L/zjH5uIjz+X7Gz47DPIzoa333bJqVJsrC/RDBoEZ58N55wDgwe794wxPlHhUZyVchZnpZzVYJv8U/lViad62ZO/h5V7VnKg4ECNe8nA9X7Sk1xi65vU1yU57+v0pPR2fy+QJRc/KCkv4b197/Fm9pu8mf0mmw5tAqBbfDe+NvBrTO43mYl9J57WTW6nKzLS9VLOPz+X2kPpFRVw8CB88gns3Om2n3wCH30Er73mJiaA682kp8OQIb4ybJhLPImJAQvdmHYvKTqJpOgkhnUdVu/7ZRVlHCw8yJ68PVVJZ0/eHvbk72F7znb+/em/OVl2ssY+CZEJ9E32JZ2+yX2rtn2T+9I5tnNIJx9LLi2UezKXf376TxbtXMRb2W9RWFJIeFg44/uM58HJD/LVAV9lRNqIkPjHr7w+06sXTJpU872SEvj0U9ixw1e2b3dDbaeqrQzSp49LNNXLkCHuPiBjTOPCw8KrhsPqo6ocOXGEPfl72J23uyrx7M7bzZ78PWTtyaoz9BYbEVuVaPom9fV9HSLJx5LLadidt5tFHy9i0c5FZO3JolzL6RbfjWuHXcu0gdOY1G8SCVGNzV8JPZGRMHSoK9WVl7vrOdu2wdatvrJkCZR6b/T3eNyQ2ogRvnLuudCjh13TMeZ0iIhbsy2uCxk96r+1JO9UXlXi2Z23uyrx7M7bzZp9a8g9lVujfe3k0y/FuwqHdzWO1JjUgCYfSy5NyD6Wzd+3/Z2/b/t71XDX0C5D+c9x/8kVg69gTM8xHXJWiMfjJgL07w+XX+6rLy1113E2b4YtW2DTJli9Gl5+2demUycYNQpGjvRtBw1yN5UaY1omOTqZkd1GMrLbyHrfzz+VX5VsqpfP8z7nvX3vkXcqr0b7ymG3fin96Jfcr8YSUP2S+7U6XvvvXo9dubuqEsqGLzcAcFHvi3js0se4fPDlDEgdEOQIgyciwnc95pvf9NXn5fmSzcaNrjz1lBt2AzchYfhwOO88GD3alaFDXc/JGNN6SdFJjIgewYi0EfW+X9nz2Z23u2oF9M/zPmdX7i7e2fUOx0uP+zUeSy5eXxZ9yd+2/I2Xt77M+oPrARjbayyPX/o43zjnGwG9GN8RJCfDxRe7Uqm0FD7+GDZscMlmwwb429/g979370dGuoRTmWwyMtzr+m4ONca0TmM9H1Xl6MmjvqST+zk/v+/nrfp+Z3RyOVF6goUfL+TFzS/y9mdvU6EVZPTI4JGpj3D1OVc3a2E+07CICJcshg+H737X1VVUwK5d8OGHbrbahx/Cq6/C3Lnu/agoN4yWkQFjxrgyeLAbpjPGBIaIVD1/aUzPMQD8HEsup62wpJCbFt3E/O3zKSwppE9SH34x/hd8Z8R3GNx5cLDD69DCwtyU6QEDfMNqqvD557BuHaxf77Z//jM8/bR7Pz7eJZsLLoDzz3fbnj2D9zMYY5p2RiaXT458whfbv+Dqc67mu+d+l4vTL+6QF+XbCxG3msBZZ/kSTnm5uydn3TpXPvgAHn/cN1OtRw+XZC64AMaOdcknLi54P4MxpqaQTC4i8gwwRFUzm9E2DLgGuBDIBtap6trG9umX0o+td27t8MsvtGcej7t585xz4PrrXd2pU27CwAcfwPvvu7Jgga/98OEu0Vx4odsOHGhToo0JlpBLLiIyFbgVyGpG207AS8AK4MeqWtGc75Eak2qJpR2Kjvb1Vm6/3dUdOeKSzdq1rlSfMNCpk0s0aWnuxrUxY6x3Y0xbCankIiJJwB1Aoz0Pb9sE4C3gn6r6YKBjM6Gpc2f42tdcATec9vHHsGaNK++9B4sXn8Uf/+h6NyNHuoQzbpwrvW0SoDEBEVLJBXgEmAM80Yy2TwEVqnpvYEMy7YnH41tx4Hvfc3WLFr1LePj4qmTz/PPwu9+593r39iWacePcKgM2M82Y1guZ5CIiM4A9qrqxqSUJROQc4HrgURF5DMgACoFfq+qagAdr2pWkpDIyM2HaNPe6tNS3ssDq1ZCVBa+84t5LSICLLoLx4909O+ef724ANcacHlFt7MG8bRSEu3byJ+BKVS0TkRUADV3QF5H7gHuBG4EXcQ+lWwBMAi5Q1Y/q2WcWMAsgLS1t9CuVv006mKKiIuLj44MdRkhp6piowqFD0WzZksjWrUls2ZLE55+79uHhFQwaVMjw4fmMGJHPsGH5JCbWfXxue2PnSV12TGqaOHHih6pa/0JnzRAqyeUF4CFV3eF9vQIaTS6vApepany1uuHAZuAfqnpVY98vIyND169f75fYQ82Kdvz42kBpyTE5dswNob37Lqxa5aZDV06DHj4cLrnEtyJBjx7+jznQ7Dypy45JTSLSquQS9GExEbkG2FyZWJopATcMVkVVt4jIXuAcf8ZnzkypqTB9uivgHra2bp0bQsvKghde8N3k2b8/TJjgEs6ECe5x08ac6YKeXIDbgEzvtZMaRESBG1X1hVpv7fXu41Gt8Xi3L4DSgEVqzlgxMS55XHKJe11W5tZKW7XKJZsFC9xEAXDPvqlMNBMmuNUI7H4bc6YJheRyC1B7oPM57/Z7uERS22ve/cYCq6vVV167MSagwsN9a5/95CduzbRt22DlSlfefhv++lfXtkcPl2QyM12xmzvNmSDoyUVVs2vXiUiR972N3td3AlcD01T1iKq+JSILgJ8DM7xtxgECPNlWsRtTKSzMt0jnf/yHmySwc6cv2axY4XvmTbduvkSTmemedWPJxnQ0QU8uzZQG9AOqTwr9FvA/IvIcrnfTG8hUVf8+lMCYFhBxT+k8+2y49VaXbD791JdoVqzwTX/u3t2XaCZOtGE00zGEZHKpPUtMVe8C7qpVdxL4cRuGZUyLibgeyqBBcMstLtlkZ8Py5S7RLF/u69n07OlLNJMmQb/WPxTQmDYXksnFmI5OxF17GTgQZs1yyeaTT1ySWb7cXbN56SXXNj3dl2gmToRevYIbuzHNYcnFmBAg4h6KNngwfP/7Ltls3+4SzbJl8PrrbvozuIQ0aZIrmZnQtWswIzemfpZcGlFcXMyxY8coLCykvLy86R1CQFJSEjt2nM4tQx1foI+Jx+MhISGB1NRUoqKi/PKZIr410v7jP9xstM2bXaJZvtyt/vyHP7i2w4f7ks2ECZCU5JcQjGkVSy4NKC4uZu/evaSkpNC3b18iIiJoas2zUFBYWEhCQkKwwwgpgTwmqkppaSkFBQXs3buXPn36+C3BVBcW5lZ0HjnSTX0uK3OPiF62zJW5c+HJJ127jAyYPNklm3HjbG00Exz2+MUGHDt2jJSUFDp37kxkZGS7SCym7YkIkZGRdO7cmZSUFI4dO9Ym3zc83D3X5he/gCVLIDfXTQyYMwciIuCRR2DqVEhJcUnmgQfcIwjK2v+yaKadsOTSgMLCQhITE4MdhmlHEhMTKSwsbLphAERFuSGx//5vtx7asWPwr3+5IbXcXLjnHrfac2oqzJgB8+f3YutWd23HmECwYbEGlJeXExEREewwTDsSERERMtfmEhLgsstcAffEzuXL4Z13XFm8eABPPw1paa5nM2WKG0pLTw9u3KbjsOTSCBsKM6cjlM+Xzp3h6qtdAZg3bw0nTlzI0qUu2VTeYzNggC/RTJzoHhVtTEs0K7mIyLnAVGAP8IaqngpoVMaYgEpLKyYzE2680TftuTLRvPQS/P73bsbaeee5RDNlinuAmk0OMM3V5DUXEckE/he4EHgI2CEi54rIJBF5S0Q+FZF3RGSWiFhPyJh2pnLa849+5O6nOXrUPcvmvvsgLg7+93/h0kvd5IDJk+E3v3GPHwiREUAToppzQT9TVSep6lWqehZuNeI/AX8DdgJvAuXAb4F1ItIzYNEaYwIuIgIuvNBNAli50k0I+Pe/3eSAo0fh7rvd4587d4arrnK9nOxsmxxgampOT6PG3yequlREVgOnvGt+ASAinYHZwGIRGauqxf4N1bS1RYsWsWnTJsrKypg4cSITJ04MdkgmCOLi4KtfdQXg8GF3b82SJa784x+uPj3dDZ9NneomCXTpEryYTfA1J7lsEpELVXVNtbrDQEX1Rqp6BPhvEdkEfB9b+r7dKi0tZc6cOVxxxRXcc889qCqTJk0iOTmZUaNGBTs8E2Rdu8K117pSuQDnkiXums38+fDHP7p2o0b5ko1drznzNDkspqqvA18XkT7Vqt8FtjXQfhFgf7O0Y48++iizZ89m3LhxgJsFlZKSwpIlS4IcmQk1lQtwzp7tejBHjsDatXD//ZCYCE884bteM2UKPPQQfPSRW87GdGzNvYlyDnCdiHxHRERVl6vqwkbah+6cTNOo999/nx49etC31oPgN2zYQLrdBGGaULlywC9/6VYMqLyZc/ZsOHQIfv5zGD3a9X6++U147jnYsyfYUZtAaNbsLlUtAx4UkfNxD+jyAOuALFU9VL2tiFwDvOf3SE2bWLBgAb/+9a9r1C1dupSKigpmzJgRpKhMexUfX/Nmzi++cNOdK6/X/P3vrn7gQN8Q2sSJkJwcvJiNf5zW1GFV/QD4QEQigDHADSLSA1DgONAHeEdV/+73SE3AlZSUUFZWhsfjYdmyZbz77rscPnyY1atXs3DhQmJjY4MdomnnuneHb3/bler31yxZAn/5CzzzjFt88/zzXaKZOhXGjnUz2Ez70qL7UlS1FNc7eQ9ARMKAUcAlwBgRGQFsxiWafX6KNWTc8eYdbPxyY7DDqGFkt5E88dUnWvUZ69atY8yYMQDk5+cjIkRFRXHgwAGOHj3qjzCNqVL9sQI/+hGUlLjrNZW9ml//2l27iY93z62pTDZnn22PgW4P/HLTo6pWAB96CwAiMhSYLiJFqvqiP76PCaxVq1Zx/fXXAzBz5kxmzpwJuB7NQw89xJQpU1r0ucuXL2fo0KF0tadamUZERsIll7hy//3u/prly33JZvFi165nT1+imTLFHpYWqgJ2R72qbqOBGWXtXWt7CKHq4MGDdO/evU69qrJr164Wf+5dd93F66+/3prQzBkoJQWuvNIVgM8/9yWaRYt8T+YcOdKXbGzKc+iwJfcNQKOr+e7cuZM+ffo0+H5jCgsLOX78OD169GhpaMYA0K8fzJoFr74KOTnwwQdu6Cw52TflOTXVJZmHH4aNG23KczDZWmAGcFON4+Li6tTn5OSwcuVKnnzS3RNbXFzM//zP/5Cens6BAweYMGECAwcO5NVXX2X//v0MHz6ckpIS9uzZw/Dhw/nHP/5BSkoKDzzwAHfccQfx8fFt/aOZDsjjgTFjXLn7bjh+HLKyfD2bn/3MlS5dfLPQpk6FXr2CHfmZw5KLASArK4u9e/fWqX/mmWcYMmQIN998MwDXXXcd3//+95k6dSrHjx9n7Nix/Od//iezZs1iwIABbN26leTkZMaOHcvs2bPZvHkzU6ZMqdrfmECIi6s75XnpUnj7bbetfKTA2Wf7Ek1mpnvujQkMGxYzAGRnZzN58mTmzZtXVbdw4UJee+013njjDSIjI1m9ejW7d+9m6tSpAOTm5nLgwAG+/vWvs2HDBjIzM0n23qBw8OBBjh49SlZWFuPHjw/Kz2TOXN27w3e+Ay++CAcPwubN8Nhj0Levu3Hz8svdENrFF7vJA2vXQnm5TUHzJ+u5GFSVyMhIbrrpJhYvXszdd98NQJcuXVi7di0x3iuk7733HpmZmVX7LVu2jEsuuYSEhATWrFlTlUR2795NcXExAwYMYNeuXQwePJicnBy62EqGJghEYPhwV37yEzh1yj1SoHII7d573QrQcXHjuPRSX8+mf3+b8twallwM27ZtY8SIEQBMnz6d6dOn19suPT2dgoICwF17mTt3Ls8//zwAa9as4Qc/+AEAv//977n33nvZsmVL1ee+8cYb3HTTTYH+UYxpUnS0W7V50iT3bJojR9wqz3/5y2E++qgHCxa4dn37+hLN5Mmup2OaLySTi4g8AwxR1cxmtj8Ld9PmdFVdEcDQOqSsrCymTZvWZLurrrqKDRs28Kc//YnPPvuMp59+mkGDBgGwfft2duzYwfbt20lNTWX27Nnk5uYSExPDc889x1e+8pVA/xjGtEjnznDNNdC16ydMmNCjapXnt9+GefPg2WddD2b0aF+yuegiiIoKduShLeSSi4hMBW4FsprZPgz38LK6U51Ms8TFxTVrUUqPx8NvfvObOvX79++na9eu3HbbbTXqU1JSmD9/vt/iNCbQKld5rlzpuazMTXmuHEJ7+GHX24mNdTd7ViabYcNsCK22kEouIpIE3AGsPY3d7gB24JaeMS1QeVd+S61Zs4aMjAw/ReAyxdoAABgbSURBVGNM6AgPd72Uiy5y12YKCtxqz5XJ5qc/de3S0mpOebbbukIsuQCP4Jb3b9Yt8CJyNnAurudyawDjMg3YunUrDz/8MImJiXz66acMHDgw2CEZEzCJiW6m2eWXu9f79vkW3nz7bXjpJVd/zjm+RDNhglsf7UwTMslFRGYAe1R1ozSjf+ld9v9h4CZgWIDDMw0YNmwY69atC3YYxgRF795w442uVFS4Kc+VvZo//AGefNL1fi680JdsMjJcXUcnqhrsGBCRTrjex5WqWiYiKwAau6AvInOALar6uohkAsuBiQ1d0BeRWcAsgLS0tNGvvPJKozElJSUxYMCA0/5Zgq28vByPxxPsMEJKWx6T7Oxs8vPz2+R7tUZRUZGtllCLv49JSUkYW7Yk8uGHqaxfn0J2djyqQlxcGaNG5XLeeblkZOTSq9fJkLxeM3HixA9VtcXj3aGSPx8DfuZ9KFmTvEv691HVXzfZ2EtV5wJzATIyMrT6/Rr12bFjBwnt8PbdwsLCdhl3ILXlMYmOjmbUqFFt8r1aY8WKFTT1f+BME4hjcumlvq8rpzwvXRrOkiVdePddd99X796+FZ4nT+44qzwHPbl4n1y5WVV3NLN9BPBr4NsBDcwYY/yocsrzNde4B6Xt2uWGz5YuhX/8A7y3jHHuub7JARdf7GamtUdBTy7AbUCmiDxW+w0RUeBGVX2hWvVFwHQgr55rM8tFZI+q9g1QrMYY02oibgWA/v3h+9+H8nL46CPf9Zrf/tYtVxMZCePGuWQzZYq716a9jHqHQnK5Bag90Pmcd/s9oPZqih/innpZXQbwrPez3vN3gMYYE0i1V3k+cQJWrXK9mqVLYc4cV5KT3coClclmwIDQvb8m6MlFVbNr14lIkfe9jd7XdwJXA9NU9QiwsVb7ZO+X2aq6PbARG2NMYMXGwle+4gq459e8807NYTSAPn18iSbUrtcEPbk0UxrQD7BnzBljzjhdusC117qiCp995ru/pvr1mhEjfMnm4ouDe39NSCaX2lOQVfUu4K5G2q8AQrRzaIwx/iPihsMGDPBdr/nwQ9ezeecdePppePxx3/01lb2a88+HiIi2i9Oe52KMMe2Yx+MSxy9+4Xozubm+pWlOnoT77oPx492qztOnu0dCb9niekCBFJI9F2OMMS0TE+MbGgM4dsyth1Y5OeCf/3T1aWm+yQGTJ0Mz1q49LZZcjDGmA0tNhSuvdAVg717fEFr1R0D37++SzOTJLum0liUXY4w5g/Tp41sPTRW2b6+ZaObO9c/3sWsuxhhzhhKBoUPhhz+E1193Q2hr1sADD7T+s63nYhq0aNEiNm3aRFlZGRMnTmTixInBDskYE0Dh4TB2rCu//GUrP8s/IZmOpLS0lDlz5nDFFVdwzz33oKpMmjSJ5OTkdrEoozEm+GxYzNTx6KOPMnv2bMaNGweAiJCSksKSJUuCHJkxpr2w5GJqeP/99+nRowd9+/atUb9hwwbS/T1X0RjTYVlyMTUsWLCAb3+75tMMli5dSkVFBTNmzAhSVMaY9sauubTAHXfAxo1Nt2tLI0e6O29bo6SkhLKyMjweD8uWLePdd9/l8OHDrF69moULFxLbXh8sYYxpc9ZzMVXWrVvHmDFjAMjPz0dEiIqK4sCBAxw9ejTI0Rlj2hPrubRAa3sIoWrVqlVcf/31AMycOZOZM2cCrkfz0EMPMaVyPQmvLVu2cPvttzN9+nTuvPPOgMU1b948fvrTn7JgwYKq5NcSp06d4oEHHqBnz54cO3aMWbNm0aVLFz9GaoypZD0XU+XgwYN07969Tr2qsmvXrjr1w4cPJzk5mYsuuiigcV199dUAjBw5slWfc/vtt/O1r32N2267jdWrV7Ns2TJ/hGeMqYf1XAwA5eXlDb63c+dO+vTpU6e+oqKixlBaoHz00UcMHTqUiFasF75nzx5WrVrFs88+C8Bzzz1nvRZjAsiSiwHcVOO4uLg69Tk5OaxcuZInn3wSgMLCQn71q1+RkZHBtm3bOOeccygtLeX5559n6dKlvPTSS2zfvp3bbruNF198kcWLF7N//34uuugiCgsLefPNN/nZz37G1q1b+eyzz0hPT+eGG25oNLYVK1bUWR1g7969/N///V+D+3g8Hn71q18RHu5O8aVLl9K1a1fmz59PYWEhu3bt4r777iM7O9svMRpjalHVM66MHj1am7J9+/Ym24SigoKCFu332GOP6XXXXVen/le/+pWOGDFCi4uLVVV1xowZunLlSlVVnTNnjt5///06b948PXXqlA4YMEALCwu1uLhYJ0+erH/+85/1iy++0C5duujhw4dVVfWb3/ymPvzww6qqumnTJp02bVqTsU2bNk3XrFnTop9L1R2TBx54QDt37qyFhYVVcTz//PN+i7FSezlvli9fHuwQQo4dk5qA9dqK37N2zcUAkJ2dzeTJk5k3b15V3cKFC3nttdd44403iIyMZN26dWRnZ3PJJZcAsG3bNiZMmMBll13Gxo0bGTJkCPHx8URGRjJ9+nSuvPJK1qxZwze+8Y2qIagdO3ZU9QK2bt3K4MGDG42rvLycjRs3kpGR0aqfLyEhgf79+xPvfe7rWWedxT//+U+/xGiMqcuGxQyqSmRkJDfddBOLFy/m7rvvBqBLly6sXbuWmJgYALKysqoSC7hrIWPGjKGkpIR//etfTJ8+HYCCggK6du1KfHw8q1at4uKLLwYgLy+P0tLSql/ir7/+OrNmzSI/P5+kpKR6Y/voo48YPnw45eXlVUNc4K6hPPPMMw3+TLWHxYYNG1bjulJYWBjl5eV+idEYU5clF8O2bdsYMWIEANOnT69KErUlJibStWtXAJYvX056ejpbt24lMjKSI0eOcP755wNuNeXKacyrVq3iJz/5CQCrV6+uWq+svLyctWvX8uKLL/L8889z66238tprr3H++efTu3fvqu+5du1aLr74Yl544QVuueUWwsJcZzs9PZ0HH3yw2T/j+PHjOXHiBMePHycuLo6PP/6Yq6666rRjNMY0jw2LGbKyspg8eXKT7a699loOHz7Myy+/TFFREenp6bz33nuMGDGCG2+8kX/961/86U9/qhp+qqioICIigl69egEuiV122WWA6zlkZGTw7LPPVv2Sf+yxx/jrX/9a43uOGTOGffv2kZqaWpVYWiIyMpK5c+fyy1/+kvvvv5+hQ4fyrW9967RjNMY0j7jrNmeWjIwMXb9+faNtduzYwZAhQ9ooIv8pLCwkISHhtPb585//XHXzZLDNnTuXWbNm+fUzW3JMWqq9nDcrVqwgMzMz2GGEFDsmNYnIh6ra4oud1nMxIZNYtmzZwrnnnhvsMIwxfmDJxYSMoqIiLrjggmCHYYzxA0suJmRceOGFwQ7BGOMnllyMMcb4XUgmFxF5RkRWNNFmvIisEZEiEflYRGaLiLRRiMYYYxoRcslFRKYCjd5QICKDgN8CLwI/APKBp4EfBzxAY4wxTQqpmyhFJAm4A1jbRNOvA1NU9ah3v78DW4AfAo8HNEhjjDFNCrWeyyPAHKCkiXYvVSYWAFU9CSwGOgUwNmOMMc0UMslFRGYAe1S1yafTq+qBeqojgCy/B2aMMea0hcSwmIh0Am4Brmzh/mHAV4Bv+TMuY4wxLRMSyQV4DPiZqpa1cP8bgH+r6vsNNRCRWcAsgLS0NFasWNHoByYlJVFYWNjCcIKnvLy8XcYdSG15TE6dOtXkuRUKioqK2kWcbcmOiX8FPbmIyDXAZlXd0cL9+wOXA1c31k5V5wJzwa0t1tQaQjt27Giz9aj8qS3X0Wov2vKYREdHM2rUqDb5Xq1h62jVZcfEv0LhmsttwGMiopUFmABM8L6+oaEdRSQV+DVwo6qWtk24xhhjmhL0ngvuWkt8rbrnvNvvAXvr20lEEoGngDtUNbdafXdV/SIQgRpjjGmeoCcXVc2uXSciRd73Nnpf34kb9pqmqkdEJAVYCPwVGOm9Md8DnAfk4W6wNK20aNEiNm3aRFlZGRMnTmTixInBDskY004EPbk0UxrQD4gRkTjclONhwCW12pUBvdo4tg6ntLSUOXPmcMUVV3DPPfegqkyaNInk5OR2cT3BGBN8IZlcVDWz1uu7gLuqVQ1v04DOMI8++iizZ8+mb9++AIgIKSkpLFmyxJKLMaZZQuGCvgkh77//Pj169KhKLJU2bNhAenp6cIIyxrQ7llxMDQsWLODb3/52jbqlS5dSUVHBjBkzghSVMaa9CclhMRMcJSUllJWV4fF4WLZsGe+++y6HDx9m9erVLFy4kNjY2GCHaIxpJyy5tMQdd8DGJpdAa1sjR8ITT7TqI9atW8eYMWMAyM/PR0SIioriwIEDHD16tIm9jTHGx5KLqbJq1Squv/56AGbOnMnMmTMB16N56KGHmDJlCgDLly9n6NChdO3aNWixGmNCmyWXlmhlDyFUHTx4kO7du9epV1V27dpV9fquu+7i9ddfb8vQjDHtjF3QN4Bb3LEhO3fupE+fPoBbp+v48eP06NGjrUIzxrRDllwM4KYax8XF1anPyclh5cqVXHPNNbz66qvMmjWLlJQUHnjgAYqKioIQqTGmPbBhMQNAVlYWe/fWXcbtmWeeYciQIdx8881ERkayefNmpkyZws033xyEKI0x7YX1XAwA2dnZTJ48mXnz5lXVLVy4kNdee4033niDyMhIwCWh8ePHBytMY0w7YT0Xg6oSGRnJTTfdxOLFi7n77rsB6NKlC2vXriUmJgaA4uJidu3axeDBg8nJyaFLly7BDNsYE8IsuRi2bdvGiBEjAJg+fTrTp0+vt92WLVuq2r3xxhvcdNNNbRajMaZ9sWExQ1ZWFpMnT26yXf/+/YmJieG5555j6tSpbRCZMaa9sp6LIS4urlmLUqakpDB//vw2iMgY095Zz8VU3ZVvjDH+YsnFGGOM31lyMcYY43eWXIwxxvidJRdjjDF+Z8mlEaoa7BBMO2LnizE+llwa4PF4KC0tDXYYph0pLS3F4/EEOwxjQoIllwYkJCRQUFAQ7DBMO1JQUEBCQkKwwzAmJFhyaUBqaiq5ubkcOXKEkpISG/Iw9VJVSkpKOHLkCLm5uaSmpgY7JGNCgt2h34CoqCj69OnDsWPH2L17d6MP0wolp06dIjo6OthhhJRAHxOPx0NCQgJ9+vQhKioqYN/HmPbEkksjoqKi6N69e72P/g1VK1asYNSoUcEOI6TYMTGm7dmwmDHGGL8LyZ6LiDwDDFHVzCba/QIYChQDh4Ffqmr7GL8yxpgOLOR6LiIyFbi1Ge0eBCao6rdV9WagF/BUoOMzxhjTtJBKLiKSBNwBrG2i3dnAncDvqlU/C8wWERtcN8aYIAup5AI8AswBSppo9/8AD/BBtbr1QAXwncCEZowxprlCJrmIyAxgj6pubEbzcYACRysrVPUEUACMCUyExhhjmiskLuiLSCfgFuDKZu7SDcir5+J9IdDVn7EZY4w5fSGRXIDHgJ+patlp7HOinjoPDQypicgsYJb3ZZGI7Dy9ENuNzsCRYAcRYuyY1GXHpC47JjUNbs3OQU8uInINsFlVd5zGbvuAfvXUJwIf17eDqs4F5p5+hO2LiKxX1YxgxxFK7JjUZcekLjsmNYnI+tbsHwrXXG4DHhMRrSzABGCC9/UN9ezzERArIomVFd6v44F32yJoY4wxDQuF5HILMKpW+dBbRgGv17PP87gL+pdUqxuNmy02L5DBGmOMaVrQh8VUNbt2nYgUed/b6H19J3A1ME1Vj6jqZyLyB+A6YLF3t1uBp1R1e9tEHrI6/NBfC9gxqcuOSV12TGpq1fGQUFxKXkRWAFQu/yIijwDXA6NVdZ+3Lhx4EEgABNgN/EZD8QcyxpgzTEgmF2OMMe1bKFxzMS0gIpki8paIXF/Pe91F5GUR+a2IzPeu19bhNXZMvO8/UX3iiIgsrq9dRyAiMSLyqIgcEJGjIrJARM6q1WaI9/x4yrs9L1jxtoXmHBNvu4W1zpPf1fd5HYE4/yUi+0UkT0QWiUjPWm1adJ4E/ZqLOX0ichlwLXAp8HKt9xKALOC/VPUV7w2qO0Rkpqqubvto20Zjx8T7fjfgbOCuatX/apvoguIJ3KSXnwAjcGvxnSci56pqnoj0BlYAX1fVNSIyGHhPRC6o7zpoB9HoMQEQkZG4iUHVz5NX2jrQNnQncAyYAUwCHsL9/7kEoFXniapaaYcFOAf3H+WGWvW/BnLwDnl6654HNgY75mAdE+97jwPDgx1jGx2HVODntepuqX5sgJeA9bXaLAMWBjv+YB2Tasela7DjbaNjIsDYWnUvAoW1jkeLzhMbFmu/TjZQ/y1gnXrPAq8PgHNFZETgwwqqeo+Jt9dyKzBXRH4pIp3bNqygeLLW61e9204iEoNbaumDWm0+AKaLSEqggwuSBo8JgIici5uV+rKI/FhE4toyuLamTu0V6HPwHpfWnieWXNqvOjMxRKQXkI47Qar7wrvt6It6NjQ75TxgCdAduB/YIiKXtllUbUxVj6lq7UQb4d1m4e4Ji6b+88SDO14dSjOOCcBI4C3cAwgfBzaIyOg2CjHoRCQV6AH80FvVqvPEkkvH0s27PVqrvtC7PSMX9VTVf6nq11W1L/A1oAxYJCJDghtZm5oG/FtV12HnSaXqxwRV/bOqzsD9EfJt3HF4S0S6BDHGgBOReBG5GffYkkygcgJQq84TSy4dU+1FPT3ebVPPyenwVPXfwESgHLg9yOG0CRGJAmZT9+c9Y8+TRo5J5XDRS7jkkwrc2MbhtbWTuJ79I7jrMH8Xkf7V3m/ReWLJpWPZ593WHgutXIPtcBvGErLUzXJ5HujfVNsO4kHgXlX9zPvazpO6x6QOdbMrF9PBzxNVLVfVvar6DDAdN4t4Oq08T2wqcsdyGDiA69ZX18u7tUU9fXYCscEOItBE5FZgm7fHVmkH7q/V+s6TU7h1/TqsBo5JQ3YC+QEOKWSo6joR2YI7D1p1nljPpQPxzhB7HhgvIlLtrfOAtar6eXAiC0mjgd8HO4hAEpHrgBhVfa5aXQJuSHA+NRd+BXeeLFT3VNcOqaFj0sjMsCHAX9okuNARAaxU1QJacZ5Ycmm/Kme61O59PoEbC50KICJpwFeBO9outKCp95iIyJ9E5KciEiUiYSLyfWCLqrbqeRWhTES+i7so/bGIfNVbrsXdIFcK/DcwQETO9rYfirtP6O5gxRxoTRwTj4i8JiI3iohHRCJF5F7gRVXdG9TAA0RE4kTkZ94bRyvrfgDMV9XK52K1/DwJ9o08Vlp089OFwLO4qbdZwOW13h+E+4vjYeBvwMXBjjmYxwSXcAtwUygXARODHW+Aj8UNuLvMtZ7y22rtLgBe854nLwHDgh17sI4J7iL1POA4bhHcvwMjgx13gI9JN2AjUAyswq2CfF097Vp0ntjClcYYY/zOhsWMMcb4nSUXY4wxfmfJxRhjjN9ZcjHGGON3llyMMcb4nSUXY4wxfmfJxRhjjN9ZcjHGGON3llyMMcb4nSUXY/xARJJF5FURyRGRpxpoc5+IHBeR+d7y1VZ8v9HVPidHRF5ocfDGBIAtuW+MH6hqnoh8E/fUvjWNNM1R1W/44ft9CHwDQERWtPbzjPE367kY4z/nAcnAsmAHYkywWXIxxn+mAltV9VCwAzEm2Cy5GOM/U4Glp7ODiIwQkSdF5KiIxIvIH0QkV0Q+9l5X6Scib4nISRFZKyK1nwpoTEiy5GKMH4hILHAR8M5p7noSOAtIBW7BPXtmHNAJ91TRa4FZwMXACOAXfgrZmICyC/rG+McluAdOZZ3OTqr6qYhsAqar6v9W1ovIMuACVf2Nt2qPt93Z/grYmECynosx/jEVWKfuueOnq6yeupP11BUDkS34fGPanCUXY/xjKqc/JGZMh2XJxZhWEpE0YDjwjoicJyJjgx2TMcFmycWY1uvj3b4PfA3Yepr7hwOIiKdaXRh1/38K7rqOMSHPkosxrbcBWAT8BvhcVYuau6OITABmel/+UkS6ishVwBSgh4jc7l1aZhYwCjhXRK73c/zG+J2oarBjMOaMICL3ATeoal8/f+4KYLeq3uDPzzWmNaznYowxxu8suRhjjPE7u4nSmLbVRUTme79+TlXfbMmHiMhofHfrDwV2+yE2Y/zGrrkYY4zxOxsWM8YY43eWXIwxxvidJRdjjDF+Z8nFGGOM31lyMcYY43eWXIwxxvidJRdjjDF+Z8nFGGOM3/1/nZ7eB/Tz5qYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "W_min = W_psf(PSF_opt,ls)\n",
    "beta_cpsf = beta_fun(W_min,ls)\n",
    "W_min_dvm = W_psf(psf_display,ls)\n",
    "beta_dpsf = beta_fun(W_min_dvm,ls)\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(ls/1e3,beta_cpsf,'-g',label=r'$\\beta_{opt}$' )\n",
    "ax.plot(ls/1e3,beta_dpsf,'-b',label=r'$\\beta_{dvm,l=6m}$' )\n",
    "ax.plot([ls[0]/1e3,ls[-1]/1e3],np.full(2,bT),'r',label=r'$\\beta_t$')\n",
    "ax.set_xlabel('$l$ [m]')\n",
    "ax.set_ylabel(r'$\\beta$')\n",
    "ax.set_xlim(ls[0]/1e3,ls[-1]/1e3)\n",
    "ax.set_ylim(4.0,5.0)\n",
    "plt.legend()\n",
    "plt.tight_layout()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative (and older) format of semi-probabilistic design equations makes use of the global safety factor ($SF$). The design equation with the global safety factor applied on the mean values reads:\n",
    "\n",
    "\\begin{equation}\\label{eq:SF}\n",
    "\\dfrac{W\\mu_{F_m}}{\\mu_G\\dfrac{l^2}{8}+\\mu_Q\\dfrac{l}{4}}>SF\n",
    "\\tag{7}\n",
    "\\end{equation}\n",
    "\n",
    "where the numerator is the mean value of the moment resistance and the denominator is the mean value of the action induced moment.\n",
    "\n",
    "The minimum cross-section plastic modulus satisfying Eq.$\\eqref{eq:SF}$ is:\n",
    "\n",
    "\\begin{equation}\\label{eq:c6e12}\n",
    "W_{min}=\\dfrac{SF}{\\mu_{F_m}}\\left(\\mu_G\\dfrac{l^2}{8}+\\mu_Q\\dfrac{l}{4}\\right)\n",
    "\\tag{8}\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "### Calibration of the global safety factor for a given span\n",
    "As before, $SF$ can be calibrated so that Eq.$\\eqref{eq:c6e12}$  provides a design with the target level of risk. The limit state function in Eq.(1) is rewritten with $W= W_{min}$:\n",
    "\n",
    "\\begin{equation}\n",
    "    g(f_m,g,q)=\\dfrac{SF}{\\mu_{F_m}}\\left(\\mu_G\\dfrac{l^2}{8}+\\mu_Q\\dfrac{l}{4}\\right)\\cdot f_m-g\\dfrac{l^2}{8}-q\\dfrac{l}{4}\\le 0\n",
    "    \\tag{9}\n",
    "\\end{equation}\n",
    "\n",
    "Consequently the reliability index is computed as:\n",
    "\n",
    "\\begin{equation}\n",
    "    \\beta = \\dfrac{\\dfrac{SF}{\\mu_{F_m}}\\left(\\mu_G\\dfrac{l^2}{8}+\\mu_Q\\dfrac{l}{4}\\right)\\mu_{F_m}-\\mu_G\\dfrac{l^2}{8}-\\mu_Q\\dfrac{l}{4}}{\\sqrt{\\left(\\dfrac{SF}{\\mu_{F_m}}\\left(\\mu_G\\dfrac{l^2}{8}+\\mu_Q\\dfrac{l}{4}\\right)\\sigma_{F_m}\\right)^2+\\left(\\sigma_G\\dfrac{l^2}{8}\\right)^2+\\left(\\sigma_Q\\dfrac{l}{4}\\right)^2}}\n",
    "    \\tag{10}\n",
    "\\end{equation}\n",
    "\n",
    "where it is clear that the reliability of the design is depending on the global safety factor, i.e. $\\beta=\\beta(SF)$.\n",
    "\\par\n",
    "The global safety factor is calibrated by setting $\\beta(SF)\\equiv \\beta_t$. As an example, considering $l=6$ m, the calibrated global safety factor is found to be $SF^*=2.51$, giving $W=3.01\\cdot 10^5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEWCAYAAACqitpwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hU5fn/8fe9S28qVbEhqIAFELArQowGwd4hKmIi2DAWMMYCghorKmBFxa4gWKMSjMIqiiJgw4IJiAWNDfQbQQSB+/fHM/tjXWcLuzPzzM58Xtd1rnHPnDN7+zjy4ZynHHN3REREUqkgdgEiIpJ7FC4iIpJyChcREUk5hYuIiKScwkVERFKuVuwCMqF58+bepk2b2GWkzYoVK2jYsGHsMrKK2iQ5tUtyapfk5s2b9527t6jKuXkRLm3atGHu3Lmxy0iboqIievbsGbuMrKI2SU7tkpzaJTkz+7Sq5+q2mIiIpJzCRUREUk7hIiIiKadwERGRlFO4iIhIyilcREQk5RQuIiKScnkRLj/+GLsCEZH8khfh8tlnsGZN7CpERPJHXoTLzz/DnXfGrkJEJH/kRbg0agSXXgrffx+7EhGR/JAX4bLllrBsGVx+eexKRETyQ16ES4MG8Oc/w7hx8NFHsasREcl9eREuAFdcEULm/PNjVyIikvvyJlxatgz9Ls8+C9Omxa5GRCS35U24AAwZAu3awbnnwi+/xK5GRCR35VW41K0Lo0fDhx/CHXfErkZEJHflVbgAHHoo7L8/jBgRRpCJiEjq5V24mMGNN8IPP8Bll8WuRkQkN+VduADsvDMMHgy33goffBC7GhGR3JOX4QIwahQ0aQJnngnusasREckteRsuzZvD1VdDURE8+GDsakREckvehguEWft77BEmVmrdMRGR1MnrcCkogNtvD6PGLroodjUiIrkjr8MFoHNnOPvsMO9l9uzY1YiI5IYo4WLBpWa2xMx+MLOnzGzzSpz3NzN70MzuNrOrzKwwFfWMHAmtW8Npp+mhYiIiqRDrymUosAw4BLgc6As8Ut4JZnY1sJ+7n+DufwK2AMamopjGjWHMGHj7bbj55lR8oohIfst4uJiZATPd/RZ3f8vdRxOCZZdyzulACKSSf/TfCZxhZmWetyGOPBL69AmLW37xRSo+UUQkf2U8XDx4vdTub4HJ5ZzWDygE3iixby6wDjgxFXWZhee9rFkD55yTik8UEclf0Tv0zawp0Bo4u5zD9gYcWFq8w91/Av4H7JqqWtq2DVcuU6bA1Kmp+lQRkfxjHml6upk1Ao4DLgYaAKe7+xNlHPse0Nrdm5ba/xmw0t3bJzlnEDAIoFWrVt0mTpxYqbp++cX485+788svBdxzzxzq1l23If9aUSxfvpxGjRrFLiOrqE2SU7skp3ZJrlevXvPcvXuVTnb3KBvhNtdWwOnA18AvQLsyjn0PWJJk/xfA/Ip+V7du3XxDzJjhDu7Dhm3QadHMmDEjdglZR22SnNolObVLcsBcr+Kf8dFui7n7Wnf/zN1vAw4GaiVek/kc2CTJ/ibAN6murWdPGDQoPPtFc19ERDZc9D4XAHefA8wHfi7jkDeBBmbWpHhH4p8bAa+ko6brroPNN4dTToFVq9LxG0REcldWhEtCbeClMt6bQOjQ71FiXzfCaLFJ6SimSRMYPz4syX/55en4DSIiuSvGPJeGZvZXM+tSYt+ZwBR3X5D4eaiZzTaz5gDuvgi4A+hf4qMGA2PdPW1PZOndG04+Oaye/Oab6fotIiK5J8aVS2PCvJXZZjbTzMYD37v7pSWOaQVsA9QvsW8I8KWZ3ZE4513gvHQXe8MN0LIlDBwIq1en+7eJiOSGWpn+he7+FdClgmOGAcNK7VtDmKWfUZtsElZOPuywcAUzfHimKxARqXmyqc8lax16KPTvD1dcAfPnx65GRCT7KVwqaezYcBUzcKBWThYRqYjCpZKaNYNbboF58+D662NXIyKS3RQuG+Doo8M2YoRuj4mIlEfhsoFuvTXcHvvjH+HnsqZ8iojkOYXLBmrRAiZMCFcuF18cuxoRkeykcKmCPn3gzDPDHJgXXohdjYhI9lG4VNG110KHDjBgACxbFrsaEZHsonCpogYN4KGH4NtvYfBgiPRYHBGRrKRwqYauXcOillOmwP33x65GRCR7KFyqaehQ2G8/OOss+Pjj2NWIiGQHhUs1FRaGq5bCQjjxRM3eFxEBhUtKbLVVmP8yaxZcdVXsakRE4lO4pEj//mEbORJmzoxdjYhIXAqXFLrtNmjbFo4/PowiExHJVwqXFGrSBB59FJYuDf0v69bFrkhEJA6FS4p16QJjxsC0aeHhYiIi+UjhkgaDBkG/fnDppfDSS7GrERHJPIVLGpjBHXfAttuGkPn669gViYhklsIlTRo3hsmT4fvv4YQTYO3a2BWJiGSOwiWNOnWCcePCysl//3vsakREMkfhkmZ/+lN4sNhll8GMGbGrERHJDIVLmpnB7bfD9tuH+S+ffx67IhGR9FO4ZECjRvD447ByJRx5ZHgVEcllCpcM6dgRHngA5s6F00/X819EJLcpXDLosMNgxAi4777Q0S8ikqsULhk2fDgceiicdx4UFcWuRkQkPRQuGVZQEG6PbbcdHHMMfPpp7IpERFJP4RJBkybw5JOwejUccQT89FPsikREUkvhEkn79vDQQ/D222EtMnXwi0guUbhEdPDBMGpUCJnRo2NXIyKSOrViF5DvLroI3nkHLrggLHR5+OGxKxIRqT5duURWUBCGJu+6a1gmZt682BWJiFSfwiULNGgATz0FzZvDIYdoiRgRqfkULlli003h2WdhxYrQF/Pjj7ErEhGpOoVLFtlpp/AMmPffh+OOgzVrYlckIlI1Cpcsc+CBcMstMHUqnHOOhiiLSM2k0WJZaPBgWLgQrr8+zOT/y19iVyQismEULlnqmmtg0SI491zYemsNURaRmkW3xbJUQQE8+GAYotyvH8ycGbsiEZHKU7hksQYNwgiyrbcOQ5TffTd2RSIilaNwyXLNm8Pzz4enWfbuDYsXx65IRKRiUcLFzOqb2fVm9oWZLTWzJ8ysbQXnbGVmd5rZCDMbbmb3mtm2mao5pq22gmnT4Oef4Q9/gG++iV2RiEj5Yl253AQ0As4Dbgf6ADPMbONkB5tZLWAqcI+7j3T3UcDdwItmVi9DNUe1447wzDOwZAn06aNJliKS3TIeLmbWFFjs7qe5+yR3vxg4C9gKKGtM1E7ADsCiEvvmJM7ZIZ31ZpO99gqTLN9+G448Elatil2RiEhysa5cxpT6eXLitVkZx38HODCyxL79gCXAh6ktLbv17Qt33w0vvAADBsDatbErEhH5rYzPc3H3ZUl21068vlzGOUvM7FbgTDNrQLitNgLo6+4r01Np9howIPS7XHABNG4chiqLiGSTbJlE2ReY6u5zyjnmbGBt4vVoYDd3f6+sg81sEDAIoFWrVhQVFaWu2iyw665w4oltuOuuNvz3v20wK8IsdlXZY/ny5Tn33zwV1C7JqV1Szzzy4lVmVheYCfRz90XlHFcPuBn4ARgC/B/hyqW8QAKge/fuPnfu3BRVnD3cw9XL9dfD+efDddehgEkoKiqiZ8+escvIOmqX5NQuyZnZPHfvXpVzs2Gey9XAiPKCJWEi8IG7DwX2B+oAz5Y1wiwfmMG118IRRyxh9GgYPjx2RSIiQdRwMbPBwPvuPrWC4/YBDgMeAHD3V4BDgebAMemuM5uZwVlnLeTUU+GKK8ImIhJbtD4XM+sP1Hf3m0rsawysc/cVpQ7fJPG6uniHu79sZm8BeX8jqKAAbr89TLK89FKoVw+GDo1dlYjks1gz9E8CTgAWmFnvxHY88Ajwi5kNNbPZZtY8ccoM4DNCR37xZzQB6gFPZ7j8rFRQABMmhIeMDRsG48bFrkhE8lnGr1zM7GRgAuGK46BSb9/s7qvNrBWwDVAfwN2Xm9kBwDVmthPwCbAFcIy7f5Wp2rNdrVrwwANhcuXZZ4d9Q4bErUlE8lOMeS73AvdWcMwwYFipff8GjkhbYTmidm2YNAmOPz4EzKpVukUmIpmXDaPFJMXq1AkBU3yL7MorY1ckIvkmWyZRSorVrh0eNlanDlxySbiCGTlS82BEJDMULjmsVi24554QNJdfHgLm6qsVMCKSfgqXHFdYCHfeCXXrhgmXq1bBjTcqYEQkvRQueaCgAG65JQTMTTeFgLnllrBfRCQdFC55wgxuuCEEzDXXwPffw/33hz4ZEZFUU7jkEbPQ59KsWVjwctkyePxxaNQodmUikmt0YyQPDRsWZvNPnw777w/ffRe7IhHJNeWGi5kdmtj2zlRBkhkDB4arlnfegX33hc8/j12RiOSSiq5cngS6AW9loBbJsEMPheefhy+/hL33hgULYlckIrmionCZ5e4j3P0nADO7xsw+NrMHzezwDNQnadajB7z0EqxeDfvsA2+8EbsiEckFFYXLrx7g5e5/BVYAJ7n7k2mrSjKqSxd49VXYaCPo2ROeeip2RSJS01UULuuS7Fvo7r/Zb2btU1OSxNCuHcyaBTvvDEccAWPGxK5IRGqyisLFK7kPSjxrRWqmVq1gxgw4/HA455ywqvLatbGrEpGaqKJ5Ln3M7OVS+zom2VcP6Axo/d0arkEDmDw5zIO54Qb45BN4+GHNhRGRDVNRuLQkXKn8XGLfcmDLUsc1qMRnSQ1RWAijR4dbZUOGwH77wTPPwGabxa5MRGqKigLhWne/sDIfZGYjU1CPZJEzzoCttw7Phdl9d3j22dAnIyJSkYr6XB4u700za2FmDRI/TklNSZJN+vaFmTND38uee8KTGiMoIpVQUbi0KTFL/+DinWa2mZkVAV8By8zsKnefn85CJZ5ddoE5c2DHHcNIslGjYF2ycYQiIgmVmaH/MNAeeB/AzAx4AugBLABuA44xs4FprFMia906TLY88UQYMQKOPRZWrIhdlYhkq8osXHmsu1/n7osTPw8EdgPmAF3d/VxgL2BAmmqULFGvHtx3X+jsf+IJ2GuvMJpMRKS0isLlQ3d/rvgHM6sDjCBMrjzD3VcBuPs3wA9pq1Kyhhmcdx489xx8+insumu4ohERKamicPlPqZ/PJgxDfsTd55V6TzMh8sgf/hDWIWvWDH7/exg3Drys6bUikncqCpcGZrYRgJltCwwnrC32t5IHmVkzYM+0VChZa/vtYfZsOOigMJu/f39Yvjx2VSKSDSoKl9uB18zsOqCIcHXyV3dfUnyAmTUE7ifM0pc8s9FGYXjyVVfBo4/CbrvBhx/GrkpEYis3XNz9ceAyYA/ga2Cwu99a/L6ZnQlMApoAs9JXpmSzggK48EJ44QVYujT0w0ycGLsqEYmpwiVb3P1R4NEy3rsFuCXVRUnN1KsXvPlmGKbcr19Yxn/0aKhTJ3ZlIpJplRmKLFJpm28ORUVw7rlw881hXbJPP41dlYhkmsJFUq527bCi8qOPwvvvh4eRPfZY7KpEJJMULpI2xxwDb70F220HRx8Np50GK1fGrkpEMkHhImnVrh288goMGwZ33BE6+997L3ZVIpJuChdJuzp14NprYdo0+PbbEDC3365JlyK5TOEiGXPggfDuu9CjB5x+erhVtnRp7KpEJB0ULpJRrVrB1KnhSuYf/4CddgrrlIlIblG4SMYVFIQ+mDfegObNwwPJTjtNS8eI5BKFi0TTpQvMnRuCZvz48PMsrfMgkhMULhJV3brhFtlLL4VHKe+7L1x0EaxeHbsyEakOhYtkhX33DZ39AweGRTB32y0sJSMiNZPCRbJG48Zw113w9NPwzTchYC66CH7+OXZlIrKhFC6SdQ45JCwbM2BAuIrZZRf1xYjUNAoXyUqbbAJ33x0mXq5cCfvsA3/5C6xYEbsyEakMhYtktQMPhPnz4YwzYOxY2Hnn8NwYEcluUcLFzOqb2fVm9oWZLTWzJ8ysbSXPLTCz481sjJkNMbM90l2vxNW4cVi+/+WXoVYtOOAAOOEE+Prr2JWJSFliXbncRHhk8nmERyn3AWaY2cblnWRmzYDngDbAue4+zt1fT3OtkiWKR5QNHw6TJ0OHDmF+zLp1sSsTkdIyHi5m1hRY7O6nufskd78YOAvYCji8nPMaA9OA2e5+tbvrj5Q8VK8ejBwJ77wTJl0OHhz6Y959N3ZlIlJSrCuXMaV+npx4bVbOOWOBde4+Ij0lSU3SoQNMnw733w8LF0LXrmGmvzr8RbJDxsPF3Ze5e+lHRtVOvL6c7Bwz2wEYABSZ2Wgze8nMnjGzPdNZq2Q3MzjxRFiwAE45Ba6/PoTOxIlazl8ktmwZLdYXmOruc8p4/1jAgA+AC4CDCIH0spl1zUyJkq2aNg19L7NmQcuW0K8fnHNOF955J3ZlIvnLPPJf8cysLjAT6Ofui8o4ZjJwkLs3KrFvZ+Bd4HF3PyrJOYOAQQCtWrXqNnHixHSUnxWWL19Oo0aNKj4wD6xdC1Onbsadd7Zh+fI6HHLIlwwcuJiNNloTu7SsoO9KcmqX5Hr16jXP3btX5dxsCJcbgefdfWo5x/wT6Ozum5Xa/ynwk7t3LO93dO/e3efOnZuSerNRUVERPXv2jF1GVvnHP17hX//ah1tvhY02giuugEGDoLAwdmVx6buSnNolOTOrcrhEvS1mZoOB98sLloTPgE3MrPQfDf8FvktLcVKjNW68hrFj4a23oFOnMAmzS5cw419E0i9auJhZf6C+u99VYl9jM2uY5PDHgLpA6QmTzYCKgkny2M47h1FlU6aEZWR69w7b/PmxKxPJbbFm6J8EnAAsMLPeie144BHgFzMbamazzaw5gLtPA54ALizxGXsTOvlLD2sW+RUzOOqosBjmDTeEJ2B26QKnngpffRW7OpHcFGMS5cnAvYQRX1NLbI8QJleuBloB2wD1S5z6R2Chmd1lZsOBk4Ge7q6ZDVIpdevCueeGeTFnnw333QfbbguXX65HLIukWox5Lve6e4G7W5JtSOKYYe7e0t0/L3HeSnc/193/7O6j3P1Ud1+S6fql5mvaFG68ET74INwiGz4c2rUL65fpCZgiqZEt81xEMm7bbUNfzKxZYfLlkCHQvj088EAY0iwiVadwkby3555QVARTp4bnyJx0UuiTefppzfQXqSqFiwih0793b5g7Nywfs2oVHHYY7LUX/OtfChmRDaVwESmhoACOOy6MLLvjDliyJDywbN994cUXFTIilaVwEUmidu0wo3/hwtDR/8kn8Pvfw377wYwZsasTyX4KF5Fy1K0LZ54ZQmbcOFi0CH73O+jZM/TT6EpGJDmFi0gl1KsHZ50VwmXsWPj3v6FXr/CgsueeU8iIlKZwEdkA9eqFIcuLFsEtt4Q+mb59w8PKJk/WEGaRYgoXkSqoXz8shrlwIdx7b1i37NhjYccdw8+//BK7QpG4FC4i1VC7NgwYEEaXTZ4cQmfgQGjbFkaPhv/9L3aFInEoXERSoLAQjj4a3nwz9MFsuy0MHQpbbgkXXBBun4nkE4WLSAqZwUEHheHKc+ZAnz7hCmabbcIVzrvvxq5QJDMULiJp0r07PPJI6Pw/4wx47DHo3BkOOACeeQbWrYtdoUj6KFxE0qxNGxgzBj77DP7+d/jwQzjkkLBI5tix6peR3KRwEcmQpk3hb3+DxYvD+mUtWsBf/gJbbBFeFy6MXaFI6ihcRDKsdu2wftmsWeGpmIceCrfdBttvH/ponnlG82Wk5lO4iES0667w4IPw6adw6aXw1lvhllm7dnDVVfDNN7ErFKkahYtIFthsMxg5MvTLTJ4c5slcdFG4Zda/P8ycqSVmpGZRuIhkkdq1w3yZ6dPDY5hPPz3Mm+nRI8z+v+EG+O672FWKVEzhIpKlOnYMo8y++ALuugs22gjOPx823xyOPx5eeEHDmSV7KVxEslzDhvCnP8Frr4VJmKefDs8/H+bLbLcdXHklfP557CpFfk3hIlKD7Lwz3HQTfPklPPwwbL01XHJJeD3wwLDvp59iVymicBGpkerVg379Qt/MokUwfDj85z/wxz+GwQGDBsGrr2oQgMSjcBGp4dq2hcsuCyEzYwYccQQ89FB4kNn224f3/vOf2FVKvlG4iOSIgoLw+OV774WvvoIJE2CrrWDUqBAyu+8Ojz++uebOSEYoXERyUOPG4bkyL74Y5s5cdx2sWgXjxm1H69ZhJYAHH9S6ZpI+CheRHLfFFuHZMm+/DRMmzGHYsPBwsxNPhJYt4cgj4dFHYcWK2JVKLlG4iOSRbbZZwVVXhcUzX301dPy/9lpY66xlyzB/5skn4eefY1cqNZ3CRSQPFRTAXnuFJf+XLAkDAU48MdxGO+KIsGJzv34wZYquaKRqFC4iea6wMAwEuP12+O9/Ydq0ECwvvgjHHBOC5qijwoPP1EcjlaVwEZH/r1atMBlz/PgwUXP6dDjllHDrrH//EDR9+8Kdd4YRaSJlUbiISFK1akGvXnDzzeHW2SuvwJlnhidpDhoErVuHW2vXXAMffRS7Wsk2ChcRqVBBAey9d1iVedGisMbZyJGwejVceCF06BC2YcPgpZdgzZrYFUtsChcR2SBmYY2zSy+FuXPDPJqbbw4TNseMCf03LVqE22gPPwzLlsWuWGJQuIhItWy5Zbhd9vzzsHQpPPZYGHH24othrbMWLWDffeHvfw9P2tR6Z/lB4SIiKdO4cZiUOWFCGHk2ezZcfHEYznzxxdC1a+irOflkmDRJVzW5rFbsAkQkNxUUwG67hW3UqDC6bNo0mDoVnn4a7rsvHLP77uHZNAceGI6tXTt25ZIKunIRkYzYdFMYMAAmToRvv4VZs8LVzLp1cMUVYRXn5s3h8MPh1lvDSs66hVZz6cpFRDKusBD23DNso0bB99+HOTXPPx+2p54Kx229Nfzud7D//uF1s83i1i2Vp3ARkeg22SSsAnDUUeFqZdGiEDIvvhjWOrvnnnBcx47rg2a//aBp07h1S9kULiKSVcxg223DdsYZsHYtvPNOCJrp08NggZtvDsd16hSGPu+3H/ToAc2axa5eiilcRCSrFRaGUWZdu4ZJmqtXwxtvQFFR2MaPD/NrIIRNcdDsuy+0ahWz8vymcBGRGqVOndD5v88+cMklIWzmzFkfNnffDePGhWO32y6ETPHWtm244pH0ixIuZlYfuBzoB9QDXgbOd/ePK3l+W+Bd4GB3L0pXnSKS/erUCUvT7L13GH22ejW8+SbMnBm2J54It9IgDAjYe++wJtree0OXLuF8Sb1YVy43AQ6cB3QChgJdzayzu/9Q3olmVgDcAzRMe5UiUuPUqQN77BG2YcPCUOcPPggLb86cGYZAT5kSjq1XL8yt2WKLbfjxx3BOixZx688VGQ8XM2sKLHb3qxO7JpnZJ8B44HDg3go+4hzgQ6BHumoUkdxRUAA77RS2004L+778MoTMq6+G10mTtuThh8N77dqtD6fdd4fOnXV1UxWxrlzGlPp5MiFcyh3rYWYdgM6EK5fB6SlNRHJd69Zw9NFhA/jnP1+hQYMezJ4Nr78eRqU99FB4r25d6NYNdt11/YoD7dqp76YiGQ8Xd0+2mlDxgg8vl3WemRUC1wKnADuloTQRyVP16q2jR48wygzCXJslS0LQvPZaGJ1WclTaJpuEsCneuncPgaXAWc88C9ZXMLOTgWPdvU85x1wMzHf3p82sJzAD6FVWh76ZDQIGAbRq1arbxIkTU1121li+fDmNGjWKXUZWUZskp3ZJrjLtsnatsXhxAxYsaMKCBY1ZsKAJixc3ZN26kCibbLKa9u1/ZPvtf6R9+7A1a7Y6E+WnTa9evea5e/eqnBt9KLKZ1QXOIIwcK+uYTsBW7n5lZT/X3ccTbrXRvXt379mzZzUrzV5FRUXk8r9fVahNklO7JFfVdlmxIkzwnDcP5s6tw7x5zXjwwWasWxfe33RT2GWXMEdnl13Cts02+XGFEz1cgKuBEe6+KNmbZlYbuBI4IaNViYhUoGHDMKx5r73W71uxAt5+OwTOm2+GZ9g8/3xYaQBg443DEOguXcJggc6dYYcdQt9OLokaLmY2GHjf3aeWc9hewMHAD/bbuJ9hZp+6e5s0lSgiskEaNlw/76bYzz/De++tD5s334Q77oCVK8P7tWqFddM6dw6h06lTeNpnq1Y19yonWriYWX+gvrvfVGJfY2Cdu68oceg8YJdSp3cH7gROBWalu1YRkeqoVy90+ncv0Xuxdi0sXBiuct55J2zTp8ODD64/pkWLEDLFYdOpUwihhjVgll+sGfonAccDY82sd2L3xoRbX0ea2VDgGKCvu38HvF3q/I0T/7jQ3T/IUNkiIilTWAjt24ftuOPW7//uO5g/H959d/3r+PHw00/hfbPQb7PTTrDjjuvn8LRvn1231mJMojwZmAAYcFCpt29299Vm1grYBqif4fJERKJq3hx69QpbsbVr4eOPQ9i8/364xfb++/Dcc7BmTTimoCCsJL3DDmHr2DG8dugADRpk/t8jxjyXe6lgFr67DwOGlfN+ESGcRERyXmFhWIRzu+3gyCPX71+9OjyxszhsPvwwLHXzzDPrQ8csPHStY8cQNMVbx44hyNLVp5MNo8VERKQK6tQJt8Z23PHX+3/5JfTnfPBBCJ0FC0LwFBWtH0QA4WFrxbfmirfttw9XQNW9xaZwERHJMbVrhyuTjh3D0z2LrVsHn38egqY4cD76CKZNg3vvXX9cQQG0aVO9GhQuIiJ5oqAg3CLbemvo3fvX7/3vf/Dvf4fto4/C9nGlHoKSnMJFRERo0uS3w6UnTar65xVUvyQREZFfU7iIiEjKKVxERCTlFC4iIpJyChcREUk5hYuIiKScwkVERFJO4SIiIiln7h67hrQzs2+BT2PXkUbNge9iF5Fl1CbJqV2SU7sk197dG1flxLyYoe/uLWLXkE5mNtfdu1d8ZP5QmySndklO7ZKcmc2t6rm6LSYiIimncBERkZRTuOSG8bELyEJqk+TULsmpXZKrcrvkRYe+iIhklq5cREQk5RQuIiKScgqXGsTMeprZNDMbkOS9zczsETMbZ2ZTzOyAGDVmWnltknj/JjPzEtszma4x08ysvpldb2ZfmNlSM3vCzBOrxXYAAAaJSURBVNqWOqZj4nsyNvHaNVa9mVCZNkkc92Sp78vNMerNFAsuNbMlZvaDmT1lZpuXOqZK35W8mOeSC8zsIOB44EDgkVLvNQZeBi5194lm1gz40MyOcPdXM19tZpTXJon3NwU6AMNK7H4uM9VFdRPgwHlAJ2Ao0NXMOrv7D2a2JVAEHO7ur5lZe2CWme3u7gujVZ1e5bYJgJl1Adbx6+/LxEwXmmFDgWXAIcDvgGsI/y/1AKjWd8XdtdWQDdiB8D/IyaX2Xwl8S2KARmLfBODt2DXHapPEezcAO8euMcPt0RS4sNS+U0u2EfAQMLfUMdOBJ2PXH6tNSrRLy9j1ZrBdDNij1L4HgB9LtUmVviu6LVazrCxj/x+BOZ74L5/wBtDZzDqlv6yokrZJ4qplMDDezC4xs+aZLSuqMaV+npx4bWZm9YEjCd+Pkt4ADjazTdJdXCRltgmAmXUGjgEeMbNzzaxhJouLwYPXS+3+lkTbVPe7onCpWX4zbtzMtgC2JnwpSvpv4nXXdBcVWVlj6bsC/wI2Ay4H5pvZgRmrKhJ3X+bupQO3duL1ZaAbUI/k35dCQrvllEq0CUAXYBqwI+GK9y0z65ahErOCmTUFWgNnJ3ZV67uicKn5Nk28Li21/8fEa8sM1pI13P05dz/c3dsAfYA1wFNm1jFuZVH0Baa6+xz0fSlWsk1w9/vc/RDCX0ZOILTDNDPL6XUJAcyskZn9CZgL9ASKBwNV67uicMkdP5X6uTDxujrThWQbd58K9ALWAkMil5NRZlYXOIPf/nvn7felnDYpvlX0ECF8mgIDM1xeDCsJV/nXEfphHjWzdiXer9J3ReFS832eeC19/7NJ4vWbDNaStTyMbJkAtKvo2BxzNTDC3Rclftb35bdt8hseRlk+Qx58X9x9rbt/5u63AQcTRhEfTDW/KxqKXPN9A3xBuJwvaYvE6yuZLSerfQQ0iF1EppjZYOD9xJVbsQ8Jf1NN9n35GZiXofKiKKNNyvIR8H9pLimruPscM5tP+C5U67uiK5caLjFCbAKwj5lZibe6Aq+7++I4lWWlbsDtsYvIBDPrD9R397tK7GtMuDU4hcQ8hhK6EoaXlr4FkjPKapNyRoZ1BO7PSHHZpTbwkrv/j2p8VxQuNUvxCJfSV5w3Ee5/HgBgZq2A3sA5mSstmqRtYmb3mNn5ZlbXzArM7DRgvrtX+eFHNYWZnUTolF5gZr0T2/GEyXG/AKOAbc2sQ+L4HQnzhS6KVXO6VdAmhWb2mJkNNLNCM6tjZiOAB9z9s6iFp5GZNTSzvyYmjxbvOxOY4u4LEruq/l2JPZFHW6UnPO0J3EkYevsycGip97cn/C3jWuBhYN/YNcdsE0Lg/o8wbPIpoFfsejPUJicTZpl7km1cieN2Bx5LfF8eAnaKXXusNiF0UE8CVgCfAI8CXWLXnYF22RR4G1gFzCQsr98/yXFV+q5oyX0REUk53RYTEZGUU7iIiEjKKVxERCTlFC4iIpJyChcREUk5hYuIiKScwkVERFJO4SIiIimncBERkZRTuIhUg5ltbGaTzexbMxtbxjGXmdkKM5uS2HpX4/d1K/E535rZvVUuXiSNtOS+SDW4+w9mdhzhaX2vlXPot+5+dAp+3zzgaAAzK6ru54mki65cRKqvK7AxMD12ISLZQuEiUn0HAO+5+9exCxHJFgoXkeo7AHhhQ04ws05mNsbMlppZIzO7w8y+N7MFiX6VbcxsmpmtNLPXzaz00wBFsprCRaQazKwBsBfw4gaeuhJoCzQFTiU8f2ZvoBnhyaLHA4OAfYFOwN9SVLJIRqhDX6R6ehAeNvXyhpzk7v8xs3eAg939xuL9ZjYd2N3dr0rs+jRxXIdUFSySCbpyEameA4A5Hp43vqHWJNm3Msm+VUCdKny+SDQKF5HqOYANvyUmkvMULiJVZGatgJ2BF82sq5ntEbsmkWyhcBGpuq0Sr7OBPsB7G3h+LQAzKyyxr4Df/n9phH4dkRpD4SJSdW8BTwFXAYvdfXllTzSz/YAjEj9eYmYtzewo4PdAazMbklhaZhCwC9DZzAakuH6RtDF3j12DSE4zs8uAk929TYo/twj4xN1PTuXniqSCrlxERCTlFC4iIpJymkQpkhktzGxK4p/vcvd/VuVDzKwb62fr7wh8koLaRFJOfS4iIpJyui0mIiIpp3AREZGUU7iIiEjKKVxERCTlFC4iIpJyChcREUk5hYuIiKTc/wNa+HPp+jnHFAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "SF = W_bT*mF/(mG*ls*ls/8+mQ*ls/4)\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(ls/1e3,SF,'-b')\n",
    "ax.set_xlabel('$l$ [m]')\n",
    "ax.set_ylabel('SF')\n",
    "ax.set_xlim(ls[0]/1e3,ls[-1]/1e3)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calibration of the global safety factor for different beam spans\n",
    "\n",
    "As before, the design performed with the calibrated $SF$ for $l=6$ m results in a different reliability for different spans. The $SF$ can be calibrated so to give a reliability as homogeneous as possible over the different design scenarios. Therefore, as previously, the $SF$ can be calibrated by solving a minimization problem. \n",
    "\n",
    "\\begin{equation}\\label{eq:c6e15}\n",
    "\\min_{SF^*}\\left\\{ {\\sum}_{i=1}^n\\left(\\beta_t-\\beta_i(SF^*,l_i)\\right)^2\\right\\}\n",
    "\\tag{11}\n",
    "\\end{equation}\n",
    "\n",
    "The optimized $SF$ obtained is $SF^*=2.58$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization terminated successfully.\n",
      "         Current function value: 6.471390\n",
      "         Iterations: 17\n",
      "         Function evaluations: 34\n",
      "The optimal global safety factor for l = [6,30] m is:\n",
      " SF = 2.58\n",
      "\n"
     ]
    }
   ],
   "source": [
    "W_sf = lambda sf,l:  sf/mF * (mG*l**2/8 + mQ*l/4)  # W as a function of SF and l\n",
    "\n",
    "# minimum W using SF for l = 6m\n",
    "W_min = W_sf(SF[0],ls)                             \n",
    "beta = beta_fun(W_min,ls)\n",
    "\n",
    "\n",
    "def betadiffSF(sf):\n",
    "    \"\"\"compute sum of square difference over lengths for given SF\"\"\"\n",
    "    beta = []\n",
    "    for l in ls:\n",
    "        W_min = W_sf(sf,l)\n",
    "        beta.append(beta_fun(W_min,l))\n",
    "        \n",
    "    return np.sum((np.array(beta) - bT)**2)\n",
    "\n",
    "res = sp.optimize.minimize(betadiffSF,3,method='nelder-mead',options={'xatol': 1e-5, 'disp': True})\n",
    "SF_opt = res.x[0]    # Optimal SF\n",
    "print(\"The optimal global safety factor for l = [6,30] m is:\")\n",
    "print(\" SF = {sf:.2f}\\n\".format(sf=res.x[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAEcCAYAAADdgl6CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdeXhTVf4/8PdJ0zZpmu6lOy0tq+zboOIoAyKOK+A4AqOojAgqWnAEFR0Vvwwqg4wo6E9RHAe3ER0UUHABKyAutWxS1rKUQqFr0rTNnnt+f9y2llJooUua9v16njzJvTm5+STDmHfPPfccIaUEERERkS/SeLsAIiIioovFIENEREQ+i0GGiIiIfBaDDBEREfksBhkiIiLyWQwyRERE5LO8GmSEEJ8KIWSt29IG2g+o094hhIhvrXqJiIiobdF6642FEAMAKABm19r9YQMvewTAHADVk98cl1Lmt0B5RERE5AO8FmSgBpjpUsrCxjSuCj75Usp/tmxZRERE5Cu8EmSEEP0B3AogVgixDsAbUsrKBl72NIChQogQAP+RUm5r6TqJiIiobfPWGJkBAL4E0BvAYgA7hBCDz9W4Krw4AJwEcA+ArUKIRUIIv9YoloiIiNom4c21loQQAsAkAMsAuAH0klIWNfCaJACvALgZwLNSyqfP0e5eAPcCgF6vH5yUlNScpbcpiqJAo+EFaHa7H8xmf5SX+0NKQK/3IDTUBaPRBSG8XV3bwH8r9eP3Ur+DBw8WSymjvV0H0fl4NcjUFCHEcABbADwmpVzYiPZ+ANYCGA4gUkrpPl/7IUOGyF9++aVZam2LMjIyMGLECG+X4RUuF/Dxx8CSJcBPPwHBwcAddwBDhmRiypSh3i6vzenI/1bOh99L/YQQWVLKId6ug+h82sSfIFLK7wGsA5DWyPYeAE8CCAHAvxY6oKIi4B//AFJSgEmTgNJS4JVXgPx84NVXgdTUhoZcERFRe+DNq5bqOgCg7ALb2wCUtEw51Bb9+qva+/Luu4DDAVxzDbB8OXDttQDPDBARdTxtKcj0AnD/BbQfBvVqJ2cL1UNthJTAxo3AokXAl18CQUHA3XcDDz0E9Orl7eqIiMibWv1vWCGETgjxiRDibiGEnxAiQAjxNICVUsrjVW0eEUL8JISIqtqeIIT4txCiS9V2fwBTADze2vVT63G5gPfeAwYNAkaPBnbtAhYsAPLygNdeY4ghIiLv9Mi4oF6htBTq3DA/A1ggpdxZq00MgC4A9FXbJQCuApAthNgO4HMAd0spXa1WNbWa8nLgzTeBl14Cjh9XA8tbbwF/+QsQGOjt6oiIqC1p9SBTNVD3tgbazEatpQuklF9DDTbUjpWUqAN2X34ZMJmAK68Eli0DrruO41+IiKh+bWmMDHVQp04Bixerp4sqK4GbbwYefxwYNszblRERUVvHIENec+wY8M9/qqeNXC5gwgTgsceAvn29XRkREfkKBhlqdUeOqHPA/Oc/gBDAnXcCjz4KdO3q7cqIiMjXMMhQqzl2DJg/H3jnHcDPD7j/fmD2bCAx0duVERGRr2KQoRZ3/LjaA7NihTpo97771FNI8fHeroyIiHwdgwy1mBMn1ADz1lvqKaR771UH8bIHhoiImguDDDW70lLg+efVS6k9HuCvfwXmzgXa8QLkRETkJQwy1GwqK9V1kBYuBCwWdRXqefPUhR2JiIhaAoMMNZnTqc7E++yzQEEBcNNN6imlPn28XRkREbV3DDJ00aQEVq1Sx70cOaLOxLt6NXDZZd6ujIiIOgpO/E4X5YcfgOHDgdtuAwwG4IsvgIwMhhgiImpdDDJ0QY4eVcPL5Zerj996C9ixA/jjH9Urk4iIiFoTTy1Ro5jNwIIF6mBePz/gqafUyeyCg71dGRERdWQMMnReHo86kd3cuerq1Hfeqc7Om5Dg7cqIiIh4aonO44cf1BWo770X6NULyMoC3n6bIYaIiNoOBhk6y6lTas/L5ZcDp08D778PfPcdMHCgtysjIiI6E4MM1XA6gUWLgB49gA8/VC+r3r8fmDiRA3mJiKht4hgZAqBeOn3ffWpwueEG4F//Arp29XZVRERE58cemQ6uqAi46y7gD38AHA5g3Tpg7VqGGCIi8g0MMh2UoqhzwPTsCbz3nnoaac8e4PrrvV0ZERFR4/HUUgeUnQ1Mnw5s3Qr8/vfAa68BvXt7uyoiIqILxx6ZDsRuB558EhgwANi7V+2RychgiCEiIt/FHpkO4ocfgClT1MG8kyerVydFR3u7KiIioqZhj0w7V1kJzJypLvBotQIbNgDvvMMQQ0RE7QN7ZNqxjRuBqVPVxR3vvx94/nnAaPR2VURERM2HPTLtUFmZuqzA1VerCzx+9x2wbBlDDBERtT8MMu3Mpk1A377qQN7Zs4Hdu4Err/R2VURERC2DQaadsNnUsTCjRgE6HbBtG7BwIaDXe7syIiKilsMg0w4cOGDEoEHAkiXAjBnAzp3qqtVERETtHQf7+jCXC1iwAHj22UGIiwO++goYPdrbVREREbUeBhkfdfAgcPvtQGYmMHp0Af7731iEh3u7KiIiotbFIONjpFTngZkxAwgMBFatAqKi9iM8PNbbpREREbU6BhkfUlYG3Hcf8MEHwIgRwLvvAgkJ6jIDRERNIRUJ52kn7Ll22I/ZYc+1e7skokZhkPERP/4ITJoEHD8OzJ8PPPaYOkcMEVFj1ASVo2pQsR21qYHlmB2OXAfsx+2QTuntMokuGINMG6cowAsvAH//O5CYCGzZAlx2mberIqK2yG1xw3bEBvsRe829/WjV41w7pOPMoOIf4w9dig7Bg4MRNT4KumQddCk66JJ1CEwOBDiJJvkABpk27PRpdUDvxo3AbbcB/+//AWFh3q6KiLxFSgnnKSdsh22wHbbBftiuPs6xwXbEBneJ+4z22nAtdF10CO4XjKiboqDrolNvKerNT89uXfJ9DDJt1HffARMmqONi3noLuPtuQAhvV0VELU1KCWe+E7YcG6yHrLAdqgoqh9TwoliV3xprAF2yDvo0PaJviYY+TQ9dqg76VD10XXTwD/dvlZq3b98+RqvVPi2ljAXnJ6PmowghzFLKXW63e8HgwYMP1deIQaaNURTgn/8E5s4FunUDvv4a6NPH21URUXNzmVywHbTBetAK6wFrzWPboTPDivAX0KXqENQtCOGjwqFL00HfVa+GlmQdNAHezQ3bt28fExgYuDQlJcWp1+tNGo2GA22oWUgp4XK5tBaL5YrTp09/npWV9dDgwYM31G3HINOGmEzA5MnAunXAn/8MvPkmF3ok8mWKW4H9mB3W/VZY91thO2Creewqdv3W0A/Qd9FD312P8D+EQ99Nr4aVbnroOusg/Npud6xWq306JSXFaTAYbN6uhdoXIQQCAgLcUVFRZTqdLujIkSNzADDItFW//ALceitw8iTwyivAAw/wVBKRr/DYPLAesMK677db5b5K2A7ZzrgSyD/aH0E9gxA1Ngr6HnoEdQ9CUI8g6Lp4v2flYkkpY/V6vcnbdVD7ZjAYrFLKLvU9xyDjZVICr78OpKcDMTHqVUlcJ4mobfLYPGpIya5EZXYlrNlWVO6thP2oHajOKxpAn6ZHUK8gRN4QiaAeQQjqqQYW/4jWGbPSyjQ8nUQtTah/2df75z2DjBc5HMD99wMrVgDXXqtOcBcZ6e2qiEhxKbAdtKHi1wpU/lqJyj1qcLEf+S2wCH+BoB5BMA4xInZyLIJ6Bam37kHQBPpm7wqRL/JqkBFCfArg5lq7lkkpZ5yn/eMAegNwACgE8KSU0tOyVbaM/HzgllvUie6efBKYNw/Q8L99RK1KSgnHSQfwI5D7Y64aWn6thHW/FdKlJhahFdB318M4yIjYO2IR1DsIht4G6LvqofHn/2mJvM1rQUYIMQCAAmB2rd0fnqf98wAGSCmvrdpeCeBlAA+0ZJ0t4ccfgfHjAYsF+PhjNdAQUctSHAoqsytRsasClbvV+4pdFXCXqnOvHMVRBHYOhKGvARHXRSC4bzAMfQ0I6sEeFqK2zJs9MrMBTJdSFjbUUAjRE8AjAMbW2r0cwHdCiDellDtaqMZm99Zb6umkxETgyy+Bvn29XRFR++Myu1C5qxLlO8pRsaMCFTsrYN1rhXSrvSwavQaGvgZE3xKN4P7BOKQcwvA7hsM/rF2OYSFq17wSZIQQ/QHcCiBWCLEOwBtSysrzvGQiAD8AP9fa9wvUHp07ALT5IONyAbNmAcuWAaNHAx9+CEREeLsqIt/nLHaiYnsFyrPKUZ5VjortFerg2yoBsQEIHhiMyBsiETwgGMH9g6FP059xSfOhjEMMMXRR3n333bCdO3fq3W63GDVqVPmNN95Y7u2aOhpv9cgMAPAlgKEARgK4TwgxUUqZdY72w6EOsSup3iGltAohLFXHaNNKSoA//UldpfpvfwOefx7Qcpg10QVzlbpQnlmO8l/Ka4KL47ij5nldmg7GIUbETY1D8MBgBA8IRmBsoBcrpvbK4XCI9PT0hPHjx5sWLVpkVhQFl112WfeIiAj38OHDOadOK/LKz6mU8h0A7wj1eqpJAJYB+FII0UtKWVTPS2IBmOsZ2FsOoFN97yGEuBfAvQAQExODjIyM5ir/guTl6fH4431RWKjD3LkHMHp0AbZubd73qKio8Nrna6v4ndTPp74XG4CDAA4A2F91n1/r+QQA3QH8seq+G2A32mGHHUWo+s/I/qpbA3zqe6E24ZlnnomZNWtWYY8ePZwAoNFoEBoa6lm/fn0Ig0zr8mq/gJRSAnhPCHEMwBYAdwNYeI7m1nr2+QFwnuPYbwB4AwCGDBkiR4wY0dRyL1hGhjo/jJ8f8O23wPDhvQD0aoH3yYA3Pl9bxu+kfm31e5Eeicp9lSj/qRyWnyyw/GhBZXalevIYQGDnQBgvNcI41IiQoSEIHhzcrKeC2ur3Qm3Tpk2bDPHx8a7qEFMtOzs76NZbby31Vl0dVZs4wSGl/L5qrEzaOZrkAahvRr8QNOrvrdb39tvAtGlAWhrw+edAaqq3KyJqO5zFTlh+tMCyTQ0t5Znl8FSoHa7acC2MvzMialwUQn4XAuMQIwJiArxcMdFvPv7447CXX375ZO19n376qVFKiQkTJpR5q66Oqk0EmSoHAJzrH8B2ANcKIUKklBYAEEKEAAgG0MwnappGUdR5YZ57Dhg1Sr28OizM21UReY/0SFTurYTlBwvKtpXB8oMFtoNqz7vQChj6GxAzOQYhw0IQcmkI9N301bN4ErU5drtduN1uodVqsWbNGuPmzZuDCwsL/TMzM4NXrVqVYzQalYaPQs2pLQWZXgDuP8dzKwA8DuBKAOuq9g2G2vH835YvrXGsVuDOO9XwMnWqeoWSPy+EoA7GY/OgPLMcZVvL1Nu2MnjK1N4W/2h/hFwegrgpcQi5TO1t8Qvy83LF1BKmfDYlaU/hniBv11Fbn059rCtuXpHXlGNs3rw5aOjQoZUAYDab/YQQCAwMVAoKCvyLiora0m9qh9HqX7oQQgfgPaiB5D9Qx7k8DmCllPJ4VZtHoF6efb2UslhKeVgI8TrUgcHVQWYagJellHtb+zPUp7gYuOEG4OefgUWLgIcf5qKP1DG4zC41sGxRg0v5L+U1CyUG9Q5CpwmdEDo8FKGXh0KXqmNvC/m0jIwM47Rp00oAYPLkyWYAZgBwOp1i4cKFsWPHjr2oy6/Xrl1rHDRokC0hIcHdjOV2CN5Ijy4AbgBLATwNdW6YBVLKnbXaxEAdE6Ovte9BAM9XBRoBYDeA51ql4gYcOwaMGQPk5qq9MePHe7siopbjLHaibEsZzN+ZUfZdGSp2VQBSXXvIONSIxFmJCL1CDS7tdJFEaoSm9ny0Vfn5+f7JycmuuvullCIvL++ir/V//PHHE9etW5fTtOo6plYPMlWXUN/WQJvZOHPpAkgp3VBn921Tdu1SF3y024GvvwZ+/3tvV0TUvFwlLpi/M8O0yQRzhhnWbPUCQo1eg5DLQpDydApCrwpFyLAQ+Ol5mojaL7f73J0lOTk5gfHx8fVeRdsQk8mksVqtmpSUlLMCEjWM5/Oa4NtvgbFjgZAQYOtWoHdvb1dE1HTuMjfMm80wf6uGl8rdlYAENEEahF4Ripi/xCDsyjAYhxqhCeAaRNRxbNu2LchgMJw1mDc/P1+bmZlpnD9/fh4A2Gw2MXfu3Ljk5GTHyZMnA0aOHFl+ySWXOFauXBl+4sSJgL59+1qdTqcmNzc3oF+/frbVq1eHhYaGeubMmRP3xBNPFISGhnLA8AVgkLlIH30E3HEH0LUrsGEDkJTk7YqILo7iVGD5wQLTNyaYvjHB8rMFUAARKBB6eShSnk1B+B/CGVyow9u0aZMxLy/vrLkAFi9eHJ2ammpPT08vBoCxY8emTp8+vWjcuHEWi8WiGTJkSM/09PTTDz/8cFGXLl36zp8//1RUVJSnf//+PR9++OHCXbt26UeOHFk+a9as4tb/VL6PQeYivPwyMHMmMHw4sGYNEB7u7YqIGk9Kico9lTB9rQYX83dmKFYF0AAhvwtB8txkhI0KQ8ilIfDT8VQRUbWcnJzAkSNHWpYvXx4+depUEwCsXLkybO3ateHr1q3L0el08quvvjKcOHEiYNy4cRYAKCoq8isoKAiYNGmSedu2bUGXXnppeVRUlAcACgsL/QsLC7Xbtm0zTpky5ZgXP5pPY5C5AFKqc8QsWACMGwe89x6g1zf8OiJvc5W4YPrGBPwb+OEvP8CZr57K1/fQI/buWESMjkDoVaFcOJHoHBRFQUBAgJw5c2bJBx98EDpjxowEAIiOjnZlZWXtCw4OlgCwZcuW4Msvv7zmyqX169eHDB06tDw8PFzZunVr8PDhwysA4MCBAwFOp1PTu3dvR15eXmD//v0d+fn52vj4eF61dIEYZBpJUdTlBpYuBe69F3j1VXXpAaK2SHokLD9bULqhFKUbSlGeWa4uu2oEQq8NRcSYCIRfEw5dks7bpRL5hKysLF2/fv2sADBx4sSyiRMn1juBa3JystNisfgB6liZt99+O+rtt98+BgA//fST4YEHHigCgCVLlkTPnj07PzMzU9+zZ08rAHz00UehM2fOLKnvuHRuDDKN4HYD99wDvPMO8MgjwMKFnCOG2h5XqQulX5ai5PMSlG4ohbvEXXO6KOXpFISPCceOyh3oPYqj0oku1MaNG43jx49vcPmBu+66y5Senh60ZMmSyMOHDwcuXbr0eL9+/RwAcOjQIX12drZuz549uoiICPdjjz1WVFRU5KfT6ZTFixdH3XTTTVze4CIwyDTA6QQmTQI++QR49ln11BJDDLUF1WNdStaVoOTzElh+UAfp+kf5I/K6SEReH4nw0eFnzuWS4bVyiXyawWBQunfv3uDl1VqtFsuWLTtZd//hw4f9IyMjXY8++mhR7f3R0dGeDRs2HGnOWjsaBpnzsFqBW25Rr0r617/UAb5E3qQ4FZg3m1GypgTFa4rhyHUAAIIHByP5iWREXh8J4xAjhB/TNlFzevDBB5t0yue7774LHjBggLW56qHfMMicg8UC3HgjsGUL8OabwF//6u2KqKNymV0o/bwUxWuKUbqhFB6LBxqdBuGjw5H8pBpeAuMuekJRImphmZmZusWLF8cajUbPr7/+Gti3b1+Ht2tqTxhk6lFaqs7Wu2MH8MEHwG3nnYeYqPk5TjlQ/Fkxiv9XDPO3Zki3hH+MPzr9uRMib4pE+KhwLrZI5COGDh1q37Nnzz5v19FeMcjUUVICXH01sG8fsHq1uhAkUWuwHbahaHURilcXq+NdJKDvpkfiw4mIGheFkN+FQGh4yoiIqDYGmVpqh5g1a4BrrvF2RdTeWXOsKFpVhKJVRajYUQEACB4YjJR5KYgeH42gS4K4WjQR0XkwyFRhiKHWYj1UK7zsVMNLyKUhSHsxDVHjoqDvwlkWiYgai0EGaogZNQrYv58hhlqGPdeOwg8LUfhh4ZnhZXEaom+Jhq4zJ6YjIroYHT7IFBerPTEHDjDEUPNyFjhRuKoQhR8UwrLNAgAwDjOq4eVP0ZxVl4ioGXToIFM7xHz2GUMMNZ3b4kbR/4pQ+H4hTBtNgAIY+hrQ5R9d0GlCJ+hTedqIiKg5ddggUz0mhiGGmkpxKzB9bULBygIUf1oMxaZAl6pD58c7I2ZiDAy9Dd4ukYio3eqQQcZiUeeJ4ZgYulhSSlTsqEDBygIUfFAAV4EL2ggtYu+KRcwdMQi5NIRXGxERtYIOF2SsVnVumJ071XliGGLoQjgLnSh4twCnVpyCNdsKESAQeUMkYu6IQeR1kdAEaLxdIhFRh9KhgozDAYwbB3z/PfD++5zszpdJKeGRHrgVNzxK1b30QJEKFKmg1FmK/PL8mm1FKme8tj4aoYEQAgKi5l4jNBBuAes3VphXmmFZbwHcgOF3BnRZ2gWdbusEXaSOvS9ERF7SYYKM2w1MnAh89RXw1ltcdqA5OD1OlDvKUe4sR6WzEpWuyrPurS4rbC4bbG7bWfd2tx0OjwMOt+Ose6fHCafHCZfigsvjOuPerbjPCCbn9EPTPl9SURKu23Edrtl1DSIqI1BqKMVXv/sKGwZsQG6nXKAYwDK1rUZooNVo4Sf8oNVoodVo4e/nD3+N/1n3AX4BZ938/dT9gX6B6k179r1OqzvnTa/Vq/f+eui1+pp7nZYhi4jatw4TZO6+Wz2VtGQJMGWKt6vxPiklLA4LSm2lZ9zMdjPMdjPKHGUos5fB7DCjzF6GMkcZLA5LTXCxOCxwehpc0f4M/hr/M35odVrdGT/WhgADIvwiEKgN/O0HXnN2GKgOCrWDg5/GD37CD34aP2iEBjmHctCzR09ohEbtVanqZakmcOaPu4SElBLSLhG0MQih/wuFfoce0k+i/Ipy5N6Yi7JLy9BD0wPdZDcoUqnpBardK1T7Vjd81d6uDmtWl7UmtDk86r66wa5Roe089Fo9gvyDEOQfBLiAyAORNdtB/kEw+BvOemwIMMDgb0BwQHDNY0OAul19q27LoERE3tQhgszx40BWFvCPfwAPPeTtalqGW3GjqLIIBZUFKKgoQGFlIYqsRSi2FqOosgjFtqp7azFKbCUw2UzwSM85j6fVaBEaGIowXRhCdaEIDQxFl7AuMAYaERIQAmOgEcYAY8197R+7uvfVPQN+mtZb5DCjMgMjBo9odPvKvZU4tfwUTr9zGm6TG/quesS9EIfYO2MREBPQcoU2gltxq71XbgfsbnvNzeFRt2v3cJ2r98vqssLqsiL3ZC6CQ4Nrtosqi2B1WWt6zyqdlXAprkbXJiDOCjjGAKN6H2hEsH9wzePa/15CAs/8NxQSGIKQwBAY/A0MRkR0QTpEkAkvOoAfk0Yg9SsAX3m7mgujSKXqr3UnnB4HHO6q+6ptp8cJrduBrcpvoSQQQFLVTQhU9WQE1OnZSIBW4w9/jVa996vu5VB7PNRejDMqAWCuurV9A8xmICzsvG2kAriKnHDkO+GxuBElgLgofwT0D4Q2TAvxBYAvWqfe89ECCK66NZXZbEZYWPW/FX8AYVW33yhSQpEeeKQHHkWpGnvkgUc5c99v29X3FnikqWrbXbPfraht6lNZdTtda5+fxg9ajR/8hFa912jPu632yNV6LNQeuguJQ43590J0MUwmk+b555/v9MILL5xuuDVdjA4RZAL8gS6p3q7ibBLqOBP1r2077DWnFKr/2nbA5Tn7r2MhBAKrTr3otXrohB5GfTD8a8Zc+NeclvHTaC/oP+gdgeJQ4Mh3wnnKAemS0Og10KXqERAbAI0/vy2NENAILbTQAs3UiSaBmjDkrgk+7jOCjke6qx5X3dfadngcv+0/Ryiq68yQc2bQ0dbZZ3c7oHFozmh/dpgnOtu7774btnPnTr3b7RajRo0qv/HGG8urn3O5XJg8eXLy0qVLT3izxrZo6tSpiU899dTppKQkd1OP1SGCTGC/HhAZGV557zJ7GQ6bDiOnNAdHTEdwzHwMR81Hccx8DLnmXDg8jjPaGwOMSApNQufQzkgKSUJiSCISQxIRb4yvuUXqI8/ofs/IyEDfESNa+ZO1bTszMjCi1nciFQnTNyacfPUkStaWAAAib4xEwgMJCB0VDqHpGD9Zdb+X1iKgZiI/AE09UadIpWasVpldHbtVPYarejxXzX3tx7XuzfaC855aBYAAvwCE6cLOvgWevS9cH67e68Jr9gVqA5v4SdsAnuY7J4fDIdLT0xPGjx9vWrRokVlRFFx22WXdIyIi3MOHD7cBwPPPP99p9OjRluTk5HOer12yZEnkpk2bQsLCwtwul0t07tzZeeONN5Z99NFH4UuWLMkHgG+++cbw4osvxmzZsiXUbreLa665xuzv7y8BwOVyiePHjwdmZ2cHzZs3L++pp54qbK7PaLfbxauvvhpZXFysLS8v99uzZ49+zJgxZY899lhRY16/Y8cO3YMPPpj05JNPnrr22msraj/37LPPnr799tuT165deyQoKKj+S0kbqUMEmZZW7ijHwZKDOFByAIdKDiHHlIOcUvVWbC0+o22kPhJdwrugX0w/3NT9JnQJ74Lk0GQkhyUjKSQJobpQL32K9sld7sbpd07j5MsnYTtkg3+0Pzo/1hnx0+K5UKOP0giNOm5LF4rEkMSLOoaUElaXFWWOMpjtZny77VukXpJas133Vr3/eNlxmGwmmOymBge767X6mpBTHXDOeKwLr9mue89B1G3fM888EzNr1qzCHj16OAFAo9EgNDTUs379+pDhw4fbSktLNe+++27U7t2799X3ekVRcMstt6QkJCS4/ve//x3181O7Pz/88MPQq6++uufy5cuPVre9+uqrK6+++uojiYmJfXv06GH94osvjtQ93sSJE5MTExMbP8CtEWbMmJHw448/GjMzM/fr9Xq5ZcuWoCuvvLJXRUWFZv78+QXnet37778f+sknn4QbjUbP999/H+J2u886rRYXF+e+4447SmbOnJnwxhtvNKnHikGmkaSUOFVxCtmF2dhXvA/7i/fjQMkBHCg+gJPlJ2vaCQgkhSaha0RXjO85HmkRaega0RVp4WlIDU+FMdDoxU/RgZwCcv6Wg1NvnoLH4oFxmCfMrB0AACAASURBVBG93u2F6D9FQxPISes6OiHUQcqGAAPijfEoDC3EiG4jLugYdrcdJpsJZY4ymGwmmO1mmOymmvBTHXiq9+eX52Nv0V6Y7CaU2csgce4/Qv01/jWhJkIf8VvQqbNd/ThCH1Gz3S56gtq4TZs2GeLj413VIaZadnZ20K233loKAK+99lrUmDFjzIGBgfX+D/36669HZGdnB61evXpv7f0TJkwomz17tuu6666z1N5/8ODBgJMnTwakp6efqu94Q4YMqezWrZujvuculkajgclk0rrdbgFADhgwwA4AmzdvDgFwziAzadKkskmTJpUdOHAg4J133ul0rnaTJ082L1y4MO7w4cP+aWlpFx3CGGTqUWwtxu6C3dhTuAfZhdnILlJvZvtvA13DdGHoEdkDo1JHoUdkD/SM6okekT2QFpEGnZZ/6XuDlBJlW8pw4qUTwGfASc1JRN8ajcT0RIQMC/F2edTO6LQ6xBnjEGeMu+DXKlJBmb0MJrupJvCc677UVorTFaexr2gfSm2lKHOUnffYeq3+jIBTO+TUfVz7FhLIZTUa6+OPPw57+eWXT9be9+mnnxqllJgwYUIZAGzYsCE0PT39nD/2q1evDk9MTKy3W+/SSy+tiIiIOGPehS+++MIIADfffHNNwFmxYkX4nXfeafLz80NCQoKrd+/e9qZ8rrqqekpqekt27dqlA4Bhw4ZVnPNFF0Cj0eCPf/yj+bXXXotatGhRvQGtMTp0kHErbhwqOYRdBbuw6/Qu9b5gF/LL82vaROgj0Du6Nyb0noDenXqjd3RvXBJ9CToZOvH/9G2E4lZQ9FER8l7MQ8X2CmgjtMAEYNgLw6BLZKiktkcjNGqvij4cCL+w13oUT00vUKmt9KzQU2orVZ+zq/eHSw/XPGd1Wc95XD/hd0b4idBHNOkzTpmCpD17ENSkgzSzPn1gXbECeU05ht1uF263W2i1WqxZs8a4efPm4MLCQv/MzMzgVatW5RiNRsXlciEzM9N41VVXnXUKqJrBYPCsXbs2Yu7cubEPPPBAcUJCQs2g1zfeeON43fbffvttSFpamr179+5OADh16pT2vffei5gyZYoJQE2AakkLFiyIHTFiRNmCBQsuOnTUNWLEiIpHH300kUGmEdyKG/uK9uGX/F+QdSoLWaeysOv0LtjcNgBqV26v6F4Y1WUU+sf0R7+YfujTqQ9ig2MZWNood4Ubp948hRP/OgHHcQeCegah+xvdEfOXGGz5eQtDDLVLfhq/mqCRhrQLem316bDak2DWDkC1bwUV5+xM6NA2b94cNHTo0EoAMJvNfkIIBAYGKgUFBf5FRUVaADhx4oS/oiiIjo4+54jyp5566vTWrVtDnnvuuYTnnnsuITk52XHLLbeUzJs373R4ePhZs2D+8MMPRqPR6LnllltSLBaL39atW0Nuu+224vqO3dzmzp0bm5OTE2iz2TQffPDBUX9//2Y79vDhwysPHjyot9lsQq/XX9Sg3w4RZPYX70fIcyE1oSU4IBiD4gZh+pDpGBA7AP1j+qNXdC8E+Hl34jNqHMdpB06+fBL5r+XDbXYj9Peh6LasGyKvi+wwVx8RXYwLPR0mpl38/5+a2vPRVmVkZBinTZtWAqhjPFA1uZbT6RQLFy6MHTt2bPmJEyf8jUbjeS+LGzhwoP3IkSN7Pvnkk5CNGzeGbN68OWTx4sXxJ0+eDPjwww9za7fdtWtXYFFRkf+CBQvy7rnnHhMAPPnkkzHnGhMzbty4lOLi4kaljcjISPenn3569HxtFixYcBoAPv/88+CePXv2feeddw7ffPPN5ed7TWNFREQoGo0G+/fvDxw4cOBFnRrrEEFGCIHpQ6ZjcNxgDI4fjO6R3aERHPDpa6w5VuQtzMPpd05DuiSixkeh8+zOHP9CRK0mPz/fv77LqaWUIi8vLxAA3G73eXvyy8rKNKGhoYrRaFTuuusu81133WVWFAXXXHNN2hdffBEO4Iwgs2HDhhAAuPrqq2vGpkgpxZgxY+oNE6tXrz52UR+uAddff31FSkqK/a677krNzc39NSQkpGnrp1QJDQ11l5SUXPSsVR3i17xHZA8sHrMYf+n3F/SM6skQ42Mqfq3A3kl78XOPn3H6P6cRNyUOvzv4O/T5uA9DDBG1Grf73HO35eTkBMbHxzsB9dJii8VS7w/zgQMHAhYtWhRdd79Go8EVV1xRYTAYzurJycjIMCYkJDhTUlJqAtS8efNOx8bGNm52yItw/PhxbadOnfpNnDgxufb+pKQkh9ls1mZlZTXbuXu32y10Ot1FzyXTIXpkyDdZfrIgd0EuStaUwC/YD0l/S0Liw4kIjOXlpUTU+rZt2xZkMBjO6oXIz8/XZmZmGufPn58HAImJiS5FUYTJZNLUHe+ybt26kKrLmc+yadOmkAkTJpTU3qcoCn788UfjH/7whzMG82q15/75bo5TS6dOnfIvKiryLy0tPeONSkpK/P38/GT1oOOmUhQF5eXl2ri4OF5+Te2HKcOE3Pm5MG80QxuuRcozKUh4MAH+Ec03wIyI6EJt2rTJmJeXd9ZgysWLF0enpqba09PTiwEgODhY9u/fv+L777833HDDDeV1jhGSnZ0ddPvtt5t69+7tAACPx4M5c+bEuVwuUfeKoJ9++klvNpu1w4cPb/Qlz81xamno0KG2K664wvLMM8/UXMZ78ODBgF9++SV46tSpBdVXWX300Uch06ZN67JixYqj48aNO2PuG49H7TBSlHOfgSosLPTTarUyMTHxopcq4DkWajPM35mxY8QO7PrDLlizrUj9Zyouzb0UKU+nMMQQkdfl5OQEjhw50rJ8+fKai+ZXrlwZtnbt2vB169bl1D49cs0115Rt2rSp3hlQN27ceOD//u//Yq+77rrU8ePHp1x77bVper1eycjIOFR95U5mZqbu+uuvT50wYUIaALz++uudrr/++tSNGzcaWvpzAuqprs8+++zwypUrI26//fbOd9xxR+d77rknecGCBceXLVtWM4eOlFK43W7hdDprepm+/PLL4Ouvvz71j3/8Y3cAuPfee7tce+21qf/5z3/OWpk1IyMjeOjQoeXnmjiwMdgjQ15n3mzGsaePwZxhRkBcALou6Yq4qXHw0zfTioVERE2kKAoCAgLkzJkzSz744IPQGTNmJABAdHS0Kysra19wcPAZP8TTp08vuf7667t6PJ786uUHAODzzz8/AgB1r0yqa+jQofbqtt4SERGhvPrqqyfP1+a2224ru+2223bW3jdmzJiKMWPGNKoHKSMjwzhmzJgmzYHDIENeY95sxrFnjsH8rRkBsQHo+lJXxN3LAENEbU9WVpauX79+VgCYOHFi2cSJE8/745ucnOwaPXp02TvvvBNePWkdnclut4uMjIyQ77///kBTjsNTS9TqLD9bsGv0Luy8aicq91Yi7V9pGHZkGBLTExliiKhN2rhxo/G66667oLlTnn322dPvvfdeRH5+PjsN6vHKK69ETpw4sTgyMrJJV1+d98sVQtxU9bBESvl9U96IqDK7Ekf/fhTFq4vhH+WPtBfTED89Hn5BDC9E1LYZDAblQq/UMRqNyooVK44/8cQTcW+//Xa7nCDwYu3Zsyfwm2++CV2zZk2TT5811CPzKYDBAHY09Y3ORQjxmhAioxHtBgghZK2bQwgR31J1UfOxHbVh3537kNk3E6ZvTEiZl4JhR4Yh6eEkhhgi8gkPPvhgScOtzpacnOxiiDnbm2++Gfnhhx8ebcog32oNdXdtk1I+Xb0hhHgBwK0AtgH4WEr5aVPeXAgxGsA0AJsb0fwRAHMAVH/o41LK/PO0Jy9zFjhx7P+O4dQbpyD8BJL+loTOj3WGfySvQCIi6sheeumlZvv9bijIHK69IaV8VAhxHYDJUsomTU0shAgFMBPAj41oOwBAvpTyn015T2odnkoP8l7Mw/GFxyEdEnH3xCH5yWQEJnAiOyIial4NBZn6wkpOfSFGCNFDSnkhI4//CeAJAC81ou3TAIYKIUIA/EdKue0C3odaifRInP73aRz9+1E4TzkRNT4Kqc+nIqhbkLdLIyKidqqhMTL1nbs61/msPzX2TYUQNwLIlVLubETbEAAOACcB3ANgqxBikRCCgyvaCCklStaX4JcBv+DAPQegS9Fh4NaB6PNJH4YYIiJqUQ31yFwnhKg7fqVXPft0APoD+EdDbyiEiAQwFcD4xhQopbQAmFD12iQArwD4G4BKqD0153qfewHcCwAxMTHIyMhozNv5pIqKCu99viMAXgWQBSAewDOA5UoLdrh2AF4qCfDyd9KG8XupH78XIt/VUJDpBLUHxl5rXwWApDrtghpxrGovAnhUSnnB6ypIKfOEELcAWAtgphDi/851HCnlGwDeAIAhQ4bIESNGXOjb+YyMjAy09udzFjtx7O/HkP9GPrRhWqQsSUH89HhoAtrG1ETe+E58Ab+X+vF7IfJdDYWPhVLKxxpzICHEvEa0+TOA3VLKfY05Zn2klB4hxJNQ+wCiAZxq4CXUjBSXgvxX83HsmWNwl7uR8EACUp7hWkhEROQdDQWZ98/3pBAiGkCllNIK4ONGvN99AEYIIV6s51gSwN1Syn834jgHANgAXNR1/XRxSjaU4PCsw7DutyJ8dDi6vtQVhktaZf0yIiKiejUUZFKEEClVjxUp5ToAEELEAfgAwO8BuIQQ/5JSPt6I95sKILjOvjer7u8BcLwxRQMYBuANKeUFzbJIF8d22IZD6YdQ+nkp9F316LOmDyJviIQQouEXExERtaCGgsynAKwA5qGqx0Wov16rAfwOwD4AXwG4VQhxUEr59vkOJqXMqbtPCFFR9dzOqu1HoE66d72UslgIMQHAtQDmSSmPCiH6A5gCNRRRC/LYPDj+wnEcf/44NP4apC5MReJDidAEto1xMERERI0ZoPtnKeUXtbbvhhpiMgFcKaV0CCGeA/ARgPMGmUaKAdAFgL5quwTAVQCyhRDbAXwO9RSUqxnei86h5IsSHHrwEOxH7Ii+LRpdX+zKCe2IiKjNaSjI7KsdYoQQAVAveVYA3C+ldACAlLJQCGG+mAKklCPqbM8GMLvW9tdQgw21AnuuHTkzc1D8aTGCegah/zf9ET4q3NtlERER1auhIHOozvZDUC+9fk9KmVXnubpjX8iHKC4FeS/mIffZXEAAXZ7rgqSHk9rM5dRERG3Ru+++G7Zz50692+0Wo0aNKr/xxhvLvV1TR9NQkAkSQoRKKcuEEF0BPAV1IrozBvZWTXJ3WQvVSC3M8pMFB6YeQOWvlYgaF4WuL3WFrrPO22UREbVZDodDpKenJ4wfP960aNEis6IouOyyy7pHRES4hw8fbvN2fR1JQ0Hm/wH4QQjxOYCJUHtdZkgpT1Q3EEIYAPwH6uy+5EPc5W4cfeIoTi49iYD4APT5rA+iborydllERG3eM888EzNr1qzCHj16OAFAo9EgNDTUs379+hAGmdZ13iAjpfyfEEIL4EEABVCvHFpe/bwQ4gEAfwQQAoALOfqQ4rXFOHT/IThOOpDwQAK6/KMLtCGNnZyZiKjj2rRpkyE+Pt5VHWKqZWdnB916662l3qqro2rwl0tK+RHUK5Lqe24ZgGXNXRS1HMcpB3LSc1C0qgiGPgZcsuoShF4a6u2yiIh8xscffxz28ssvn6y979NPPzVKKTFhwoQyb9XVUfFP8A5CSomCdwuQ81AOPDYPuvyjC5Ie4WBeIqILYbfbhdvtFlqtFmvWrDFu3rw5uLCw0D8zMzN41apVOUajUfF2jR0Ng0wH4Mh34OD0gyhZW4KQy0PQc0VPBPUI8nZZRNSO7Z+yP6lyT2Wb+g+NoY/B2nNFz7ymHGPz5s1BQ4cOrQQAs9nsJ4RAYGCgUlBQ4F9UVMTfVC/gl96OSSlRsLIAOek5UOwK0hanIfGhRAg/Li1ARHQxMjIyjNOmTSsBgMmTJ5sBmAHA6XSKhQsXxo4dO/aiLr9eu3atcdCgQbaEhAR3U+pbuXJl2A8//GAoKSnR/ve//81tyrF8BYNMO+XId+DgtIMoWVeCkOEh6Pl2TwR1a1N/HBFRO9bUno+2Kj8/3z85OfmsmeWllCIvL++ipz9//PHHE9etW3fWMj4X6u9//3viunXrDmVnZ3eYK4kZZNoZKSUK3itAzoM5UBwK0v6VhsQH2QtDRNRUbve5O0tycnIC4+PjL2ohY5PJpLFarZqUlJQmLb2zf//+gMjISFefPn0cffr0cTTlWL6EQaYdcZW6cHD6QRStKmIvDBFRM9u2bVuQwWA4azBvfn6+NjMz0zh//vw8ALDZbGLu3LlxycnJjpMnTwaMHDmy/JJLLnGsXLky/MSJEwF9+/a1Op1OTW5ubkC/fv1sq1evDgsNDfXMmTMn7oknnigIDQ2td8Cwx+PBCy+8EB0VFeXOzc0N/Nvf/lYYERGhAMA333xjWLx4cYwQAo8++mjsP/7xj9Nabcf4ie8Yn7IDKP2qFPvv2g9XsQtdnuuCzrM7sxeGiKgZbdq0yZiXlxdQd//ixYujU1NT7enp6cUAMHbs2NTp06cXjRs3zmKxWDRDhgzpmZ6efvrhhx8u6tKlS9/58+efioqK8vTv37/nww8/XLhr1y79yJEjy2fNmlV8vvdfvXp1SGVlpd/cuXOLpkyZknTs2LGAiIgIOwBcffXVlV9//bU1MjLSM2fOnKKW+QbaJl576+M8Vg/wMrB7zG5ow7UY9NMgJD+WzBBDRNTMcnJyAkeOHGlZvnx5zUq6K1euDFu7dm34unXrcnQ6nfzqq68MJ06cCBg3bpwFAIqKivwKCgoCJk2aZN62bVvQpZdeWh4VFeUBgMLCQv/CwkLttm3bjCNHjmxwkHCnTp3cr7zySuzIkSO7jh071jxo0CB77ef37t2r79evX4ebVZhBxoeVZ5Uja3AWsBpInJmIwb8MhnGg0dtlERG1O4qiICAgQM6cObMkODhYmTFjRsKMGTMSjhw5EpCVlbWvW7duTgDYsmVL8OWXX14TStavXx8ydOjQ8vDwcGXr1q3Bw4cPrwCAAwcOBDidTk3v3r0deXl5gf3793fk5+ef9yzJlVdeaf3ll1/2XnXVVZZ77723S93nDx06pB84cCCDDLV9UpE4/sJxbL90O9zlbmAR0PVfXeGn9/N2aURE7VJWVpauX79+VgCYOHFi2dKlS08uXbr05NNPP10YHBwsq9slJyc7g4KCFEAdK/P2229HLV68+AQA/PTTT4ZLLrnEDgBLliyJnj17dn5mZqa+Z8+eVgD46KOPaqZZ//e//x2Wk5PjX729fv364N///vfdunfv7nz66acLhw4dekYPjsvlQnl5uV9MTIynJb+HtohjZHyM47QD++/YD9M3JkTfGo3ur3fH97u+93ZZRETt2saNG43jx49vcPmBu+66y5Senh60ZMmSyMOHDwcuXbr0eL9+/RyA2mOSnZ2t27Nnjy4iIsL92GOPFRUVFfnpdDpl8eLFUTfddFPN8ZcsWRJ74MAB83PPPXcaALp16+YYPXp02dKlSyMLCwu18+bNO1X7fXfv3q1LTU21owNikPEhJRtKsH/yfngqPOi+vDvi/hoHITgWhoiopRkMBqV79+4NXl6t1WqxbNmyk3X3Hz582D8yMtL16KOPnjEQNzo62rNhw4Yjddvv2LFj/6JFi6Kqt7t27ep66qmnCs/1vtu3b9dfccUVFzUZn6/jqSUfoDgVHJ59GL/+8VcExARg8C+DEX9PPEMMEVErefDBB0ua8vrvvvsueMCAAdbGtv/555/1gwYNanC8y8aNGw1XXHFFt127dgXdfvvtpqbU6KvYI9PG2Q7bsHfiXpRnliP+/nikLUrjWBgiIh+SmZmpW7x4cazRaPT8+uuvgX379m1wsjqLxaK5+uqrKxtqFxwcrPTv398aFRXl7t27d4eZBK82Bpk2rHBVIQ789QCEn0Dv//VG9Lhob5dEREQXaOjQofY9e/bsu5DXNCbEAMCwYcNsw4YNO+tUVkfCINMGKU4Fh+ccxsklJxFyaQgu+e8l0HXuMMtmEBERNRqDTBtjz7Nj75/3wvKjBYkzE5H6Qio0ARzKREREVB8GmTak9KtS7PvLPigOBZesugSd/tTJ2yURERG1aQwybYD0SOTOz8Wxecdg6GNA7497I6g7F3skIiJqCIOMl7lKXdg7aS9MX5oQMzkG3V/rDr8gXpVERETUGAwyXlSxuwJ7xu6B46QD3V/vjripnOCOiIjoQjDIeEnhR4XYf/d+aMO0GLh5IEKGhXi7JCIiIp/DINPKpEfiyBNHkPdCHkIuD0HvT3ojMDbQ22URERH5JAaZVuQyubB3ojoeJn56PLou6cpLq4mIiJqAQaaVVGZX4tebf4XjuDoeJv7eeG+XRERE5PMYZFpByecl2DthL/yC/TAgYwBCLw/1dklERETtAoNMC5JS4sRLJ3D4kcMIHhCMvmv6IjCB42GIiIiaC4NMC1FcCg7NOIRTb5xC1Pgo9PpPL/gZOD8MEVF78u6774bt3LlT73a7xahRo8pvvPHGcm/X1NEwyLQAV6kL2bdmw7zJjM6Pd0aX+V0gNJwfhoiovXA4HCI9PT1h/PjxpkWLFpkVRcFll13WPSIiwj18+HCbt+vrSHjJTDOzHrRi+6XbUbalDD3f6YnUBakMMURE7cwzzzwTM2vWrMJrrrmmEgA0Gg1CQ0M969ev56RgrYxBphmZMkzYful2uE1u9N/UH7GTY71dEhERNbNNmzYZ4uPjXT169HDW3p+dnR2UnJzsPNfrqGUwyDSTgvcLsPua3QiIDcCgnwYh7Iowb5dEREQt4OOPPw677777Smrv+/TTT41SSkyYMKHMW3V1VBwj00RSSuQtzMORx44gbEQYeq/uDf8wf2+XRUTkXVOmJGHPniBvl3GGPn2sWLEirymHsNvtwu12C61WizVr1hg3b94cXFhY6J+ZmRm8atWqHKPRqDRXudQ4DDJNID0Shx46hPxX89FpQif0/HdPaALZyUVE1F5t3rw5aOjQoZUAYDab/YQQCAwMVAoKCvyLior4m+oF/NIvksfqwd5Je1HyWQmS5iQh9TkO6iUiqtHEno+2KiMjwzht2rQSAJg8ebIZgBkAnE6nWLhwYezYsWPLAWDt2rXGQYMG2RISEtxeLLdDYPfBRXAWObFz5E6UrClB11e6Iu2FNIYYIqIOID8/3z85OdlVd7+UUuTl5dXMePr4448nulwu/jC0Aq/3yAghXgPQS0o5ooF2jwPoDcABoBDAk1JKT8tXeCbbYRt2X7sbjhMO9P6kN6LHRbd2CURE5AVu97k7V3JycgLj4+OdAGAymTRWq1WTkpJyVuCh5ufVHhkhxGgA0xrR7nkAV0kpb5dS/hVAIoCXW7q+uip2VWD78O1wmVzov6k/QwwRUQeybdu2IIPBcNZg3vz8fG1mZqZx/PjxphUrVoRPnjw5OTQ01DNnzpy4srIynvloYV7rkRFChAKYCeDHBtr1BPAIgLG1di8H8J0Q4k0p5Y6Wq/I35q1m/HrDr9AateiX0Q+GnobWeFsiImojNm3aZMzLywuou3/x4sXRqamp9vT09GKdTid37dqlHzlyZPmsWbOKvVFnR+PNpPhPAE8AaGjyoIkA/AD8XGvfLwAUAHe0TGlnKvmipGaOmIHfD2SIISLqgHJycgJHjhxpWb58eXj1vpUrV4atXbs2fN26dTk6nU4CwLZt24wjR47kmkutxCs9MkKIGwHkSil3CtHgWKjhACSAmsmHpJRWIYQFwNCWq1JV8H4B9t+5H4b+BvRb3w8B0WeFcSIiaucURUFAQICcOXNmyQcffBA6Y8aMBACIjo52ZWVl7QsODpYAYLPZRF5eXmD//v0d+fn52vj4eF611MJaPcgIISIBTAUwvpEviQVgrmdgbzmATud5n3sB3AsAMTExyMjIuPBiVwN4BUB/oGJeBbZlb7vwY7SCioqKi/t87Ri/k/rxe6kfvxdqSFZWlq5fv35WAJg4cWLZxIkT653BNzMzU9+zZ08rAHz00UehM2fOLKmvHTUfb/TIvAjgUSnlhaRUaz37/HCe01JSyjcAvAEAQ4YMkSNGjGj0m0kpkft/uTj28jFE3hyJSz68BH46vwsot3VlZGTgQj5fR8DvpH78XurH74UasnHjRuP48eMbXH6gV69eDp1OpyxevDjqpptu4nIFraBVg4wQ4s8Adksp913Ay/IAdKlnfwiA/c1SWC1SShz+22Gc+NcJxN4Vi+7Lu0Oj5aBzIqKOzGAwKN27d29wQcjo6GjPhg0bjrRGTaRq7V/o+wC8KISQ1TcAVwG4qmr7rnpesx1AkBCiZmn0qsfBALY2Z3FSkTj0wCGc+NcJJDyUgB5v9WCIISIiPPjggzxF1Ea19q/0VAAD69yyqm4DAayp5zUroA72vbLWvsFQr1r6b3MVJj0SB6YeQP5r+Uiak4SuL3XlbL1ERERtXKueWpJS5tTdJ4SoqHpuZ9X2IwBuBXC9lLJYSnlYCPE6gEkA1lW9bBqAl6WUe5ujLsWtYP9d+1H4XiGSn0pGyjMpaMTVVERERORlXl+ioB4xUMfE6GvtexDA81WBRgDYDeC55ngzxalg31/2oejjInRZ0AXJjyc3x2GJiIioFXg9yNRdY0lKORvA7Dr73FBn921WikNB9q3ZKFlbgrTFaUialdTcb0FEREQtyOtBxls8Ng/2jNsD05cmdFvWDQn3J3i7JCIiX6RIKXk6nlqUlPKcz3XIS3I8dg/2jN0D01cm9HizB0MMEdFFEkKYXS5Xh/2jmFqH0+n0F0KY6nuuwwUZxaEge3x2TYiJ+2uct0siIvJZUspdFouFC9BRizKbzUZFUdbX91yHCjKKQ0H2n7JRur4U3d/ojrgpDDFERE3hdrsXnD59WlZUVASdr/uf6EJJKeFwOPwLCgoiQPZuegAACs9JREFUCgoKzB6P54362nWY7kDFqSD7tmyUrCtBt9e6IX5qvLdLIiLyeYMHDz6UlZX10JEjR+ZIKbtAvbKUqFkIIUyKonzo8XjeGDx48LH62nSMICOBvRP3ouSzEnRb2g0J0zkmhoiouQwePHgDgA3eroM6pg4RZGxHbSjeXoyuL3VFwgMMMURERO1Fhxgj4za5kfZiGhLTE71dChERETWjDhFkAhMCkfQwJ7sjIiJqbzpEkAmIDfB2CURERNQCOkSQISIiovaJQYaIiIh8FoMMERER+SwGGSIiIvJZDDJERETksxhkiIiIyGcxyBAREZHPYpAhIiIin8UgQ0RERD6LQYaIiIh8FoMMERER+SwGGSIiIvJZDDJERETksxhkiIiIyGcxyBAREZHPYpAhIiIin8UgQ0RERD6LQYaIiIh8FoMMERER+SwGGSIiIvJZDDJERETksxhkiIiIyGcxyBAREZHPYpAhIiIin8UgQ0RERD6LQYaIiIh8FoMMERER+SwGGSIiIvJZDDJERETksxhkiIiIyGcxyBAREZHPYpAhIiIin+WVICNUfxdCnBBCmIUQnwkhEhrxugFCCFnr5hBCxLdGzURERNT2aL30vo8AKAVwI4CRAF4A8AGAKxvxujkAZNX2cSllfksVSURERG1bqwcZIYQAsEVK+WPVrh1CiAEAxjbwugEA8uX/b+/uY66s6ziOvz/eoYHADJQHtfKhTYgJdZPDslKbt6HQpoWTHqZQgS2zVdZ6wkE2tc1WqLiKOSqd0yVNLaexkpBIeXLS0AlLh9lYCyQNpVsN+/bHdd1xOpz78dznus7vnM9rO7s419P95bsf9778fr/r+kXc2OgYzczMLA2FFzIREcDGqt17gXv6uXQpcIakscDtEfFoI+IzMzOzdJQ+2VfSOOB44It9nDMWeA3YDXwW2CDp+5I6ionSzMzMmlFZc2SQNBq4FPg2MAroAu6tdW5E7Afm59e9FbgFuBo4QNZTU+v+i4HF+ddXJO0czvibzLHAC2UH0WSck9qcl9qcl9pOKzsAs/4oG+kp4QdnvSknAHOAZcA4YEpEPDvAa38NnAWMj4iDDQy16UnaGhHvKTuOZuKc1Oa81Oa81Oa8WApKG1qKiDci4vmI+BEwl6x3aO5ArwWWAGOB4xoXpZmZmTWz0ufIAETEFmA78OogLtsJdAP7GhKUmZmZNb2mKGRyI4BHBnH+LGBlRLzeoHhSsrLsAJqQc1Kb81Kb81Kb82JNr/A5MpKOBr4ArImIbfm+K4FJEXFN/v2rwCXAnIh4QdJ8YDbwnYjYJWkG8DVgUUR0F/oXMDMzs6ZRRo/MGODjwCZJf5C0Enixp4jJTQROBkbm3/cBZwNPSdoAXAgsdBFjZmbW3kp7asnMzMysXs00R8b6IekcSWskXV7j2GRJd0m6RdJqSV1lxFi0vnKSH19etdDoA0XHWCRJI/OXRe6WtE/SvZJOqTpnat5Gbs63nWXFW5SB5CU/776q9rKijHiLMpAFfNuxvVhaSnshng2OpAvIXgp4PtkCm5XHxgDrgWsi4m5J44GnJV0cEX8sPtpi9JWT/PgkYArZfKoeDxYTXWmWky2q+hVgOtlCq52SZkTES/kLJdcBF0XEY5JOAx6VNCsinikt6sbrMy/wv/Xc/sP/t5e7iw60YH0u4NvG7cVSEhH+JPIB3kn2y3hB1f7ryNarUsW+VcC2smMuKyf5sR8Ap5cdY4G5GAd8o2rfosr8AHcCW6vOWQvcV3b8ZealIjcTyo63wLwIOLNq3x3Ay1U5aav24k96Hw8tpaW3yc2fBLZEROWEp83ADEnTGx9WqWrmJO+NuQJYKWmJpGOLDas0N1V971mMdbykkcBHydpGpc3AXElvaXRwJeo1LwD5k5CXAHdJ+nL+dGVLi0yvC/i2eXuxhLiQScthM7MlnQi8newXUKW/5dszGh1UyXqbrd4J/BaYDHwX2C7p/MKiKkFE/CMOf5JvRL5dD8wE3kztttJBlrOWM4C8ALwLWANMI+vJe0LSzIJCbAo1FvBty/Zi6XEhk75J+bb6Dccv59sJBcbSNCLiwYi4KCJOIntc/yBwv6Sp5UZWuDnAQ5G9Pdtt5ZDKvBARP4+Ij5AVvp8iy8UaSS2/BIqk0ZI+A2wFziFbwBfcXiwRLmRax7+qvnfk27Z/83FEPAScC7wBXFVyOIWRdBTweQ7/O7d1W+kjLz3DLXeSFTrjgIUFh1eGbrLeyxvJ5s38QtKpFcfbur1Y83Mhk76/5tvq8eqx+XZPgbE0rciesFgFnNrfuS3ke8DSOLSivNtKpjovh4nsab8HaIP2Er0v4Ov2Yknw49fp2wPsJusSr3Rivt1QbDhNbScwquwgiiDpCuCpvDeqx9Nk//uu1VZeBR4vKLzS9JKX3uwE/tngkJpKRGyR1LOAb9u3F0uDe2QSlz+ptAp4vyRVHOoENkbErnIia0ozgR+XHUSjSfoEMDIibqvYN4ZsaG01+TtCKnSSPU5bPYTQUnrLSx9PKE0Fbi8kuOYyAngkIvbTxu3F0uFCJi09T1pU96QtJxuv7gKQNJFskc0vFRdaaWrmRNJPJV0t6ShJR0j6HLA9IrYWH2JxJF1GNll1h6TZ+Wc+2UvO/g1cC7xD0pT8/Glk7+L5VlkxF6GfvHRI+qWkhZI6JB0paSlwR0Q8X2rgDSTpaElfz18E2LPvSmB1ROzId7Vle7G0eGgpEZLeC3w6/3qZpD0R8SvIHi+V9CHgeknnkXX9zouITSWFW4i+ckI2JLCU7M2lm4HlEfH7EsIsjKQFZL1zAi6oOrwiIl4HnpH0YeA6Sc8CJwDntnLPXX95AQ6QPdW2gqzNbAauj4htBYZZhp4FfK+VtJlsKGldRNzac0JEtF17sfR40UgzMzNLloeWzMzMLFkuZMzMzCxZLmTMzMwsWS5kzMzMLFkuZMzMzCxZLmTMzMwsWS5kzMzMLFkuZMzMzCxZLmTMzMwsWS5kzOog6RhJ90jaK+nmXs5ZJumApNX5Z3YdP29mxX32SvrZkIM3M2sBXmvJrA4R8ZKkS4F9wGN9nLo3IuYNw897HJgHIGldvfczM0ude2TM6tcJHAOsLTsQM7N240LGrH5dwJMR8feyAzEzazcuZMzq1wX8bjAXSJou6SZJ+ySNlvQTSS9K2pHPgzlZ0hpJ3ZI2SprcoNjNzJLmQsasDpJGAe8DHh7kpd3AKcA4YBGwHDgLGA+sAuYDi4EPANOBbw5TyGZmLcWTfc3q80GgA1g/mIsi4s+S/gTMjYgf9uyXtBaYFRE35Lv+kp83ZbgCNjNrJe6RMatPF7AlIvYP4dqDNfZ119j3GnDkEO5vZtbyXMiY1aeLwQ8rmZnZMHEhYzZEkiYCpwMPS+qUdGbZMZmZtRsXMmZD97Z8uwm4EHhykNe/CUBSR8W+Izj836XI5uGYmVkVFzJmQ/cEcD9wA7ArIl4Z6IWSzgYuzr8ukTRB0seA84DjJV2VL3+wGHg3MEPS5cMcv5lZ8hQRZcdg1tIkLQMWRMRJw3zfdcBzEbFgOO9rZpYS98iYmZlZslzImJmZWbL8QjyzYhwnaXX+59si4jdDuYmkmRx6y+804LlhiM3MLFmeI2NmZmbJ8tCSmZmZJcuFjJmZmSXLhYyZmZkly4WMmZmZJcuFjJmZmSXLhYyZmZkly4WMmZmZJeu/IbFvpqzTBjQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "W_min = W_sf(SF_opt,ls)\n",
    "beta_SFopt = beta_fun(W_min,ls)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(ls/1e3,beta_cpsf,'-g',label=r'$\\beta_{opt}$' )\n",
    "ax.plot(ls/1e3,beta,'-b',label=r'$\\beta(SF=3.1)$')\n",
    "ax.plot(ls/1e3,beta_SFopt,'-m',label=r'$\\beta_{opt,sf}$')\n",
    "ax.plot([ls[0]/1e3,ls[-1]/1e3],np.full(2,bT),'r',label=r'$\\beta_t$')\n",
    "ax.set_xlabel('$l$ [m]')\n",
    "ax.set_ylabel('SF')\n",
    "ax.set_xlim(ls[0]/1e3,ls[-1]/1e3)\n",
    "ax.set_ylim(3.5,5.5)\n",
    "ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}