coating-delam-mechanics.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import sin,cos,tan,pi,sqrt,exp\n",
    "from scipy import special\n",
    "import scipy.integrate as integrate\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve\n",
    "eps=np.finfo(float).eps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evans-Hutchinson delamination model\n",
    "\n",
    "## Part 1 thermo-mechanical stress in semi-infinite solid\n",
    "\n",
    "$T(y,t)=T_{surface}^{i}-(T_{surface}^{i}-T_{surface}^{f} ) erfc \\left( \\frac{-y}{\\sqrt{\\kappa t}} \\right)$\n",
    "\n",
    "$\\sigma(y,t)=\\frac{E\\alpha (T_{surface}^{i}-T_{surface}^{f} )}{1-\\nu} erfc\\left( \\frac{-y}{\\sqrt{\\kappa t}} \\right)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "y=np.linspace(0,-10,100000)\n",
    "T=lambda t,y: 1-1*(special.erfc(-y/np.sqrt(0.01*t)))\n",
    "sigma=lambda t,y: 120e3*1e-6/(1-0.3)*(1-T(t,y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x737d165310>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8W9WZv59z79UueZH3LbETOyFrE7KRQiA0BAi0ZGg7FJgOUJrSoQvDUtq0M9BCp5NQOrSlUPgxME3oDKUFZkho2Qt0odlJIAuQffG+yvsi6Z7fH7IUL7Itx7a85DyfOLLuPTrnyJJevfc97/m+QkqJQqFQKCYW2mhPQKFQKBTDjzLuCoVCMQFRxl2hUCgmIMq4KxQKxQREGXeFQqGYgCjjrlAoFBMQZdwVCoViAqKMu0KhUExAlHFXKBSKCYgxWgOnpqbK/Pz80RpeoVAoxiW7du2qllKmDdRu1Ix7fn4+O3fuHK3hFQqFYlwihDgRSzsVllEoFIoJiDLuCoVCMQFRxl2hUCgmIMq4KxQKxQREGXeFQqGYgAxo3IUQ/yWEqBRC7OvjvBBCPCyEOCyE+EAIce7wT1OhUCgUgyEWz30DcHk/51cBRZ0/twCPDX1aCoVCoRgKA+a5Syn/LITI76fJauBpGarXt1UIkSSEyJJSlg3THLvxm/vvofbQiHTdAxGHMcKjiB7D9Rxb9Pgtyvle0+3ZUpy+J0SXs+HfBIge9zvbhttFzgjt9DERn7+TYiSRUX4/fSuitgnfNxEyCARDv3e7DSBkO4L207e0IzBH7JnETOfbVorQT/gzZGoCv1UQsGq0OnUak600JRqdn4OBugx9Rgw00qSTOTIVG3poHARVriKOpn2Kv5uXQ36qa8SeWpjh2MSUA5zqcr+481gvCyyEuIWQd8+kSZPOaLC6o2W0dZSf0WMVI4EG6CB0QEd03nY/ZoCwABaEsPT6XQgLYO383QrCjhA2ELbO/hQjghidJTfD34StvR5rRwO2dh+O1gqcLRW4mstxtFYhGNm6zoN91nUu2DFN8PvFgork6EZeRjH+kzv8/Kq0khTTRBOh5/TFvd/lYMVKfvkPCwY77UEzHMY92rON+upIKZ8AngBYuHDhGb2CX9vw5KAf02cR8D6Oy77eXP3MuO9C432MIUP/SWToVoZvO2dg9jgHSNM8PVbn8dA5Gfonzaj9htubpokZDCK73ErTxDSDmKaJDIZuTbN7m8itGSQYCBD0+0M/Af/p+wF/5Higx7FARzv+tjb87U10tLXhb28j0N7e9x+zCxa7A7vLjd3lwuZ2Y3e5cSYk4UxKxpWYhCspGWf4NikJq90RU7+K00Teu7Lz3So7PwGRt1bP86cfYwZl548Z+T0YNDEDkmDAxN8WpKMtQEf4tjVAS30HzfXttDR00FDbRnl9R2QuNqdBRkEi2UWJ5M9JxZvtGrGrw/DnAtMM3UqJ9PsxW1sxm5vxV1TQcew4Cdu24n3jTS7bL0hf+x2811/fZ3+mNGkPtrOjfAe3v3M7j116O/cuvRc6mmFdHtekHuPJutYReT49GQ7jXgzkdbmfC5QOQ7/DRp9vjj6Oq0DDyCNNE3/Y6Le10dHWir+tjfaWZtqam2hraqK9uYm25tO3bU1N+MrLKD34ES0N9VG/nC02O57UNBJS00hITSchLZ2E1LTQsbR0PN5UhKaSxLoS+XwIuoTr4kdHa4C6ihZqS5soP9pA2ZF6tr5Yw9YXj5KQ5mDaogxmLcvGnWwf1nGF6AxFdnk/CKsVzeWC1FSskyfjWryY5C9cg7+ikrJ7/pWK+3+IbGsn5eYvRe1PFzpOzclFeRdxZcGVvHLsFdYuXovV6oKUqUxqK6WxLTCsz6MvhsO4bwa+IYR4FlgC1I9UvF0xcRCahtXuOGNP2wwGaW1soNlXR4uvjuZ6H82+Opp9dTRWV9FQXUnF0cO0NjZ0e5xhsZKUlY03K4fk7Fy8Obl4s3Lw5uRidTiH46kpBonVYZCRn0BGfgIzPpkNQHN9O8fer+bo7kp2vnKcXa+eYMq8VBZ9uoCUbHfc52jJSCfvl7+k5K5vUfnggzjmzsG5cGG/j1kxaQWbjmzi/ar3WZS5CJwpOFsaCY5w2CnMgMZdCPEbYDmQKoQoBr4PWACklI8DLwNXAIeBFqD3V5pCMcxouo4rKRlXUnK/7fxtbTRUV9FYXUl9VSV1ZSXUlZVQeeIoh7ZvQcrTi3tJmVmkT55CWv4U0gumkD55Cq5kr1o0HgVciTZmX5jD7AtzaKhuZd+fS9j35xKO7K7inCWZfPJzhTg81rjOSRgGWT/6EW0HDlB2z71MeWkzwujbhM5JmwPAR7UfhYy7PRGnWT12jLuU8roBzkvg68M2I4ViGLHY7aTk5pGSm9frXDDgx1deTm1ZMTUnT1B54iiVx49ycNu7kTaupGSyp80ge9o5ZE+fQXpBIYbFEs+ncNaTkOrgk58t5NxLJ7PrtRN88NYpju+t4cJrp1G0KCOuc9HdLtK/fTcl37yN+k2bSfrcZ/tsm+pIJdWRyke1H4UO2BNxBJswtTFi3BWKiYpuWCKGv2jR0sjx9pZmqk4co/L4USqOHKL04Ecc2v63zscYZEwpIm/WHCbPmUfWtBnK2McJu9vC+Z8r5Jylmbz19Ee8/tR+Sg75WHZNEboRv3UUzyWXYJsxg9oNG0j87NX9XtlNTZzK8YbjoTv2JBxm09jx3BWKsw2b00XujNnkzpgdOdbsq6P00EeUfvwhJR8fYPum59n2f7/DsNnImzGbSXPmUTBvYdQrBMXwkpLt5nN3n8vWTUfZ/fpJakuauOJrc7G74vMlK4Qg+QtfoPwHP6Bt714cc+f22TbDlcH28u2hO1YXtmALplDGXaEYM7iSkilatDTi4be3NHPqwD5OfLCbE3v3cOzXT/GnXz9FcnYuhYvOo3DheWQVTlOZOSOEpmt88rOFpE3y8OaGA7z4091cdds8nAnxicMnfPpKKtavp/7FF/s37s4MqlqqCJpBdMOGTpBgMD6buJRxVyjOAJvTReHCJRQuXAJAQ3UlR3ft4PDOrez6/f+xY9PzoS+EJZ9kxgUXk1U0XS3MjgBFCzOwOQxeeXwvL/50N5/91rlx8eB1txvXBefT+PY7ZNxzT5+vbaYrk6AMUt1aTYYempeQ4ycVUqE460lITWfeZVcy77IraWtu4tjunRzevoV9b73Bntf+QFJGFjOWLWfGsotJzswe7elOKCbNSuHKr8/lpUfe55XH9/KZ2z6BYRn5nc3uiy6i6c0/0n7wEPbp06K2yXRlAlDeUk6GHrqqMKR/xOcGyrgrFMOO3eVmxgXLmXHBctpbWji0/W98+Je32PLCs2x5/jfkzpjNJ1auomjJJ9ENtRg7HOSe4+WSG2fy+lP7efu/P+KSm2aO+JWS+8KLAGj605/6NO4p9hQA6trqoNO466Yy7grFuMfmdDJ7+SXMXn4JjTXVfPjXd/jgzVf4w8MP4kxMYvbFK/nEJatISEsf7amOe4oWZeCrbGH7S8fILkxi1rKcER3PkpGOtXAqrbt2AV+J2ibBmgBAQ0cDdIZldOW5KxQTC09KKotXf55Fn/ksxz/YzftvvMyOTS+wY9MLTDvvfBat/jwZBVNHe5rjmoWr8ik/Us9ffnuIjIIEUnM9Izqec/58Gl9/A2maURfPE2wh417fXn/ac4+TcVdL+QpFnBGaRsG8Bfzd3few5pEnWfDpv+PYnp3899p/5vkf3cOJvXv6EaJT9IfQBJd8aSY2l8GbGz4c8cwUx7x5BOvr6Th+POp5t8WNQHR67sq4KxRnDQmp6Vz0xZu55ZcbWHb9TVSfPM7z//av/Oaeb3Fy3/ujPb1xicNj5aLrplNT3MTu106O7Fjz5wPQunt31PO6puO2ujs993BYJhiXL29l3BWKMYDN6WLx6s+z5hdPccmar9NYW8NzP/wXnvvh9yg9+NFoT2/cMWVeGoUL0tnx8jHqyptHbBxrQQGa203b/v19tkmwJnTz3K0EMONwYaaMu0IxhjCsVj6xchVf/tkTLL/hK1SdPMFv7vkWLz74Q+rKSkZ7euOKZV+YhmFovPvC4REbQwiBraiItoMH+2yTaEukob2rcfcTjIN1V8ZdoRiDGFYrC65czZpfPMkF197Aqf0fsOGur/On//4v2ltaRnt64wJngpUFq/I5sbeGUx/Vjtg4tmnTaD94qM9QS4I1gfqO0wuqFgKYKiyjUJzdWO0Ollx9DTf/7AlmXngxO3//f/zX7bew9+3XI5W5FH0z91O5eLx23n3+MOYIecu2aUWYDQ0EKiqinvdYPTR2NJ427iKgPHeFQhHClZTMZf/0z/zDjx4iKSOL1x9/mN/e911qS4tHe2pjGsOis/TqqdQUN3HkvcoRGcM+LbSBqb2P0IzDcNAaaFWeu0Kh6JvMqUVce/+PufSfbqP61HGe/vY32fq/vyUYiI9eyXikcEE6yZlOdr58HDkCHrMtVuOuhSQRLASHfQ7RUMZdoRhnCCGYc/GlfOmhx5m68Dze/e2v+e/v3k7FsSOjPbUxidAEC1blU1vazLEPqoe9fz0xET0lhY4TJ6KedxpOWv2njbuGGRdFd2XcFYpxiispmc/c/h1W330PbY0NPPMvd7F90/OYZnw8w/FE0cJ0EtIcIe99BEIi1rw8Ok6einrOYTjoMDsIdhYe1+JUrEMZd4VinFO4cAk3PPgIUxcu5i/PbOC5H/4LDdUjE18er2i6xrmXTqLqZCNlh33D3r9lUh4dp6JvmHIYoSLwrcEOAHRM4rEBWRl3hWIC4PAk8Jk7vstlt95OxdEjPH33Nzm49a+jPa0xxbQlmdicBh+8Pfz7Bax5kwiUlWN2dPQ6FzHunbIDGvHJclLGXaGYIAghmL38Em748S/wZufy0k/X887TT6rF1k4sVp2Z52dzdE8VTXVtw9q3dfIkkBJ/ce/sJafFCUCrGTbuknhEZpRxVygmGEkZmXzhvvXMv/wz7PrDizz3w+/RVFsz2tMaE8y+KAcpJfv+PLzeuyUvVDu342Tv0Ey0sEw8UMZdoZiA6IaFT33pq1xx291UHDvCr9f+M8UH9o32tEadhFQH+XNSOfDX0mFVjLROmgSAvz/jLkPGXRMmMg6uuzLuCsUEZsb5F/EPP3oIm9PFc//2r+x7+43RntKoM+uCbFob/ZzcN3xXM7rXi+Z00hElLBM27i3B02GZeCyoitHSjV64cKHcuXPnoB/XUdaMv7gR27RkjETbCMxMoZh4tDU18dLP1nNy7x4WXfU5ll13Y7fiEn6/n+LiYtrahjcWPRaRUtLsa0c3NBwe67D166+sRBgGhtfb/XjQT1VrFV5bIvbmGupwk5iQjKb1XwbQbreTm5uLxdK9FKMQYpeUcuFA8xl3lZja9lfT8Gbo0sdId2IvSsI+LRlrQSKadeSL4ioU4xG7281n1/6At371ODs2v0BdWQlXfONbWOx2AIqLi/F4POTn54947dGxQGNtG62NHaTmutH04QlgtDscEAxim9q9mlZ7oB3dp5PjyiSpTqdEppCRmYvRz7hSSmpqaiguLqagoOCM5jPujLtnxSQcs1NpO1RH28E6mraV0/RuKRgCW34i9mnJ2KcnY6Q7z4o3qUIRK7phcMmar5OSk8c7Tz/Fb+9by2e/ex/OhETa2trOGsMOYHdbaG3soK05gDNheLx3YVgw29p7HddEyIgPRk9GCEFKSgpVVVVnPJ9xZ9yFEFgyXVgyXXiW5SL9QdqPNdB2sI62Q3XUv3yM+pePoSfbsJ/jDf1MSURYlFevUAghOPeK1SRmZPL7nz7As/d+m8//6w8j584WLFYdw6rT1uwfPuNuMZABP1LKbn/LiHHvzJKJ9a881Ndj3C+oCouOfVoySZ+eQuYdC8hcu5ikqwuxZLpo2VlBza/2U3r/Vqo37KdpaykB38SPKSoUAzF1wRI+968/pKXex2/uuXvUc+F9Ph+//OUvY25/7NgxlixZQlFREV/4whfoiLJ5CGDdunUUFhYyffp0Xnvttcjxm2++melzC/jk8kUE/MOTNSM6Y+Oyx98ybKRP++1KfuCMMJJsuJdkkXrjLLLvXUrql2bhXJiBv7IF34tHKF+/g/Kf7qL+lWO0H6tHBlUhYsXZSe45s/jCD9YjTZOW+jo6RnExdbDG/Tvf+Q533HEHhw4dIjk5maeeeqpXmwMHDvDss8+yf/9+Xn31Vb72ta8RDIZ0d2666SZe/sPLALS3DE/B6ohx93fvT3T66uHkFUF8zHtMxl0IcbkQ4mMhxGEhxNoo5ycJId4WQuwWQnwghLhi+Kc6eIRFwz7dS/LqQjLvXkjGnQtIvKIA3WWh8S8lVP2/Dyj7963UPneQ1gM1SL8SXFKcXaRNLuDa+x9ECI26shI6WltHZR5r167lyJEjzJs3j7vvvrvftlJK3nrrLT7/+c8DcOONN/Liiy/2ardp0yauvfZabDYbBQUFFBYWsn37dgAuvPBC0tJTEQLaW4bnqiVs3Olp3IVACBGX3PauDBhzF0LowKPASqAY2CGE2CylPNCl2b8Cv5NSPiaEmAm8DOSPwHzPGCEElnQnlnQnngtzMdsCtB2so/VADa37q2nZVRH6MpiWjH1WCo5zvGhOy8AdKxTjnKSMTJzV1ei6QV15KcmZ2VgdjrjOYf369ezbt489e/bQ2NjIvHnzorZ75plnSE9PJykpCcMIma/c3FxKSnrvOC0pKeG8886L3I/WTmiCQEeQYMBEN4YWyBCd8+npuUPIez/tucfHyMeyoLoYOCylPAoghHgWWA10Ne4SSOj8PREoHc5JjgSa3cA5Nw3n3DRkwKT9WD2t+2s6jX0NdRrYpiThmJmCfWYKRpLKqVdMXDRNJzk7h7qyEu55YTfHGmW3PPihMjM7ge9/ZlZMbT0eD3v27OnzfLQMkmiLj9H28PRsp3Xeb2/x40wY4mdc10HToht3IeKk4n6aWIx7DtBVqLgYWNKjzQ+A14UQ3wRcwCXDMrs4IQwNe1Ey9qJkkq6air+kqdPQV+PbfAQ2H8GS68YxOxXnnFSMlPh6NQpFPNANg+SsHDStjEBHC4bVNqwGPlYaGxtZtmxZ1HPPPPMMM2bMwOfzEQgEMAyD4uJisrOze7XNzc3l1KnTpitqOwGGVae9JTBk4y6EQBhGrwXV0DCiW8w9HsRi3KPNpedX0HXABinlfwghlgK/FkLMllJ2W4YWQtwC3AIwqVOLYawhNIE1z4M1z0Pi5fn4q1po3V9D2/4aGl49TsOrx7HkuHHMUYZeMfHQDYN/u2YhdWUlBAMBkrNysHZudBpJPB4PjY2Nkd/789wBLr74Yp5//nmuvfZaNm7cyOrVq3u1ueqqq7j++uu58847KS0t5dChQyxevLhXO5vDoLm+HTNoDnlDU5/GvVvMfexkyxQDeV3u59I77PJl4HcAUsotgB1I7dmRlPIJKeVCKeXCtLS0M5txnLGkOUlYnkf61+eRuXYRiVcWIDRBw6vHKX9wJxW/2E3jn04RqFUploqJQcSD13V85SX423tvzBluUlJSOP/885k9e/aAC6oADzzwAA899BCFhYXU1NTw5S9/GYDNmzdz7733AjBr1iyuueYaZs6cyeWXX86jjz6Krof2u1x33XUsXbqUjz/+mOmzp/I/v32ajtahL6z2Zdy1kKJMXAMzA2rLCCEM4CCwAigBdgDXSyn3d2nzCvBbKeUGIcQM4I9Ajuyn8zPVlhkrBOraaN1bTcveavynQh6HJdeNc04ajjmpGN6R93YUiuHiww8/ZMaMGd2OBfx+aktDQlje7FwMy8RMMJBSUlPShMVmkJg2tCvxjtJSzPp67D3+lkd8R7BoFvKaaqmSCSRn5mOJ4Soh2usybNoyUsqAEOIbwGuADvyXlHK/EOJ+YKeUcjNwF/CfQog7CF1z3NSfYZ8IGMl2PBfm4rkwl0BtG637qmn5oIr6V45R/8oxrJM8OOen45iTiu4ePnEihSJeGBYLyVnZ1JWWUFdWgjc7F90Yd5vaB0QIgdVu0N4a6LW7dNB9GQYyGESaZrf1CiEEpoyPjnuYmF4pKeXLhNIbux67t8vvB4Dzh3dq4wfD28PQ762iZXcVvk1H8L10BHtRMs756dhnpihxM8W4wmK1kZSZTV3ZaQOv6RPvPWx1GLQ1+/G3B7Haz/wLLJIOGQx2M+7hsMzpLU0jz8T7Gh5lDK8dz0V5eC7Kw1/eTMueSlp2V1H77McIq4ZjZgqO+enYC5MR+tmj5aEYv1jtdpIysvCVl+KrKCc5KwshJtbmdqsjZAo7WgNDM+7hL75AALqEsYQQmGY4GTI+ZfaUcR9BLJkuEi8vIOHSfDqON4QM/d5qWvZUobksOOam4pyfjjXPc1aJNinGHzank4S0dOorK2ioqiIhLX1CvWc1TWCx6aFF1eQhdBT23HvqyyAGLRw2VJRxjwNCE9imJGKbkkjSVVNp+7iWlj1VNO8op3lLGUaaA+eCDFznpqMPdSOFQjFCODwJBP1+mupq0S0W3MnegR80jrA6DJp97QSDJvoZpkSKvoy76MxzFyJejrsy7vFGGBqOWak4ZqVitgVo3VtN866KUA79a8exT0vGuSADx8wUxBC3QysUw40r2Usw4KeptgbdMHB4EgZ+0DjBajdoph1/WxDdNTTjTg/jfjrmHj/5AWU9RhHNbuBalEn6P32CjG8txLM8FKevfeYjSn+0jbpNh+koboy6jVqhGA2EECSkpWN1OGioqqSjbXiExkZC8rempoaLL74Yt9vNN77xjQH7NKwaQhN0tA0h313TQAhksLsIYcRzj1tQRhn3MYMl1UHiZflkfmcxqTfPxj4tmeYd5VQ+sofKn79H41+KCTZF16xWKOKJEBpJGVlohoGvooxgYOiSuSMh+Wu32/nhD3/IT37yk5j6FEJgten4285cHbYvCYKw/EA8zbsy7mMMoQns05JJue4csv/lPJL+rhAsOvV/OEbZv2+n5n8+pO1QHdJU3rxi9NB0naTMLKQp8VWUI82h5XCPhOSvy+XiggsuwD4I+QSL3SAYMAkOoYBHVOMuwguq8VJzVzH3MY3mMHCfl4X7vCz8Fc0076ig5b0KWvdWo3vtuBZn4lqQgT6MFdwVilixWG0kpmfgKy+joXpoGTQjIfl7JljtoVTGjvYADsuZfa6EYUQt2DEWhcMUYwBLhoukT08h8bJ8WvdX07y9PLQI+/oJHLNScC3OxDY1CaFNnPQ0xSjxyloo3xtTUzuQFvATDPgxLVZ0vQ+TkjkHVq2Pqc/hkvw9E3SLhqYLOlqDONxn2omO7FHVKlJqLzLNkffelXEfZwiLhnNeOs556firWmjeXk7LLuXNK0YPzTAwTZOgvwMhBJo2tB2swyX5eyYIIbDYDPztZy5FIHS914JquEi2DPnwwzLXgVDGfRxjSXOSdOUUEi+N7s27l2ZhLUicUJtNFHEgRg87jAD0YJDakmJMGSQlZ9KgNWhGQvL3TLHaddpb/AQDJoblDL6odANMs5u+TKSO6rDNcmDUguoEIOzNp90yl4y7FuD+ZDbtR3xUPbGXip+9R9O2MswOVR9WMXJouk5iRiYyaFJfWTHo9N2RkPwFyM/P584772TDhg3k5uZy4MCBvrqMYLGFDLq//cw+M8IIPb6r9x52sMw4assMKPk7Uox3yd+xjtkRpPX9Kpr+Voq/rBlh13EtyMC9NBsjVRUYUXQnmrTsmdDSUE9DVSXuZC9ub8owzCz+SCmpLm7C5jBIOIPPSrC+no5Tp7AVFqJ1ZurUtdVR2lTKlKCgI6DjyCjEagx8VTCikr+K8Ylm1XEtysS5MIOOk400/a2Upi1lNL1bin16Mq6l2dinJasFWMWw4vAk4G9rpamuFovdgc3pHO0pDZpQ3F0/Y8+dTvEwGejtucczz10Z9wmOEALb5ARskxMIXtlB8/YymraV07ZhP7rXjntpFq4FGWjOiVmIQRFfhBB4UtPxt7dTX1lOSu7g4+9jgbCI2JnozESUIYOnc91Px9zjF3tXMfezCD3BSsIlk8lauwjv9eegJ1hDm6PWbadu02EC1cOzlVxxdqNpGkkZmUh5ZvH3sUA47h44E+9djxJz7+avq2wZxQghdA3n3DScc9PoKG2i6d1SmreX07y1DPuMFDwX5GAtSFBZNoozxrDa8KSk0VBVSUu9D1fSUHR044/FqoMQ+NuC2AZ5VSuiGXcR/2wZZdzPcqzZbrx/P43Ey/Jp2lpK89Yyqg7UYMlx41mWg2NOKmKIFeEVZycOTwIdLc001dZgdTix2MaPnLXQBBarRseZeO6d4mFE8dyVtowi7ugJVhIvzSdz7WKSri5EdgSpffZjyh/YQcM7pzBbhi4OpTi7EELgSUtH6Br1leWYQ9SfiTcWm06gIzhoHSchRNSNTNBlh2ocXHhl3BXd0Kw67iVZZNyxgJSbZmGkO2l49fjpuHyNissrYkfXDRLTMgh0dNBUW9Nnu8GqQj7yyCMUFhYihKC6urrPdhs3bqSoqIiioiI2btw4qLlH8t3PYI+I0PXunvsohGWUcVdERWgCxzle0tbMIf22+TjmpNK8vZzyn+yk5pkP6ShpGu0pKsYJNqcLZ2ISLfU+2luao7YZrHE///zzefPNN5k8eXKfbWpra7nvvvvYtm0b27dv57777qOuri7mMYa0mUnvrgypdqgqxiTWbDfea6aT9Z3FeC7Mpe3jOip/sZuqp/bSdtg3LrMhFPHF7U3BsFqpr6rAjBKuGIzkL8D8+fPJz8/vt81rr73GypUr8Xq9JCcns3LlSl599dWY56zpGrqhETgTz93Qoy6oQvwqMakFVUXM6AlWElcV4FmeR9O2Mpr+WkL1k3ux5LpJWJ6HfWaK2hSliIqmaSSmZ1JbcorGmioS0zO7nR+M5O/MmTNjGrOkpIS8vLzI/TORBjasZ7iZqWdYZhQ8d2XcFYNGcxgkLM/Dc34Oze9V0PinYmr++0OMVAeei3Jxzk9X9V/HMQ9sf4CPaj8a1j7P8Z7DdxZ/B1dSMk11tdhcbuyu6Jq6sQiHxUK0K8rBpvdabFpIRGyQm5nCC6phZUll3BWuZOY5AAAgAElEQVTjCmHRcC/JwrUwk9Z91TT+6RR1Lxyi/o0TeC7IwbUkE82m3mKK07iSk2lrbqahugqr3YGm99ZXGUjyN1bPPTc3l3feeSdyv7i4mOXLlw9qvob19GYm3Tk4446UYJqg610WVEOVmOJh5NUnTzFkhC5wfiINx9xU2g/5aPzTKepfPkbD26fwnJ+N+/wcNId6q40XvrP4OyPWtxAaienp1PQIzwxW8jcWLrvsMr73ve9FFlFff/111q1bN6g+wsbd32FiG4xMTpeNTKLLF5iMY9RSXTsrhg0hQvVf074yl/Svz8OWn0DDmycpW7+d+tePE2xWufIKsNjsuJO8tDY20tYcyp4ZrOTvww8/TG5uLsXFxcydO5c1a9YAsHPnzsjvXq+Xe+65h0WLFrFo0SLuvfdevF7voOaqaQLDog1ahiBi0Dtz+0djE5OS/FWMKB0lTTS+dZLW/TUIq457aRbuZTnoblUpaiwxXJK/sSKlSU1xMdIMkpI7KWp4ZqzQUN1Ke2uA1Fx3zDH7YFMTHcePYy0oQHe5CJgBPq79mHSp4+ww0dOnY4+hEMhQJH+V564YUaw5blL+cSYZt5+LfYaXxj8XU/7ADny/P0qwoWO0p6cYJcLhmWAwQGNN35uQxgIWm440JcFA7I5wuAJTOGNGLagqJiyWTBcp152Df8UkGt8+RdPfSmjaWop7cRbui3IxEseP7ohieLDY7LgSk2n21eHwJGB1jM0iMpFF1Y4ghiVGfzgccw+HZdQOVcVEx5LuxPuF6WTetRDnvHSatpZR/uMd1L14mEB9+2hPTxFnXMledIuFhupKpByb2jOGNSQENph899Oa7r099zElHCaEuFwI8bEQ4rAQYm0fba4RQhwQQuwXQjwzvNNUTDSMFAfez08j81sLcS3MoHlHOeUP7sD30hGCjSpcc7agaRoJqWkEOjpo9vlGezpREaJzUXUwO1U7wzI9xcPGVFhGCKEDjwIrgWJghxBis5TyQJc2RcB3gfOllHVCiPSRmrBiYmF47SRfXYTnojwa3jpJ05aQtrxraTaei3LRXapC1ETH5nRhd7tprqvF7nZjWMbeYrvFqtPW4o9sShoIEZb97RKWEULE1brH4rkvBg5LKY9KKTuAZ4HVPdp8BXhUSlkHIKWsHN5pKiY6htce8uTvXIhjdipNfwktvNa/dlzJDZ8FeFLSQAgaqqrGpFaRYdWQpsQczKJqD9lfEVGVic/zi8W45wCnutwv7jzWlWnANCHEu0KIrUKIy6N1JIS4RQixUwixs6qq6sxmrJjQGKkOvF+YTsYdC7Cfk0zj26co+/EOGv54ErMtMHAHinFJY1MTv/nf/6OjtYW2psYB2/cl+Sul5LbbbqOwsJC5c+fy3nvvRX38rl27mDNnDoWFhdx2220DfqFENjP5BxOa0SOee2R+jK2Ye7S59PxLGEARsBy4DnhSCJHU60FSPiGlXCilXJiWljbYuSrOIizpTlKun0H6bfOxTUmi4Y0TlP+4s3DIGaj0KcY2Pp+PJ3+1AYvdTmNNdVTlyK70Jfn7yiuvcOjQIQ4dOsQTTzzBrbfeGvXxt956K0888USk7UBqkeEsmUBH7Iu+Qtd6KUOGDedYKZBdDOR1uZ8LlEZps0lK6ZdSHgM+JmTsFYohYc12k3rDTNK/MQ9rnoeGV49T/uMdNG0pRQbGZnaFYvCEJX8/dcWn+cGP/p2mutp+2/cl+btp0yZuuOEGhBCcd955+Hw+ysrKurUpKyujoaGBpUuXIoTghhtu4MUXX+x3PKEJ9MEuqvahDBkvYslz3wEUCSEKgBLgWuD6Hm1eJOSxbxBCpBIK0xwdzokqzm6suR5SvzSb9hMN1L96HN+mIzT+pYTESyfjmJumpIbHOWHJ3/fff5+SY0f55IUXolusvRYvBxIO60vmNysrq1ub3NzcXm0GwmLV6WgbRDqkpmN2nF4vEgik6PTZ4+C6D2jcpZQBIcQ3gNcAHfgvKeV+IcT9wE4p5ebOc5cKIQ4AQeBuKWXfNbUUijPENjmBtFvm0HawjoZXj1P77MdY/lxM4uUF2IqSBi3pquhN+b//O+0fDq/kr23GOWR+73sxtc2aNJm3Xv4DhtVKclbOoF7TWGR+z1QK2LBqtDX7MYMmWizyv7oGZs+wTPwWi2PaoSqlfBl4ucexe7v8LoE7O38UihFFCIFjuhd7UTKt71dR/8YJqv9rH7YpiSSuKsCa5xntKSqGQHNLC5dctRozEEAzLGjaaUM6kOeem5vLqVOn8z+Ki4vJzs7u1aa4uLjfNtEwLOGdqiZWx8DGva9smXiZeCU/oBi3CE3gnJ8equ+6rYyGt05R+egeHLNTSLg0H0v6YDRaFWFi9bCHk56Sv++//wG1Jacwg0FS8iZ3M/D9cdVVV/HII49w7bXXsm3bNhITE7uFZACysrLweDxs3bqVJUuW8PTTT/PNb35zwL4Na2gO/o4g1lgkrDUNpESaZmfee0xPYdhQ8gOKcY8wNNzn55D57YUkXDKJtoM+Kn62i7oXDhFUkgbjgp6Sv0IIPKlpBAMBmqMsrvYl+XvFFVcwZcoUCgsL+cpXvtKt6HbX0n2PPfYYa9asobCwkKlTp7Jq1aoB56jpGpquxZwxE03297THPvK+u5L8VUw4gk0dNL51iqZtZSAEnmU5eJbnqqpQ/RBvyd9Yqa+soK2pkZTcSRjW0d+56qtsIRgwScmOXiKwKwGfD39xMbaiIjSbjaP1RxH+drLag8i0c3BYB34/KslfhaILuttK0lVTybxrIY7ZKTS+fYryB3fStLUMGRx7ux8VfeP2poAQY0YW2GLVCfpNTHPg91E02V+lCqlQDAOG107KteeQ/vV5GKkOfC8epuLnu2j9qHZMbnFX9EY3DNzJXtpbmmlvaR7t6UTi7sFYdqr2lP3tsqAaD5RxV0x4rHke0r46l5R/nAEm1GzYT/VT++gobRrtqSliwJmQiG6x0FhTM+pfynp4p6p/4Lh7L9nfOKfpKuOuOCsQQuCYlUrG7eeS9Jkp+EubqPzFbmqfO6gWXcc4QtPweFMIdLTT2tgwqnPRDQ0hRGyLqj1kf+MdllErTIqzinBmjfPcDBrePkXTuyW0flCFe1kOnovy0Gxjt5bn2YzN5cZit9NUV4Pd7Yk5NXK4ESIkQxAcjOfeLSwTP11I5bkrzko0h0HSFQVk3rUwVNv1rVOUP7iD5u3lyBgWyxTxRQiBJyUNMxCkxVc3qnMxLBqBWGLuPQt2iLGnCqlQTFgMr52U62eQ9rVPYKQ4qPvfQ1Q+uof24/WjPbWzCp/P1y0nPRpWux2720Ozr46Hf/7zQUv+bty4kaKiIoqKiti4cWPUMWpra1m5ciVFRUWsXLmSurreXySGVcMMSsxg/967ECKUMdPFc48nyrgrFIBtUgJp/zQX73XTMZs6qHr8A2p+8xEBX9toT+2sIBbjDp2pkcC82TMHJflbW1vLfffdx7Zt29i+fTv33XdfVMO9fv16VqxYwaFDh1ixYgXr16/v1SYiQxBDaAZN71YkW6VCKhSjgBAC5yfSybhrIZ4Vk2jdX0PFf+yi4c0TSkN+hAlL/s6bN4+77767z3aGxYIzMYlp+fnkZGX2Ot+X5O9rr73GypUr8Xq9JCcns3Llyqga7ps2beLGG28E4MYbb4wqBawPRttd1yDYc4dqfEy8WlBVKHqgWXUSV07GtTCD+leO0fDmSZp3VpB4RQGOOalKeXIECEv+7tmzh8bGxm5SAV155plnOGf6dFobG6JubOpL8rev4z2pqKiIaNFkZWVRWdm7YqimC4QmYoq7C01Dml2zZWTcgjPKuCsUfWAkh+Lx7efV43vpCLXPfIQ1P4Gkz0zFmjPw9vPxyl9+d5DqU8O7ByA1z82ya6bF1Nbj8bBnz55+27iTvTRU96632pec75nK/EZDCIERY8YMXWLuSjhMoRhj2KYkkv7N+SRdXUigqoXKR3ZT97+HCDZ1jPbUJiRhzz3az4EDBwBwJCSgWyxI0+xmuPuS/I1FChggIyMjUrmprKyM9PT0qHM0LDqBDnPATVVC16OEZeKD8twVihgQmsC9JAvn3DQa/niSpr+V0vJ+FQkrJ+Nemo3QJ06oJlYPezjpKfk7kOcuhIY72YuUkrbmZuisydyX5O9ll13G9773vcgi6uuvv866det69XvVVVexceNG1q5dy8aNG1m9enXU8Q2rhmySmEGJbvTz2vcIy0SIg5VXnrtCMQg0h0HSp6eQcfu5WCcnUP/7o1T+4j3aj6rUyaHQU/J3IB5++GGKZsykrLycRUuW8OUvfxnoW/LX6/Vyzz33sGjRIhYtWsS9996L1+sFYM2aNYQVateuXcsbb7xBUVERb7zxBmvXro06/ulF1f7j7l1TISFeS6mdYyvJX4XizJBS0nagBt9LRwn62nHOTyfxigJ0z+hL0w6WsSr5OxDtzc3UlZeSkJaOMyExbuOaQZPq4iZcSTZcibY+2/nLywnU1OCYNYvKlkqqWqqY3h6gPXUWrhgkqIci+avCMgrFGRLWq7EVJdP49ika/1xM64GaCRmqGatYnU6sdkfcZQnChTsGXFTV9Eg1pjDxcqdVWEahGCKaVSfxsnwy7ljQPVRzTIVqRhohBG5vSkiWoN4X17ENy8BVmUS4kLZpqh2qCsV4xZLqIPVLs0j54gzMtiBV/+8Dan/7McFGlVUzklgdDmxOFy31dZjB+G02M6wagcAAGTNhfRnTPJ16GScbr4y7QjGMCCFwzE4l484FeC7Oo+WDKsp/spPGd0tUFagRxO1NwQyaNMdRVEy3hApgBwN9e++RakwqLKNQTAwioZrbz8U6yUP9S0ep/MVuJUg2QlhsNhweDy31PoKBQFzGNCzhqkz9hGa0zmpMwaAKyygUEwlLmpPUm2fj/YcZmK0Bqh7/gLoXDmG2+Ed7ahMOV7IXiYyb964bMVRlihJzV567QjFBEELgnJNKxl0LcF+YQ/Oucsr/YxfN71WMetm4sUKsqpBhHnnkkV6Sv4bFit3l4c677hqU5O+uXbuYM2cOhYWF3HbbbVFfk2hSwqGMGdGv594tLKPkBxSKiYlm1Um6Ygrp35iP4bVT97uDVD+5F39Vy2hPbdQZrHE///zzo0r+/mXbdo4eP8Guv70bs+TvrbfeyhNPPBGRCo6mFtmXlLBuaP3G3LsW7DgdllEFshWKCYk1203arZ8g6e8K6ShpouJn79Hw5glkf0ZighOr5G+Y+fPnk5+f3+v47//wB754/XW0NTWyaMGCASV/y8rKaGhoYOnSpQghuOGGG6LK/PYlJWxYdAL+vjNmunru8Q7LqE1MCsUoIDSB+7wsHLNS8P3+KA1vnqRlTxVJVxdin5o02tOLO4OR/J05c2af/ZSUlDD1uusAaPbVDSj5W1JSQm5ubq/j0fqN9viZ01KQZj8aM511VOUohGWUcVcoRhHdYyXlunNoW5BB3YuHqf7PvSEZgysL0N2jI2Pw9oYnqDxxdFj7TJ88hYtvuiWmtrEIh/WFlBLdMHB4EmhtbEBK2a/kb6xSwH21i2TMBMzIAmuvvoRQm5gUirMV+7RkMu8493Ru/H/sonnH2VmsOxbJ374IS/u6kpMBOHXyZL+Sv7m5uRQXF/c63le/PdtFBMT6XVTVkcHuYZl4vKoxee5CiMuBnwM68KSUsndhwVC7zwPPAYuklEoVTKEYBMISyo13zkuj7v8OU/fCIZp3VZB8dSGWDFfc5hGrhz2cDFbyty+6Sv7uO3gIt8tJWmpKn5K/Xq8Xj8fD1q1bWbJkCU8//TTf/OY3++23q5Rw+Mqg31x3XQOzy87ZsbJDVQihA48Cq4CZwHVCiF5BLyGEB7gN2Dbck1QoziYsGS7SbplL8ueKCFS2UPHw7gm/4Homkr9hr3vu3LmsWbMG6C75e/u37mb9fffRXFfXr+TvY489xpo1aygsLGTq1KmsWrUKgMcff5zHH3+8V79dpYSFEOgDVGUKy/6Gwz3xuhYbUPJXCLEU+IGU8rLO+98FkFKu69HuZ8CbwLeAbw3kuSvJX4ViYIJNHfheOkrr+1UYGU6SP1eEbVLCsI8zXiV/B6KhuoqWBh+puZMxrCOzhtFQ3UpHW5DU3OilF9uPHgUh6MhJ5WTDSQr8fkzvbNwjLPkbS8w9BzjV5X5x57Gug80H8qSUv4+hP4VCESO6O7TgmnLjTGRbgKrH3se3+Qhme/wEssYzrqRkBIJmX+2IjaFbNMygidnX+khnqb1uqZBx2LwWS8w9WoQoMjMhhAb8FLhpwI6EuAW4BWDSpEmxzVChUOCYkYKtIJH6147TtKWU1gM1JF9diH26d7SnNqbRDQNnQiLNDT5cySkYFsuwj9FVY0az6b3OC03DNMemtkwxkNflfi5Q2uW+B5gNvCOEOA6cB2wWQvS6bJBSPiGlXCilXJjWWfNQoVDEhmY3SF5dSNpX5yIsGtW/2k/tsx8RbFY6Nf3h7PTeW0ZIc+Z0xkwfV1PhUnudtn0sacvsAIqEEAVCCCtwLbA5fFJKWS+lTJVS5ksp84GtwFUqW0ahGBls+Ylk/PO5eD6VR8sH1VQ8tJOWPZVKp6YPuua9BwPD/0UYzm/va1FVaFpIz32see5SygDwDeA14EPgd1LK/UKI+4UQV430BBUKRW+EoZF4aT4Zt83H8DqoffZjajbsJ1DXNtpTG5M4k0J5782+4a/WNGDGTDhbpvOu7PL/SBJTnruU8mXg5R7H7u2j7fKhT0uhUMSCJdNF2q2foOlvpTS8dpyKn+4i8bJ8XEuzEZqq4RrGsFiwuz20NtTjSkpGN4Z3c35/JfdEpwQBnQuuUsTHh1c7VBWKcY7QBJ4LckI1XPMT8b10lKrH3x9XapPDIfkL0aV5w7zw0u9Z+qlLmDZtWkySv7W1taxcuZKioiJWrlwZ2QDVk40bN7LgvLksuuATbNiwoXeDiHhYfMNmyrgrFBMEw2sn9UuzSL5mGv6qVip+/h6NfyoeFxIGwyX525c0b21tLf/2ox/xx1df4eX/fT4myd/169ezYsUKDh06xIoVK1i/vvfG/LCU8J/e/iuvbnqL+++/v/eXQKdxF7JLQCYOL4ky7grFBEIIgevcDDLvXIB9mpf6V45R+dj7+CuaR3tq/TJckr99SfOGJX9zp0wh0eNh+YXLBpT83bRpEzfeeCMAN954Y1Qp4HC/6RmpJCUm86mLV/TSgw/L/orOOqpjKuauUCjGF7rHSso/zqD1gyp8m45Q8fBuElZOxrMsF6GPvVj8cEr+9iXtm5eXh8Vqw+5yk+71UnzqVL+SvxUVFWRlZQGQlZVFZWVln+OFM2ays7J7Swb3iLnHi3Fn3FtbWwFwOByjPBOFYmwjhMD5iXRsU5LwbTpMw6vHad1Xjffz07Bk9i1E5nvpCB2lw+vpW7NdJH1makxthyr525Oe0r6u5GRMKfG3t8cs+TvQeJquITSBGZS9Hh/a5wnIrp77yDPujPuePXt47bXX8Hq95ObmkpOTQ05ODpmZmRjDvAKuUEwEdI+VlC/OpOWDKnybDlPxi90krJiE56JchD72IrONjY0sW7Ys6rmBPPf+pH3feecdACw2O5XV1UwpyCc7O7tPyd+MjAzKysrIysqirKyM9PT0qOOF+zUsGsUlxcya20OjR+++oKqMex8UFBSwYsUKiouLOXr0KB988AEAmqaRmZlJTk5OxOh7vV40bey9eRWK0cA5Nw3blER8m4/Q8PoJWvdVk/z307FmdffiY/Wwh5ORkPztKs3bU/L3nb/8le/c/s8kuZx9Sv5eddVVbNy4kbVr17Jx40ZWr17da7yu/TY0tfL2O2/x0M9+0r1RpycvuoVlVMy9F5mZmWRmZgKhS6KGhoZITK2kpIQ9e/awY8cOAOx2O9nZ2RHvPjc3F7c7unKbQnE2oLutpFw/g5Y51fg2Habykd0kXJwH2aObUdNV8nfVqlU8+OCD/bZ/+OGH+fGPf0x5eTlz587liiuu4Mknn+SKK67g5ZdfprCwEKfTya9+9SuAbpK/AN///vfJyMqipd7HL3/5KF/60s20trayatWqiOTv2rVrueaaa3jqqaeYNGkSzz33HAA7d+7k8ccf58knn+zWrzQld37z2yR1bpgK07WOKmIMSf6OFCMl+WuaJlVVVd0MfkVFRSQ2lpiYGDH2OTk5ZGdnYx0hKVCFYiwTbPbje+kIrXuq8F2dwMz5s9GsvYWvJirtLc3UlZWSmJaBI2HoMsptzX4aqltJznJh6fJ3lKZJ24ED6OlpHNKqyQgEsCdMx+20D9jnUCR/x53nPhCappGRkUFGRgbnnnsuAB0dHZSVlXUz+OFyXUII0tPTuxn8tLQ0dP3seZMrzk50l4WUa8+hdU4qvsZiApUtaAlWdI91UIuK4xWrw4lhs9FcX4fd4xnyc9a7qEN2Ne5C0yJ1VNHoFBBTYZlhwWq1Mnny5G4bHpqamigtLaW4uDhi7MO72SwWC1lZWd0WbBMTE8+KN7zi7MMxKxVtfyWa04LZ0IFsDaB77WiWie3gCCFwJSVTX1FOe0szdtfQQrYRAbEoFbNC1ZjCC6rxsSNnhXGPhtvtZtq0aUybNg0Ixe9ra2spKSmJGPxt27YRDIZkPF0uV7fYfXZ2tkrHVEwYhCYwvHZMh06grp1AZQt6gg3NbZnQTo3d5abJYqHZV4fN6RrSc9U0gab3ISAmOmV/OxkzBbLPBoQQpKSkkJKSwty5cwEIBAJUVFREQjnFxcUcPHgw8piUlJRu4RyVjqkY72gOCxarTtDXTrC+HbM1gJFsR1gmZtaZEAJXYhIN1VX429qwDtFh0y0agWjGvZvnHp8a2coS9YNhGBHDHaa1tZXS0tKIwT9y5EgkHVPXdTIzM8nNzY38JCUlTWjPRzHxELqG7rUjWgMEfe34K1vQE61oronpxds9CTTV1dLsqxuycTcMjbaW3prxQhNItYlpbONwOJg6dSpTp4ZygaWU1NfXd/Pud+3axbZt24BQOCccuw/f2my20XwKCsWACCHQnRY0q07A107Q18WLNyaWF69pGs7EJJpqa/C3t2MZwudTtwikKTGDJlrXDWKa1llHNczIm/iJ9SqNAkIIkpKSmDVrFpdeeik333wz3/3ud/nqV7/KlVdeSWFhIdXV1bz11ls8/fTTrFu3jl/+8pds3ryZXbt2UVFRgWn2IfKvUIwywtAwUuzoSTZkh4m/ooVgs3/Yqz7FQ/J348aNFBUVUVRU1Evyd+mFF7F0xSV84+tfH1Dyt78xui6qdpUSvvP++5FmEBBxU4VESjkqPwsWLJBnE83NzfLgwYPy7bfflr/+9a/lunXr5Pe//335/e9/X/7oRz+SGzZskG+++ab86KOPZGNj42hPV3GWceDAgQHbmP6g7Kholu2nGmRHVYs0A8FhG//YsWNy1qxZMbd/77335LFjx+TkyZNlVVVV5Pgf/vAHefnll0vTNOWWLVvk4sWLpZRS1tTUyIKCAllTUyNra2tlQUGBrK2tlVJKuWjRIvm3v/1N1ldWyIsvXCZf2rxZSinl3XffLdetWyellHLdunXy29/+dr9jSCmlvyMgK47Xy5bG9ki/pmnKSy9aLl/8z/+UB6r2y9Ly92VjU3NMzzPa6wLslDHYWBWWiRNOpzPiNUDoS7WmpiaSmVNcXMxf//rXiNeQlJTULXavFmsVo40wNIw0B2aTn2BDO/6KIHqSDd1pGXLfXSV/V65cOeAO1fnz50c93pfk7zvvvMPKlSvxer0ArFy5kldffZXly5dHJH+DAT/XXH01Lzz/HJ/+zGfYtGlTRDfmxhtvZPny5TzwwAN9jpGVlRXx3EuKSyL9Anzx85/jpTfe4K6rPwmobJkJjRCC1NRUUlNTI/Km4c1WxcXFFBcXc/LkSfbt2weoxVrF2EAIEdrkZNcJ1rYTrG1DtgXQE21DEiGLl+RvtONhyV/dsDB5yhRe2LyZYDDQp+RvX31lZWWF/j6GRvGpHlLCOTmUVlQgOsMyQm1iOruIttmqoaEhYuzVYq0iHrzyyiuUl5fH1jgokYGQZoowNOijbmtmZmZEs2UgRlryt7/jDrcHgaC1vn7QY4SJmg7ZwwlTnruChIQEZs6cGfFWgsEglZWV3Qz+xx9/HGmfnp7ezeCnpaUpZUzFyKELhKYjAybSb4buDzGbZqQlf8PHly9fTm5ubjfJ37KKCrJzsmlpqO9T8revMcLohkZmWlZ3KeHycrLS0hASpFB57ooo6LpOVlYWWVlZEYW71tbWSNy+uLi4m5SC1WqNGPqw0VfKmIr+iNXD7oo0JcGGDsymDoQRypMfjAhZvCV/X3/9ddatW4fX6+0l+ftPX70FMxhk1WWXRpX87WuMMIZFIz09A4/7dL//89xzfPVzn0OTIfmBeOg1KuM+AXA4HBQWFlJYWAicXqztavDffffdSMqlWqxVDDdCExhJNkx7V/kCK1qMImTxlvy99957I4urjz32GDfddFNE8veq1X9HbWkxt375Zr5257d6Sf72NQbAvHnz2L4lpHb78M9+wZo1a2htbeWyiy/msmXLKOs06vGIuU84yV9FdLou1oaNfkNDA3D6aiAvL4/c3Fzy8vJIGAYJVMX4IZq07JkiTUnQ14bZEkBYdQzv+Nv41NbUhK+ijKSMLOyDvNINBkxqSprweO04PCE58YDPh7+4mPI0A0P4SXIW4ImhXyX5qxiQgRZrT506xfbt29myZQsQivWHDX1ubi5ZWVnKu1fEREiEzEHQ7g/JF1S0oCdZ0ZzjR77A5nKhWywhOeBBGndNFwghui2qhgt2aPESlkEZ97Oanou1gUCA8vLyiLEPx+9BefeKwRORL6hrI1jXjtkaxEgeWspkvOgqKNbR1orVHrvmTCgdUnSX/u007kIqbRnFKGAYRiQOf9555wEDe/ddjdlVe8cAABqASURBVL2K3St6IgwNI7X7xicj2Y7mGPvvk5CgWA0tPh/WzMEJiumGRiAQ3XMPiYirPHfFKBOLd79//35AefeK6HTf+NRGoKYVzWUJbXzqIy9+LKBpGg5PIs2+OgJ+P4Yl9p24ukWjvS2IlDIUiuriuccLZdwVg0J594ozRbPoiHRnKGWysQPZHkRPtqPZxm7FJ2diEi31PlrqfSSkpsX8ON3QQErMoEQ3Tht3TYVlFOMJ5d0rYkUIgZHYmTJZ206gamzXbdUNA7vbTWtjA+5kL1qMtZW7qkPqhhYJy0Ri7nGw8GN/ZUMx7gh79+eddx5///d/zx133MGdd97JNddcw+LFixFCsH37dp577jkeeughHnroIZ577jm2bNlCcXExgUBgtJ+CYoTRbAaWDGekbmv14TIe/cWjMT9+uCV/w9K8t912Wy/J34WfvIBr/vEGSk6eiHmMWXNn8Nvnnzm9qNrFc/f56lm9+qpeUsLDTizSkcDlwMfAYWBtlPN3AgeAD4A/ApMH6vNsk/xVdMfv98tTp07JLVu2yN/97nfyoYceikgg33///fLJJ5+Ur776qty/f7+sr68f7elOeGKR/B0pAs0d8uOt++TM6TNkoLFdmqY54GNGQvLXNE15+eWXy5dffllK2V3y957vrpXf+Oot0jSDMY1RU1MjJ+VNliePlkoppTRNU7bs3Ssrj38ov/L1m+QPfvB9KWV3KeFojKjkrxBCBx4FVgLFwA4hxGYp5YEuzXYDC6WULUKIW4EfA18Yri8gxcSjv9h9OJSjYvdnB7rTwr0P3c/RE8c4d8lCLrl4BT/5+X/0mzI5EpK/ADfccAMvvvgiq1at6ib5e/OX13DJypW0NTXFPMbyCz/F66+/zpe/eiNCCISmoUnJm6++zeuvvA50lxIebmL5dCwGDkspjwIIIZ4FVhPy1AGQUr7dpf1W4IvDOUnF2cFgY/dhQaiw0Vex+/HL+gceYN/+/by3ZSe+0mrmzZ2HMEQvNcWRlPztehzoJvk7ecoUqmtqaPb5Yh8j53RfAGgawpRUV9WSkZkJdJcSHm5iMe45wKku94uBJf20/zLwSrQTQohbgFsAJk2aFOMUFWcrg/XuExMTycvLi/xkZGSgx7gApjjNwYM/pLHpw2Ht0+OewbRp9wzYTvdYSS7IYOcftyADJprbip4Y+2KrHKLkb/h4tGNCCAId7QSjrAlF60voAmnSLR1SyFCW+1jRc4/2V406MyHEF4GFwEXRzkv5/9s77+A4yzuPf37vu011V12yZLlINjayje3QQgnEgWAIwZAEsC+X4RKnUi65OSaXHDepcxlfkpk7SJuheFKoIXGCATtAgNiMKTZxoQTbsmWQ5CbbklaypNW25/54V6tdaSWtZe2utHo+Mx62vOX3Cun3PM/3+RV1P3A/WLVlkrRRo4ky1uw+tsGJ3W6nuro66uxramrIzc3NpPmaJOjp7+XylZejQgrCyqoVbxogqS/5G1u+d1jJ34oKDNOkoqw0qXscPXqYC5ZdggorxLRkGVEhSstKOHrsGPMK3HGlhCeaZJx7KzAz5n0NcGToQSJyFXAPcIVSqn9izNNoRifR7N7r9dLS0hL9F9u+sLS0NG52X1JSouvdDyGZGfZEM1LJ33BfkGCHDxSYbgdG3uiJRBNR8veuu+6KXmtoyd9ct5urrvgIv/n1r8e8x4svv8g3v/FfhIJhDNOwZu4hxYprruTRRx/ju9/9Xlwp4QlnrB1XrAGgCZgDOIA9QMOQY5YBB4F5yeziKh0to0kj/f396tChQ2rr1q3qkUceUevWrYtG5qxbt049/PDDasuWLaqpqUn19/dn2tyMkMlomQHWrFmjGhoa1N133x33eTgYUv62+Mbc9957r6qurlamaaqqqiq1du1a69hwWN1+++1q7ty5atGiRWrHjh3R6zz00EOqrq5O1dXVqfXr10c/37Fjh2poaFBz585Vd9xxRzRa5+TJk2rFihWqvr5erVixQp06dUoFgwF19MB+9cXP/8uY93jggQfV8fe9qq/br9auXau2PfWU6tr3rnp97yvqiis+EnfdkTibaJmkSv6KyHXA/wEmsF4p9d8i8oPITTaKyF+BxcDRyCnNSqkbRrumLvmryRQqUu8+dnZ/4sQJwNJOKysr42b3brd7UibYTCQTWfI3FSilovVpMMSqT+PKTLRU14k2+rq7KK2djTlKxJYKK060dJPndpLnceJvbsbfd5ojxWFmOGvIL/SMea+Ul/xVSm0CNg357Dsxr69K5joazWQgtjn5QEhdb28vhw8fjur2u3btYvv27YAlE8Q6ex2GmX6i9Wmckfo0J/swChyYhenPbM11e+jt8tLX5SW/uGTE48QQDNOIS2SSsEKhe6hqNGkjNzc3mr0IVq/a48ePx83uB8of22w2ZsyYEefw8/LyMmn+tMFwROrTePut+jS+EGaxE8Oevqgom8OBMzeP3i4veZ6iaGmBRJj2mLruYoDuoarRZJaBOPoZM2Zw0UVW5G9XV1ecs3/ttdfYtm0bYLVxq62tjTr70tJSvVGbIiQiywy29OtLezOQXLeHjqOH8fWcJqdg5PwK02bg7w1G7DaQiAyuZ+4azSSisLCQhoYGGhoaAAgEAhw5ciTq7Pfv3x+N8nA6nXEz++rqapxOZybNHxM1EI89RTBy7NgdJsH2SDMQXwibJz3NQBw5OdgcDnq8nbjyC0b8uZk2g3BYEQ4rMCRaOCyZOPdk9kNHQzt3jWac2O32uNaFSina29vjZvcvv2wlb4sIFRUVcQ7f4/FMGmfqcrk4deoUJSUlk8amZBBzSDOQtlBaNltFhFy3h64TbQR8Phw5iZt5RKtDBsLDK0OOwsCmv8vlGr+NZzs6jBcdLaOZDvT19UU3ageyav1+PwD5+flxzj6TfWoDgQCtra34fL6M3H8iUMEw4d4AKqQwXDYMlzmsfMGE3k8pTrefxLQ7yC10JzwmFAzT6/XjKrBjBnyEvF7aC8Bj8+DMGT2hzuVyUVNTg31IkxDdIFujmQTk5ORQX19PfX09AOFwmLa2trjZ/XvvWan+Azp/rMPPP8PmzOPFbrczZ86ctNwrlYT9IbzPNNGz/Rj26nyKV5+DvSx1WcmvPPprdmzcwNr7HsBdXjHse78vyAPf2MrFN86l1neAo9/6Nj/6isF3F32fpVd/JmV2gXbuGk1aMQyDyspKKisrueCCCwDo7u6Olk9oaWnhjTfe4NVXXwWgqKgoztmXl5frjdpRMBwmRZ+ah2t+ER0bGmm7bxeeVXXkfqgiJXLT0muuZ8fTG9j13DNc+bm1w753uGzkFNjxnujDKLcGGUdAIBweduxEo527RpNhCgoKWLhwYTRZJRgMcvToUZqbm2lpaeHgwYO89dZbADgcjmglzIF6OWejy2YrOYtKccwsoP2JfXT8oRFfYydFN9VPuBZfUFLK/Isu5Z2XnueSm/8Jh2u49u4uy6XrRB/GbMu5OwMQSoMcrp27RjPJsNlsUecNlrbb0dERJ+Vs3bo1Gt1SXl4eDcOsra3F4xk783E6YLqdlH5xMd1/a6Hrrx/gb+mmZM0CHDMLJvQ+y6+7gX2vvcK7W15k2TXXD/veXZ7D4X0dGLmW43cGFIrUdxvTzl2jmeSICMXFxRQXF3PeeecB4PP54iph7t69mx07dgCDjU0GHP50Ln0shlC4ohZnnYf2x/bS9qs9uK+ZRf7lNYgxMTJN1bwFVNbPZ9fmp1l69XXDkprcZTnse/0YYUcRAC4/BMPauQ8jFOpFqRCmmT+lQrY0monE5XLFbdTGZtQOyDkDjU3sdjs1NTVRZz8dpRznrEIq/nUZHRsa8W5+H9+BTopvOQezwHHW1xYRll97A5t+9lMO7fk7c5ddEPe9u8yasff0W+7WFQDQmvswDh9+jMYDP0LEhs3mxm4vwm73xP+zFWG3D35ni/ncNCd3IolGMx4SZdR2dnbGOfvRpJzpUBzNyLVT/NmF9Gw/RufTTRy/dyfFt5yDa37RWV97/sWXsvXh9ezctHGYcy+MOPeuHmtG7/JDKNK0I5VMOefu8VxIff23CAS8BAIdBAKdBAId+HytdHe/QyDQQTg8cjl5w3BFBoHIoGAbHBSig4C9CHvMwGGzuTGMKfej0kxzPB4PHo+HxYsXA5aUc/jwYZqbm4dJOQUFBdTW1ma9lCMi5F9UhXN2Iace3cvJ9e+Q/5Fq3B+fjdjGH4Vk2uyc9/FPsO2J33GqtZmSmsFOc55IKGZ3D+RgbaiGtXMfTmHhYgoLF496TCjkG3T8wc7oABAMDLwe/Px0z37ru6AXNcoP3GYriKwIPNhiVgWDg8PQFUSRlo40kwqXy0VdXR11dXXA9JZy7BV5VNy5lM5nD3F662H6m7yUrF6ArTRxpmkyLLlqJa9veJydmzdy9ZfujH7uzLPhyLHR3Rkkh8jMPayd+7gwTRemWYXLVZX0OUopQqHTMauBwUEhEBxcJQwMEH29HxAIdhIMdo14TREzRjqK/HfYSiF+lWC3ezDN8f+CaTTJMt2lHLGbFN1Yj6veQ/sfGjl+3y6Kbqond9n42t7lFrpZeNlH+cfWl7lszW3k5FtROSKCuyyHrpM+Su0mzkCYMNq5pw0RwWYrwGYrICcn+ebd4XCQYNCbYJUQIxtFBgef7wjdgXcJBLyEw30jXtMwnFGnb0u0Kki0p2BzYxijtyDTaMZiNCmnpaWFPXv2ZJ2Uk7OolIqafNof30f7E/vwNXbgWVWP4TzzZ1l+3Q288/LzvP3ic1y4ajAD1V2WQ1tzNyGnDZffT1jpDdVJj2HYcDhKcDhGLtqfiFCon0BwyCAQlYzi3/f0HCAYGTiUGjmEyjTz41cJo+4pRD6zFSCiMx41iZkuUo7N46LsS0voeqmZ7pea8Td3U7xmAY7qMyv/UFY7m5kNS9j93LOcf/1NGJHBzl2WQ9OuEwSdLlwBv565ZzOm6cQ0K8FZmfQ5g9JR/GZyYGDlMLBiiAwOfX3NlpQU7GLkOnRGxOG7Iw6/KH5PISobxa8gDCNnSi/JNeMjm6UcMQX31bNw1blpf3wfbb/cjfvaOeRfOuOMbF5+3Sqe+skPadz+Gud8+DLAipgJhxX9+aW4Al06WkYTT7x0VJP0eUqFCAa7BgeDUfYU+vuPc/r0XgLBTkKh3hGvaRgObHGbyYODQ6yEZIvbU3BjGGcfV6yZXIxHyhlw9pNRynHO9VD+9eV0/LER7zNN9Dd2UHTzfMz85H535y4/H3dFJTs3PRV17p5yax+tN68cZ18TYb2hqpkIRMyIcy0Ckq/8Fw73x68SgoODQjB2cAh00tt7KPqdUoERr2maeXF7B/GrhHgZaWDgsNkKtXQ0hUgk5bS1tUWdfXNz86SXcsw8OyWfW0jP60fpfLaJ4/fuovjWc3DVj13awTBMlq/8JC//5gGOHWyksm4ehaVWOKTPVYazS+HXSUyaTGIYTpzOcpzO5KMHLOmoN7p3EEwoIQ2899Lna41IR15Glo4kEnUUs0pIGGkUGSwiIaummTtpJYDphGmaVFVVUVVVNaWkHBEh/8MzcMwqpP2xvZx86G0KrphJ4dW1Y3Z7arjyarb9/mF2bt7IdXf+O3luBza7gc9ZQk4AfHpDVTPVsKSjPGy2PHKoTvo8pcIR6Wjo5nK8hBQMePH7T9DT00gg4CUUOj2KLY7BgSBGQhp5tTAgHeks5lQzlaQcx4x8yu9ahvfpJrr/1kJ/UyfFqxdgKx55heHMzaXhyqvY8/xmPvLZz5NfVExhWQ4+XwkeP4R14TDNdEHEiDrbMyEc9g9KR0EvwRH3FDrp7fuAQNeeiHTkH/GappmL3RYbXZQ4H2FQNnJjt7sRmVza8VRisks5hsOk6NPzcNZ76NjQyPH7dlL0qXnkLikb8ZxlKz/Jrr88w54XNnHpLf+MuyyHY0fcuALoUEiNZiwMw4HTWYbTOfIf2VCUUoTDffEDQEw+wtBchdP9RyMDSCcjF3wSbLbCYfkI8dnMQ6KQ7G6dxTwCk1XKyT2vzKoT//he2h/dS39jJ+5PzsVwDB/YiypnMHf5Bex5YTMX3XgLhWU5vE8hzoAQ1pq7RjPxiAimmYtp5uJyzUj6PEs66o6Ri+KzlmNlJL//FD09TQQCHWNIR3bsdndM5NHQZLUhq4TInsJ0LIA3VMrp7++ntbU17VKOrdhF2VeW0PVCM91bWuj/wEvxmoU4qvKGHfuh61bx5A/vYe+2LXjKFqKwYVBISJf81WgmD5Z0ZEkwMCvp88LhwGAWc3Sl4B2WrHZmBfBy4pPVhuQjJKp/lG0F8JxO5xlLOQPO/mylHDEN3Ctn46x30/7EPtp+sQvPJ+aSd3FV3IphZsMSSmfOYufmjVxx2zLAipghqJ27RjPlMQw7DkcpDkfpGZ0XCvUl2DtIVACvg9On9yVZAK8wqUqoU7EAXiIpx+v1xjn7V155ZUKlHFd9ERVfX07Hk/vpfOqg1c7v0/Mw86xSICLCsmtv4IX7f0Zv5/sA9OWUQf/I4cIThag09PJLxPnnn6/efPPNjNxbo8lWRi2Al7D+USRkNdg94jXjC+ANz0eIq39kix0UMh+vPpShUk5rayt+v7W5fjZSjgorTm87gvcvhzDz7RTfugDnXDcAAX8/99/+earPOZcjBy+m9oMXaLuxiy985VfjegYR+btS6vyxjtMzd40mi5joAniJ9hTOuABeRBIaupk8fE9hQDpKXQG8VEk5YggFl1fjnGPFxJ944C0KVtRSuKIWu8PJeVet5I0/P0lBSR19OWUY/lMpe8aoTXrmrtFoxkso5IsOBsEx9xQ6Y6SjZArgDd1gTpCrMFAPyZY/YVnMQ6Wc48ePn5GUE+4P0vnUQXp3tuGYXUjx6gX0hbt58M61uPIXktezALXsZW67+zfjsm9CZ+4ishK4FzCBB5VS64Z87wR+C3wIOAXcqpR6/0yN1mg0Uwurd0IlrnEVwEucjxDdV4gMGn19HxAIeCNZzCNhxAwCCVpsnkEBPLfbzeLFi8cfleO0We375hXR8acDVju/T89j/sWXse+1VyH3InL9qdfcx3TuYmVm/AK4GmgFdojIRqXUP2IOWwt0KKXqRWQ18D/ArakwWKPRTG3ipaOZSZ+nVCjq5ONKWkST2AYHC1//saQL4CWTj1BcUkRFZQ2XXLIIw8jn5MnOpKScqq+ey+kNhzj18HssXfBRGtU2/KEm8nypT3hLZuZ+IXBAKdUEICKPA6uAWOe+Cvhe5PUfgJ+LiKhMaT4ajSbrEDFxOIpxOIo58wJ4Q+ShERrr9PY2JVkALx+73U3tLA91dUVAHn19Bt2nFR3t77N3Xz9vv+0kGHTgdldTMrsCz/58Lq39LLuPbSHcNzwmfqJJxrlXAy0x71uBi0Y6RikVFBEvUAKcnAgjNRqNZrxYBfAqcDorkj7Hko564vYOhlZCja971Eog2InN5qWsXFGWoNZeaJ5wKuhgRsBJz96lE/iEiUnGuScK/Bw6I0/mGETky8CXAWprk9/J12g0mnRiSUf52Gz5Z1gALxTJYo7fTPb72+nsPMzxw/vxdZ0gTOozjJNx7q1ArDBWAxwZ4ZhWEbEBbqB96IWUUvcD94MVLTMegzUajWayYvVOSFwAb9Ys4Lz02ZJM7NAOYJ6IzBERB7Aa2DjkmI3AbZHXnwFe0nq7RqPRZI4xZ+4RDf1O4DmsUMj1Sql3ReQHwJtKqY3AQ8DvROQA1ox9dSqN1mg0Gs3oJBXnrpTaBGwa8tl3Yl77gJsn1jSNRqPRjBfdmFKj0WiyEO3cNRqNJgvRzl2j0WiyEO3cNRqNJgvRzl2j0WiykIyV/BWRE8AH4zy9lOlX2kA/8/RAP/P04GyeeZZSasyO8Blz7meDiLyZTD3jbEI/8/RAP/P0IB3PrGUZjUajyUK0c9doNJosZKo69/szbUAG0M88PdDPPD1I+TNPSc1do9FoNKMzVWfuGo1GoxmFKeXcReRmEXlXRMIicv6Q774tIgdEZJ+IXJMpG1OJiCwVkddFZLeIvCkiF2bapnQgIndF/r++KyI/zrQ96UJE7hYRJSKlmbYl1YjIT0Rkr4i8JSJ/EpHhBdGzABFZGfldPiAi30rlvaaUcwfeAT4FbI39UETOxSoz3ACsBH4ZaeydbfwY+L5Sainwncj7rEZEPorVo3eJUqoB+GmGTUoLIjITqyl9c6ZtSRMvAIuUUkuA/cC3M2zPhBPxSb8ArgXOBdZEfFdKmFLOXSn1nlJqX4KvVgGPK6X6lVKHgANYjb2zDQUURl67Gd4RKxv5GrBOKdUPoJRqy7A96eJ/gW+SoF1lNqKUel4pFYy8fR2r41u2cSFwQCnVpJTyA49j+a6UMKWc+ygkauKdfOPDqcM3gJ+ISAvWDDbrZjcJmA9cLiJviMgWEbkg0walGhG5ATislNqTaVsyxBeAzZk2IgWk1U8l1awjnYjIX4HKBF/do5R6aqTTEnw2JWc8oz0/8DHg35RSfxSRW7A6YF2VTvtSwRjPbAOKgIuBC4Dfi8jcqd7GcYxn/k/g4+m1KPUk87ctIvcAQeCRdNqWJtLqpyadc1dKjcdZJdPEe0ow2vOLyG+Br0fePgk8mBajUswYz/w1YEPEmW8XkTBWXY4T6bIvFYz0zCKyGJgD7BERsH6Xd4rIhUqpY2k0ccIZ629bRG4Drgc+NtUH7xFIq5/KFllmI7BaRJwiMgeYB2zPsE2p4AhwReT1CqAxg7akiz9jPSsiMh9wkMVFppRSbyulypVSs5VSs7EcwvKp7tjHQkRWAv8B3KCU6s20PSliBzBPROaIiAMrCGRjqm426WbuoyEiNwE/A8qAZ0Vkt1LqmkjD7t8D/8Ba0t2hlApl0tYU8SXgXhGxAT7gyxm2Jx2sB9aLyDuAH7gtS2d1052fA07ghciK5XWl1Fcza9LEopQKisidwHOACaxXSr2bqvvpDFWNRqPJQrJFltFoNBpNDNq5azQaTRainbtGo9FkIdq5azQaTRainbtGo9FkIdq5azQaTRainbtGo9FkIdq5azQaTRby/0UlyatI90u+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x737d14ea10>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for t in np.logspace(-2,6,9):\n",
    "    plt.plot(y,T(t,y),label='t='+str(t))\n",
    "\n",
    "#plt.plot(y,T(100000*t,y),label='t='+str(100000*t))\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x737b111bd0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl8nNV19793Fu2jfbVG8iYTsLHBYBsDYcc1JolJG5qYpMG0IbxAE/ombyg0TdJCmwBt2hRe0rS8kNSkITRNm0ASs4Y4W1lsgwPYgDeMNdq3kUaa9Xme8/4xi0fSjDTaRtv9fj7zseZZ7nNnLN1z7znn/o4SETQajUajsc12BzQajUYzN9AGQaPRaDSANggajUajiaENgkaj0WgAbRA0Go1GE0MbBI1Go9EA2iBoNBqNJoY2CBqNRqMBtEHQaDQaTQzHbHdgIlRWVsqyZctmuxsajUYzr9i/f3+3iFSNd928MgjLli1j3759s90NjUajmVcopd7L5LqMXEZKqauUUu8opY4qpe5Mcf5ipdSrSilDKXVt0vHLlFIHkl5BpdSHY+f+TSn1btK5szP9cBqNRqOZfsZdISil7MA3gS2AB9irlHpSRA4lXXYSuAH4QvK9IvIL4OxYO+XAUeDZpEtuF5EfTuUDaDQajWZ6yMRltAk4KiLHAZRSjwPXAAmDICInYuesMdq5FnhKRPyT7q1Go9FoZoxMXEb1QHPSe0/s2ETZAXx/xLGvKqVeV0p9QymVO4k2NRqNRjNNZGIQVIpjEyqioJSqA9YCzyQd/gvgdGAjUA7ckebem5RS+5RS+7q6uibyWI1Go9FMgEwMggdoSHrvBlon+JyPAj8SkUj8gIi0SZQQ8B2irqlRiMhDIrJBRDZUVY2bNaXRaDSaSZKJQdgLrFJKLVdK5RB1/Tw5wedcxwh3UWzVgFJKAR8G3pxgmxqNRqOZRsYNKouIoZT6DFF3jx34togcVErdDewTkSeVUhuBHwFlwIeUUneJyBoApdQyoiuMX45o+ntKqSqiLqkDwM3T9Jk0Go1m3mFZBoYxQCTiJWL0EYl4GfR18LvfPI0Z9NPQuINNW68dv6EpkNHGNBHZDewecewrST/vJepKSnXvCVIEoUXk8ol0VKPRaOYDIoJpDkYH9khf7F8vEcObOGZE+k+dix03jIGU7bnqQUTRfWz9jPd9Xu1U1mg0mmximsHowG2cGsCNyKmBPXp8+MBvGP2IGGnbdDhcOB1l2B3FmGY+weASBgfr8HojBAN2IkYOdlsxpdRQ3FVC/lAeR9pOsG7naTP+ebVB0Gg0C56oO8Z7arY+1qw96TrLCqZt02bLw+ksjb4cpRQVnobTWYojfix2PPpzGeGwg9bWAZqbWzjc3ExbWxuWFd26VVlZSWNjIw1nNrCksArZ3YHR7ifYGOEnb95PgfMyHAX5M/49aYOg0WjmDSKCYfiiA3jSrD36b3/C9578MgwvhuFL26ZSDhyOEpzOMpzOEvLy3LhcZ+JMHIsO6A5n0ntHKXZ7Xto2Lcuiu7ub9947SXNzMydP/pa+vj4AHA4HS5Ys4YILLqChoYGGhgYKCgoQEfz7OvB+9xgqx0bFztX8x8NfBnsxxSEbKl8bBI1GswARESwrmNLHPswlM2wWH3fHmGnbdTiKEwN4jrOMwoIVaWbspYnB3W4vIprsOHnC4TCtra2cPBk1AM3NzQSD0dVFQUEBjY2NbNy4kYaGBurq6nA4hg+9VtCg70dHCfyui9yVJZR/7H14TrxFj+ckzvzLyR3sxu5Mb4CmC20QNBrNlLCsyDg+9iT/emLw78OywmnbtNnykwbtEoqKTk85oDucJTgd0YHd4SjGZsvOkObz+WIz/6gBGOn+Wb16NQ0NDTQ2NlJeXj6mwQmdHKD38XcwvUGKty7FdUkDyqZ49aknyHeVYNnPxBH+EXabfcY/lzYIGo0GABEr4Y5JNaCnyoyJRLyY5mDaNpVynBrYHaXk5TfgKl6XdCzZLVMaG9hLsdvnjpJN3P0TH/xPnjw5rvsnE8QSfL/yMPDse9iLc6j6X2eRu7QYgL62Fo6/to+1l32YI685sBndOGzOGfuMcbRB0GgWGNG0R/8IH3vyrH20jz36cz+QTp9Sxfzs0QE8J6eSwoKmU4P4CDdMfCZvtxdO2R2Tbabq/skEcyBM7w/eIXTUS/7aSsr+YBW2/FPtvPb0T7HZ7NSsej9HXmsBqwubXiFoNIsbywqPyIoZkcM+IlMmEommQYqkd8fY7YXDZuZ5eUuG+9VTztqLiSrhLzym0/2TCYF3eun7wWEkbFL6B00Ubqwd1mbIP8Sbe57n9AsuIjTkBLGw6MGONggazYJAxDy1CzUjH3vcHTOUtk2lcmIDdnQALyhYPqaPPX6tzTZ33DHZJtn9EzcA0+H+yQQxLPqfPsHgb1pw1hZQft1anDWFo6578xfPEQkGOOfqa3h9TwCn0U/IYWLPQnxEGwSNZgJE3TFDKWbmycHU0bN2w+gnvUiwLTaoR/3nObnVFBadluR6GT5bdzrLcDhKsNsL5p07JtuM5f4pLCykoaFhyu6fTIh0B+j9/ttEWgYp3FxH6QeWo5yjZ/yWZfLa0z+h/vTV1Kxoov8/95Eb7iGUgw4qazQziWWFRm9UGjVrj7lojFODe5Jo7yjs9qJhKY75ee6kHPZULpkyHA4XSmVUzVYzDtl2/2TC0KsdeH98DOyKik+eQf6ayrTXHt+/l/7ODi7+xB8D0N8VoCjYxaCTrGRQaYOgmfeImAnf+dg+9uEzessKpG3TZss5lc7oLKWgYEXKAX1YMNVRjM2Wk8VPvrixLIuurq5hBiBb7p+M+hcy8P74GP7XOslZXkz5x07HUTq2u+7V3U/gqqyiaeP5hAIGwcEIVYFuegogM3HqqaENgmbOML4oWNwtMyIdMo0oWBTbMHdLXm4trqLTRw/oI3zvNluedsfMMZLdPydPnsTj8cyK+yejvnp89H7/bYzeIMVXNuK6vBFlG/v3qfPEcZoPvcHFn/hjbHY7A13R3dWFQx2ESwSy8PuoDYJmRjDNYHp1x1GiYP0JKYLxRMEcSQN3fn5jbDAvSwRWk33s0X+LtDtmnjIX3T/jIZYw+NsW+p8+gb3ISdVN68hdXpLRva8+9SSO3FzWXr4ViLqLAAp9XYRzyMrn0wZBMyanRMFSb0oavit1sqJgq0bpxAwPopbicJRgy8LGHM3sMNfdP5lgDobp+8/DBN/pI291BeXXrsJWkNnvrH+gn7d/+0vOvHQLeUVFAPR3+QEoDHQRcaBXCJrpIy4KZozIjEnOb08e9I2YUNjYomD2YTP2VKJgw33s0eNjiYJpFgfhcJiWlpaEAZjL7p9MCB7po/cH72AFDEqvWUnh5roJzehff+4pzEiE9ds+lDjW3xUgv8iBwwwRyYHU5e2nl7n9LWtSYpqBYa6WlKJgI/PbMxUFc5QmiYKlm7VPnyiYZnGQifunsbGRhoaGOeP+yQQxLQaeew/fLz04qvKp+tRanLWj9xaMhWlEOPDsz1h21jlU1J8qX9/fGaC4NLrCiDhAsuD61AZhFkkWBUutE5PkX08KsFpWKG2b6UXBUszaZ0EUTLPwWQjun0wweoP0fv9tws0+CjfVUvLBFdhyJr5X4PCLv2HI28fWW/73sOMD3QHq6qJ/l5Ecxg1KTwd6FJgGRomCjQqmpp61T1oUzFE6bNBPDPRzTBRMszhYaO6fTPD/rou+/z4CCso/fjoF66om1Y6IsH/3k5QtcbNs3akSmUbYZLAvRNGqqIExHKBdRlkmqtEeSO1jT1VGL2NRsCSN9pyK9KJgSbP4+SgKplkc+Hy+Ycqf7e3tCfdPVVXVvHX/ZIIVNvE+eQz/vg5yGl2U7zgdR/nkY2Kth9+m4/gRrviTW1C2Uy6hge6oQXXlRd28hlNQc2WnslLqKuB+wA48LCL3jjh/MfBPwDpgh4j8MOmcCbwRe3tSRLbHji8HHgfKgVeBT8pYilxTwOc7RDDoGeWGGb1haTxRsIJhUgK5eXUpfOyLRxRMs/BZLO6fTAi3Dkb3FnQHcF3WQPGVjSj71Pz6rz71JLmFhay+5PJhxxMZRjkhQkDECXNihaCio9k3gS2AB9irlHpSRA4lXXYSuAH4QoomAiJydorj9wHfEJHHlVL/AnwK+NYE+58Rx47/Az09exLvlXIOC5AWFCwbFjQd7mM/VVpvMYuCaRYHi9H9Mx4iwtCLbXh3H8eW76TyU2vJayqdcrsD3V0cefm3nPuBD5OTN7w8ZnwPQpE9QAiwnJKNrNOMVgibgKMichxAKfU4cA2QMAgiciJ2Lp3fZBgquoa8HPh47NAu4K+ZIYPQtPIOVq74fCJFUouCaTRRFrP7JxPMoQh9PzxM8K1e8k4vp+zaVdiLpkee5MCzPwOB9Vs/OOpcf1eAnHwHTiNqGOZSDKEeaE567wHOm8Az8pRS+wADuFdEfgxUAF45tS3VE3vOjFBUdNpMNa3RzBtGun9OnjyJ1+sFhrt/GhsbcbvdC9r9kwnBY176/uMdzKEIJR9cQdGFS6bNIEZCQd54/mmaNm2muKp61PmBrgAlVflIIOo6Mp0CtrmRdprqG0in45uKRhFpVUqtAF5QSr0BpBKfSdmmUuom4CaAxsbGCTxWo1nchEKhhPunubk5pftn06ZNi8r9kwliCgM/fw/fL5pxVORTvXMNOfVF0/qMQ7/6BcGhQc7Ztj3leW9XgOpGF5Y/ukKwHKDmiLidB2hIeu8GWjN9gIi0xv49rpTaA6wH/gsoVUo5YquEtG2KyEPAQwAbNmyYiCHSaBYNIoLX600M/s3NzXR0dCAS/ZPR7p/MMLxBeh9/h/CJAQrOqab0miZsudObFCIivPb0T6hevpL609eMOm+aFoM9QVadW4113I+lQOzMGemKvcCqWFZQC7CDU77/MVFKlQF+EQkppSqBC4G/ExFRSv0CuJZoptFO4InJfACNZjFiGAZtbW3DDMDgYHRfi9PpxO12c9FFF9HQ0IDb7SY/P3+cFjX+17vo+++jIEL5x95HwfrRrpzp4L03DtDjOclVt34upVEe7A1iWUJxVT7WwQChHIUdGZaWOlOMaxBExFBKfQZ4hmja6bdF5KBS6m5gn4g8qZTaCPwIKAM+pJS6S0TWAGcA/xoLNtuIxhDiweg7gMeVUn8LvAY8Mu2fTqNZIAwODuLxeBL+/9bWVkwzmqNeWlrK8uXLE6mf1dXV2O061TlTrFBsb8H+DpwNLip2vA9HxcwZ0Fd3P0FBSSnvu+DilOfjGUal1flYAT8Rp4oN1HNjhYCI7AZ2jzj2laSf9xJ1+4y873+AtWnaPE40g0mj0SSRHPyNv3p7ewGw2+3U1dUlfP8NDQ24XK5Z7vH8Jdzso/fxaN2C6dpbMBa9rS28+9o+zr/24zicqZVQ+zujBqG4soB+vz+6QhDJiotPxX2M84ENGzbIvn37ZrsbGs20EgwGRwV/Q6GoXlU8+Bt/1dXV4UwzkGSbSCQyLFA9r5BoRTMrYKBsCluhE+WYeZdMcNBHOBigqLwSWxoXUNAfIRI0cZXlYfT2YERCDBUKpYVLsI+jZ5SXl4fb7R71O6KU2i8iG8brn04r0GiyiIjQ19c3Kvgbp6amhrVr1yYMQFlZ2ZwN/no8HlwuF8uWLZuzfUyFGBZGXxAJmdjyHdhLc2d0VRDHMk26Tp4gr7CIkuqatNd5O/2YhkXFkiJC775LMOynv8yiruJ0nGP0U0To6enB4/GwfPnySfVRGwSNZgaJRCKjgr9DQ0MA5Obm4na7OeOMM2hoaKC+vp68vPlTKyIYDM47Y2D6I5jeEAjYy/KwFTiy1v+AbwCxLApKxq6gZhoWjvhqxbKwMrRVSikqKiro6uqadB+1QdBoppG47n/81dbWlgj+lpeX09TUlJj9V1VVpXUbzBfmizEQSzC9ISx/BJVjx1GWh3Jm77sXEfwD/eTk5ePMTW/0RQTTEHLzon0Ty0JU5uHkqf5/aIOg0UwS0zTp7OwcZgDiO3/tdjv19fVs3rw5kfpZVDS9m5sWO16vl8cee4xbb711zOussInZG+T48eN88rZP0dffxznnnMN3v/tdcnJGy1Dcc889PPLII9jtdh544AG2bo3WOP6TP/kTfvrTn1JdXc2bb745ob6G/EOYkQiu8sqx+2oKiGB3Jq0QsjhKa4Og0WRIIBBIpH42NzfT0tJCOBxVxy0qKqKxsZHzzjuPhoYGamtr9c7fGcbr9fLP//zPaQ2CiGANRjD7Q2C38aV/vJvPf+Hz7Nixg5tvvplHHnmEW265Zdg9hw4d4vHHH+fgwYO0trZy5ZVXcvjwYex2OzfccAOf+cxnuP766yfc1yGvF7vTSW7h2NXUTCOqI2VPchlNZIUwVfRvrEaTAhGht7c3IfzW3Nyc8M0qpaipqeGss86ioaGBxsZGSkpK5o37ZKFw5513cuzYMc4++2y2bNnC3//93yfOjQwc20py+MWeX/D9x78PwM6dO/nrv/7rUQbhiSeeYMeOHeTm5rJ8+XKampp45ZVXOP/887n44os5ceLEhPsZDgaJBAO4KqvG/R2JGwSb45TLyMrir5U2CBoN0eBva2vrMPeP3x8VFoun8p155pmJ4G9urpZCn23uvfde3nzzTQ4cOIDP5+Pss2Mq+5YgZjSdXtkV3/v+Y9RIDaWlpYlVm9vtpqWlZVSbLS0tbN68OfE+3XUTwd/fh81mIz+D/SJmJL5CUFHZEREspVCSnVWCNgiaRUl/f/8w909y0feKigpOO+20RPC3sjJ9zrgmyl0/Ocih1lSalZNn9ZJi/upDo7V+UuFyuXjt1deGBY7tZXnYYr74VJk3qWbrqfZlTWXlZ0QiBAcHKSwtw5ZBxTPTEOwOG0opJJaMIHqFoNFMH8m6P3Ej4PP5gKjsc319faLql9vtpnAcP69m7tHf4+XiSy6OaibbFMp+ahR97LHHOOOMM/B6vRiGgcPhwOPxsGTJklHtuN1umptPqf2nuy5T/P1elFLjpprGMQ0rET+Q2ATFUmRB5zSKNgiaBUfy7N/j8QxL/SwtLWXp0qWJwb+2tlbr/kwDmc7kpxOXy4XP58P0hckP2Nn73P/gKMvDlpd6WLvsssv44Q9/yI4dO9i1axfXXHPNqGu2b9/Oxz/+cT7/+c/T2trKkSNH2LRpcgo7lmkS8A2QV1SE3TH+7nIRwYxYOAtj12qDoNFMjPFm/0uWLElk/rjdbq37s4AoLynj/A3nsW7D2Vy1ZStfv/8fxtxxfN9997Fjxw6+9KUvsX79ej71qU8B8OSTT7Jv3z7uvvtu1qxZw0c/+lFWr16Nw+Hgm9/8ZmLCcN1117Fnzx66u7txu93cddddiTZSEfD1xzailWX0ecQSZETKKZDVLCOtZaSZV4w3+3e73YnBv6amRqd+ziBvvfUWZ5xxxqw8e9iO49LcrO44zgQRi+6T72F3OilfMkr3MyWRkEFfu5+SqnxyC5yYQ0OE332X9jJFnsOiqnINjgwkNlL9v2gtI828Z+Ts3+PxMDAQDVzq2f/iREwrGjgOGNEdx+V5WRGlmyjBwUFMw0hZHjMdZiQ6OU/egwBg2fQ+BM0iZLzZf7zWb0NDg579L0KsoIHRFwRTsBfnYHPlzKlVQRwRwd/vxZGTQ05+5nWpU21KA/Q+BM3CR8/+NZkilmAOhLEGwyiHDXt1PracuZsIEA4GiIRCFFdVT8hgmYaFzW5DxSSuk7OM9ApBs6DQs3/NZLDCJmZfEIlY2Iqc2ItzEwPmXMXv9WKz28kvmtgkJppymvTZkoLKZCnUq//qNNOOnv1rpkpCh2ggBDaFozI/bTrpXMIIhwn5hygqK59wDWQzYpFTkPQZk1cIFllZJsz9b1gz5xlr9l9SUqJn/5oJMVsFbKaDof4+lFLkZ7gRLY5lCZYlp+IHnHIZ6Z3KmjlLfPafbABSzf7jBkDP/jWZIiJYfiOaTsr4BWwylb+O8+6777Jjxw56e3vTyl/39PRw7bXXsnfvXm644QYefPDBjPtvGgZBn498VzF2+8SG1lMaRkmGz7LAZgMExRzSMlJKXQXcD9iBh0Xk3hHnLwb+CVgH7BCRH8aOnw18CygGTOCrIvIfsXP/BlwC9MeauUFEDkz1A2mmj3i5R4/HQ0tLCx6Ph/b2dj3710w7k0knHU/+eiR33HEHn/vc58aUv87Ly+Nv/uZvePPNNydc8yAw0I+IUFBSOqH7ICnDyJm8QhBQWQwgkIFBUErZgW8CWwAPsFcp9aSIHEq67CRwA/CFEbf7getF5IhSagmwXyn1jIh4Y+dvjxsPzewTDAZpbW1N+P09Hk9C8dPpdLJkyRI2b95MfX29nv1rpo3h6aS52FzOjLJzxpK/HomI8MILL/DYY48B6eWvCwsLef/738/Ro0cn9hksC/9AP7kFhThSFN0Zj1Epp9FGY3EIa8LtTZZMpnObgKMichxAKfU4cA2QMAgiciJ2bljPReRw0s+tSqlOoArwoplVLMuiq6tr2OCfrAhZWVnJaaedhtvtpr6+nurqaq35o5lWoumkIazByKTSSdPKX4/gscceo7q6OiP568kSHPRhmSaFpRNfHUDUZWSz27AlZ1BJ3GU0t9JO64HmpPce4LyJPkgptQnIAY4lHf6qUuorwM+BO0UkNNF2NZnh8/kSbh+Px0Nra2ui2ld+fj5ut5s1a9YkDEB+fv4s91gzr3jqTmh/I+PLRSSaSiqC3W4Dh0KNHPZq18K2e1M3MAKXy8WBA+k9zpnKX08GEWHI24czNw9n3uT+boyIhWNEjWexLBiWYjvzZiETg5CqFxNyaiml6oDvAjtFJL6K+AugnaiReAi4A7g7xb03ATcBNDY2TuSxi5ZIJEJ7e3ti8G9paUnU+rXZbIlqX263G7fbTXl5+Zzc8alZeAgCpiCGBUqhnPZp2Vfg8/m46KKLUp6biPz1ZAgNDUbrJddUTurvSEQwDQtnwQhFVPPUCiFbZGIQPEBD0ns30JrpA5RSxcDPgC+JyEvx4yLSFvsxpJT6DqPjD/HrHiJqMNiwYcP8UeLLEsmB3/jgn1zspbi4GLfbzaZNm3C73dTV1eF0ji/Fq9FMiAxm8lYkWuxeIha2Amc0nXQKxiAufx3/eawVAmQmfz1R4qsDRwb1ktO2YQliybCAMhDNMoq5aeeSy2gvsEoptRxoAXYAH8+kcaVUDvAj4FER+c8R5+pEpE1FTeqHgYmF9BcpwWCQlpaWYe6fkYHf888/P+H6KS4unuUeaxY7wzaZKYWjIg9b/tQnJRUVFVx44YWceeaZbNu2bcygMmQmfw2wbNkyBgYGCIfD/PjHP+bZZ59l9erVKdsMByYnU5GMkVQ2MxkRC7HFvqe5slNZRAyl1GeAZ4imnX5bRA4qpe4G9onIk0qpjUQH/jLgQ0qpu0RkDfBR4GKgQil1Q6zJeHrp95RSVUSN3wHg5un+cPMdHfjVzHckEttkFp6ZTWbxrKFMWLFiBa+88sqo49u3b2f79u2J9ydOnMi4zaH+PmyOictUJJMq5RSI7UOIGolR8ZUZIqOkcRHZDewecewrST/vJepKGnnfvwP/nqbNyyfU00WADvxqFgoigjUUwewPgwJ7eR62/LlVs2CqREJBwn4/roqKCctUJJNyUxqxoLKaezEEzQwwMvDr8Xjo74/u0dOBX818Jll6QuU5cJTmzsmaBVNlyNuHstnId01MpmIk8TrKyX/fIjJshQBzaKeyZmpYlkV3d3fC99/S0kJHR8eowG9c8kEHfjXzkVHSE6W52Aoz22Q23zAiYYKDgxSWlmGbopvWjKQJKAOScBllB20QZoCBgYFhg39LS0vC9ZOTk0N9fT0XXHAB9fX1OvCrWRCIaWH0hZDg3K5kNl34vV6UUpOSqUgmkXKa5xx5Ivpvlo2pNghTJBQK0draOizzJ54KF3f9rFu3jvr6etxuNxUVFdiynFus0cwkp+obC/aSXGxFC3NVEMc0DAK+AfJcLuxT1O6yTEFEUm9KQ68Q5jSmadLZ2ZkY+FtaWoZl/ZSVlbF06dJE0Le2tla7fjQLFrEEozeA5TdQTjv28lxszoWf5ebv9yIiFJaUTbmtREA5jcsIpbKpbacNQjpEBK/Xmxj44xu+DMMAhmf9xF0/BQWZ10/VaOYzgTe7sQbDWCXGrNU3nqj89YMPPsg//dM/cezYMbq6uqisrEx53a5du/jbv/1bAL70pS+xc+fOxDnLNAkM9JNXWDQpEbuRpBS1Y8QKwZzyYzJGG4QYfr9/lN8/vuHL4XBQV1fHhg0bEoN/WVnZgl4WazSpMAfDeJ88RuD1bvj9YhzVBbNW33ii8tcXXnghH/zgB7n00kvTXtPb28tdd93Fvn37UEpx7rnnsn37dsrKoquBgG8Ay7IoLJ366gCim9KUUtjsI8aSRHGcmMtIslNYeVEahHjKZ/Lg39vbmzhfVVXFaaedlhj8a2pq9IYvzaLH/3oX3ieOYQUNircsxecanNVi9xORvwZYv379uG0+88wzbNmyhfLycgC2bNnC008/zXXXXYdYFv5+Lzn5+Tjz8qblM6RKOQVOuYyyXD96URiE7u7uYUVeklM+XS4X9fX1rF+/PpHymTdN/9kazULAHAzjfeIYgTe6cdYXUfWHa3HWFtLy1luz2q+JyF+nk54YSUtLCw0Np6TbkmWyA4M+TMOguKp66p2PYUYsHCmM6myUz4RFYhB2797N8ePHycnJYcmSJTrlU6PJABEh8Ho33ieOYoVMircuw3WxGzXSvQHc98p9vN379rQ+//Ty07lj0x0ZXZuJuF0miIyO4CqlTklc5+WRkz89scJ4ymluYYphWGcZzRxXXnklDoeDyspKnfKp0WSA6QvT9+OjBA/24GxwUXXtKpw1k1PzzAbjyV9nukJwu93s2bMn8d7j8XDppZcSHPRFJa4rJidxnYp0AWVgVAwhWywKgzBduucazUJHRAgc6ML7k2NYYZOSbcsoen/qVUEymc7kp5OJyl9nwtatW/niF79IX18fAM8++yxf+9rXohLXubnkFkyfUYynnI7cgwCz5zLS02WNRgOAORCm59FD9P7HOziBNl9FAAAgAElEQVQq86m57RxclzSMawxmi2T569tvv33c6x944AHcbjcej4d169Zx4403ArBv377Ez+Xl5Xz5y19m48aNbNy4ka985SsU5uZghMMUlU5vZuG4KwSlhhmEbPwvqFQ+s7nKhg0bZN++fbPdDY1mQSEi+F/txPuT44hhUbJ1KUUX1o9bvOatt97ijDPOyFIvZwcRocfTDAgV7sZpNQi+ngBBv0FVw2jp7HBrK1Z/P/5lNbQNtdEYhqK61Rk9P9X/i1Jqv4hsGO/eReEy0mg0qTF6g/T96AihI15ylhZTdu0qnFV6g2WckH8IIxyipLpm2vcdGcboOsoJTAts9mjJUSBbYWVtEDSaRYhYwuD/tDLwzAlQitLtKyncXDct9Y0XCiLCUF8fdqeTvCkUwEmHGbHIyUszBFsmym5LMgjZQRsEjWaREWkfove/jhBp9pH3vjJKf78JR6neezOScMBPJBScUnnMdFiWYJkpZK9jiGVBUkakmislNDUazcJADIuBF07i2+PBlm+nfMf7yD+rSkuwpGGorw+7wzGl8pjpMCNRgaIxXUYOe8p9ETOJNggazSIgdKKfvv86gtEVoGB9NSUfXIG9UCvxpiMcCBAOBnBVVk2pPGY6jHQqpzHEMrHZpy6eN1G0QdBoFjBW0KD/6RMMvdSGvTSXyj9eQ977yme7W3OeQW8vNrudfNfMKBmYEStaazpdEaGYyyg5qJyNlVxGpk8pdZVS6h2l1FGl1J0pzl+slHpVKWUopa4dcW6nUupI7LUz6fi5Sqk3Ym0+oPS6VaOZVgJv9dDxjf0MvdxG0YVLqPncuQvKGMTVTjPlwQcfpKmpCaUU3d3dieMiwm233UZTUxPr1q3j5Rf/h7DfHy2PmbQ62L9/P2vXrqWpqYnbbrttSu4cI5JG1C6OZaFstqy7jMY1CEopO/BNYBuwGrhOKTVyH/hJ4AbgsRH3lgN/BZwHbAL+SikV1439FnATsCr2umrSn0Kj0SQwB8P0fP9tenYdQuU5qLrlLEo/tBJb7sJS7J2oQbjwwgt5/vnnWbp06bDjTz31FEeOHOHIkSM89NBD3PqnfxpdHRSXDLvulltu4aGHHkpc+/TTT0+672YkfcqpiAwLKmdzppzJCmETcFREjotIGHgcuCb5AhE5ISKvA9aIe7cCz4lIr4j0Ac8BVyml6oBiEXlRoibwUeDDU/0wGs1iRixhaG87Hf+4n8Cb3RRvWUrNZ9eT27gwBRyT5a8z2am8fv16li1bNur4E088wfXXX49SinPOOot+bz++YGjY6qCtrY2BgQHOP/98lFJcf/31/PjHP55Uv8WKitrZ01WXi8lWqGH7ELJDJjGEeqA56b2H6Iw/E1LdWx97eVIc12g0kyDSMUTfj44SPjEQ3WD2B01zWoxuOpgu+etkyevBvh7q6mrpGfCxcsQ1brc78T5ZFnuixCUr0q4Q4rUQYvsQsrlCyMQgpOpPpmYr3b0Zt6mUuomoa4nGxsYMH6vRLA4kYjLwQjO+X3lQOXbKPrKKgnNrsr7BrP1rXyP01vTKX+eecTq1X/xiRtdORdwu7qcPB/yEAwHsDueogljpZLEnw3gZRqeK49gSo2K2zEImBsEDNCS9dwOtGbbvAS4dce+e2HH3iOMp2xSRh4CHIKpllOFzNZoFT/BwH30/PorZG4ymkn5gOfai7KcqzgWmIn/tdrs5efIkq1csw+Zw0NbePkohOS6KF8fj8UxaRTmhcpouw8iM7lFQtrm5QtgLrFJKLQdagB3AxzNs/xnga0mB5N8D/kJEepVSPqXUZuBl4Hrg/06s6xrN4sT0hfH+9DiB33XhqMyn8sa15DWVzmqfMp3JTyfTJX+9fft2Hrj/fi7btIG3T7xHSUkJdXV1w66pq6vD5XLx0ksvcd555/Hoo4/y2c9+dlLPS2QYpVnFnXIZ2TP3xUwT4waVRcQAPkN0cH8L+IGIHFRK3a2U2g6glNqolPIAfwj8q1LqYOzeXuBviBqVvcDdsWMAtwAPA0eBY8BT0/rJNJoFhljC4EtttP/DvmjQ+MpGav7snFk3BrPFdMlfb9u2DfeSOs6/Ygu3fe7zwzKXkuMS3/rWt7jxxhtpampi5cqVbNu2bVL9jtdRTksiqJz9FYKWv9Zo5gHhtiG8PzpC+KSP3BUllP5+06yrki4U+euQf4i+tlaKK6spKCkZ/4YpICJ0Nw+SV+TEVZ5aP8rweol4POSuWkVruAt/qJ+6kB1X/ekZPUPLX2s0CxQrbDLw/EkGf+PBlu+g7KOnUbB++sXWFisiwmBfb1SzyDX9mkUjsQxBRNJrGMGwGEK0k9nzHGmDoNHMQUSE4MEevD89jukNUbixlpJty7AVaP2h6SQc8BMJxhRNs1Bv3YiJ2qXNMCIphjDMZTR3sow0Gk0WiXQH8D55jNDhPpy1BZTfvI7cZTPryliMiAiDvT3Ync6srA5g/D0IQMq002yhDYJGM0ewwia+X8T2FDhslHxwBUXnL5mzNY3nO6GhQSKheDW07JSXNyIWNpvCZh/bICibHaXUnEw71Wg0M4iIEDzUg/cnUfdQwfpqSrYtx168OPcUZIP46sCRkzMj1dDSYUasMd1FAGJaEDMYc1G6QqPRzBBGdwDvT44RfKcPR00BVTetI3eFdg/NNAHfAEYkQmltXdYC9CKCEbHIGy8OZJnDq6WRPc9RdtZJGo1mGFbYpP/ZE7R/Yz+hEwOUfGAFNbet18ZgAkxF/vrksaM48/LILSgcJX/96quvJu7ZtWsXq1atYtWqVezatStlu729vWzZsoVVq1axZcsW+vr6Ul5nmYJYgj1nnBVCTPoaokYkm0FlbRA0miwiIgQO9tDxj/vxvdBMwdpKav/PBlwX1aPG8itrRjFZ+evGhgZMw8BVXolSapT89S233AJEB/q77rqLl19+mVdeeYW77ror5WB/7733csUVV3DkyBGuuOIK7r333pTPj2sYjRlQhkRxHMi+y0j/Bmo0WSLS5afn3w7S891DqFw7VTetpXzH6TpWMEkmI3/d2NCAWBa5+YXk5OcDw+WvN2/ejNfrpa2tjWeeeYYtW7ZQXl5OWVkZW7ZsSVkD4YknnmDnzmjtr507d6aVxTbD49RRTlxoopLE9XRQWaNZQFhBg4EXTjL429Zo9tDVyym6cIleEUyRychfD3n7EKCwrCxxPln+Gk5JW6c7PpKOjo6E9lFdXR2dnZ0p+2FELGz2cTKMGO0yyuZviTYIGs0MIZbg399B/zMnsIYiFJxbQ8nWZdhdC29F8OsfHKa7eXBa26xsKOKij56W0bWZiNuZhoG/34tSCmdubuJ4Omnr6ZS8hqhBcKQripNMskGIuYz0TmWNZh4TOtGP9yfHibQMkrO0mNIb1pDjzl5642IjE/nrob6orqZtRK0Dt9tNc/OpOl5xaWu3282ePXuGHb/00ktHtV9TU0NbWxt1dXW0tbVRXV096pp4hlF+0fg7zZPLZwKoLIYRtEHQaKYRoz9E/+53CfyuC3txDuU73kf+WVULXnso05n8dDIR+WsjHCbgGyDfNbqc6Pbt23nwwQfZsWMHL7/8ckL+euvWrXzxi19MBJKfffZZ7rnnnpT379q1izvvvJNdu3ZxzTXXjLrGNCwQwZFBhhEiYIsarVNBZZ1lpNHMGyRiMvDzk3R8fR+Bgz24Lm+g5gsbKDhbC9HNFBORv/b1dvPwrkdZvf6cUfLXV199NStWrKCpqYlPf/rTicyl8vJyvvzlL7Nx40Y2btzIV77yFcrLywG48cYbiSsv33nnnTz33HOsWrWK5557jjvvvHPU809lGI3jMopLX9uT006zt0TQ8tcazRQQEQJvdNO/+11Mb4j8tZWUbFuOI4208UJivshfhwJ++lpbKCqvoKisfFb6MOQNMdQforLBhW2M8qZWOEzo8GGc9fU4yso43HeY/EgYVySX0vpVGT1Ly19rNLNA2OPD+7PjhN8dwFlXSPlHTyN3xeIsVjNXEREGe7qxOxwUlsze/40RMbE7bGMaAyCl9HU2dyprg6DRTBDDG2Tgmffwv9aJrdBJ6YebKNxUm/XC9prxCQ76ogJ2NbVZkbdORzTDaPznDyufidYy0mjmLFbQwLfHg+83LYDgurQB16VubHn6z2guYlkWg709OPPyyCssmrV+iCWYEYvcTGpZjFghaLVTjWaOIaYwtLeNgedORvcTrK+meOtSHKULP04wn/F7+zANg5Lq2lkN7BuZ1ECIMWqFkGUtI20QNJo0iAjBt3vp3/0uRleAnOXFlH5A7yeYD5iGwVC/l7zCooRExaz1JZy5QUisEJJcRnNuhaCUugq4H7ADD4vIvSPO5wKPAucCPcDHROSEUuoTQHI+2DrgHBE5oJTaA9QBgdi53xOR1Hu+NZosE24ZpH/3cULH+nFU5lPxydXkrS7XKaTzhMHeHhChqKJitruSUdnMOBIzCCQHlbNYU3ncHiql7MA3gW3AauA6pdTqEZd9CugTkSbgG8B9ACLyPRE5W0TOBj4JnBCR5N0jn4if18ZAMxcw+kP0/uAdOh98jUjbEKXbV1LzuXPIX1OhjcEcI53aaSQYJOAboKC4BIfzlExIsvx1d3d34vhk5K/379/P2rVraWpq4rbbbkspcxFvd905a7hs2wW89tpr438oywKlUDYbIpL1FUImYfdNwFEROS4iYeBxYORWvGuA+Lf1Q+AKNfqv5zrg+1PprEYzU1j+CN7dx2n/+734f9dF0cVuam/fSNEFWoRurpLKIIgIAz1d2Ox2CkfsOYjLXy9dunTY8cnIX99yyy089NBDiftSqaDG2335l6/xwD9+M9HuWEiS0mk8w0jFc0+zQCa/6fVAc9J7T+xYymtExAD6gZFrtY8x2iB8Ryl1QCn15RQGBACl1E1KqX1KqX1dXV0ZdFejyRwrbDLwi2ba/m4vg79uoWBdFbX/ZwOl25Zjy9chtrlMKvnr4KCPSDCIq6JylGbR+vXrWbZs2ah2Jip/3dbWxsDAAOeffz5KKa6//vqUktdPPPEEf/SJP0IEzr/g/ES7Y2JZp2Qrhq065k5QOVVPRq6PxrxGKXUe4BeRN5POf0JEWpRSLuC/iLqUHh3ViMhDwEMQ3amcQX81mnER02JoXwcDz5/E8oXJO6Ockq3LcNYWznbXNBmSSv7aCIdRSmF3nkrxTJa/TsVE5a9bWlpwu92jjqdqt642Ond25NgS18WlslMRXSEMn6fPtaCyB2hIeu8GWtNc41FKOYASoDfp/A5GrA5EpCX2r08p9RhR19Qog6DRTCdxqYmBZ9/D6A6Qs7SYkk+cTu4yXbpyKvzi3x6i873j09pm9dIVXHbDTRld63K5+PXzzzHU76XC3YAzN/OU4InKX2cqiy0iUVE7wJFjT3vdMExz1Ka0ubZTeS+wSim1HGghOrh/fMQ1TwI7gReBa4EXJPatKaVswB8CF8cvjhmNUhHpVko5gQ8Cz0/xs2g0YxI82kf/0yeIeAZx1BRQcf1q8s7QmUMLgd6ebi6++GJsNjs2x/BhbbwVwkTlr91uNx6PZ9T1qdo9+d5Jzlm7CZtNpb0uGTEtbLFAeNzwzKkVgogYSqnPAM8QTTv9togcVErdDewTkSeBR4DvKqWOEl0Z7Ehq4mLAIyLJ04dc4JmYMbATNQb/b1o+kUYzgrDHR/8zJwgd8WIvzaXsD0+jYH21lpqYRjKdyU8ncflrEYFQkBd2/5SKhqXY7ROL/UxU/rq8vByXy8VLL73Eeeedx6OPPspnP/vZlO1+4+v385Hf/0NeeumlRLtjYplgH14cJ5tppxl9cyKyG9g94thXkn4OEl0FpLp3D7B5xLEhonsWNJoZI9I+RP+z7xE81IOtwEHJB1ZQtLkOlckGIc2cJyF/vWYNl1x4AV//+tfHNAYPPPAAf/d3f0d7ezvr1q3j6quv5uGHH+bqq69m9+7dNDU1UVBQwHe+8x1guPw1MEz++lvf+hY33HADgUCAbdu2sW3bNgD+5V/+BYCbb76Zq67axn/9x48557y1FBYVJtodi3RZRtlKPtXy15oFR6TTz8Dz7xF4oxuVY8d1UT1F76/XmkPTzFyQv7Ysk57mkyi7nYr6hjnl/gsHDbwdfkqq8jPSMRLLInjoEI7qapzV1QSNIMe8x6g3TEyriIolKzJ6rpa/1mgAoyfAwPMn8R/oRDltUfG5i+qxZSIqppmXDPb2YhoG5TWzq1eUCiM8PKA8LoniOMOvn1MxBI1mrmP0BRn4+Un8r3ag7DaKLqrHdbEbe9HCK2avOUUkGMTf76WguIScvNnVK0qFETax2RQ2e2ZDekK2wj58H8JcyzLSaOYkZn+IgV80M7S3HYCizUtwXdaA3aUNwUJHRBjo7sTmsFNUPvt6RakwIhaOHHvmK5cU0tdAVktoaoOgmXeYA2F8v2xm8OU2sKBwYw2uyxpxlObOdtc0WcLf7yUSClFaUztqR/JcQEQwwiYFxZlPTtIWx4km8E9vB9OgDYJm3mB4Q/h+GVsRWELB+hqKr2hcFPWLNacwIxEG+3rJLSgkdxYL34zFhOMHMFr6ei7uQ9BoZhujN4hvTzND+ztAoPDcGlyXunFUzD2/sWZmiYvXgeCqrJpzgeQ4Rjg6uDsnYBBGSl8n71TOFjohWzNnMboD9P7nYdq/vpeh/R0Ubqyl9vYNlH1klTYGi5TQ0CChoSGKyioYHBpKKX+djmzIX/f29rJlyxbOPGs1f/hHH6bf503Zl5TtjsgyOrVCyF4MQRsEzZwj0umn9/G3af+HfVEp6s1LqPvzjZR9uAlHmXYPLVZM02Cguwtnbi4FJaVp6yGkIxvy1/feey9XXHEFr/zmAJdecin33XffqH6ka3dUcZwYeoWgWZRE2ofoeewtOr6xn8DBHoouqqfujo2Ubl+JvUQHjBc7vu5uxLIorqpBKZVS/nossiF//cQTT/DJT16PGbG4/o9Sy2KnaxfTQtlOZSUNF7fTQWXNIiF0oh/fHg/Bt3ujO4svaaDo/Uv0PgJNguDgIMFBH0XlFThzo5ODVPLXqcim/HVHRwdVFdV4O/y4G+vp7BxdCDJdu5KkYwRJLqMsiklog6CZFcQSgu/04tvjIfzeALYCB8VXNlJ0wRK9s3ge4v3JMcKtQ9PaZs6SQko/tBLLNBno7sSZm0thaVnKa10uFwcOHEh5bjymW/46EooFlHNTB5TT3p+kYwR6H4JmESCmhf93Xfh+6cHo8GMvzaX0Qyso2FiLbSIpeppFw0B31zBXUSp8Ph8XXXRRynPZlL+uqanB09xCZVk1HZ0dVFdXp3xeqnbFOFULAWYny0gbBE1WsMImQ3vbGfx1C6Y3hKOmgLKPnkbBWVW6ZvECoPRDK2ek3eBQzFVUVp5wFcWJy1/Hf57sCmE65a+3b9/Ov3/vu3zh83/O/3tkF9dcM7L8PGnble4ebM5TSROzsQ9B/yVqZhRzKMLA8+/Rfu8r9P/kOPbSXCp2rqbmz86h8JwabQw0aTENg4GuThy5uRSWlY86n5C/PvPMjILKDzzwQGJ2v27dOm688UYArr76alasWEFTUxOf/vSnE5lLyfLXGzduHCV/feONN9LU1MTKlSsT8te33/7n7PnVC5y7eS3PPfccd955JwD79u1LPC9tu6YxxgpBy1+PQstfzx8i3QEGf9uCf18HErHIO70c16VuXapyATGT8tcigre9jXDAT7m7AWfO/MgyC/kj9HcFKK0pIGcCcusiQvDgQRxVVThragDo8nfR6e/kjHCYbsqpXrJ0nFaiaPlrzZxARAi/24/v1y0E3+4Fm6Lg7GpcF9Xr4vWaCREYGCDkH8JVWTVvjAFAZDKSFZBW+jrbaIOgmTJiWgRe78b3mxYiLYPYChy4Lmug6PwlWnlUM2GMcBhfTxe5BQUUFM+vFWUkaODIsWObYHnWkdLXoIPKmnmG5Y8w+Eo7Q//TijkQxlGVT+nvN1F4TjXKqTOGNBNHxKK/sx1ls42ZVTQXiSqcWuQVTSJt2jCA4SsEkeh2NG0QNHMaozuALyk+kNtUSulHVpG3qkwXrtdMicHe3qisdW0ddsf8Gp6MsIWIpN1/MBZpVwgJgziHdiorpa4C7gfswMMicu+I87nAo8C5QA/wMRE5oZRaBrwFvBO79CURuTl2z7nAvwH5wG7gz2Q+RbgXGWIJwSN9DP1PK8HDfYn4QNGFS8hZMjcliDXzi5Dfz5C3j/ziYvLmqKz1WERC0Vn+VAyCSjKCgqCyuj7IwCAopezAN4EtgAfYq5R6UkQOJV32KaBPRJqUUjuA+4CPxc4dE5FUe8q/BdwEvETUIFwFPDXpT6KZESx/hKH9nQy+1IrZE8TmclJ8RSOFm+qwT6D4h0YzFqZh0N/ZjiMnB1dF1Wx3Z1JEQiZ2hw27YxKp1CNqIUDcZRTXNcoOmfR8E3BURI6LSBh4HBi52+IaIK4N+0PgCjWG808pVQcUi8iLsVXBo8CHJ9x7zYwRbhui77+P0HbPK/T/7Dh2Vw7l151O3R2bKL5yqTYGmmlDROjvbEfEilZAs2U2oE5U7XQm5a8tyyISMvH5+9myZQurVq1iy5Ytic1nYz0DTq0Q9h84kGj3y7d/OePPNl1k8s3XA81J7z2xYymvERED6AfihU6XK6VeU0r9Uil1UdL1nqT7U7UJgFLqJqXUPqXUvq6urgy6q5ksYlr4X++i819/R+f9r+J/rZOCs6upvm091TefFd1VPJnZj0YzBoN9vYQDAYorq3FMIMV0Lslf7/7ZbixTeOCb/8gVV1zBkSNHuOKKK7j33nvHfEYCw0TZ7Nx6662Jdo8fPc6vf/7r2AXZcR1l8tedqicjVzDprmkDGkVkPfB54DGlVHGGbUYPijwkIhtEZENV1fxcSs51zP4QA8+/R9t9e+l97G3M/jAlVy+n7i82UfaRVTpGoJkxQv4hhvp6yXcVk+8qntC9c0n++kc/ispc/2z3T9m5cycAO3fuHCaLneoZccQ0aevpGdbuR677CD/f/fMJfSdTJZOgsgdoSHrvBlrTXONRSjmAEqA35g4KAYjIfqXUMeC02PXupPtTtamZQcQSgof7GHq5LbqJTCB3VSlFv99E3vvKdbaQZsYxjQj9nR3RuEHlxCd7c0r+2tOCsik6Ojuoq6sDoK6uLiF/na6t+LViGrT1dA9rt3ZJLR1tHdHzGX8rUyMTg7AXWKWUWg60ADuAj4+45klgJ/AicC3wgoiIUqqKqGEwlVIrgFXAcRHpVUr5lFKbgZeB64H/Oz0fSTMWZn+Iob3tDO3twOwPYSty4rrETeHGWl2WUjNpnnrqKdrb2yd0jxEKISI4cnNT7jeora1NaASNx2zLX4uMnV00nmQ2pomM+A4EOXVNlvZjjGsQRMRQSn0GeIZo2um3ReSgUupuYJ+IPAk8AnxXKXUU6CVqNAAuBu5WShmACdwsIr2xc7dwKu30KXSG0YyRbjVQ8sEV5J9RruMCmqxjRiJYYuHIyZmWzWezKX998r2T1FTVkpPnoKamhra2Nurq6mhra0vIX6d7RhwxTdz19cPabWtpo6Y2qmuUtXx8EZk3r3PPPVc0mRPxBqX/uRPS+rWXpfmOX0nL37wo3qeOS6TbP9td0ywADh06NKn7Br190nb0sPh6uqf0/O7ubmlsbJzwfUuXLpWurq7E+5/+9Kdy1VVXiWVZ8uKLL8rGjRtFRKSnp0eWLVsmvb290tvbK8uWLZOenh4REdmwYYO8+OKLYlmWbNnye/K97/ynREKGfOELX5B77rlHRETuueceuf3228d8RpzAoUMSbmkZ1u4lV14i3/6Ph0VaXpXW1uaMP1+q/xeik/dxx9hZH+Qn8tIGYXyssCFDBzqk8+HXpfnOX0nzHb+Szodfl6HXu8SKmLPdPc0CYjIGIegfkrZjh6W3rUUsy5pyH6677jpZs2aNfOELXxj32vvvv1/q6+vFbrdLXV2dfOpTnxIREcuy5NZbb5UVK1bImWeeKXv37k3c88gjj8jKlStl5cqV8u1vfztxfO/evbJmzRpZsWKFfPpP/pd0vtcvlmVJd3e3XH755dLU1CSXX355woCM9YyzzjpL/G+8IeGOjmHtfvLGT8q7fcdiBsGT8XcyFYOg5a8XACJCxDPI0L52/L/rQoIm9tJcCs6ppvDcGh0b0MwIE5W/NiIReluasdntlC9xY5tlZc/pQEToaRnEmeugpGpyf2dWOEzo8GGcS5bgKD9V9+Go9yg5ykHjUB+tqpoldSkz80eh5a8XKaYvjP/VTob2d2B0+sFho+DMCgo21JC7olRnCmnmDJZp4m1vBRFKa+oWhDEAMA0LyxRy8qbweVLIVkBsp/Jc1DLSzB3EsAi81Yt/fwfBw71gQU6ji9I/aKJgXRW2CRTl0GiygYjg7WjHjEQorV2CI2fh7HIPB6ODuXMKBkHiSqcpDIIty2qvevSYB4glhE/04z/Qhf/1biRoYCvOwXWxm4JzanBWF8x2FzWalIgIA12dhAN+SqpryC1YWL+rkaCBzT5J/aIYcYPASIOAJNYFc2kfgmaWiLQPMfRaJ4EDXZj9IVSOjfzVFeSvryavqQxl1y4hzdxmyNtHwDdAUVn5hHciz3VEhHDQJDffMaXUWUlRCwHAEitJ7VS7jBYlhjeI/0AXgQOdRNr9YIO8VWWUbFtG3uoKbBMtzafRzBIBn4/B3h7yXS4Ky8rHv2GeEQmZiCXk5E9xGDUMsNlGGYRhG9OyhDYIcwBzKELgYDf+1zoJvzsAxOIC16wkf20l9qKF43PVLA5C/iEGujrIyc+nuLJ6XlU+y5RwIDqzn6pBEMNA2UfHD2SY/PXcEbfTzACWP8LQ3na6vv0mbV99Ge9/H8UajFC8ZSm1t2+g+tazozWJtTHQzDPCgQDe9jYcOTmU1tShMpSzniizLX+96YJz2Xzpev73//6zhDvNgZgAABiZSURBVDRFb2/vhOWvH/3BDzhz6+8Ne0a8nnL8m+vr60vZ7rSTyWaFufKa7xvTzKGwDO5tk85H3pDmv/i1NN/xK2m97xXx7j4uIY9vWjbqaDTZItUGqHAwKB3Hj0rXyRNiGJEZff67774ra9asyfj6V199Vd59991RO5V/9rOfDdtFvGnTJhGJ7lRevny59PT0SG9vryxfvlx6e3tFRGTDho3y0/96Tga9Abnqqqtk9+7dIiJy++23D9up/Od//ufjPmOZ2y1tB3437BmGacibXW9Kl69FpOVVufmWW1O2m4qpbEzTK4QZxvJHGNrXQfd33qT1qy/T98MjGN0Bii6qp/ozZ1N7+wZKti0np75oQS6rNYsHIxymr60FZbNRVrsEu31mPdKzKn/d38/GczeRm+/k+uuvHyZzPRH562eeeYYrLriA8sqKYc+IrxDiLqNnn3kmZbvTjY4hzADmYJjgW70E3uwmeNQLpmAvy6XownoK1lXi1IO/ZoFhRKLGAKCsrh670znjz5xN+eu62iXRcplOW+I4QEfHxOSvPR4P9dXViT0I8eMicYMQpbu7K2W70402CNOE0RMgcKiHwMEewu8NgBAzAksoWFuF062NgGZhYkTCHHz9LwmEjuLIyaW5d3ocD66iMzjttMzKSGZT/to0LSwrGkyO/02P97ed9hmWFf05aQ+CUiophqCzjOYFIkKkdYjAoR6CB7ujKaKAs64Q1+WN5K+pwFlXqI2AZkFjRML0tUZntI6c3BkLII9HNuWvqytraWtrIbfAMex6YMLy1/W1tbywbx8qtqKKP8OSmKGIXV9ZWZWy3elGG4QJIKZF6MQAwYM9BA71YHpDoCBnWTElH1hB/upyLSSnWTRYppEwBmvWfhVnbub1kKcDl8uFz+dL/DzZFcL27dt58MEH2bFjBy+//DIlJSXU1dWxdetWvvjFLyYyep599lnuuece7GY+RUVF7H9tL5s3b+bRRx/ls5/9bKKtXbt2ceedd7Jr1y6uueaaMZ/xe5dcwl/+5V/i9fmwmWbiGSNjCFdu3Zqy3Wknk8jzXHnNRpaR4QvJ4L526f7eIfH81W+l+Y5fSfNf/lq6/u1NGXylTQxfKOt90mhmm6733pWXf/0r6Xj3mISDwVnrR7blr03Tks73BmTPc79JyFT/6Z/+aSJDcKLy15HubvnW3XePktgeCg/JH3ziD2TPr54WaXlVDhx8O2W7qdDy19OIWEKkbYjg270E3+4l7PGBgM3lJO995dHXaWXYxiiXp9EsZFreeYsf3ffXbPqTz7B+0yacOdldGcwmQX+Ega4ApdUFU9+hDEQ6OjC6u8lbvXqYe3koMsSJ/hMsy6+isL+F92xultZmVnday19PEStkEDriJfB2L8F3+rB8YVDgdLsovqKRvNPLcS4p0nLSmkXP8df28pN/vBdXRQWFpWWLyhgAhIYMlE1NSd00GYlEUI7RWkgjYwhay2gGEUuItA8RPNxH6HAfofcGwBRUnp280/5/e3ceJVV1J3D8+6uqrup9oWmWpmloNpVFiI1gosbRCOJKjEbJppkkkknEzHhmNFFPcpIcM2Ni5sxozCRxIRNNPIxLFiY6QRAhLkFBAVltmr27gW6a3qur6tWr3/zxXkPTVNMNdFX1cj/n1Km33Pfevb2833v33ndfgXsnUGCeEjaMTna8+QZ/+cV/Mrx0PLc88AP21xxKdZaSKhZTIu1RAlnnNphdZ2pZxxuUT1repdtpsvQqIIjIAuAxwAs8raqPdFkfAJ4FyoF64HZV3Sci84BHAD8QAe5T1dXuNmuA0UC7u5v5qpqYzrU4L5MJVTYSrmggtKuBWKsFQNqoLLIvHUPG+QX4x+UiXvOsnmF0pqq8+/v/4e0XfsvYqTNYeN93nWGsh1hACActVJX0rL57xkItC0/GqR1RunY77TfDX4uIF/g5MA+oAtaLyHJV3d4p2VeBBlWdJCKLgB8DtwNHgRtVtUZEpgMrgM7vgfuCqib8nZj1z++g/UNn7BJPlo/A5ALS3Y8319wFGEZ3opEIr/3qcXa8tYYLLr+S+YvvOekFN6rJH5EzVUJtFl6fh7Q+aj9UVWdgO9+pp+GOKiPPGQ5/fa5twr25Q5gDVKrqHgARWQYsBDoHhIXA993pl4AnRERUdWOnNNuAdBEJqGr4nHJ9hgIT8kkbnUX65ALTFmAYvRRsauSPP32YQxU7ufT2LzH35ttOOvmnp6dTX19PYWHhoA8KthXDCtlk5QX6rqy2DbFY3Cqj420IZ3AsVaW+vp709PSzzlJvAsIY4GCn+SpgbndpVDUqIk1AIc4dQodbgI1dgsGvRcQGXgYe1gR1ecq+ZHQidmsYg9aRvbtZ/u//SrCxgRvv/Q5TLrnslDQlJSVUVVVRV1eXghwmV7g9SiQYJaslgOdQ37UfROvq8EYieLr8DFsiLbREWhB/OxKs56jHpr2h5xr19PR0SkpKzjpPvQkI8Urf9cR92jQiMg2nGml+p/VfUNVqEcnBCQhfwmmHOHnHIouBxQClpaW9yK5hGOdiyxuv8fozvyAjJ5fbv/8IoyZNiZsuLS2NsrKyJOcu+WIx5Xff+xs5hel8+t7un3g+Uy2rVlG15B7Gv/giGV26iT7+weMs3bGUjTPuQ1Z8ky/nPMV///NtfXbs7vSmBbUKGNtpvgSo6S6NiPiAPOCYO18C/AG4Q1V3d2ygqtXudwvwPE7V1ClU9UlVna2qs4uKetcP1zCMMxeNRFjxy8d57ZePM+a8qXzpx493GwyGkgNb62k+GmL6J8/+yjueyAGn4sVfOvaUdcFokAxfBuJeV8ckOZ1denOHsB6YLCJlQDWwCPh8lzTLgTuBvwG3AqtVVUUkH3gFeEBV3+5I7AaNfFU9KiJpwA3AqnMujWEYZ6XhUDV/fuwn1O7dzdybb+cTt30ej8c8fAnw4ZoqsvIDlM0a3qf7jRzYjycvD29e3inr2qPtZPgyIGYDEOsvzyG4bQJLcHoIeYGlqrpNRH6I8zj0cuAZ4DkRqcS5M1jkbr4EmAR8V0Q6hi2cD7QBK9xg4MUJBk/1YbkMw+gFVWXL6hW88Zun8PnS+PT932ViedcmwqGr4XAbB7cfY+5NZXj7uEu6deAg/m6qwdstNyBoR0DoP3cIqOqrwKtdln2v03QI+Gyc7R4GHu5mt+W9z6ZhGH0t2NzEa7/6Gbs3rKN0+kwW3H0vOcP69ip4oNuyphqPT5h62ZieE5+hyMGDZMyYEXfd8TsEt7eRSnLu1obkk8qGMdRVbniXVU89Qai1hb+742tcdO1NKRu6ur8KNkfY8XYNUy4eSWYfP6+kloVVU0Pu9dfFXX+iysgNCP3pDsEwjMGhrbGB1b/+FRXr3qKodDy3PPhDisYN/p5CZ2Pz6weIRmOULxjf5/u2amrAtvGXjou7vj3aTlZaVv+sMjIMY2BTVbauWcna554hGolw2aI7mH3jZ/DGeUrWcJ5K3rKmmsnlI8gfmdnn+4/s2weAf1z8NoRgNEhRZtHxKqP+1MvIMIwB7MieSt74zZNU79xOyQXTmbd4CcOK+7YL5WCz+fWDWGGb8mvHJ2T/oYoKAAKTJsVdf6LKKAqAmoBgGMa5CDY18tayZ9nyxkoysnOYt/geZlw5z7QV9KC1IcymlQeYVD6CwjHZCTlGuGIXvlGj4nY5hU4BwYoAYEtyxlwzAcEwBhkrEmbTildY9/IyopEw5dfdxCW3fI70rMSc3Aabd5fvJqbKx2+emLBjhCsqCEyZHHedqtIaaSU7LRtCIWeZ6WVkGMaZsKMWW1av5N3fL6O14Rhls8q54o6vUTjm1CdhjfjqDrSwc91hPnZ1KbnDE/N+dLUsInv2kH35qeNDAYTsEJFYhNxALjQ0YZGGJ0nD8puAYBgDXMy22fHWGv720vM01R6h+LypXPet+xg7NX4fdyO+mB3jjd/uJCM7jfJr4/f+6QuR/ftRyyIwJf6wIM3hZgBy/blgW1iShjdJo8magGAYA5QVCrHljZW8/8ofaK6rZUTZRD7znW8wflb5oB+OOhE2v15F3YEWrrlrOoHMvnsJTlftW7cCkN5lQLsOzREnIOQF8sCOEMWHJ0lD9puAYBgDTLC5iU0rXmHjij8TammmeMoFXPnlrzOxfI4JBGep4XAb7/7vHspmDmfiRYkdRLN94yY8OTn4J8Zvo+gICM4dQoSo+MwdgmEYJ6gqNRU72bzyVSrWvYVtWUwon8Ocm25lzPl9NyTzUGRFbFY8tZW0gJcrPndewoNq+8aNZMya1W1vr6ZwE9Bxh2ARxYfX3CEYhhFqa+Wjd/7K5tdepe7APvwZGcy4aj6z5l9PYYl5P8i5UlXeXFZBfU0bNy6ZSVZ+IKHHs1taCO/aRc4187tN0/UOwRJTZWQYQ5Ydtdi78X22v7maPR+sx7YsisZPYN5dSzj/sivwpyem98tQtGnVQXa8c4jya8dROq0w4ccLrt8AqmRedFG3aY43KgfcKiPS8CapJtAEBMPoB+yoxcFtW9j13jtUrHubUGsLGbl5XHj1AqZefhUjJ0wy7QN9bNf6I7zzciUTLypi7o0TknLM1rVr8WRmklHe/WDPDeEGfB6f8xyCbWGZKiPDGPzCwSD7Nr9P5fp17PlgPZH2IL5AgEmzL2Hq5VdSOmOWGWsoQXZtOMKqX29n9KQ8rv77qUgSTriqSuvatWRd+gk8/u6fPD7cdpiRmSPxiMepMsKHxzQqG8bgErNtDu/exYEtm9i/ZRM1FTuJ2VEycvOYcsllTLr4EkpnzCTNn9h67KFu25vVrHn+I0ZPzOOGu2fiS0vOU8ChrduIHj5M9pK7T5vuSPAIIzNHOjNuQDB3CIYxwNlRi9q9e6ip2EHVjm0c3PYh4WAbACPKJlJ+/UImXHQxxeddYF5XmQR2NMbbL+5iy9pqSqcNY8HXZ5DmT97PvfHll5BAgJz53Tcog3OHML1wujNjW0TMcwiGMbCoKs11R6jdu4dDlR9RU7GDI7sribqDk+UWjWTy3EsZd+EsSqfPJDM3/qBmRmLUV7ey+tkd1O5vYdbVY/n4zROTNhwEgN3aRvOfXyHnmvl4c3O7TaeqHGk7wtWlVzsLIm2ESTfPIRhGfxUJtdNQU03dgX3U7dtD7f491O3be/zq3+P1MXLCRGbOv47i8y6gePL5ZA9LfA8W41ThoMUHrx1g08oD+DN8XHPXdCaVj0h6Phqee5ZYayvDvvjF06araashEotQkuMOTx5qoplCU2VkGKlkRcK0HK2jqfYIDTVVHKuppuGQ8916rP54Ol8gQFHpeM6/9JOMGD+RovFlFJWW4TtNo6GReMHmCNverGbz6wcJB6NMmTuSy26dTEZO8n8v1qFD1D+zlOyrriLjwgtPm3bnsZ0AnDfsPGdBqIkWsklWB7NeBQQRWQA8BniBp1X1kS7rA8CzQDlQD9yuqvvcdQ8AXwVs4FuquqI3+zSMRFBVIu3tBJsaaGtsINjUSGvDMZrramk+WktzXR3NR2tpb246abtAZhYFxWMonXYhBcUlDCseQ+HYcRSMLjb1//2EFbY5sK2eyvdr2bOpjpitjJtRyNybJlA0NicleYqFQlTfdx8aizHy2/f3mH57/XYEYXL+ZIjZEG6i2ZvZf6qMRMQL/ByYB1QB60Vkuapu75Tsq0CDqk4SkUXAj4HbRWQqsAiYBhQDq0SkY4i/nvZpGN1SVaJWBKu9nVBbK6HWVsJtrc50Wyvh1hPTodZWNwA0EmxqJBoJn7I/nz9AzvAicocXMaJsArmFReQWjSC3aAQFo8eQmZdvngPoR1SVYHOEY9VtHNrdyKHdTRza3YRtxUjPSmPGFSVM+2QxBaOyUpZHq7aWmvvup/39Dyh+9FH843oeQXXtwbXMLJpJZlomBI8B0KKZ/arKaA5Qqap7AERkGbAQ6HzyXgh8351+CXhCnP+ehcAyVQ0De0Wk0t0fvdinkWIaixGL2cRiMdR2vmO2fXy5xmLE7Bgas098x2Lu8hPb2dEodtTCtiyi7rdtWceXdUxHOy2LhsNEQiGscAir47vTdCQUAtXT5j8tkE4gO5v0rGyy8gvIH1VMVn4BmXn5ZOXlk5lfQFZePlkFw8jIyTUn/BRRVWIxJWYrdiRGJBQlErKJhKJYIZtQm0WwKUJbc5hgY5jm+hCNR4KEg87rJUWgsCSbaZcXM2FmEaMn5SW1wfh4OWyb6OHDhHbtonXtWpr+tBxiMYof+Tfybri+x+3/svcvfNTwEQ/OfdBZcNR5zWZltIiR/uTU7vfmKGOAg53mq4C53aVR1aiINAGF7vJ1XbYd4073tM8+87MvfIWoHYqz5vQnlDMTf19nd4SetnLXn3RC1C7bqTvXNU287+6WJZsX8CB4QHwIzqdjGvEjZCL48Hvd9eJD8DvrxI/Q8Z2GiBeCEAlCpA4ajh+n3f0cSlVBBzA3aJ4UPAU9ZbkcT6sd8+KmFe+JD17w9O5kJ7EwPquZNKuBzPBhCsKH8YcOk9G+H++HIaJAhfvpU+qWxv3XEHe+Y7kvCgEL/BHoCEOWDzZPEVbOFY7WPgRPPRR31zGBGNDugUavMCmszH7hR+znRxRoI5l4eC8ynsV56X1dqrh685uId9nU9azRXZrulscL33HPRCKyGFgMUFp6doN5eX0eNJb4H6jELW7fH+XEpMRdflI+pMs2cabllOWd58Xdn/uRLvPx0nTMHz+2e5LH4y7vMo/H2eacr9CjzkeDzmwqY9ug1tMFhRsCTrlgcb5FbcCO+y0aRTSEJxbq9N2O125CNHzKX6kdgNYAJLp/jMqJb4VOwQ2iXrD8QiQNWnI8HCvwcGi0F9sn+IBRp/k79MSc/wJ/VCiN+rkslE1LpnN6rMXH5uzL+eSwGXz6Y2O630kf6s1PsQro/A6+EqCmmzRVIuID8oBjPWzb0z4BUNUngScBZs+efVb/4t/8zdNns5lhGEZKXZzk4/Wmom09MFlEykTEj9NIvLxLmuXAne70rcBqVVV3+SIRCYhIGTAZeK+X+zQMwzCSqMc7BLdNYAmwAqeSd6mqbhORHwIbVHU58AzwnNtofAznBI+b7gWcxuIocLeq2gDx9tn3xTMMwzB6S7SHnhr9yezZs3XDhg2pzoZhGMaAIiLvq+rsntIlv2+WYRiG0S+ZgGAYhmEAJiAYhmEYLhMQDMMwDMAEBMMwDMM1oHoZiUgdsP8sNx8OHO3D7AwEpsxDgynz4Heu5R2nqkU9JRpQAeFciMiG3nS7GkxMmYcGU+bBL1nlNVVGhmEYBmACgmEYhuEaSgHhyVRnIAVMmYcGU+bBLynlHTJtCIZhGMbpDaU7BMMwDOM0Bn1AEJHPisg2EYmJyOwu6x4QkUoR+UhErklVHhNJRGaJyDoR2SQiG0RkTs9bDXwico/7e90mIj9JdX6SQUT+RURURIanOi+JJiKPishOEflQRP4gIvmpzlOiiMgC92+5UkS+k8hjDfqAAGwFPgP8tfNCEZmKM0z3NGAB8F8i4k1+9hLuJ8APVHUW8D13flATkStx3tF9oapOA36a4iwlnIiMBeYBB1KdlyRZCUxX1Qtx3pr5QIrzkxDuOennwLXAVOBz7rkrIQZ9QFDVHar6UZxVC4FlqhpW1b1AJTAYr54VyHWn8+jmzXSDzDeAR1Q1DKCqtSnOTzL8B3A/Q+TFoar6mqpG3dl1OG9dHIzmAJWqukdVI8AynHNXQgz6gHAaY4CDnear3GWDzT8Bj4rIQZwr5UF5JdXFFOByEXlXRNaKSLLfRJhUInITUK2qm1OdlxT5CvB/qc5EgiT1PJXYN1MniYisAkbFWfWQqv6pu83iLBuQV1enKz/wKeBeVX1ZRG7Debvd1cnMXyL0UGYfUABcgvNa2hdEZIIO4C51PZT3QWB+cnOUeL35vxaRh3Dexvi7ZOYtiZJ6nhoUAUFVz+YEVwWM7TRfwgCtTjld+UXkWeAf3dkXgaeTkqkE66HM3wB+7waA90QkhjMWTF2y8tfXuiuviMwAyoDNIgLO3/EHIjJHVQ8nMYt9rqf/axG5E7gB+NRADvY9SOp5aihXGS0HFolIQETKgMnAeynOUyLUAFe401cBu1KYl2T5I05ZEZEpgJ9BOhCaqm5R1RGqOl5Vx+OcQC4a6MGgJyKyAPg2cJOqBlOdnwRaD0wWkTIR8eN0hFmeqIMNijuE0xGRm4GfAUXAKyKySVWvUdVtIvICsB3nlvNuVbVTmdcEuQt4TER8QAhYnOL8JMNSYKmIbAUiwJ2D+ApyqHoCCAAr3Tujdar6D6nNUt9T1aiILAFWAF5gqapuS9TxzJPKhmEYBjC0q4wMwzCMTkxAMAzDMAATEAzDMAyXCQiGYRgGYAKCYRiG4TIBwTAMwwBMQDAMwzBcJiAYhmEYAPw/irLwtdAwRV8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x737d14e9d0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(2)\n",
    "for t in np.logspace(-2,6,9):\n",
    "    plt.plot(y,sigma(t,y),label='t='+str(t))\n",
    "\n",
    "#plt.plot(y,sigma(1000000*t,y),label='t='+str(1000000*t))\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2 Bilayer on a substrate - Heat Transfer\n",
    "\n",
    "$T_{,yy}=\\frac{1}{\\kappa}T_{,t}$ **(10)**\n",
    "\n",
    "$\\lim_{\\epsilon \\rightarrow 0} k_1 T_{,y}(-h+\\epsilon,t)= k_2 T_{,y}(-h-\\epsilon,t)$ **(11)**\n",
    "\n",
    "$T_{equil}(y,t)=A(\\eta)T_{surface}(t)+B(\\eta)T_{substrate}(t)$ **(12)**\n",
    "\n",
    "where \n",
    "\n",
    "$A(\\eta)=1+c_1 \\eta,~\\eta=-h/H...0$                       \n",
    "\n",
    "$A(\\eta)=c_2(1+\\eta),~\\eta=-1...-h/H$                     \n",
    "\n",
    "$B(\\eta)=c_1 \\eta,~\\eta=-h/H...0$\n",
    "\n",
    "$B(\\eta)=1-c_2(1+\\eta),~\\eta=-1...-h/H$\n",
    "\n",
    "$c_1=\\frac{k_2/k_1}{1-h/H(1-k_2/k_1)}$\n",
    "\n",
    "$c_2=\\frac{1}{1-h/H(1-k_2/k_1)}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$T(y,t)=T_{equil}(\\eta,t)+\\delta T(\\eta,t)$ **(12)**\n",
    "\n",
    "where\n",
    "\n",
    "$\\delta T (\\eta,\\tau) = \\sum_{n=1}^{M} a_n(\\tau)f_n (\\eta)$ **(A.5)**\n",
    "\n",
    "note: $\\tau=\\frac{\\kappa_2 t}{H^2}$\n",
    "\n",
    "eigenvalues:\n",
    "\n",
    "$\\sqrt{\\frac{k_1}{k_2}}\n",
    "\\tan\\left(\\sqrt{\\lambda_n}\\left(1-\\frac{h}{H}\\right)\\right)\n",
    "+\\tan\\left(\\sqrt{\\frac{k_2 \\lambda_n}{k_1}}\\frac{h}{H}\\right)=0~~(n=1,\\infty)$\n",
    "\n",
    "eigenfunctions:\n",
    "\n",
    "$f_n(\\eta) = -\\mu_n \\sin(\\sqrt{k_2 \\lambda_n / k_1}\\eta),~ \\eta=-h/H ... 0$ **(A.4)**\n",
    "\n",
    "$f_n(\\eta) = -\\sin(\\sqrt{\\lambda_n}(1+\\eta)),~ \\eta=-1 ... -h/H$ **(A.4)**\n",
    "\n",
    "$\\mu_n = \\frac{\\sin(\\sqrt{\\lambda_n}(1-h/H))}{\\sin(\\sqrt{k_2 \\lambda_n / k_1} h/H)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "62 eigenvalues range between 13.983058561 and 48613.7195126\n"
     ]
    }
   ],
   "source": [
    "k1k2=2 # ratio of k1/k2 the conductivity of the dense layer and TBC layer\n",
    "k2k1=1/k1k2\n",
    "hH=0.4 # ratio of dense thickness to total thickness 0 is no dense layer, 1 is all dense layer\n",
    "eigprob=lambda Ln: np.sqrt(k1k2)*np.tan(np.sqrt(Ln)*(1-hH))+np.tan(np.sqrt(Ln/k1k2)*hH)\n",
    "Lmax=50000\n",
    "l=np.arange(1,Lmax,1)\n",
    "plt.figure(3)\n",
    "plt.plot(l,eigprob(l),'.')\n",
    "plt.plot(l,np.zeros(np.shape(l)))\n",
    "plt.ylim((-10,10))\n",
    "plt.title('zeros are eigenvalues')\n",
    "i_cross=np.where(np.diff(np.sign(eigprob(l))))\n",
    "\n",
    "# Solve for eigenvalues between 1 and Lmax\n",
    "Ln=fsolve(eigprob,l[i_cross[0][1::2]],xtol=1e-11)\n",
    "print(str(len(Ln))+' eigenvalues range between '+str(Ln[0])+' and '+str(Ln[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fn(Ln,eta):\n",
    "    eta=np.atleast_1d(eta)\n",
    "    mun=sin(Ln**0.5*(1-hH))/sin(Ln**0.5/k1k2**0.5*(hH+eps))\n",
    "    fn=-mun*sin(Ln**0.5/k1k2**0.5*eta)\n",
    "    fn[eta<= -hH]=sin(Ln**0.5*(1+eta[eta <= -hH]))\n",
    "    return fn\n",
    "    #if eta>=-hH:\n",
    "    #    return -mun*sin(Ln**0.5/k1k2*eta)\n",
    "    #elif eta<-hH:\n",
    "    #    return sin(Ln**0.5*(1+eta))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5637808662ad479f9192924ac7cd3ffa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f59967d4be0>]"
      ]
     },
     "execution_count": 194,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau=np.linspace(0,1)\n",
    "eta=np.linspace(-1,0,1000)\n",
    "plt.figure()\n",
    "plt.plot(eta[:-1],fn(Ln[0],eta[:-1]))\n",
    "plt.plot(eta[:-1],fn(Ln[1],eta[:-1]))\n",
    "plt.plot(eta[:-1],fn(Ln[20],eta[:-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\delta T_{,yy}-\\frac{1}{\\kappa}\\delta T_{,t}=\n",
    "\\frac{1}{\\kappa}\\left(A(\\eta)\\frac{dT_{surface}}{dt}+B(\\eta)\\frac{dT_{substrate}}{dt}\\right)$\n",
    "\n",
    "$a_n(\\tau)= \\frac{e^{-\\lambda_n \\tau}}{C_{n}}\n",
    "\\left(A_n\\int_0^{\\tau}e^{\\lambda_n \\xi}\\frac{dT_{surface}(\\xi)}{d\\tau} d\\xi+\n",
    "B_n\\int_0^{\\tau}e^{\\lambda_n \\xi}\\frac{dT_{substrate}(\\xi)}{d\\tau} d\\xi\\right)$\n",
    "\n",
    "$A_n=\\int_{-1}^{0} A(\\eta)f_{n}(\\eta)d\\eta$\n",
    "\n",
    "$B_n=\\int_{-1}^{0} B(\\eta)f_{n}(\\eta)d\\eta$\n",
    "\n",
    "$C_n=\\int_{-1}^{0} f_{n}^2(\\eta)d\\eta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [],
   "source": [
    "def AB(eta):\n",
    "    '''Define eta-functions A(eta) and B(eta) for T_equil'''\n",
    "    eta=np.atleast_1d(eta)\n",
    "    c1=k2k1/(1-hH*(1-k2k1))\n",
    "    c2=1/(1-hH*(1-k2k1))\n",
    "    A=1+c1*eta\n",
    "    B=-c1*eta\n",
    "    A[eta< -hH]=c2*(1+eta[eta < -hH])\n",
    "    B[eta< -hH]=1-c2*(1+eta[eta < -hH])\n",
    "    return A,B\n",
    "\n",
    "def dAn(Ln,eta):\n",
    "    '''Define dA_n/deta for integration'''\n",
    "    A,B=AB(eta)\n",
    "    return A*fn(Ln,eta)\n",
    "\n",
    "def dBn(Ln,eta):\n",
    "    '''Define dB_n/deta for integration'''\n",
    "    A,B=AB(eta)\n",
    "    return B*fn(Ln,eta)\n",
    "\n",
    "def dCn(Ln,eta):\n",
    "    '''Define dC_n/deta for integration'''\n",
    "    return fn(Ln,eta)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).\n",
      "  \n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "26b01c727198440f9fbc5894d9b03c58",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5996551c50>]"
      ]
     },
     "execution_count": 199,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A,B=AB(eta)\n",
    "plt.figure()\n",
    "plt.plot(eta,A*10+B*0)\n",
    "plt.plot(eta,1*(1+eta)*10+(1-(1+eta))*0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [],
   "source": [
    "def An(Ln):\n",
    "    An=integrate.quadrature(lambda x: dAn(Ln,x),-hH,0,maxiter=200)\n",
    "    An+=integrate.quadrature(lambda x: dAn(Ln,x),-1,-hH,maxiter=200)\n",
    "    return An[0]\n",
    "\n",
    "def Bn(Ln):\n",
    "    Bn=integrate.quadrature(lambda x: dBn(Ln,x),-hH,0,maxiter=200)\n",
    "    Bn=integrate.quadrature(lambda x: dBn(Ln,x),-1,-hH,maxiter=200)\n",
    "    return Bn[0]\n",
    "\n",
    "def Cn(Ln):\n",
    "    Cn=integrate.quadrature(lambda x: dCn(Ln,x),-hH,0,maxiter=200)\n",
    "    Cn=integrate.quadrature(lambda x: dCn(Ln,x),-1,-hH,maxiter=200)\n",
    "    return Cn[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [],
   "source": [
    "def T_eqbm(eta,t):\n",
    "    A,B=AB(eta)\n",
    "    eta=np.atleast_1d(eta)\n",
    "    t=np.atleast_1d(t)\n",
    "    Ti_surf=10\n",
    "    Tf_surf=0\n",
    "    t0_surf=0\n",
    "    T_surf=Ti_surf+(Tf_surf-Ti_surf)*t/t0_surf\n",
    "    T_sub=Tf_surf\n",
    "    T_surf[t>t0_surf]=Ti_surf\n",
    "    T_eqbm=A*T_surf+B*T_sub\n",
    "    return T_eqbm\n",
    "\n",
    "def deltaT(eta,tau,Ln):\n",
    "    ETA,TAU=np.meshgrid(eta,tau)\n",
    "    dT=np.zeros(ETA.shape)\n",
    "    for ln in Ln:\n",
    "        A=An(ln); B=Bn(ln); C=Cn(ln)\n",
    "        dT+=(A+B)/C*(exp(-ln*TAU)-1)*fn(ln,ETA)\n",
    "    return dT,ETA,TAU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.47476359996224465"
      ]
     },
     "execution_count": 185,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Cn(Ln[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e83a83617a884d3eb9c708e3e0ad0bc4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f5996adbb38>]"
      ]
     },
     "execution_count": 188,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tau=np.linspace(0,1)\n",
    "eta=np.linspace(-1,0,1000)\n",
    "dT,ETA,TAU=deltaT(eta,tau,Ln)\n",
    "\n",
    "plt.close(5)\n",
    "plt.figure(5)\n",
    "plt.plot(ETA[-5,:].T,dT[-5,:].T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 177,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "04ce0b444f2942bea3cbc4d06b455ad1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f5996ee14e0>"
      ]
     },
     "execution_count": 177,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib import cm\n",
    "\n",
    "plt.close(6)\n",
    "fig = plt.figure(6)\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot_surface(ETA,TAU,dT,cmap=cm.coolwarm,\n",
    "                       linewidth=0, antialiased=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: \n",
    "\n",
    "Solving a PDE with initial condition $u(x,0)=f(x)$\n",
    "\n",
    "for \n",
    "\n",
    "$\\begin{align*}& \\frac{{\\partial u}}{{\\partial t}} = k\\frac{{{\\partial ^2}u}}{{\\partial {x^2}}}\\\\ & u\\left( {x,0} \\right) = f\\left( x \\right)\\hspace{0.25in}u\\left( { - L,t} \\right) = u\\left( {L,t} \\right)\\hspace{0.25in}\\frac{{\\partial u}}{{\\partial x}}\\left( { - L,t} \\right) = \\frac{{\\partial u}}{{\\partial x}}\\left( {L,t} \\right)\\end{align*}$\n",
    "\n",
    "$\\begin{align*}{A_0} & = \\frac{1}{{2L}}\\int_{{\\, - L}}^{L}{{f\\left( x \\right)\\,dx}}\\\\ {A_n} & = \\frac{1}{L}\\int_{{\\, - L}}^{L}{{f\\left( x \\right)\\cos \\left( {\\frac{{n\\,\\pi x}}{L}} \\right)\\,dx}}\\hspace{0.25in}n = 1,2,3, \\ldots \\\\  {B_n} & = \\frac{1}{L}\\int_{{\\, - L}}^{L}{{f\\left( x \\right)\\sin \\left( {\\frac{{n\\,\\pi x}}{L}} \\right)\\,dx}}\\hspace{0.25in}n = 1,2,3, \\ldots \\end{align*}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mechanism maps for deep delaminations\n",
    "\n",
    "Impose $G=\\Gamma_{C}^{(i)}$ and use (1) and (25)\n",
    "\n",
    "$Y^2-3YX+3X^2= 6\\left(1+\\tan^2\\left(\\left(1-\\lambda\\right)\\psi\\right)\\right)$ (35)\n",
    "\n",
    "$X=\\frac{\\Delta\\alpha \\Delta T_{substrate}}{\\sqrt{\\zeta}}$\n",
    "\n",
    "$Y=\\frac{\\Delta\\alpha \\Delta T_{sur/sub}}{\\sqrt{\\zeta}}$\n",
    "\n",
    "$\\zeta=\\frac{\\left(1-\\nu_2\\right)\\Gamma_{IC}^{(i)}}{\\left(1+\\nu_2\\right)E_2 H}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: invalid value encountered in sqrt\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n",
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: invalid value encountered in sqrt\n",
      "  after removing the cwd from sys.path.\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "2eaf5cfcae8f46eba4e7c3ba2a9caa0f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f599cabb0f0>]"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x=np.linspace(0,10)\n",
    "LL=np.array([.2,0.3,0.5,1])\n",
    "Yp=3*x/2+np.sqrt(9*x**2-4*6*(1+tan((1-LL[-1])*psi)**2))/2\n",
    "Yn=3*x/2-np.sqrt(9*x**2-4*6*(1+tan((1-LL[-1])*psi)**2))/2\n",
    "\n",
    "Yp\n",
    "Yn\n",
    "\n",
    "plt.close(7)\n",
    "plt.figure(7)\n",
    "plt.plot(x,Yp)\n",
    "plt.plot(x,Yn)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in sqrt\n",
      "  import sys\n",
      "/opt/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:14: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.70710678,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [ 0.91287093,  0.91287093,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [ 0.57735027,  0.91287093,  0.57735027,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,  0.70710678,  0.70710678,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,         nan,  0.57735027,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,         nan,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,         nan,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,         nan,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,         nan,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,         nan,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan],\n",
       "       [        nan,         nan,         nan,         nan,         nan,\n",
       "                nan,         nan,         nan,         nan,         nan,\n",
       "                nan]])"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#X,Y=np.meshgrid(np.linspace(0,10,11),np.linspace(0,10,11))\n",
    "LL=np.zeros(X.shape)\n",
    "sx,sy=X.shape\n",
    "psi=-7.9*pi/180\n",
    "#fl= lambda l,x,y: y**2-3*y*x+3*x**2-6*(1+tan((1-l)*psi)**2)\n",
    "\n",
    "LL=-np.arctan(np.sqrt((Y**2-3*Y*X+3*X**2)/6-1))/psi+1\n",
    "\n",
    "#for i in range(0,sx):\n",
    "#    for j in range(0,sy):\n",
    "#        LL[i,j]=fsolve(lambda l: fl(l,X[i,j],Y[i,j]),10)\n",
    "        \n",
    "\n",
    "np.sqrt(1-(Y**2-3*Y*X+3*X**2)/6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "FigureCanvasNbAgg()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7f599ca8a7b8>"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.close(7)\n",
    "plt.figure(7)\n",
    "plt.contourf(X,Y,LL)#,levels=np.array([0.2,0.3,1,2,3,4]))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[         nan,          nan,   6.69620253,   8.83221447,\n",
       "          9.7715429 ,  10.3126658 ,  10.66670605,  10.91699711,\n",
       "         11.10353577,  11.24802139,  11.36327818],\n",
       "       [         nan,          nan,   3.81109548,   8.06444762,\n",
       "          9.38649635,  10.07875599,  10.50900449,  10.80331037,\n",
       "         11.01762974,  11.18079391,  11.30922102],\n",
       "       [         nan,          nan,          nan,   6.97542568,\n",
       "          8.90189269,   9.7996384 ,  10.32677864,  10.67479245,\n",
       "         10.92205971,  11.10691465,  11.25038846],\n",
       "       [  5.4638468 ,          nan,          nan,   5.4638468 ,\n",
       "          8.30233756,   9.46918111,  10.11705268,  10.52982866,\n",
       "         10.81587641,  11.02579012,  11.18639088],\n",
       "       [  7.61250077,   3.81109548,          nan,   3.81109548,\n",
       "          7.61250077,   9.08859003,   9.87878878,  10.36745248,\n",
       "         10.69838604,  10.93694042,  11.11689403],\n",
       "       [  8.67925642,   6.97542568,   3.81109548,   3.81109548,\n",
       "          6.97542568,   8.67925642,   9.61590732,  10.18826201,\n",
       "         10.56948461,  10.84013917,  11.04168414],\n",
       "       [  9.34242499,   8.30233756,   6.69620253,   5.4638468 ,\n",
       "          6.69620253,   8.30233756,   9.34242499,   9.99597175,\n",
       "         10.43020765,  10.73563594,  10.96076804],\n",
       "       [  9.7996384 ,   9.08859003,   8.06444762,   6.97542568,\n",
       "          6.97542568,   8.06444762,   9.08859003,   9.7996384 ,\n",
       "         10.28354941,  10.62448725,  10.87452013],\n",
       "       [ 10.13550078,   9.61590732,   8.90189269,   8.06444762,\n",
       "          7.61250077,   8.06444762,   8.90189269,   9.61590732,\n",
       "         10.13550078,  10.50900449,  10.78388475],\n",
       "       [ 10.39326974,   9.99597175,   9.46918111,   8.83221447,\n",
       "          8.30233756,   8.30233756,   8.83221447,   9.46918111,\n",
       "          9.99597175,  10.39326974,  10.6906299 ],\n",
       "       [ 10.59761794,  10.28354941,   9.87878878,   9.38649635,\n",
       "          8.90189269,   8.67925642,   8.90189269,   9.38649635,\n",
       "          9.87878878,  10.28354941,  10.59761794]])"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}