{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#Snippets and Programs from Chapter 7: Solving Calculus Problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P182: Simple example of finding the limit\n",
    ">>> from sympy import Limit, Symbol, S \n",
    ">>> x = Symbol('x')\n",
    ">>> Limit(1/x, x, S.Infinity).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-oo"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P182: Specifying the limit direction: negative\n",
    ">>> from sympy import Limit\n",
    ">>> Limit(1/x, x, 0, dir='-').doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "oo"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P182: Specifying the limit direction: positive\n",
    ">>> from sympy import Limit\n",
    ">>> Limit(1/x, x, 0, dir='+').doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P183: Indeterminate limit example\n",
    ">>> from sympy import Symbol, sin \n",
    ">>> Limit(sin(x)/x, x, 0).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "p*exp(r*t)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P183: Continuous Compound Interest\n",
    ">>> from sympy import Symbol, Limit, S\n",
    ">>> n = Symbol('n')\n",
    ">>> p = Symbol('p', positive=True)\n",
    ">>> r = Symbol('r', positive=True)\n",
    ">>> t = Symbol('t', positive=True)\n",
    ">>> Limit(p*(1+r/n)**(n*t), n, S.Infinity).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10*t1 + 2"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P184: Instantaneous Rate of Change\n",
    ">>> from sympy import Symbol, Limit\n",
    ">>> t = Symbol('t') \n",
    ">>> St = 5*t**2 + 2*t + 8\n",
    ">>> t1 = Symbol('t1')\n",
    ">>> delta_t = Symbol('delta_t')\n",
    ">>> St1 = St.subs({t: t1})\n",
    ">>> St1_delta = St.subs({t: t1 + delta_t})\n",
    ">>> Limit((St1_delta-St1)/delta_t, delta_t, 0).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10*t1 + 2"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P185: Finding the derivative\n",
    ">>> from sympy import Symbol, Derivative \n",
    ">>> t = Symbol('t')\n",
    ">>> St = 5*t**2 + 2*t + 8\n",
    ">>> d = Derivative(St, t)\n",
    ">>> d.doit().subs({t:t1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2*x + 1)*(x**3 + x**2 + x) + (x**2 + x)*(3*x**2 + 2*x + 1)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P186: Derivative of a complicated arbitrary function\n",
    ">>> from sympy import Derivative, Symbol \n",
    ">>> x = Symbol('x')\n",
    ">>> f = (x**3 + x**2 + x)*(x**2+x)\n",
    ">>> Derivative(f, x).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter a function: 2*x**2 + 3*x + 1\n",
      "Enter the variable to differentiate with respect to: x\n",
      "4⋅x + 3\n"
     ]
    }
   ],
   "source": [
    "#P186: Derivative calculator\n",
    "'''\n",
    "Derivative Calculator\n",
    "'''\n",
    "from sympy import Symbol, Derivative, sympify, pprint\n",
    "from sympy.core.sympify import SympifyError\n",
    "\n",
    "def derivative(f, var):\n",
    "    var = Symbol(var)\n",
    "    d = Derivative(f, var).doit()\n",
    "    pprint(d)\n",
    "\n",
    "if __name__=='__main__':\n",
    "    f = input('Enter a function: ')\n",
    "    var = input('Enter the variable to differentiate with respect to: ')\n",
    "    try:\n",
    "        f = sympify(f)\n",
    "    except SympifyError:\n",
    "        print('Invalid input')\n",
    "    else:\n",
    "        derivative(f, var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Theta: 44.997815081691805\n",
      "Maximum Range: 63.7755100185965\n"
     ]
    }
   ],
   "source": [
    "#P193: Use gradient ascent to find the angle at which a projectile has a maximum range\n",
    "\n",
    "'''\n",
    "Use gradient ascent to find the angle at which the projectile\n",
    "has maximum range for a fixed velocity, 25 m/s\n",
    "'''\n",
    "import math\n",
    "from sympy import Derivative, Symbol, sin\n",
    "\n",
    "def grad_ascent(x0, f1x, x):\n",
    "    epsilon =  1e-6\n",
    "    step_size = 1e-4\n",
    "    x_old = x0\n",
    "    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()\n",
    "    while abs(x_old - x_new) > epsilon:\n",
    "        x_old = x_new\n",
    "        x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()\n",
    "    return x_new\n",
    "\n",
    "def find_max_theta(R, theta):\n",
    "    # Calculate the first derivative\n",
    "    R1theta = Derivative(R, theta).doit()\n",
    "    theta0 = 1e-3\n",
    "    theta_max = grad_ascent(theta0, R1theta, theta) \n",
    "    return theta_max\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    g = 9.8\n",
    "    # Assume initial velocity\n",
    "    u = 25\n",
    "    # Expression for range\n",
    "    theta = Symbol('theta')\n",
    "    R = u**2*sin(2*theta)/g\n",
    "    theta_max = find_max_theta(R, theta)\n",
    "    print('Theta: {0}'.format(math.degrees(theta_max)))\n",
    "    print('Maximum Range: {0}'.format(R.subs({theta:theta_max})))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEACAYAAABVtcpZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucTfX+x/HXZwYZt1wLUWmSpERKSSdTJ9HN5XRDF5Uu\nSqhUpNSIjuhIFOkkRVGoODQlSpN+pTu5Ry5FhRBF7vP5/bHWZpvmsmZmr732nv15Ph77MXutmb33\n22B/9/p+vhdRVYwxxiSmpKADGGOMCY41AsYYk8CsETDGmARmjYAxxiQwawSMMSaBWSNgjDEJzNdG\nQETqicj8sNt2EekhIpVFZLaIrBCRWSJS0c8cxhhjcibRmicgIknAz0BToDuwWVWHiEhvoJKq9olK\nEGOMMQdFszvoIuAHVV0HtAHGuefHAe2imMMYY4wrmo1AB+B19/7RqrrRvb8RODqKOYwxxrii0giI\nSCngCmBK9u+p0x9la1cYY0wASkTpdS4BvlHV39zjjSJSXVU3iEgNYFP2B4iINQzGGFMIqipefzZa\n3UEdOdQVBDAd6Oze7wxMy+lBqhpTt8ceeyzwDJapeOWyTJYp0reC8r0REJGyOEXht8NOPwm0FJEV\nwIXusTHGmCjzvTtIVXcCVbOd24rTMBhjjAmQzRgugLS0tKAj/I1l8i4Wc1kmbyyTf6I2WaygRERj\nNZsxxsQqEUFjsDBsjDEmBuVbE3DX9WkGHI8znn8tME9Vt/uazBhjjO9y7Q4SkX8AD+C8+c8HfgEE\nqAE0xmkMhqjq//kSzLqDjDGmwAraHZTXlUB7oJeqrszlhU4CugK+NALGGGP8Z4VhY4wpRiJeGBaR\ne0TkSHG85O4L0KpoMY0xxsQCL6ODbnGLwBcDlYEbsBm+xhhTLHhpBEKXFZcBr6rqYh/zGGOMiSIv\njcA3IjILuBSYKSIVgCx/YxljjImGPAvDIiJAbaAasEpVt4lIFeAYVV3oazArDBtjTIEVtDDspRFY\npKqnRiJcQVgjYGLFqXLo/1NtoBywA1gHLLZ/oybGRHKeAKqqIvKNiDRV1S+LHs+Y2DI3I4MnunRh\nw8aNJAP73PM5vbU3A14MO74Np4GwhsDEs3znCYjI98CJwI/ATve0qmpDX4PZlYDxQehNf+vGjRwA\nygB7cD4N7XF/Zh+HGoHaYY99L4fnuwR4z/6dmhgS0SsBl80JMHFtVHo6UwYPZvvu3ZQBknHe/I9w\nb/txhsDl9FZeLp/nLgeUuP94UvbXpEJSDaqWOoZ6VetxTuoptGrcgAbHHxXJP4oxEedpxrC7jtCJ\nqvqyiFQDyqnqGl+D2ZWAKYK5GRkM79GDNatXUw3nU345Do13BigFlAQOkHsjsCPsfm5XAr0XrGLp\nul/4YeOvrP5tPd9vXc4ve5fyZ8oSJKsUR+9tRpOjzqNt4+Z0SmtCmdIlI/OHNCYHES0Mu0+YDjQB\n6qnqSSJyDDBZVZsXKWl+wawRMAU0NyOD8f36sWzRIsrs308JnDHQ5YC/cGY67gn7+VDfZl5XAuvC\n7mevCdwKfE7uxeGsLOWL5euY+MmnfLzmU1bu+YQ9pX/imD0tuezEK7j3ikuoV7tqjo81prD8aAS+\nw1k19BtVbeyeW2g1ARNLenfqxIrXX0eAvUBZ93yK+3W3+3Vn2GN24vSH/k7uNYHywJ9hjynq6KCF\nqzcw7J13eX/NO/ya8iGV/2rKlXVv4PGO/6J65fw6n4zJnx+NwJeq2lRE5qtqY3fj+HnWCJighX/y\nr7h/P42BlUBpDnX77HK//gUcDfzE4VcD1XDe5DdArqODSrmPDb35Vz//fF7++OMi59/6xy4GTJ7B\nxCWvsqn0J9TZewV9W3bj1tbnFPm5TeLyoxF4AGd00MXAIOAWYKKqjihK0HyDWSNgchHq789avZrK\nwCacrp46OJtc7ObQm3gLYAZOY1AVOBJYD2wNe76qOG/wO0U4rUMHBk+cGJU/R7glazfR943XeHfz\nSI44UIUbTurB4BuvpkLZI6KexcS3iDcC7pNejNMIALyvqrMLmc8zawRMTnp36sTK118nCwjNYFyP\n88Z/Is6VQHmcN/ldwHFAQ+BNnE/xZdzvlwRKVK/O3WPGcP5ll0X3D5GHvfsOMOCNdxn1zQi2lVxG\nu2oP8mLX26hcISX/BxuDP1cCg1W1d37nIs0aARMuNL6/zMaNNMDpj6/jfm8Zzhu7cqh/vzywDaeb\n50jc+QClSnHJQw9xV3p6tOMXyoQ53/LAjAFsLPkFV1R+gDF33kHVI8sEHcvEOD8agfmhgnDYuUWq\neprHQBWBMUADnP+nN+N8YJuE80FtLXCNqm7L9jhrBMzBrh9dvZp9OMPUAJbjfPIHqAlk4rzxC7Da\nPV8GOEKEciecQJfhw2PqE39BTJ77HfdOe5yNyV9y8/EDeb7rDZRI9rL2o0lEEWsERORO4C4gFVgV\n9q3ywKeqep3HQOOAj1V1rIiUwBm48TCwWVWHiEhvoJKq9sn2OGsEElyo66ckUBen26eW+70NHOr3\nrw4cA2TgXAWUAJJjsKunqP773jx6zbqPA+xh4PlDua/9BUFHMjEoko3AkUAlnA1kenNowMWfqrrF\nY5gjgfmqekK288uBFqq6UUSqA5mqenK2n7FGIIGFhnw24NC09mVAN2Bc2M8J8CuH3vzLpabG9af+\n/GRlKb3GTuG55b2pkXU207sNo1FqjaBjmRjiV2G4UDOGRaQR8AKwFDgd+Aa4B1ivqpXcnxFga+g4\n7LHWCCSgUPfP7tWrOcs9t9/9WhNYCHQCxuO8+QuQUrIkR512Gtc+/nixffPPbvP2v2jzn4F8vu9F\nrqqazms9u1KqZHLQsUwMiNaM4Smqeq6HMGcC84BzVfUrEXkGZ1j23eFv+iKyVVUrZ3usPvbYYweP\n09LSSEtL8/rnMnFoVHo6mYMGIXv3ksKhrp+LcT79h7p95uLMBdgsQoOAhnTGiv99toQbJ3flgOxj\ncsdXuLTpyfk/yBQrmZmZZGZmHjzu379/7MwYdrt65qlqHff4POAh4ATgAlXdICI1gI+sOyixjUpP\nZ/aAATTIyqIEh3f9VMdZxXA8zizf/UCZYt7tUxD7D2RxwzMvMOm3R2lX+REm9+puheMEFnMzhkVk\nLnCrqq5wrypCY9y2qOpgEekDVLTCcOKZm5HBrBEj+H75cvb99BON3PP7+XvXz06c2bxJxbDgGykf\nzv+B9uM6k0wp3rntFZo3OC7oSCYAMTdjWEROxxkiWgpnlNHNODP0JwPHYkNEE9Ko9HQWDhlCp127\nGAnU51Df/8XA+xze9fMbcGrHjgnd9ePF3n0HaDdkKDP/+A8P1H+BwTe1DzqSiTKbMWxiXqjrZ2pW\nFo/gjOrZz+F9/62A2TifFhaXKsWFcTTJKxa89P4XdJ3dgVNKXMEn/Z6y5ScSiC+NgPvER+L8f1UA\nVd2a9yOKxhqB4inUAJyelUU6kM6hBuB9Du/73weUtb7/Qlvz6+80H9KF7fzEezdP4fyGdfJ/kIl7\nBW0E8q0eicgdIrIBp4v2a5xhnl8XPqJJVOENQKjrJ3sDMBtnVNDupCQufOwxxv3wgzUAhVSnRiXW\nD32LVjVu4IIJzXjqrQ+DjmRikJeawA/AOaq6OTqRDr6uXQkUI9kbgPA3/vAGIBlYlpJCiwcftO6f\nCHpmWia9PutImyq9eeuBniQlef6gaOKMH4XhWUB7Vd2Z5w9GmDUCxUNozf+tCxbQUDXHT/6bcJaE\nqHnccVQ/+WRadu9un/598OmSH2k5ph01kk9jfv8XrU5QTPnRCJwBvIIz6Wuve1pVtUdhQ3oKZo1A\n3AuNAKq6a9dhxd/sn/y/S0qiZb9+9sk/CjZv/4vG/W/ij6wNfH3/VOrWqhJ0JBNhEa8JAP8FPsDZ\nTjVUE/imcPFMogh1/4y2BiCmVD2yDGueeoOTyzWjwbBzmbNgVf4PMsVaifx/hGRVvc/3JKbYmJuR\nwcdDhnB6VhZgDUCsKZGcxBcDB3PdsBNoOfE8nv/1bW6/pFnQsUxAvFwJvOeOEKohIpVDN9+Tmbg0\nNyODkZ07U3/Xrr9N/spe/LUGIFgT7r2DxxqPpWtmW/pPfDfoOCYgXmoCazl8321wagIn5PDjEWM1\ngfgTqgFUdxuA0Jv/Ezgzf2cDP4hQuXHjhFrxM9aNmfk5d3zUjq51nmZk105BxzFF5NtksWizRiC+\nZJ8FbN0/8eV/ny3hyqmtaX9Ub6Y8cHfQcUwRRKwwLCJpHl7MtjYyf6sBWPdP/Gl7bgPm3DCXaRue\noeWAJ4KOY6Ior53F/gOcjzMy6GucPTyScJZ2ORO4CGcJ6Ad9CWZXAnHhYA1gyxb2AwND5znUACyv\nUoW7xo2z7p84sGDVr5w98p80q3ANcx59zCaVxaGIdgeJSHmgLdAcZ1N4gB+B/wP+p6o7ipA172DW\nCMS8vGoAIXekpHDdlCnWAMSRxWs2cuazF3FmubbMTR9gDUGcsZqAiYq5GRmMvPpqJu3aZTWAYmjZ\nT79xxrCLaFi2NfMef9Iagjjix2QxYw4TPgwUrAZQHNU/thoLe81h0c7ZnPfYw2Rl2Qey4squBEyB\nZO8CshpA8fb9us2cPiyN8yp24INHHwk6jvHArgSMb0KjgEaH1QAedr93PjAA+DUlxRqAYqRe7ap8\nefcHzN0+nisGDQ06jvGBl/0EyopIPxF50T2uKyKX+x/NxJK8uoD64WwOc21KCqc/+KA1AMVMwxOq\n88mtHzJzy0g6DB0VdBwTYV6uBF7GWT30XPf4Fw4fAGKKubkZGbzfs+fBYaDgfPIPrwEsr1KFblOm\nWA2gmDq7fm0+vPFD3twwiG6jbZ/n4sRLI5CqqoNxl5GO9r4CJnhv9OvHE6tWWRdQgju/YR2mtHuP\n59fcy8A3ZgYdx0SIl1VE94hISuhARFKBPf5FMrFkbkYGO5YtA/7eBRS+C5g1AImhffNTGbntbbp9\n0o4ald6hS6uzg45kishLI5AOzARqichEnIljN/mYycSIUB2g7u7dgPPJH7J1AdkVQMK587LmbNz2\nCrfPacvRFT/i8rPrBx3JFIGnIaIiUhU4xz38PBr7DdsQ0WCF6gAlV63iQv4+E7hr6dJ0evNNawAS\n2O0jxzN2zaN8fec8GqXWCDqOcUVsxrCINOHwJaRDT6oAqvqtx0BrgT+AA8A+VW3q7kcwCWcpirXA\nNaq6LdvjrBEISOgKYNKWLTyCMxcgfB7AAeDXxo0Z862nfwKmGLvo8YF89vtUVvf7mOqVywUdxxDZ\nRiCTv+8jcJCqelpBVETWAE1UdWvYuSHAZlUdIiK9gUqq2ifb46wRCED4FUA6zpt/9quAvqmptB4+\n3K4CDFlZSv3et7Ft/wZ+HDyN0qW89DAbP8Xc2kFuI3Cmqm4JO7ccaKGqG0WkOpCpqidne5w1AgG4\n64wzGDV//sErALDZwCZvf+3ex7F9Lqf6EaksHDTS1hkKWCT3E7jQ/XqliPwr+60AmRT4QES+FpHb\n3HNHq+pG9/5G4OgCPJ/xSfaRQNmHgu5NTbUGwPxNmdIlWfDwFFbt/ZT2Q4YFHccUUF7Xbi2AOcAV\n5Nwt9LbH12iuqr+KSDVgtnsVcJCqqojk+JE/PWziUVpaGmlpaR5f0hRUbiOBQkNBl1epwl3WBWRy\nUataBebcNoPzXjmH9An1SL/O/p1ES2ZmJpmZmYV+vJc9hk9Q1dX5nfP0YiKPATuA24A0Vd0gIjVw\nNqex7qCA2EggEyn/fW8eXTPbMrXtR7Q9t0HQcRKSHwvIvZnDuSkew5RxN6ZBRMri9DIsAqYDnd0f\n6wxM8/J8JvJCVwChGcGh5SBC6wH1A/bXr28NgPHk9kuaccfxT3PV21fw/TrfR5KbCMi1O0hE6gOn\nABXdGoDgdAtVAEp7fP6jgakiEnqtCao6S0S+BiaLSBfcIaKF/hOYQgtfEwgO1QGe4FB3UN/UVG4c\nMCCghCYePX/n9Szot4SmQ//Fz4M+oFxKqaAjmTzkNUS0LdAepyYwPexbfwJvqOpnvgaz7iBf5TQX\nAGwkkImM/QeyqH1/e6qWqsWiwSODjpNQIj5EVESaqeq8IicrIGsE/GNzAUw0/LRpO3UHn02n4x7k\n5R63BB0nYfhRE7hTRCqGvUAlERlbqHQmJoSvCgp/rwN0qFLFGgBTZMcedSRTr53GuJ/78PKsL4OO\nY3LhpRFoGL6kg6r+DpzhXyTjJ5sLYKLp0qYn06fBi9w2+0oWrt4QdByTAy+NgLhr/YQOKuN0GZs4\nNGvECI4NmwtgVwDGb/++sS3Ny97CP0Zcw1+79wUdx2TjpREYCswTkQEiMhCYBzzlbyzjh7kZGaz7\n8sscrwA2lC5tVwDGNx/2e4wjpBwtBvQNOorJxutS0g2A0IJxc1R1qa+psMJwpIWKwbJqla0KagKx\ncv0W6j/ThF6nDmPwTe2DjlNs+VEYBqgM7FTV54DfRKROodKZQIRPCAtdBYSuANKBAzYXwERB3VpV\nePHiyTy17A5mf7My6DjG5WWIaDrQBKinqieJyDHAZFVt7mswuxKIiOzDQeHwq4DvK1XizldftW4g\nEzUdho7if+tf4Of+n1O5Qkr+DzAF4seVQHugLbATQFV/BsoXLp6JtlkjRhw2HBQOvwo4sWlTawBM\nVE28906OTmpAswHdg45i8NYI7FHVrNCBuwaQiQOhQjAcPhw0pG9qKi27239EE11JScJnfV9gTdYn\ndBs9Meg4Cc/LNkBTROQFnDWEbgduAcb4G8sUVagbqPY2Z4qHLQ1tYknNKuV5tc0kOr7XknbfnEXL\nJnWDjpSwvI4OuhjnwyTA+6o629dUWE2gKMLXBbIlIUws6zB0FNPXj2HDwHlUKHtE0HGKhZjbXrKw\nrBEoHCsEm3iSlaUce//VVClVk++eHBF0nGIhkttLfup+3SEif2a7/SEia0WkWyRCm8ixQrCJJ0lJ\nwqe9x7B0/zv0Hf+/oOMkpFwbgdAQUFUtp6rls90q4Awb7RGtoCZ/Vgg28ei4oysy8p8TGLzkDr5e\n8XPQcRKO15rA6TgfKBX4RFW/c8/XVNVffAlm3UEFkn1GMNjeACa+/PPxAczf+jGbhs6iRLLXeawm\nu4jPExCRnsAEoBrOTmGviUgPAL8aAFNwoW4gWxnUxKv3HurLAfbQdvDQoKMkFC8zhhcB56jqTve4\nLPC5qp7mazC7EvBsbkYGL11/PePc4aBWCDbx6tMlP/KP8Wcx/uL3uP6fTYKOE5f8WjsoK5f7JmDZ\n5wOAFYJN/Gre4Di6n/gst7zbkU2/7ww6TkLw0gi8DHwhIuki0h/4HLCdxWJETt1AIVYINvFo+G3X\nUlvO4cInHwg6SkLIsztIRJKAZsBu4DwOFYbn+x7MuoPyZd1Aprj6adN2TniqIY+eMZpHO14SdJy4\n4sdG8wtUtVGRkxWQNQJ5y2k0ULh+rVoxYObMqOcyJlKenvoRD8y7nuU9FlK3VpWg48QNP2oCH4jI\nVSLi+UmN/6wbyBR397W/gMYlOnDB03eQlWUfCP3i5UpgB1AGZwOq3e5pdSeM5f8CIsnA18B6Vb3C\n3aN4EnAcsBa4Jnwj+7DH2ZVALqwbyCSKbTt2U/3RM7npxN6MvuuGoOPEhYhfCbgzhpNUtWS2GcNe\n9QSW4tQTAPoAs1X1JOBD99h4ZKOBTCKpWK4049q+xn9/7MVX368POk6x5GWymIjIlSIyTESGiojn\nzUFFpBZwKc7S06GWqQ0wzr0/DmhXwMwJzbqBTKK5tkUj0srczaWjb7NuIR94qQmMAu4AFgJLgK4i\nMsrj8w8DHuDwuQVHq+pG9/5GnFnIxoPwtYHOB1rh7A+QDnSsVMmWhzbF1ju9H2Inm7j52ZeCjlLs\neNlU5gLglNDuYiLyCk73Tp5E5HJgk6rOF5G0nH5GVVVEcm3a09PTD95PS0sjLS3Hp0kIuXUDHdws\nxrqBTDFWpnRJJlw1jitnXMDtS1rSvMFxQUeKGZmZmWRmZhb68V4Kw+8Ad6vqWvf4eOA5Vb08n8f9\nG7gB2A+UBioAbwNnAWmqukFEagAfqerJOTzeCsNhHmnVioGzZtkmMSahtR74JF9u/sAWmcuDH0NE\nKwDLRORjEcnEuQooLyIzRGR6bg9S1b6qWltV6wAdgDmqegMwHejs/lhnYJrXsImsxJ49gHUDmcQ2\n7cH72Ss7uOGZF4KOUmx46Q56NIdzilPoLchH9dDPPglMFpEuuENEC/AcCWduRgazRozgh+++O3jO\nuoFMoipdqgRvdHiZNlPP5+4ll1q3UATY9pIxLFQHeGLVKusGMiZMq4GD+GZzJpuenklSks1jDWd7\nDBcjoTpASGhS2E+VKnFs06a07N7dGgCTkP7avY+qD51DhxPuZmz3m4OOE1P8WkraRFn4cNCQ0KSw\nOg0bMmDmTGsATMIqU7okr7Qfyyvre/PtStvbqii8TBZr464maqIkp+Gg4Q6ULh3lRMbEnmvOP51/\npNzJpaO62iSyIvDy5n4t8IOIDBGRvw3lNJFns4KN8Saj98NskzXcM2ZS0FHilpe1g64DGgOrgVdE\nZJ6I3C4i5X1Pl6BsOKgx3pRLKcVzF7/IyB/uZeX6LUHHiUueunlUdTvwJs7qnzWB9sD80IbzJnLm\nZmSwbPHig8e2OJwxebu19TmclnwNlz5zf9BR4pKXmkBbEZkKZAIlgbNU9RKgIXCfv/ESS6gW0G3L\nFusGMqYA3r1/IGtkDkPe/CDoKHHHy2SxfwHDVHVu+ElV/UtEbvUnVmIK1QJC+uHsEbC8ShXusm4g\nY3JVs0p5Hmk0ioc/v4NbWi6i6pFlgo4UNzzNE3DX+GmKsxroV6q6wfdgCTZPIPtGMeHSW7QgvQgL\nRBmTKI7r1ZGjStfiqyeeCjpKYCI+T8D9tP8FzhXB1cAX7pIPJkJsSKgxkfFu9+F8s38ckz5eEHSU\nuOGlMPwg0FhVO6vqjcAZQG9/YyUWGxJqTGQ0OP4objxmEF2m3cHefQeCjhMXvDQCm4EdYcc73HMm\nQmxIqDGRM6bbzSTrEVw/fHTQUeJCroVhEenl3v0BpwsotORzW5xdxkwR2QqhxkReieQkXu0wmnbT\nWvDtyvacUbdm0JFiWl5XAuWBcsAqnDX/1b39D2fimCmCUB1g4KxZ3LVtm3UDGRNBbc45heZHdKXN\n8z2DjhLzbBXRgNgKocb4a+sfu6jevyF9Gg/j8evz3AixWLFVRONEqA4QYiuEGhNZlSukMLDZ8wxa\n0J3N2/8KOk7MskYgANmXhghnw0GNiZwHr7qImnoObf4zMOgoMcsagSizpSGMia6pXYfy+b7/8s4X\ny4KOEpNyrQmIyLNhh6E9hQ8eq6qvi8cV15pAeC0gVAc4uDTEuHHWDWSMD64cMoI5v0xly9Nziv12\nlJGsCXzj3o7AmSC2AlgJNAJKFSVkIguvBYSvEHryqadaA2CMTybccxd7ZDt3/3di0FFiTq6NgKq+\noqqvAKcDF6jqs6o6ArgQZ38BU0BWCzAmGKVLleDZ1s/zwuoH+HFjzsuzJCovNYGKQIWw4/LuOVMA\nVgswJlhdWp3NSVxBm2GPBh0lpuQ7T0BEbsbpsch0T7UA0t2rBP+CFbOagNUCjAneyvVbqPdsfV6/\ndBbXtmgUdBxfRHyegKq+DJwDvO3ezvHSAIhIaRH5QkQWiMhSERnknq8sIrNFZIWIzBKRhLiqsFqA\nMcGrW6sKnWoM5Lapd9vm9C4vS0knARcBp6vq/4BSItI0v8ep6m6cWkIjnF3ILhCR84A+wGxVPQn4\n0D0utuZmZPBIq1YsmZ/z0rZWCzAmusZ268IB2UPX518NOkpM8FITGAU0Azq6xzvcc/lS1dA0vVI4\nvR+/A22Ace75cUA7r2HjTfj6QN3/2G61AGNiQKmSyTx3yUhe+rG3FYnxVhOYr6qNQ1/dc9+p6un5\nPrlzFfEtkAo8r6oPisjvqlrJ/b4AW0PH2R4b9zUBWx/ImNhV/8HbOSI5hQWDhgcdJaIKWhPwssfw\nXhFJDnuBajjbTOZLVbOARiJyJPC+iFyQ7fsqIrm+06enpx+8n5aWRlpampeXjRk5rQ90PpDesCHp\nM2cGkskY45jW/d/Uf+4U3vykC1f9o2HQcQotMzOTzCJsP+vlSuB64BqgCU73zVXAI6o6uUAvJNIP\n2AXcCqSp6gZ37+KPVPXkHH4+rq8E5mZkMLJzZyZt2fK37/Vr1YoB1ggYE7gOQ0cxc91ktj79UbGZ\nSezH6KDXcLaTHAT8ArT10gCISNXQyB8RSQFaAvOB6UBn98c64+xVUKzYnABj4sMr3W9nj/xOr7FT\ngo4SmLzWDqqc/ZT7VQFUdWueTyxyGs6VQ5J7e1VVn3KfdzJwLLAWuEZV/1adiecrAZsTYEz8eHb6\nXO795Hp+6buMoyqVDTpOkRX0SiCvRmAt7ht+TlS1ToHTFUA8NwLpaWmkf/zx38+3aEF6EfrujDH+\nOL5XJ2qXS+WT/gOCjlJkEesOUtXjVbVObrfIxC1eQnMCloftGRzO5gQYE5um3D6ET/c8T+Z3ibdz\nrqfJYiJyg4g86h4f62WyWKKxPYONiV9n1atFy3K96PjKvUFHiTovo4NG4wwJvVBVT3b79Gep6pm+\nBouz7iCbE2BMfNu2YzfV0hsw4OzR9Lm6ZdBxCs2PPYbPVtW7cIZ3hgrCJQuZr9iyPYONiW8Vy5Xm\n/tOGkj7vHnbv3R90nKjx0ggUerJYorB9AowpHp64oS1ls2pyw/DRQUeJGi+NwLPAVOAoEfk38CnO\nnAGDzQkwpjhJShLGXD2MtzY/zsr1f5/oWRzlWxMAEJH6wD/dww9V1fcdm+OlJmBzAowpfhr26Q4o\nC598LugoBRbxmoC74XwlVX3OvfneAMQT2yfAmOJnao/+LNbJvPV/i4KO4jsv3UHfAI+IyGoR+Y+I\n+DoqKN7sP+KIHM9bLcCY+JVaszJXVn2U2968t9hvPuNl7aBXVPVS4Czge2CIiPzge7IYF5oY9vNP\n67gx21qsVgswJv6N634HO5N+4bEJ7wQdxVdelpIOORE4GTgOWOpPnPgQKgY/sWqVcwxcm5JCjdRU\nyh9zDK2/2wL6AAATfElEQVRtToAxca9M6ZL0a/o0j3/Vg95XtaJcSqmgI/nCy2SxIUB7YDXwBjA1\npwXfIh4shgvD2SeGhdgS0cYUP0fdeynNq1/M1N73BB3FEz82lVkNNFPVzYWPVbxknxgWkrx7d5ST\nGGP8NrbDUNpMPZ/v111PvdpVg44TcV5qAqOBAyLSVETOD92ikC0m2cQwYxLL5WfX51Q6cPXI9KCj\n+MLLENHbcLq9ZwH9gfdxRkEmHJsYZkxieuvudBbrZKZ/XvzKoV5qAotxRgbNU9VGInIyMEhV2/sa\nLAZrAjYxzJjE1e7JYXy28QM2DcsIOkqe/FhAbreq7nKfvLSqLgfqFTZgPLOJYcYkrtd6dmNb0goG\nTf77oJB45qURWC8ilXD2Ap4tItNxtoVMODYxzJjEVS6lFPeeNoT+n/di774DQceJGC+F4Xaq+ruq\npgP9gDFAO7+DxZLQxLAfV6+1iWHGJLBBN7bjiKxK3DpybNBRIibPmoCIlAAWq+rJ0Yt08LVjoiaQ\n08SwkWETw2yzGGMSy/gPvubmWW1Y1/t7alYpH3Scv4loTUBV9wPfi8hxRU4Wp2aNGHGwAQCnFjBp\n1y7KH3OMbRZjTAK68aIzOT7rIq4ZPjjoKBHhpSZQGVgiInNEZIZ7m+53sFhhE8OMMdlN7PIEn+19\nni+WrQs6SpF5mTH8CJD90iL4fpoosIlhxpicnF2/Ns1KdaXTS4+w6j/jgo5TJF6uBC5T1czwG3Cp\nlycXkdoi8pGILBGRxSLSwz1fWURmi8gKEZklIhWL8GfwhU0MM8bkZVL33qxNfp/XM+cHHaVIvEwW\nm6+qjbOdW6Sqp+X75CLVgeqqukBEyuHsTdAOuBnYrKpDRKQ3zqY1fbI9NtDCsE0MM8bkp+PQ53l/\n3ZtsfvoDkpI812J9FbHCsIjcKSKLgHoisijsthZY6OXJVXWDqi5w7+8AlgHHAG2A0DXUOGJwyKlN\nDDPG5OelbreyM+kXBrzxXtBRCi2v7qCJwBXAdOBy9/4VQBNVva6gLyQixwONgS+Ao1V1o/utjcDR\nBX0+v9nEMGNMfsqULsmDjYcw6OsH2L13f9BxCiXXwrCqbge2Ax2K+iJuV9BbQE9V/VPk0JWKqqqI\n5Njvk56efvB+WloaaWlpRY2Sr7kZGcwaMeLgxLDxYX+vfVNTaW21AGNMmP7XXc5z3zzNraPG8to9\nt0f99TMzM8nMzCz04/OtCRSViJQE3gHeU9Vn3HPLgTRV3SAiNYCPsk9IC6ImYBPDjDGFEZpA9nOf\nFVSvXC7QLH4sIFeUMAK8BCwNNQCu6UBn935nnHWJAmcTw4wxhXHjRWdS+0AaHZ8dGnSUAvPUCIjI\n8SJykXu/jIhU8Pj8zYHrgQtEZL57aw08CbQUkRXAhe5x4GximDGmsMZ3foKPd41g4eoNQUcpEC+b\nytwOTAFecE/VAqZ6eXJV/T9VTVLVRqra2L3NVNWtqnqRqp6kqhdHY89iL6wYbIwprPMb1uGMpJvo\n+EL/oKMUiJcrgW7AecAfAKq6AjjKz1BBmJuRwS+bNnFT8uFdaTYxzBjj1aS7H2aZvMm7Xy4POopn\nXpaN2KOqe0IjetyVRYvVshGhgvDYVauYi7Ne9o+lS1P+lFO49vHHrRZgjPEktWZlLj2yN11e78Ov\nTWOi1JkvL1cCH4vIw0AZEWmJ0zU0w99Y0RVeEA5NDBu/ezeVq1WzBsAYUyCvdb+b35IX8NyMT4KO\n4omXRqAP8BuwCLgDeBdnUbliwwrCxphIqViuNLeeMICH5jxIVlbsd5p4aQRKAy+p6lWqehUwFkjx\nN1Z0/ZHLr8EKwsaYwnjujuvYzy76jPM0hiZQXhqBORz+pl8G+MCfONEV2jZywZKFdC6RfNj3rCBs\njCmsEslJPHLOYJ5Z/BB/7d4XdJw8eVlFdIGqNsrvXMSD+Txj2GYHG2P8lJWlVL2vJa2PvYqJ93WN\n2uv6MWN4p4g0CXuBM4FdhQkXS2x2sDHGT0lJwogrBjNpw+Ns2Loj6Di58tII3ANMFpH/E5H/AyYB\ncd9PYsVgY4zfrv9nE2odaEGnZ58OOkqu8m0EVPUroD5wJ9AVOFlVv/Y7mF9CdYDlC3PeEsGKwcaY\nSBp34xNk7hrOkrWbgo6SI0+riIrIuUAdnMllCqCq430N5kNNILwOMBd4H3gi7Pt9U1NpPXy4dQUZ\nYyKq0UM9UVW+e3KE769V0JqAl8Lwa8AJwALgQOi8qvraJeRHIxC+ZSQc2jbyp0qVOLZpUysGG2N8\nsWTtJk4bXZ85Hb8i7fQTfH2tgjYCXpaNaAKcEuiGvxGSvQ5wvntLb9iQ9JkzA8lkjCn+Ghx/FGkp\nPblpfD/WDp0QdJzDeCkMLwZq+B0kGmyVUGNMUCZ2v491yXN4PXN+0FEO46URqAYsFZFZIjLDvU33\nO1ikzc3IYMNvv9GlxOEXPzYpzBgTDdUrl+PKox6h+7Q+QUc5jJeaQFpO51U104c84a8bsR6o7AXh\n2dgqocaY6Nuxay+VHjmFf5/7Ag9c+U9fXiPiheGgRLIRyF4QDunXqhUDrBZgjImini9OYszSp/hz\n6FckJXl+r/Ys4jOGRaSZiHwlIjtEZJ+IZInIH0WLGV02McwYEyuG3nI1inL/y28GHQXwVhN4DugE\nrMRZUbQLMMrPUJFmBWFjTKwokZxEv2aDeG7pwzGxuJynjeZVdSWQrKoHVPVloLW/sSIjNDv4x9Vr\nuTHbYFgrCBtjgtL7qpaUO1Cb20e/HHQUT/MEdorIEcB3IjIE2ABEviMrwnJaJfTasFVCW9vEMGNM\nQJKShGGXPckts9vxzPbrqXpkmcCyeBkddBywCSgF3AtUAEap6g++BitiYdiKwcaYWFfrvqs4rcpZ\nvPdw74g9px9LSbdT1V2qul1V01X1PiDmP0JbMdgYE+te7PgE7//5H9b8+ntgGbw0AjflcO5mL08u\nImNFZKOILAo7V1lEZovICncCWkWPWQvEisHGmFh3yVn1qJfVng4jnwwsQ66NgIh0FJEZQJ2wmcIz\nRCQT2OLx+XMqIvcBZqvqScCH7nHEhIrBm37+mRuTD78ismKwMSbWjO/yKF/tH8PXK34O5PVzrQm4\ntYA6wJNAbw4Vg/8AFqrqfk8vIHI8MENVT3OPlwMtVHWjiFQHMlX15BweV+CagG0ZaYyJR00ffpAd\n+/5g6ZDRRX4uP5aSLgfsUtUDIlIPqAe8p6qeBrjm0Aj8rqqV3PsCbA0dZ3tcgRsBKwYbY+LRql+2\nUnfESbx/9TxaNqlbpOfyYynpj4F/iEglnH1YvgKuBa4rXMRDVFVFJNd3+vT09IP309LSSEtLy/P5\nrBhsjIlHqTUrc1HZ++gyoR8/NXmjQI/NzMwkMzOz0K/tpRFIUtW/RKQLztDQISLyXaFfETaKSHVV\n3SAiNXCGn+YovBHwworBxph49drdPakxqC4T5nzLdRee4flx2T8g9+/fv0Cv62nGsIg0w/nkn1GQ\nx+ViOtDZvd8ZmFaE5zrMxT16cGuVaoeds2KwMSYeHFWpLFcf/Qj3zHg4qq/rpSbQAugFfKqqg0Uk\nFeipqj3yfXKR14EWQFVgI/Ao8D9gMnAssBa4RlW35fDYAtUE5mZk8M7Tz/DJdx9x1N5y1DmuthWD\njTFxZceuvVTqV5+nznuJe9qlFeo5EnIp6eyjggAeTk2llW0ab4yJM3eNnsCr3z/H9qGfFWqp6YjN\nGBaR4e7XGTncYmpnsVkjRhzWAAA8sWoVs599NqBExhhTOCNu68h+2Um/CTOi8np5FYbHu1+H5vC9\nmLp8sFFBxpjiokRyEvc3foKn5vflsQ6XUapksq+vl+uVgKp+437NBJYAS1Q107197GuqAvo9K+fz\nNirIGBOP+l93OSWzytNzzOu+v1Ze3UEiIukishlYAawQkc0i8pjvqTwKLRGxaNkibipxeGtpo4KM\nMfEqKUkYeMEgxqx6lB279vr6WnktG3EfcAlwu6qucc+dAIwGZqrq074Gy6cwbEtEGGOKu6r3tuai\nWm14o9ddnh8TsdFBIrIAaKmqv2U7Xw1nAbhGnlMVQn6NgC0RYYwp7ibM+ZYbZ17Orw+t5KhKZT09\nJpL7CZTI3gAAuOe8zDT2lRWDjTHF3XUXnkHN/edx40j/Rjrm1QjktUBc4Lsj2xIRxphE8N+OA5i1\nY6hvG8/k1Qg0FJE/c7oBp/mSxoNQMXjF9z/Y5vHGmGLvkrPqUfdAW64b9ZQvzx9XM4atGGyMSUTz\nlv5E8/GNWXD7EhqeUD3Pny3Wy0ZYMdgYk6jOeOhe9ut+Fj6Zd33Aj43mY4YVg40xiWrCnX1ZzETm\nLlwT0eeNq0bAisHGmERV/9hqnHdEN24enx7R542LRsA2jzfGGHitWy/WJL/H9M+XRuw5PdUERORF\nVb0tt2M/hGoCVgw2xphDLvv3UyzY/Dk/P/1Wjt/3pTAsImeq6tdhx01CC8z5JdQIWDHYGGMO2frH\nLqoNPJGxLafRueVZf/u+L4XhbA2AAHW9vkBRWTHYGGMOqVwhhWur9+O+jMhsQ5nXKqLlRKSXiIwS\nkbtEJElE2uMsK90pIq/ugRWDjTHmcP+98xb+KLGKp6d+VOTnyutKYDzOzODvgH8CnwP3Ap1UtU2R\nX9mDuRkZbPjtN24pcfjUYCsGG2MSWbmUUtx24uM8NrcvWVlFm+uV1yqiC1W1oXs/GfgVOE5VdxXp\nFb0GE9G+qak8sWoVc4HZwI+lS1P+lFO49vHHrRhsjEloe/cdoELvRvRp8iTp1x16P4zkUtLzVbVx\nbsd+E5Eck1lB2BhjHLO/WUmTurWoXCHl4LmCNgJ5LQnd0F0sLiQl7FhVtULB4kaGFYSNMcbRsknR\nx+jktcdwsqqWD7uVCLtf5AZARFqLyHIRWSkivb0+zgrCxhgTOYHMGHZrDM8BrYFTgI4iUj/7z11b\n8fB4QReEMzMzA3vt3Fgm72Ixl2XyxjL5J6hlI5oCP6jqWlXdB7wBtM3+Q6su/hf9WrUivUUL+rVq\nRevhwwMtCMfiX7pl8i4Wc1kmbyyTf4LaJvIYYF3Y8Xrg7Ow/NGnYC6TWrBy1UMYYk2iCuhLwNLDV\nGgBjjPFXIJvKiMg5QLqqtnaPHwKyVHVw2M/E5m43xhgT42J+ZzERKQF8jzMT+RfgS6Cjqi6Lehhj\njElggdQEVHW/iNwNvA8kAy9ZA2CMMdEXs3sMG2OM8V/M7SxW2ElkPuQYKyIbRWRR2LnKIjJbRFaI\nyCwRqRjlTLVF5CMRWSIii0WkR9C5RKS0iHwhIgtEZKmIDAo6U1i2ZBGZLyIzYiGTiKwVkYVupi9j\nJFNFEXlTRJa5f39nx0Cmeu7vKHTbLiI9YiDXQ+7/vUUiMlFEjoiBTD3dPItFpKd7rkCZYqoR8DqJ\nLEpednOE6wPMVtWTgA/d42jaB9yrqg2Ac4Bu7u8nsFyquhu4QFUbAQ2BC0TkvCAzhekJLOXQaLSg\nMymQpqqNVbVpjGQaDryrqvVx/v6WB51JVb93f0eNgSbAX8DUIHOJyPHAbcAZqnoaTjd2h4AznQrc\nCpwFnA5cLiKpBc6kqjFzA5oBM8OO+wB9AsxzPLAo7Hg5cLR7vzqwPODf1zTgoljJBZQBvgIaBJ0J\nqAV8AFwAzIiFvz9gDVAl27nAMgFHAqtzOB8T/57c178Y+CToXEBlnMEslXBqqTOAlgFnugoYE3b8\nCPBgQTPF1JUAOU8iOyagLDk5WlU3uvc3AkcHFcT9ZNIY+IKAc7kbDi1wX/sjVV0SdCZgGPAAkBV2\nLuhMCnwgIl+LSGiP7iAz1QF+E5GXReRbEXlRRMoGnCm7DsDr7v3AcqnqVmAo8BPOiMZtqjo7yEzA\nYuAfbvdPGeBSnA8/BcoUa41A3FSp1WlmA8krIuWAt4Ceqhq+0msguVQ1S53uoFrA+SJyQZCZRORy\nYJOqzgdyHC8d0N9fc3W6OC7B6cr7R8CZSgBnAKNU9QxgJ9m6DgL+d14KuAKYkv17AfybSgXuwekd\nqAmUE5Hrg8ykqsuBwcAs4D1gAXCgoJlirRH4Gagddlwb52ogVmwUkeoAIlID2BTtACJSEqcBeFVV\np8VKLgBV3Q5k4PTjBpnpXKCNiKzB+RR5oYi8GnAmVPVX9+tvOH3cTQPOtB5Yr6pfucdv4jQKG2Lh\n3xNOY/mN+/uCYH9XZwKfqeoWVd0PvI3TfR3o70pVx6rqmaraAvgdWEEBf0+x1gh8DdQVkePdTwHX\nAtMDzhRuOtDZvd8Zp08+akREgJeApar6TCzkEpGqodEHIpKC0086P8hMqtpXVWurah2c7oQ5qnpD\nkJlEpIyIlHfvl8Xp614UZCZV3QCsE5GT3FMX4ewhPiOoTNl05FBXEAT7/285cI6IpLj/Dy/CGXQQ\n6O9KRI5yvx4L/AuYSEF/T9EqYhSg2HEJTgHmB+ChAHO8jtP3txenTnEzTnHoA5zWdhZQMcqZzsPp\n416A80Y7H2cEU2C5cPah/tbNtBB4wD0f6O8qLF8LYHrQmXD63xe4t8Whf9tB/55wRpV8hbOX+Ns4\nxeLA/+6AssBmoHzYuaB/Vw/iNJKLgHFAyRjINNfNtABnlF6Bf082WcwYYxJYrHUHGWOMiSJrBIwx\nJoFZI2CMMQnMGgFjjElg1ggYY0wCs0bAGGMSmDUCxhiTwKwRMMaYBPb/3bPCxYdH1WQAAAAASUVO\nRK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107360160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Theta: 44.997815081691805, Max. Range: 57.1394215709871\n"
     ]
    }
   ],
   "source": [
    "#P193: This is the program corresponding to Figure 7-5\n",
    "'''\n",
    "Find the angle at which the projectile has maximum range for\n",
    "a fixed velocty, u. This also shows the intermediate points.\n",
    "'''\n",
    "\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "from sympy import Derivative, Symbol, sin, cos, solve\n",
    "\n",
    "def plot_range_theta(u):\n",
    "    g = 9.8\n",
    "    angles = range(0, 90, 1)\n",
    "    R = [u**2*math.sin(math.radians(2*angle))/g for angle in angles]\n",
    "    plt.plot(angles, R)\n",
    "    # Use LaTex for the X-axis label\n",
    "    #plt.rc('text', usetex=True)\n",
    "    plt.plot(angles, R)\n",
    "    #plt.xlabel(r'$\\theta$ : Angle of projection (degrees)')\n",
    "    plt.ylabel('R: Distance traveled by projectile (meters)')\n",
    "\n",
    "def grad_ascent(x0, f1x):\n",
    "    theta = Symbol('theta')\n",
    "    epsilon =  1e-6\n",
    "    step_size = 1e-4\n",
    "    x_old = x0\n",
    "    x_new = x_old + step_size*f1x.subs({theta:x_old}).evalf()\n",
    "\n",
    "    X = []\n",
    "    while abs(x_old - x_new) > epsilon:\n",
    "        X.append(x_new)\n",
    "        x_old = x_new\n",
    "        x_new = x_old + step_size*f1x.subs({theta:x_old}).evalf()\n",
    "\n",
    "    return x_new, X\n",
    "\n",
    "def find_max_theta(R, theta):\n",
    "    # Calculate the first derivative\n",
    "    R1theta = Derivative(R, theta).doit()\n",
    "    theta0 = 1e-3\n",
    "    theta_max, X = grad_ascent(0.001, R1theta)\n",
    "    return math.degrees(theta_max.evalf()), X\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    g = 9.8\n",
    "    # Assume initial velocity\n",
    "    u = 25\n",
    "    plot_range_theta(u)\n",
    "    theta = Symbol('theta')\n",
    "    # Expression for range\n",
    "    R = u**2*sin(2*theta)/g\n",
    "    theta_max, X = find_max_theta(R, theta)\n",
    "\n",
    "    # calculate R for all theta's traversed\n",
    "    Y = [u**2*math.sin(2*angle)/g for angle in X]\n",
    "    X = [math.degrees(angle) for angle in X]\n",
    "    plt.plot(X, Y, 'ro')\n",
    "    plt.show()\n",
    "    print('Theta: {0}, Max. Range: {1}'.format(theta_max, R.subs({theta:theta_max})))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter a function in one variable: 25*25*sin(2*theta)/9.8\n",
      "Enter the variable to differentiate with respect to: theta\n",
      "Enter the initial value of the variable: 0.001\n",
      "theta: 0.785360029379083\n",
      "Maximum value: 63.7755100185965\n"
     ]
    }
   ],
   "source": [
    "#P195: Generic program for gradient ascent\n",
    "'''\n",
    "Use gradient ascent to find the maximum value of a\n",
    "single variable function\n",
    "'''\n",
    "from sympy import Derivative, Symbol, sympify\n",
    "def grad_ascent(x0, f1x, x):\n",
    "    epsilon =  1e-6\n",
    "    step_size = 1e-4\n",
    "    x_old = x0\n",
    "    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()\n",
    "    while abs(x_old - x_new) > epsilon:\n",
    "        x_old = x_new\n",
    "        x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()\n",
    "    return x_new\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    f = input('Enter a function in one variable: ')\n",
    "    var = input('Enter the variable to differentiate with respect to: ')\n",
    "    var0 = float(input('Enter the initial value of the variable: '))\n",
    "    try:\n",
    "        f = sympify(f)\n",
    "    except SympifyError:\n",
    "        print('Invalid function entered')\n",
    "    else:\n",
    "        var = Symbol(var)\n",
    "        d = Derivative(f, var).doit()\n",
    "        var_max = grad_ascent(var0, d, var)\n",
    "        print('{0}: {1}'.format(var.name, var_max))\n",
    "        print('Maximum value: {0}'.format(f.subs({var:var_max})))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter a function in one variable: log(x)\n",
      "Enter the variable to differentiate with respect to: x\n",
      "Enter the initial value of the variable: 0.1\n",
      "Cannot continue, solution for 1/x=0 does not exist\n"
     ]
    }
   ],
   "source": [
    "#P198: \n",
    "\n",
    "'''\n",
    "Use gradient ascent to find the maximum value of a single-variable function. \n",
    "This also checks for the existence of a solution for the equation f'(x)=0.\n",
    "'''\n",
    "from sympy import Derivative, Symbol, sympify, solve\n",
    "def grad_ascent(x0, f1x, x):\n",
    "    # check if f1x=0 has a solution \n",
    "    if not solve(f1x):\n",
    "        print('Cannot continue, solution for {0}=0 does not exist'.format(f1x))\n",
    "        return \n",
    "    epsilon = 1e-6\n",
    "    step_size = 1e-4\n",
    "    x_old = x0\n",
    "    x_new = x_old + step_size*f1x.subs({x:x_old}).evalf() \n",
    "    while abs(x_old - x_new) > epsilon:\n",
    "        x_old = x_new\n",
    "        x_new = x_old + step_size*f1x.subs({x:x_old}).evalf()\n",
    "    return x_new\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    f = input('Enter a function in one variable: ')\n",
    "    var = input('Enter the variable to differentiate with respect to: ') \n",
    "    var0 = float(input('Enter the initial value of the variable: '))\n",
    "    try:\n",
    "        f = sympify(f) \n",
    "    except SympifyError:\n",
    "        print('Invalid function entered') \n",
    "    else:\n",
    "        var = Symbol(var)\n",
    "        d = Derivative(f, var).doit() \n",
    "        var_max = grad_ascent(var0, d, var) \n",
    "        if var_max:\n",
    "            print('{0}: {1}'.format(var.name, var_max)) \n",
    "            print('Maximum value: {0}'.format(f.subs({var:var_max})))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "k*x**2/2"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P200: Basic example of finding the indefinite Integral\n",
    ">>> from sympy import Integral, Symbol \n",
    ">>> x = Symbol('x')\n",
    ">>> k = Symbol('k')\n",
    ">>> Integral(k*x, x).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2*k"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P200: Basic example of finding the definite integral\n",
    ">>> from sympy import Integral, Symbol \n",
    ">>> x = Symbol('x')\n",
    ">>> k = Symbol('k')\n",
    ">>> Integral(k*x, (x, 0, 2)).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.135905121983278"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P203: Probability of the grade falling between  11 and 12\n",
    ">>> from sympy import Symbol, exp, sqrt, pi, Integral\n",
    ">>> x = Symbol('x')\n",
    ">>> p = exp(-(x - 10)**2/2)/sqrt(2*pi)\n",
    ">>> Integral(p, (x, 11, 12)).doit().evalf()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.00000000000000"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#P204: The probability density function is such that the integral between -infinity and infinity is 1\n",
    ">>> from sympy import Symbol, exp, sqrt, pi, Integral, S\n",
    ">>> x = Symbol('x')\n",
    ">>> p = exp(-(x - 10)**2/2)/sqrt(2*pi)\n",
    ">>> Integral(p, (x, S.NegativeInfinity, S.Infinity)).doit().evalf()"
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}