Lagrangian.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Rungie Kutta that is totally typed out that at least runs.\n",
    "#\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib as mpl\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.constants as sci\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Truncate( number, digits = 2 ):\n",
    "    #radixPosition = str( number ).find( '.' )\n",
    "    #exponentialPosition = str( number ).find( 'e' )\n",
    "    #if radixPosition != ( -1 ):\n",
    "    #    numberBuffer = str( number )[ 0 : radixPosition + digits ]\n",
    "    #    #print( numberBuffer )\n",
    "    #    if exponentialPosition != ( -1 ):\n",
    "    #        numberBuffer += str( number )[ exponentialPosition : ]\n",
    "    #    number = np.float64( numberBuffer )\n",
    "    #if number == 0.0:\n",
    "    #    number = .1\n",
    "    #try:\n",
    "    #    number = np.float64( int( number ) )\n",
    "    #except:\n",
    "    #    number = 1.0\n",
    "    #print( number )\n",
    "    return number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def LengthVelocity( previousLengthVelocity ):\n",
    "    return Truncate( previousLengthVelocity )\n",
    "\n",
    "def AngleVelocity( previousAngleVelocity ):\n",
    "    return Truncate( previousAngleVelocity )\n",
    "\n",
    "def LengthAcceleration( length, angle, angleVelocity, springConstant, mass, gravity, initialLength ):\n",
    "    return Truncate( length ) * Truncate( angleVelocity ) ** 2 - Truncate( springConstant / mass ) * ( Truncate( length ) - Truncate( initialLength ) ) + ( Truncate( gravity ) * Truncate( np.cos( Truncate( angle ) ) ) )\n",
    "\n",
    "def AngleAcceleration( length, angle, lengthVelocity, angleVelocity, gravity ):\n",
    "    return ( ( -gravity ) / Truncate( length ) ) * Truncate( np.sin( Truncate( angle ) ) ) - ( 2.0 / Truncate( length ) ) * Truncate( lengthVelocity ) * Truncate( angleVelocity )\n",
    "\n",
    "def Spring( step, timeEnd, mass = .04, springConstant = 100, initialLength = .3, initialAngle = .4, initialLengthVelocity = 0.0, initialAngleVelocity = 0.0 ):\n",
    "    time = 0\n",
    "    lengthZero = 0.2\n",
    "    gravity = 9.8\n",
    "    cells = int( np.ceil( timeEnd / step ) )\n",
    "    length = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    angle = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    lengthVelocity = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    angleVelocity = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    time = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    x = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    y = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    length[ 0 ] = initialLength\n",
    "    angle[ 0 ] = initialAngle\n",
    "    lengthVelocity[ 0 ] = initialLengthVelocity\n",
    "    angleVelocity[ 0 ] = initialAngleVelocity\n",
    "    x[ 0 ] = length[ 0 ] * np.sin( angle[ 0 ] )\n",
    "    y[ 0 ] = 1.0 - length[ 0 ] * np.cos( angle[ 0 ] )\n",
    "    for i in range( cells ):\n",
    "        #print( '<k1>' )\n",
    "        k1Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] ) )\n",
    "        k1Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] ) )\n",
    "        k1LengthVelocity = Truncate( step * LengthAcceleration( length[ i ], angle[ i ], angleVelocity[ i ], springConstant, mass, gravity, initialLength ) )\n",
    "        k1AngleVelocity = Truncate( step * AngleAcceleration( length[ i ], angle[ i ], lengthVelocity[ i ], angleVelocity[ i ], gravity ) )\n",
    "        #print( '</k1>' )\n",
    "        #########################################\n",
    "        #print( '<k2>' )\n",
    "        k2Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] + k1Length / 2.0 ) )\n",
    "        k2Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] + k1AngleVelocity  / 2.0 ) )\n",
    "        k2LengthVelocity = Truncate( step * LengthAcceleration( ( length[ i ] + k1LengthVelocity ) / 2.0, ( angle[ i ] + k1AngleVelocity ) / 2.0, ( angleVelocity[ i ] ), springConstant, mass, gravity, initialLength ) )\n",
    "        k2AngleVelocity = Truncate( step * AngleAcceleration( ( length[ i ] + k1LengthVelocity ) / 2.0, ( angle[ i ] + k1AngleVelocity ) / 2.0, ( lengthVelocity[ i ] ), ( angleVelocity[ i ] ), gravity ) )\n",
    "        #print( '</k2>' )\n",
    "        #########################################\n",
    "        #print( '<k3>' )\n",
    "        k3Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] + k2Length / 2.0 ) )\n",
    "        k3Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] + k2AngleVelocity / 2.0 ) )\n",
    "        k3LengthVelocity = Truncate( step * LengthAcceleration( ( length[ i ] + k2LengthVelocity ) / 2.0, ( angle[ i ] + k2AngleVelocity ) / 2.0, ( angleVelocity[ i ] ), springConstant, mass, gravity, initialLength ) )\n",
    "        k3AngleVelocity = Truncate( step * AngleAcceleration( ( length[ i ] + k2LengthVelocity ) / 2.0, ( angle[ i ] + k2AngleVelocity ) / 2.0, ( lengthVelocity[ i ] ), ( angleVelocity[ i ] ), gravity ) )\n",
    "        #print( '</k3>' )\n",
    "        #########################################\n",
    "        #print( '<k4>' )\n",
    "        k4Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] + k3Length ) )\n",
    "        k4Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] + k3AngleVelocity ) )\n",
    "        k4LengthVelocity = Truncate( step * LengthAcceleration( length[ i ] + k3LengthVelocity, angle[ i ] + k3AngleVelocity, angleVelocity[ i ], springConstant, mass, gravity, initialLength ) )\n",
    "        k4AngleVelocity = Truncate( step * AngleAcceleration( length[ i ] + k3LengthVelocity, angle[ i ] + k3AngleVelocity, lengthVelocity[ i ], angleVelocity[ i ], gravity ) )\n",
    "        #print( '</k4>' )\n",
    "        #########################################\n",
    "        #print( '<taylor>' )\n",
    "        length[ i + 1 ] = Truncate( ( length[ i ] + ( k1Length + 2.0 * ( k2Length + k3Length ) + k4Length ) / 6.0 ) )\n",
    "        angle[ i + 1 ] = Truncate( ( ( angle[ i ] + ( k1Angle + 2.0 * ( k2Angle + k3Angle ) + k4Angle ) / 6.0 ) ) )\n",
    "        lengthVelocity[ i + 1 ] = Truncate( lengthVelocity[ i ] + ( k1LengthVelocity + 2.0 * ( k2LengthVelocity + k3LengthVelocity ) + k4LengthVelocity ) / 6.0 )\n",
    "        angleVelocity[ i + 1 ] = Truncate( angleVelocity[ i ] + ( k1AngleVelocity + 2.0 * ( k2AngleVelocity + k3AngleVelocity ) + k4AngleVelocity ) / 6.0 )\n",
    "        #print( '</taylor>' )\n",
    "        #########################################\n",
    "        #print( '<xyt>' )\n",
    "        x[ i + 1 ] = Truncate( length[ i + 1 ] * np.sin( angle[ i + 1 ] ) )\n",
    "        y[ i + 1 ] = Truncate( 1.0 - length[ i + 1 ] * np.cos( angle[ i + 1 ] ) )\n",
    "        time[ i + 1 ] = step + time[ i ]\n",
    "        #print( '</xyt>' )\n",
    "    fig = plt.figure()\n",
    "    ax = fig.gca( projection = '3d' )\n",
    "    plt.plot( x, y, time )\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\ipykernel_launcher.py:9: RuntimeWarning: overflow encountered in double_scalars\n",
      "  if __name__ == '__main__':\n",
      "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\ipykernel_launcher.py:9: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsvXmQXOV5PXx636ZnXzSafR9tSJqRkOQQB4NtXCaxK16wXC7LDjgO5GO1XWV5S8B2QMLYGAN2hbIChBgExqkoPywoxzhKYmwkwKwS0sz0Mj090zM9ve/rvd8f4n11u/t2Ty+3l+m5p2pKRTPT+z33uc9zznkkLMtChAgRIkRUH9JqPwERIkSIEHERIiGLECFCRI1AJGQRIkSIqBGIhCxChAgRNQKRkEWIECGiRiASsggRIkTUCERCFiFChIgagUjIIkSIEFEjEAlZhAgRImoE8gJ/X7T1iRAhQkThkOTzS2KFLEKECBE1ApGQRYgQIaJGIBKyCBEiRNQIREIWIUKEiBqBSMgiRIgQUSMQCVmECBEiagQiIYsQIUJEjUAkZBEiRIioEYiELEKECBE1ApGQRYgQIaJGIBKyCBEiRNQIREIWIUKEiBpBoeFCIkTkBMuySCaTAACZTAaJJK9MFREiREAkZBECgWEYJJNJJBIJRKNRertEIoFMJqM/UqkUUqkUEolEJGsRItIgErKIksAwDBKJBK2KASASiUCj0VDCJUSd/nfhcBgtLS2Qy+UiUYsQAZGQRRQBlmXBsizi8TgYhqG3rayswGw2Q6FQIBaLgWVZaDQa6HQ6aLVa+q9MJkMsFoPZbIZOp0MsFqP3LZFIIJVKIZPJRKIWseEgYdmCMufFgPoNDJZlaUVMiJhhGCwtLWFhYQHt7e0YGBig5MmyLCKRCILBIP0JhUJgGAYKhQKRSASDg4PQ6XTQ6XSQyWSU7Pm+l9zWB2l/iEQtYp0gry+pSMgi1kQ6EUskEiQSCSwsLGBpaQnd3d3o7++HQqEAy7KIxWI5SZJlWfj9fly4cAFdXV2UqJPJJFQqVUpFrdPpIJfLU4iaZdmU++cSdHqfWoSIGkFeX0axZSEiK4hiIpFIUBKMx+OwWCyw2+3o7e3FgQMHIJPJCrpfiUQCtVoNhUKB/v7+lMeLxWK0mrbZbAgGg0gmk1AqlZSgCWGTEwAZKKYTNSFovoGiCBG1CJGQRWSAj4hjsRhMJhPcbjf6+/tx4MABSKXFy9hJSyP9NpVKBZVKhdbW1pTnE4/HKVEvLy8jGAwikUhAoVCkkLROp4NSqUzpc3Mr9mg0ilgshtbWVpGoRdQcREIWQUF6vh6PB83NzZBIJAiHwzCZTPD7/RgaGsLk5OSaxCU0sUkkEiiVSiiVSrS0tKT8v1gshlAohGAwiNXVVZjNZsTjccjlcl6iDofDcLvdaGxsRDweT7kvsaIWUW2IhCwCLMsikUggkUggEonAZDJhbGwMRqMRsVgMQ0ND2LZtm6DExFchFwNC1M3NzSm3x+NxStROpxMWi4UqP0glTghbqVTSv0skEhlETbTUXNWHaHoRUQ6IhLyBwTVzABeJx+/3w+PxYHZ2FsPDwxkVqVAQipCzQaFQoKmpCU1NTSm32+12rK6uQi6Xw+VyYWFhAdFoFFKpNKNHrVar6d+lt3BYls1ZUYtkLaIYiIS8AcFn5nA6nTCZTJDL5dBoNJieni76/tOHa7l+r9KQyWRQq9Xo7u5OuT2ZTNKK2uPxYHFxEZFIBFKpNEXxodVqodFoUv4u3fRCtNRyuVwkahEFQSTkDQI+MwdwsWI0mUzQ6/XYtm0b1Go1Xn311bI/n1ojJplMBr1eD71en3I7l6h9Ph9sNhsikQgApBC1TqeDWq2mr4thmBQLOXDxve7u7hZNLyKyQiTkOgefmYNlWdhsNlgsFrS0tGD37t308pxhmBTCLuVxc5FMuVsWQiEbUTMMg1AohFAohEAggJWVFYTDYQCg7kTyo9FoIJVKYbVa0dXVhWQymeJOJI8jml5EiIRcp+AzczAMA6vVisXFRXR2dmLPnj0pAy2gdKLkuvTqGVKpFA0NDWhoaEi5nWR0cJUfoVCIKliIXZy0P6RSaU4tNV+PWlR+1C9EQq4z8GmIiavOZrOhp6cH+/btg1zO/9FX6kCvFqGU+0TBHQ52dHSkPO6ZM2eg1+up8oPYyNVqdcZAkdjIRdPLxoJIyHWCbK46s9kMh8OBvr6+olx19YhqkBYZ9HV0dGQQNcn7CIVCcLvdCAaDYBgmRZpHiJprI083vQAiUa93iIS8zkE0xGazGT09PZBKpfTS2OPxYHBwEGNjYyW56kSUDxKJBBqNJkW5AVz8XKPRKHUnLi4u0ryPdBs5N+8DyNRSsywLt9uNrq4ukahrHCIhr1NwzRwAsLi4iObmZpjNZoRCIQwNDWHLli3iAVcjKLRVQvI+1Go12traUu6Hm/eRzUZOfhQKBeLxOBYWFtDW1sbrTswm0RNReYiEvM7AZ+bw+XwIhUI4f/48RkZG0NraKh5QNYZ8tdlrIVveB4AUorbb7QgGg4jH45DJZIhEIlhaWkohagLR9FI7EAl5nYDPzOF2u2E0GiGTyaDVarFz506oVKoqPsvaR7XUH0IRci5ky/vw+/2YmZmBRCJJyfuQyWQZFTVXdSOaXioPkZBrGNnMHA6HAyaTCRqNBpOTk9Dr9Xj99ddTyFpEdlSDOCpByNkgk8mgUqnQ09OTcnsikaAVtdPppDZycoLnEjX3RM9nehE3vQgDkZBrENnMHCsrK5ifn0djYyN27NgBrVZL/0YqlQpi6OA+h2IPpI2gQy4U1SRkhmF4h7pyuZw37yORSFAdtdvtzmoj57oTuRK9dNNLJBKh1bdoeskNkZBrCNnMHNwVSVNTU7xtCSEJmRxghR4wDMPAZrNhfn4+I6ynoaEhpW+50UD6stUAwzAFyR3lcjkaGxvR2NiYcjvXRu71erG0tIRIJAKJRJJB1GTJLcuyMBqNGBwcpAUE+W6JppdMiIRcA+DTECeTSeqq6+7uxuWXX56T0GQymWCETMg9XwIhJ435+Xm0t7dj586dYFk2xa1mMplSlAANDQ30NWczqdQTqlkhJ5NJQU4Ga9nIg8Eg/H4/lpeXU/I+AoEA/H4/pFIpNBqNaHrJgfo/EmoY2Vx18/PzsNvt6Onpwf79+/MiLKlUKlgPmdh51wLDMFhcXITFYkFHRwf27t0LpVJJe97Nzc0ZOcVECRAIBBCPx/Hmm2+m7NLj/pTDxFKtVgq54qnWY5ezOl/LRv72228jHA7D5XKl5H2kp+iR791GNr2IhFwFEC2p3W5He3s7XZFkNpvhcrmKWpEkdMsi130lk0ksLi5iYWEBnZ2da1bvXHCVAMvLy5ienqbvRyAQoCYI4lbj2oobGhrogVvq66s0arGHXG6QtpVMJsPw8HBKEh53G7nD4aB5H3xETSpqgH+BQD0RtUjIFQTXzMEwDObm5tDQ0ACj0Qi/34/BwUFMTEwU9UUqR8siHclkEgsLC7BarXm1UfIFV1ubboLgHrgulwvBYJAeuISkuYlqtYqNSMgEfC0JrVYLrVab1UZOPu/0vA8uWedD1A888AAOHz68buYXIiFXAHxmDpII9vbbb2N4eLjkFUnlbFmQcKLFxUVs3rw57zZKqeDaitvb2+ntLMsiHA7T1ofdbk+5FCYkzR0uVRsbmZDzRa7PO5uNnLS6uETNtZH/6le/wre+9a1qvaSCIRJyGZFu5pBIJPB4PDAajWBZFkqlEvv27RPkscrRskgkErBYLFhaWiqon11ukKl+eoVFepbc4VI4HE6Ra5Grk0oT5EYm5FJfd742cpvNhmAwSK/kfvOb34BlWfzxj3/E1q1bs64ju/766/Hcc8+hs7MT77zzDgDA5XLhM5/5DMxmMwYHB/HMM8+gpaUFLMvitttuw8mTJ6HVavHYY49hamoKAPD444/j+9//PgDg29/+Nr7whS8U/Fqrf3TVGbKZOVwuF4xGI5RKJcbGxtDY2Ig//OEPgj2u0Dpki8UCl8uF3t7eglPiqkU8XKldZ2cnvT1drkVS1fikeUqlsizPv9qyt2qdSMs5RM1mI2dZFiMjI1CpVHj99dfx5JNP4uzZs/jABz6Af/zHf8y4ny9+8Yu4+eabcejQIXrbkSNHcPXVV+Pw4cM4cuQIjhw5gqNHj+L555/H7OwsZmdncfr0adx00004ffo0XC4X7rrrLrz66quQSCSYnp7Gxz72sYJ3UoqELBD4zBzAxbU9JJR869atGZNooSCTyTLcU4UiHo9jfn4eDoejKCKuVXDlWhKJBJFIBAMDA0gmk7S6crvdsFqt1KnG7U+nW4qLQbVlb6U+/2JRjepcIpGgra0N11xzDe677z48/PDDOX///e9/P8xmc8ptJ06cwKlTpwAAX/jCF3DllVfi6NGjOHHiBA4dOgSJRIL9+/fD4/HAZrPh1KlT+NCHPkRPDB/60Ifwwgsv4LOf/WxBz10k5BLBZ+ZgWRbLy8uYn59HS0sLdu7cmRGvyP17IQ7UUirkWCxGpXb9/f3YtGkTOjs7K0LG1bYU8xkguJZibvaDXC5PIen0kJ5cqLbsrVon1mQyWbXH9nq9GS7EfLGyskIX4XZ3d8NutwO4mKrY19dHf6+3txeLi4tZby8UIiEXCT4NMdHlWq1WdHR08K5I4oKQqBBf2GIIORaLwWQyweFwYGBggErtzp8/X3L7o5DN07UwdOMim6U4Ho9Tad7KygqNvcyWT8zFRu0hV5OQfT5fhg6+VPC1YLJFBRTzeYuEXCCyueoWFhawtLSEzZs351yRxIVMJhPsC0vuKx9Eo1GYTCa4XC4MDAxkBNhXKouiWnrgYqFQKNDS0pLSFyTzAkLU3MES1+xCvi/VwEYlZI/HU3SF3NXVBZvNhu7ubthsNjqT6O3txcLCAv09q9WKzZs3o7e3l7Y4yO1XXnllwY8rEnKe4A7quCuS5ufnsbq6WlTPtRASXQv5VMiRSAQmkwlutzun5lnoAWE2VCuESMgTgUQigVKpRGtra8ZgiSvVcjqdtFdNNLWk/SGE2SUXNiohl9Ky+NjHPobHH38chw8fxuOPP46Pf/zj9PaHHnoIBw8exOnTp9HU1ITu7m5cc801+OY3vwm32w0A+M1vfoN77rmn4McVCXkNcM0cL7/8Mg4cOEArTI/Hg4GBAYyOjhb1ha8UIYfDYZhMJni9XgwNDWFycjInKVWKkIHqWZnLjXSplkajQSAQwODgIDU/BAKBFJcaX0CPEES6kQk5n5bFZz/7WZw6dYoOs++66y4cPnwY1113HY4dO4b+/n788pe/BAB89KMfxcmTJzE6OgqtVotHH30UANDa2orvfOc72Lt3LwDgH/7hHzIWCOQDkZCzgM/MkUwm8c477yAYDAqyIklIQua7r1AoBJPJBJ/Ph+Hh4byfb75ZFqWi1nrH5QS5qspmfiB2YtL64JpdsiWp5QuRkHPjqaee4r39xRdfzLhNIpFkVW1cf/31uP766wt7kmkQCTkNfGYOn88Ho9GIaDSKzs5OdHZ2CkIm5aqQQ6EQDAYDgsEghoeHsXXr1oKe71pZFvn8fb6/V68VcjrWIkWunTj979KT1NLNLqT1oVKpeN/7ahJyIpGoag+Zq3xYDxAJGdnNHMRVJ5VKMTw8DIZh0NTUJFhlJzQhx2IxvPXWWwiHwxgeHqbBRcXcV732kNfbCqdsSWpcs4vH40nRUGu12hR5XiKR2JAVcjlUFuXGhibkbGYOsiJJrVZjYmKC5r/KZLKMHWOlQChCDgQCmJ2dhd/vx8jICNra2krOxagUIVcD9ZD2li2bmGz7CAQCcDqdsFgs8Pl8ePvttzM01JUwixC1STXg8/mKHupVCxuSkIl0LZlMpgj2V1ZWYDabeVckARf1qULurSuVkP1+PwwGA+LxOB0WcfuSxaLUlkUh2Cgti0rpkPm2fbzyyivYvn077VFzzS5kYQC39SGkzTqZTFbNtp1vD7mWsKEImU9DzLIszfZtbW3F7t27oVaref9eLpfXRIXs8/lgMBiQSCQwMjKC1tZWsCyL2dlZQZ5XPbcsqoVqZlkAoJkP6RUjN5xneXmZ1+xCiLqY1kO1h3qFZklUGxuCkLO56siKpK6urjVddYCwPV9yf4XkT3i9XszNzQEARkZGUs7+QlZfxaosWJaF0+mEwWBALBZL0dvyhctvNEKuRVUJd2EAQXqKWnrcJbf1QQLks0Ek5MJQ14TM1RADoCuSLBYLVlZWCnLVAdWrkN1uNwwGA6RSKUZHR8veFyumZeFyuTA3Nwe1Wo0tW7ZAoVCkONicTidCoRCAS5nF0WgUkUiEbi6uBKo51KtWhVzoa86VokbMLoFAICVAniwMSF/JVE1CDgaD0Ol0VXnsYlGXhMyyLAKBAKLRKHQ6HSQSCaLRKE0yK2ZFElB5Qna5XDAYDJDL5RgfH88IwSkXCmlZuN1uzM3NQalU0jQ7hmFohaxWqzP0tuFwGIFAAIlEgvYy0xPWyrmluhqVajXDhYRCrlxikkPNXckEXLTpK5VKNDU1oaGhoWILA8hJaD0E83NRV4TMNXN4PB643W709fVRcwRfbkMhkMlkdJuuEOAjZJZlKRErlUpMTk5mTNLLjXwI2ePxYG5uDnK5vKDnyM0gdjgcGBwcpNIsUnmRLdXxeJzmQXDbHusxErSaLYtyP26uhQF/+tOf0NDQkNXsQk7A5bhKIkac9YS6IGQ+M0c8HsfKygo8Hg+GhoYKNkfwoZwqC27/Va1WlzU7eS3k6u2SPrZUKk2RBBYL8jh8CWvcXmYgEMDCwkLKJTIh6Vpa1ZQNtdpDLidI4dPV1ZVSBHHNLl6vF0tLS4hEIhkLA3KZXdbCer0iWdeETA5YUs1JJBJ4vV4YjUZ6ybx3717BPphytCwSiQRWV1dhNBqh0Wiwffv2kvpeQhz4fBWyz+fD3NwcWJYVrI+91vPM1cskbY909xpXGUA2gNQCqkXIlZIvZgNf75xrdunq6qK3E7NLIBDgXRjA/WwVCkXO99Pn81X8ylIIrGtC5n4gbrcbRqMRcrmcrm85e/asoAeBkITMsizcbjecTicUCgWv7rlQkMpWSEL2+/2Ym5sDwzAZyo5SUazKgnuJnL6qiZuuNj8/n6K1bWhoQCgUqgpJV5OQ10sfNZfZhfu5WiwWxGIxyOXyDGkemTt4PJ51p0EG1jkhA8Dq6ioMBgN0Oh22bNlCL/OTyWTGWvBSIYTsjWVZrKyswGQy0S/R9u3bBX1+pR6AEokEsVgMb7zxBuLxOEZHR8siHxJa9pZtAwi37eH1ehGNRrG8vAy1Wp0yRBQqXY0P1STk9dhz5yLXwgBC1Ha7nQ6Kf/Ob3+DChQvweDx4+eWXsXXr1rwH4vfffz9+/vOfQyKRYMeOHXj00Udhs9lw8OBBuFwuTE1N4YknnoBSqUQ0GsWhQ4fw2muvoa2tDU8//TQGBwdLe60l/XUNgGEY3hVJ5TA3lFIhk7VOZrMZTU1N2LVrF1QqFU6fPi3Y8xPiNQcCAVy4cAE+nw+7du0qKkIwX1SKoLhaW0K4mzdvTklXW11dpQMnvrZHqc+1WpXqeqqQC4VCoUBzc3NGJdzf349nnnkGv/3tb/HYY4/h3LlzuPvuu3HFFVfkvL/FxUX85Cc/wblz56DRaHDdddfh+PHjOHnyJO644w4cPHgQN954I44dO4abbroJx44dQ0tLC+bm5nD8+HF8/etfx9NPP13Sa1r3hNzd3c1LQuU42IshZJZlYbPZYDab0dLSkuEEFLpCLLaCDwaDMBgMCIfD6O/vB4CSyDjf979aumBuDGa6MoBv8Sm5POZW1IVYgjdiy6Jag7VNmzZhaGgIf/7nf47vfve7Bf1tIpFAOByGQqFAKBRCd3c3fve73+HJJ58EcHHh6Z133ombbroJJ06cwJ133gkA+NSnPoWbb7655M953RNyJT/wQipQhmFgs9kwPz+P1tZWTE9Plz1kpZgKmUR1hkIhGkwUi8WKWtBYKGrRqSeVSnn7mOTyOBAIYGVlBQaDIcO5xudGJNiIhFxtl16hg+eenh587WtfQ39/PzQaDT784Q9jenoazc3N9OTLXV7KXWxK2ipOp7OkPJl1T8hrQcgvZD4HFMMwWFpawvz8PNrb2/OyZAuFQgg5HA7DYDAgEAhgZGQkJaqznrMsiiVGvstj4lzjcyOm62yrVS1uVEL2eDwpCo584Ha7ceLECZhMJjQ3N+PTn/40nn/++YzfI5+jUItNuVj3hJzrDVAoFDQopdwg2RgLCwvo7OzE3r17Kz7Nl8lkee3VMxqN8Hq9GBkZwbZt2zLeQyE2hmwE3S3XuZbuRkzX2Xo8HgSDQTQ2NlbEjch9LhuRkH0+H8bHxwv6m9/+9rcYGhqiLaxPfOIT+MMf/gCPx4NEIgG5XE6XmgKXFp729vYikUjA6/WWPHNZ94ScC6TnW05iTCaTsFqtsFqt6OrqwuWXX17wQSYUeZHsAD6QBacej2fNdU5CbAzJ5/XUw5JTPvDpbN955x0MDAyAYZgUNyL5fqa3PYQiso1KyMUEC/X39+Pll19GKBSCRqPBiy++iD179uADH/gAnn32WRw8eDBj4enjjz+OAwcO4Nlnn8VVV10lVshrVchCS98IWbEsi4WFBVitVnR3dxcUUsQFaQ8I8cXlazWQhawulwvDw8NrLjgFKkeUtdhDLhdYloVCoYBareZ1I5K2h9BuxI26vqmYLOR9+/bhU5/6FKampiCXy7F79258+ctfxrXXXouDBw/i29/+Nnbv3o0bbrgBAHDDDTfg85//PEZHR9Ha2orjx4+X/LzXPSHngtDOOuBiW8BoNGJ5eRmbN2/G/v37SwrgJsoIoQk5FovBZDLB6XRicHAQExMTBe26qxQ2CiFn6yFz3YjpgT3Zdumlqz1yXQEKoUsvFuutQgaAu+66C3fddVfKbcPDwzhz5kzG76rVarqNWijUNSGTHrIQILGdZE9XqURMIPTm6Wg0ipmZGRrcMz4+XrMri+q9x8xFoW0piURC9dB8bkSyoonPjcgNlK+mMaSa20LW4/omoA4IOdeXXC6Xl9yyiMfjsFgsWF5eRm9vLzo6OtDd3S3YF00oQo7H43A6nQgEAhgdHcX+/furagio1R7yeltymo5cbkTS9lhcXEQwGKRXS6QCJ22PSn0vqlkh+/3+isXVCol1T8i5UAohx+NxmM1m2O129PX1Yf/+/ZDJZHTFjVAolZBJpvDKygoaGhowNDSE3t5ewZ5fOVGvQz0+lFt1olQq0dramhHCRAaH6fGXfG0PoZ9fMpmsWrgTy7Lr0jK+7gl5raEe0YXmi1gsBrPZjNXVVd4g+1rZPJ1IJDA/P4/l5WX09fXhwIEDWF5eLmglVLkQiURgt9vpJfd6PDCERjVkgBKJhCal5ZOqxnUjCrHwtFoV8nqeS6x7Qs6FQoZ60WgUZrMZDocDAwMDWTeKVHvzdCKRwMLCApaWltDb20srd6Byho5s4A4S29vbqfaWYRhotdqUA50Ekq/ng6cQ1JJTL1uqGteNuLy8jEAgUJAbMR3VImQyQF2PM4q6IORsB3Y+sjeuLGxwcHDNjSLV2quXTCapzC6diAmqRcjctgl5D+PxeIqjiVRk3EDyRCIBjUaDRCJBybpeq+laIuRsWMuNSAaJ2dyI6Vs/qkXIgUBgXWYhA3VCyNmQizyJUcLtdmNoaChvWVilWxZc48laMjuht2KvRSLJZBIWiwWLi4sp7Z30kyNXMcC9dDabzUgmkzSAKRAIZOhvGxoaBF3vU82KvNYJmQ+FuBEjkUjKbsRwOFyV99vr9a7LgR5QJ4ScrULmI+RwOAyTyQSv14uhoaG8jBLp9ylknzYbiXKt2Js2bcrLeCJkhZwr7J5hGCwuLsJisaC7uxsHDhwoqhKSyWRQqVTo7u6mt6VvA7HZbPRA55J0Kf3N9XgpWyzKJXvLtvWDuxsxFArR5QbldCOmo5hgoVpBXRByNnAJORQKwWg0wu/3r2kdzvc+hQDRDhNwya6rq6sgB6CQhEzui1tdcaNEOzo6irKJc8F3Is22DYQc6MQkQfqbJGSe/NT6br1Ko9LGEG6YvN1ux/bt2yGXy7O6EdPbHkJ8fut1WwhQJ4S8VibD22+/jWAwiOHhYd4wnUJQrpYFNyWus7OzKLITsmXBbT2wLEs3szQ3NwsWJVrIUC/bElQSMk9iMYmbjUvShWYX1xNqIcui3G7EdBRjm64V1O23NBAIwGg0IhwOY2JiIiVeshQIrbKQSqXwer14+eWX0d7eXlJKnNAtC4Zh4HQ6MTc3B51Oh127dmVsZikVpfQYs4XMcy+budnFarUayWQSOp0OjY2NG6KariYhrzWDKMaNyCXpbJJKkZBrCH6/H0ajEdFoFCMjI/D7/YKRMSBcy4KsdJqdnYVMJhMkrlNIQk4mk3jjjTegUqlK3oSdDeWSveWqpk0mE6LRKIxGI0KhUEo1Rn7KHYlZSVR7hVMxx12hbkSNRkM/Q7/fD5fLlULw6wl1QcgSiYRuR04kEhgZGaGOJVLRCnXJWiohkyWnRqMRra2t2Lp1K2w2myCOpnzykNcCeR9DoRC2bdtW1i92JXXIpJom8iwyiOJWY9xITKK95fam1+NuumoTspDI5kYMh8P0M7znnntw+vRpyOVy/N///R927NiBO+64Y80BosfjwZe+9CW88847kEgk+Jd/+RdMTEzgM5/5DMxmMwYHB/HMM8+gpaUFLMvitttuw8mTJ6HVavHYY49hampKkNdYF4RMtjiPjIxkJDwRLXK1sydYloXdbofRaERzczOmpqagVqsRDAYF7fsWe19kIh6JRDA2NgaLxQKtVlv0c6nVLIt08FVj6dpbsgCVXGKvp2q62i2LcoM7BO7o6MDPf/5zfPOb38RHPvIR9Pb24vz583mpOW677TZ85CMfwbPPPotYLIZQKIS7774bV199NQ4fPowjR47gyJEjOHr0KJ5//nnMzs5idnYWp0+fxk033STYsuK6IOSurq6UYQEXQqsiCt2mwR2INTU1ZSw5FXoQV2iFHIlEYDAY4Pf7MTo6ira2Nkgecy1ZAAAgAElEQVQkElit1qq6/qqJbNpbbjXtcDhSqmkuUefrZKsEWJat2rbravXnvV4vOjs7sXPnTuzcuXPN3/f5fPjf//1fPPbYYwAubSk/ceIETp06BeDictMrr7wSR48exYkTJ3Do0CFIJBLs378fHo8HNpstRb5ZLOqCkCsdUp8PWJaFw+GAwWBAQ0ND1oGYkIRcyAHAtTkPDw9j69atKX9f6hqnfOM311PaW7ZqmvQ2051sXJLeaCe3asZ+FpqFbDQa0dHRgb/5m7/Bm2++ienpaTzwwANYWVmhJNvd3Q273Q4gdbkpcGnxqUjI72GtCE6hQ+pzgWVZOJ1OGAwGaLVaXHbZZTkv/YV2162FdJtztrzkUtc45YN6SHvLJuliGCZFKRCJRHDmzBlqkKjFalpIVDucvhCVRSKRwJ/+9Cc8+OCD2LdvH2677TYcOXIk6++XY7kpQV0Qci6Ui5D5KkBCxGq1Om9lghALRfNBNptzrudVr4RcCUilUhrgw7IsPB4P9u7dm7Klen5+PiUXQq/X06paCJ03UD2rOFkKWg34/f6CnHq9vb3o7e3Fvn37AACf+tSncOTIEXR1ddFWhM1mowNustyUgLv4tFTUPSGXo2VBqlryhXO5XDAYDFAqldi6dSsaGhoEfbxSUKzNWQiy3Aibp/MBd7CWrZrmxmEuLCwgFotR3S3XLr5equlqVsgMwxR0Mti0aRP6+vpw4cIFTExM4MUXX8TWrVuxdetWPP744zh8+HDGctOHHnoIBw8exOnTp9HU1CRIuwKoE0Jeq2UhdEYwqbqJREyhUGBycrKmEqZKtTmXUiGT6MO1CL2eK2Qu1joxcZ2FXHB70wsLCwgGgwCQEmW6Vrh8tU6I1cxCLuY79eCDD+Jzn/scYrEYhoeH8eijj4JhGFx33XU4duwY+vv76f68j370ozh58iRGR0eh1Wrx6KOPCvb864KQgdIiOAsFwzB48803oVQqa5KIic65FJtzvbYsqnECKPZKgU93m15Nc8PluS2PalfT1Q6nL/T93rVrF1599dWM21988cWM2yQSCR5++OHinuAaqBtCzgYhe8g+nw+zs7MIBoMYGxtDT0+PIPcLCHN5Tyb8dru9ZJtzpXrb632olw+EbN3kqqbJEJHrYiOLF8jfqFSqirz+ahFyKBQqi6u0UqgbQs5VIZdKyKQ1wTAMRkdHsbS0lKIlLhWkGi32C+z1ejEzM0N7jpOTkyUbFkpVWQSDQSgUipwOxI3SX65EL51oZ7lyr1gshrfeegsajQZerxeLi4u0mk7vTQtNntUiZI/Hs26jN4E6IuRsKGXRqd/vh8FgQDwex+joKP2y2+32sqxxKvQLzD1RTExMoLGxEa+++qogrYZiWxZ+vx8zMzM0wS6RSGSNyBR7yOV/XKVSia6urpTM4ng8TnvT3GqaG4VZ6mIAIeMKCsF6DhYCNgghF1ohBwIBGAwGxGKxlFyMUu4zFwrVIqfbnLlVkRB5FkDhhBwOhzE7O4toNIqxsTFotVpKtnwRmTKZDEqlEolEAj6fr66XoVbLtZYtC1mhUKClpSXle5PvYoB812yRXXyVxnoOpwfqiJCzfeEL6YUGg0EYDAZEIhGMjIxktWNXa/N0NpszF6XkWXCRb8siFovBaDTC7XZjdHSUJuvFYjF6P9kiMpeXl2G322mVxrJsioJAr9cLvka+WkO9atmX833cXIsByMnUZrPR7JW11mxVq2UhEnIdIBQKwWAwIBQKUSKupJRuLUJey+bMhVDqiLXuJ5lMYn5+HjabDYODg3nvJCQgqoBgMIiJiQkAmQoCi8WCeDwOpVIJvV6f4m4rpeJcz0O9QiBEsJBcLuddfMqtpkmwPKmmdTodQqFQQfZloSC2LGoExXzhw+EwDAYDAoEARkZG8s5NlsvlVBMqBLIRcr42Zy6EJGS+apJrNOnp6eHdfg3kJ2lL/x0+BQE3K8Lv92N1dZXmGHMr6VJ27JUb65mQ+bDWmi1SUc/NzQFARddsiYS8DkAu4wlxkEt/n8+HkZGRgtc6lbtlkUwmsbCwAKvVmpfNOdd9FYv0lgWJDzUYDGhvby95nx55jHxIm8/dlkwm6YHP3bGn0WhoNU1aHtVWc9QbIWcDdzGAy+XC6Ogo1Gp1Rddseb3elOCf9Ya6IeR8Aobi8TiMRiO8Xu+al/65IPQaJ+5evVK3OZejZeF2uzEzMwOdTkdznKsNmUzGuxUkHA7D7/enyLy4FuRYLFYVM8pGIGQuuPv0+GYI3BMqyQlPV+TodLqC21NihbwOIJVKceHCBQQCgZI2ThOUI2PZ6XTSGMBSqk8hCTkSieBPf/oTJBIJtm3bJnhGh9CyN+6lNFfmxbUge71euN1uzM/PVyxsfiMTcjZkO6FyFwPY7faCFwP4fD6RkGsVZBjmdrsxMDCAHTt2CLboVKi9equrq7BYLNBoNIJscxaiZREOh2EymeDz+bBr166yDWcqpUPmWpDj8TgdUoVCIdqXJmHzpEIjbY9StLgE1SLGbLK3Wn3sfBYD5FqzpVAo4PP5RJVFLYB70MRiMZjNZjgcDgwODqKnpwdNTU2CVSlCkB53m/PAwABYlhVEt1lKhUxaOk6nE5s3b85wfpUD1TKGyGQyGo/JfS6k38nV4nKdbXq9HlqttqBWUjUr5Gpqu4U83tZas2WxWHDjjTciGAzizjvvxP79+3HgwAEaqbkWkskk9uzZg56eHjz33HMwmUw4ePAgXC4Xpqam8MQTT0CpVCIajeLQoUN47bXX0NbWhqeffhqDg4OCvE6gjggZuEgoZrMZdrsdAwMD2L9/P730rqaRgwuuzZlkJq+srMDv9wvy3IohZJKVvLS0hIGBAYyPjyMQCMDn8wnynLKh2sO2dGTrd3KdbVarFYFAAMClDGNu6hofqknI1VKelPv1plfTg4ODeO211/AXf/EXuP3223H+/HlcuHAhb0J+4IEHsGXLFvqd//rXv4477rgDBw8exI033ohjx47hpptuwrFjx9DS0oK5uTkcP34cX//61/H0008L9rrqhpAZhsGrr76K3t7eDFWC0D3fYr5sfDZnAiG3hshksrw10izLYnFxEfPz89i8eXOKhK3UdsJ6WXKaD/icbUQz7ff74XQ6MT8/j3g8nnIZrdfrodFoNmQPuVpgGAb79u3D/v378/4bq9WKX//61/jWt76FH/3oR2BZFr/73e/w5JNPAri4T+/OO+/ETTfdhBMnTuDOO+8EcDHI/uabbxb0860bQpbJZDhw4ADv/1MoFAiHwxV+RheRy+ZMUOlFp6R3PTc3h7a2Nt4hohi/mRvZNNN8QynSOuCub6pEK6FahFytk2yxj3v77bfj3nvvpVepTqcTzc3N9OqC7MwDUvfpEZmf0+lM6XmXgrohZCD7AV7pvXpAfjZngkoSstvtxuzsLLRabU4JW6UIuRoo1+NmG0otLS3B7/fTpQGBQICG+XAHiEJrpqs5TFwvvevnnnsOnZ2dmJ6ephumc+3MK+c+PaDOCDkbhIjgTAcxTqR/4QuxORMI3bLgu69AIIDZ2VmwLIstW7asGapfz3nIlQbpTff29tLbWJalNnGPx5OytindJl4sqW40Qo5EIgVngL/00kv4z//8T5w8eRKRSAQ+nw+33347PB4P3QvI3ZlH9un19vYikUjA6/VmhI+Vgroi5FwVcjn26iUSCTrIKcbmzL2vclXIkUgEc3NzCAaDGB8fz1s1IcTW6bV6a+ulh1wq+MKFiLZWp9NlaKb9fj/dVh0Khejvcok6n2HdRiNkj8eTMpvJB/fccw/uueceAMCpU6dw33334Re/+AU+/elP49lnn8XBgwcz9uk9/vjjOHDgAJ599llcddVVYoVcKMrRsiBuPWJzJr2lQmzOBOUg5Hg8DpPJBIfDUZQ9vF57yNVAIUMfpVKJtra2DJs40eGurKzAYDDklbi20QhZSJfe0aNHcfDgQXz729/G7t27ccMNNwAAbrjhBnz+85/H6OgoWltbcfz4cUEej6CuCDnbl75cm6etVivsdju6u7uzhuzke19CEbJEIoHX68WZM2dSpH/F3E89kuV6jN/MpsPlyy8mmmm9Xo9IJCLE0y8Y6zV688orr8SVV14JABgeHsaZM2cyfketVtNlp+VAXRFyNghZ7ZHBjMPhKNnmTCAE+bEsi6WlJRiNRkgkkqJyMIR8Th6PB9FoFI2NjVkNL9Ui/XqI38yWuEY0036/H6FQCOfOnYNUKk2xHuv1+rLZxIHqtizWs20a2CCELMTBQKRiBoMBzc3N6O7uRmdnpyBf7FKeH8uycDgcmJubQ0tLC3bt2oWZmZmSD4hin1MgEMCFCxeo4oAE/KhUKuqM40Yw1mMVno5K6pC5mmmXy4WtW7fSuFjSlzabzdR6zO1LCxWLWQ8ti2qhrgi5XF96rs2ZbHM2GAwVl9Klw+PxYGZmBmq1mj6vWCxW9t4vH6LRKObm5hAIBDA2NkaT1ciBmUgk4Pf74ff7U9Y4hcNh2Gy2mlhdXy5U2xgilUp5beJEM53+maTnTBdKriIhF4+6IuRcyCZTywU+mzOB0JnIhSAYDGJmZgYMw2RI2ITsR+eDZDIJs9mM5eVlmqTHMAw9KFmWpT3UpqYmNDY20ljGaDSKN998E/F4HBaLJUNRQKq39b5rr1rDtVyPm00zzV3ZVOxqrWoRss/nw6ZNmyr+uEKirgg5n0zkfHa05bI5p9+fkFirkiJmE1KF8ukfK6GOAC71rM1mM90cAlw8GCUSCWQyGZVmkeeTTCYpQTMMQyuyTZs2YfPmzXSRAMk0TjdRcEm62F171RrqVWvrdKGPy7eyiXxWfr8fbrebaqa5zkMSuiSRSJBIJKDVaoV+OWtCrJDXEfIh5Hxsztz7E3KvHqng+SqLRCIBk8mE1dXVNc0mlTjwnU4nZmZm0NLSgr1796ZsupZIJLyLV7n/RiIRzM7OIhaLYXJyEnK5nFbVAGjAT1dXF+0zkxQ2bg+UG5Wp1+uhUqlqLrAIqB4hA8J8H8hQkHuFCCCl5eFwOGh2McMwaGpqooRdqYAjkZBrDLm+fLmkb4XYnAmEblkQXTOXkBmGgcViweLiIvr7+4uWsAkFMrCTyWS47LLLoFarwTAMXXO/1ntGzDPExci3w5BhGFpBk38BQKVSQalUpvwNd9fe0tIS3Q7CXeHEt3GiHlQWtYBsq7XeffddyGQyXs00+WzKcfIUCXkdga/FUIzNmXt/5VjjBFyS1plMJmzatAn79u2r6gJP7sBufHwcjY2NlCzzIWLS3rBYLOjt7cXevXuznljI7eknJj6S5l5eS6VSerlMnG4OhwOhUChlUFWNoWe9EjIfZDIZZDIZurq66GyDq5nOtlqLnDxLKTi8Xm9VNl0LiQ1DyNw8i1JszgTlWOOUTCbhcDgwOzuL5uZm7N27t+heqRBgWRZzc3Ow2+0pAztCxPkcPE6nEwaDAS0tLdizZ09RMsF8SJp7gmhsbIRer0dvby8laZIbEQgEqCyQa0fW6/VlO+ltJEIGMod62TTT3NVaFouFbnJPt4nn+51Z7+ubgDoj5LWGetFoFGazGVartWibM/f+hCRkhmHw9ttvQ6fTYefOnVUZihCQnORgMAiZTEZDvsnALp+qmChBZDIZduzYUXDoy1rIRtLkXy5hE+WGVqtFIBCg1RsZVK2ursJoNPJurRZii8tGJ+Rs4K7WImAYBsFgsKjVWrFYrCYW8JaCuiJkgN/9xTAMfD4fHA4HtROXWg0JJS8LhUKYmZmBz+fDyMiIYCvMiyUBUqG3trZCr9ejp6cn58AuHbFYDEajEX6/H2NjYxWtWNKHhwQk28NqtcLn86Gvrw8sy1LZV2dnJ31d0WiUbq22Wq1UTcBnaskX1ZC9VdNwU4rsLZtmmgx1A4EAtYmTNVzJZBJOp7Og93hhYQGHDh3C8vIypFIpvvzlL+O2226Dy+XCZz7zGZjNZgwODuKZZ55BS0sLWJbFbbfdhpMnT0Kr1eKxxx7D1NRUUa8xF+qOkLlgWRbLy8swmUzQarXo7u7G8PCwIPddaoUcjUZhMBjg8/kwNjYGjUYj2NmdnJQKIQ2/34+ZmRnI5XI6sDt//jxeeeUVNDQ00DZAtkt7MoC02WwYHBzExMRETVSFxMloNBppP55IA/n60sTpxh3s8plauJkRa5laqlEh1/LG6UKx1mqts2fP4mc/+xksFgv27NmDyclJ/O3f/i3NpeCDXC7HD3/4Q0xNTcHv92N6ehof+tCH8Nhjj+Hqq6/G4cOHceTIERw5cgRHjx7F888/j9nZWczOzuL06dO46aabcPr0acFeI31egt9jlUFkN1yb8/T0NILBIFZWVgR7nGLzgkn/2m63Y2hoCFu2bIFEIoHH4xE8EzmfAzIajWJ2dhahUAhjY2MpA7uJiQmwLItgMAifz4eVlRXMzc0hmUxCq9WisbERDQ0NiEajWFhYQFdXFy6//PKaMXIQVYharcbu3btT2g+FDA+JqaWpqYkOD5PJJFV4LCwsIBgMZjW1bDRCBiqjZCEnzyuuuAJ79uzBtddei5deegkXLlzIkOilo7u7G93d3QAAvV6PLVu2YHFxESdOnKBB9V/4whdw5ZVX4ujRozhx4gQOHToEiUSC/fv3w+PxwGaz0fsQCnVHyC6Xi34gxE4MXLyUFjrxrRAwDIOFhQXav06XsJGhnhDIxxzCPTGMjIygo6ODd2AnkUh4LyHJCe7cuXNgWZauYJ+fn4der88ZKlRuxONxKmMcHx/POwGskOEhAFoZE+ch19SyvLyMQCCAZDKJeDwOm82G1tbWkkwthaDahFxpkKQ34qotBGazGa+//jr27duHlZUVSrLd3d2w2+0AUlc3AZfWOomEvAai0Si2bduWsusMqM4aJ+BS28RoNKKrqyurhK1Sa5y4i017e3tTBnZAfn1iMhyNxWKYmppCQ0MDlTb5/X66ASMajUKtVlOC1uv1vMMYocAwDBYXF2G1WgVrmxQyPAT4TS1vv/02lEplVlNLtiFVKdiIhFzMvCIQCOCTn/wkfvzjH+cMty/36iaCuiPkzZs38xJbOdY4Abn7g06nE7Ozs2hsbMSePXtyVoxCOv+ykTt3YJePwy4dXGPHyMhISgYCV9pENmCQABufzwe/34/FxUVEIhEolUpK0I2NjYKkjLlcLszOzqKtrQ179+4tq2471/AwveURj8cRi8XQ1taWYobIZmrhOg/5TC35YqMtOPX5fAVvC4nH4/jkJz+Jz33uc/jEJz4BAOjq6qKtCJvNRmV6ZHUTAXetk5CoO0LOhnKEAZE2Q/rB7/P5aCjRZZddlpeErZwtC7/fjwsXLkChUGDnzp1QqVQFOewKMXZwwQ2w4epPiZIhfUjGraR1Ol1eZBQOhzEzMwMAuOyyywSX1xUCbjXNfc8GBgag0Wjoew5c2ljc0tJCPwM+UwtRHRSavlaLgUblRKFZyCzL4oYbbsCWLVvwla98hd5OVjQdPnw4Y3XTQw89hIMHD+L06dNoamoSvF0B1CEhVzLjgbj1CCGHQiGa0VBI7xIoT8uC7NMLhUIYHx+HXq8vyGEHXIoebW1tLdrYkQ5iueVW2PF4HH6/n8oTiQaaVIuNjY0pSoZkMkldltmClqoFn8+HCxcuoLGxkbdazzY8JKaWxsZGOjxkGIY3fW0tU0u2XJRyY71Eb7700kt44oknsGPHDuzatQsAcPfdd+Pw4cO47rrrcOzYMfT399PtIB/96Edx8uRJjI6OQqvV4tFHHy3L66g7Qq4kSF9aIpHAYDDA4/FgbGwshWgKuS8hrdgWiwWBQACjo6Nob28v2GFHtlST3IpyV54KhSLDJMCVm83Pz1Mnl0wmQygUwqZNmzA1NVXW7ReFIB6P04WyW7ZsyZhjEBQ6PCSmlu7ubkrSa5la8lXZCI31QshXXHFF1vbKiy++mHGbRCLBww8/XPTzyxd1R8j5XH4LVS1LpVKYzWZ4PB4MDQ1hcnKy6PsWomVBBnYrKyt0gAgUNrCrprEjHXK5nG6/AC4edKT10t/fj1AohNdffx0Mw+SllS4XuO2JwcHBor4HhQ4P1zK1OJ1OmsZWiqmlUFSTkIXyGFQTdUfIuUDaAqUerAzDwGq1YnV1Fd3d3SVZsNOfW7FYXV2lrYWenh7o9fqCBnZcY8fQ0FDNGDuAS+FGkUgEW7duzag8id02m1aatDzKUUmv1Z4oBYUMD4FLVxltbW3Q6XQIhULo6ekpydRSKNZLhVyr2FCETCI4iz1oWJbFysoKjEYjOjo60Nvbi5aWFkG+zMUSMhnYKZVKOrBzOBwwGAwwm830wCPElH6wsCwLu90Ok8lUc8YOcpIg20g6Ojp4TxLZ7LYkE8HhcMBkMiEej9Oge/J+FKuVJu2JUCiUsz1RDqzV8ggGg1QjS0wtzc3N9MRciKmlUFRzW4hIyDWIfLaGFAMiq9Lr9ZienoZKpaLBJ0KAK0HLByTkPRwOZwzs2traaN84EAjA5/NhaWmJbt8gl/dSqZTus5uamqpqshwXxO5sMBjQ1dVFJXqFQCKR0LQwMg0nWmmfzwe32435+XkaSMOtpHNl9ZK20MLCQtHtiXKADHItFgscDgfddMNXSbMsS98brvOQz9RS6KaWalbIhQzRaxV1R8i5UIwWmWQ8yGSysu7Vy/dkwd0estbATiqV0qk9AcMwcDqdMBqNVPtKdgeS3+WrpCsFkhKnUCiwa9cuQdO7uFppsnuNBNeQy3qilSYbmcn7odFoqJyxqamp7FrnQuFyuTAzM4NNmzZlyBL5+tLp67QA0DyVrq4uSvCkL+10OjE/P494PJ7T1JJIJKryvtRDFjJQh4S8VoWcr306HA5jdnYW0Wg0q4RNSDMHX0odF8SFZrFY0NfXV9TALpuxg1tJ22w2ukA1vd1RzgMtHo/DaDTC6/VifHy8Ypef3OAaPq00ubrwer1gWRYdHR3Q6/WIRCJ5a6XLiVgshpmZGSQSCezcuXNNNQxfX5pveEgIW6lU0isuPlMLSV4jppZ4PA6NRlPx/A6RkGsY2cgtnwqZqAzcbjetQLN9seRyOZViCfGc+UAu32dnZ9He3o7LL788xfiRr7GDXGb39fVlVFDZKmkyKFteXsbs7CwYhoFOp0uppEslae5zGxgYKGpZQDlAVkZFIhFEIhFMTk6itbWVVtKrq6vUuMFtdwg5IMsF7vs2MjKScjIpFIUOD8mmlnRTCwmaD4VCWF1dLdrUUgwikUhVTUFCoS4JORtyVcjcdfb5qgzK4f7jgkzwVSoVdu3aVbDDDkg1dhRymc0dlPX09ADIVDOkkzQhpXwfw+12Y3Z2li5LraUWAJHYkc0t5Lm1tbWl7JBL10oHAoGUQCaSiCckEfn9fpw/f77srZO1hodc5yF5zTqdDk1NTSkzDD5TC7flUar6hRRftXAiLxW1cwRUAHwtBm4roKenpyAJm9BmDgIysItEIhgfH6fhPYUQMTF2kHxjIaoHPjUDd8OD3W6nSy3JxJ4QNfegi0QitC2yffv2qm5HSUcsFsPc3BzC4TCvxC4d6VppAFTF4PP5sLi4CL/fTwdppWilk8kkDAYDvF4vJicnUz6HSiEbSSeTSczPz8Pj8aTkyRRiaknfIF4oREKuUeRqWZAKmci9DAYDbQUUeqYWOkGOZVnMzMzA4XAU7bCLxWIwGAwIBAIVMXZwSZqErXAzlLkkrdFokEgkEA6HMTo6SgdrtQBuC2BoaIimtRUDmUxG85MJSLVIVAzcq4tsJy4uiM68r68PY2NjNUU+pGLv7OzMaKnxDQ/VajVUKhU6OztprjgJoSpmU0u1BonlQH28ijxBCNTtdmNmZgY6nQ5TU1NFT/KFyp8gRpNgMFh0JGYymcTCwgI1dlRTjsWVnG3evDklgpT0WRcWFmA2myti3lgLpD1RztYJt0/PbQGFQiH4fD66UCGRSKS8J0qlEiaTCVKpFFNTU1XLmOZDMpnE3Nwc/H5/hgKJFg8SCTyhOByBOOz+CByBOFYDUTgCUTgCMTgCMTiDcTiDcXx4Sxu++eGRgje1eL3egpPeahV1ScjZiCgWi2F1dRXxeJw3M7lQlFohsyxLK5+Ojg40NTWhq6ur4IHdysoKTCYTuru7a8rYAVwyrmi12owt2izLUkIi5g1CSISgy0nShbYnhIZUKqUnLgLue2I2m+Hz+aBUKqHT6WC1WvPSSpcbLMtiYXkVfzo3B2VTByTqPrz7jhuOwDIcQUKy7/0biCPJc7WqUUjR3qBEu06J4Q4d9g4qsKevibY1gEsnsVymlv/4j//Av//7v4NlWTz00EPYvXs3pqeniy6yXnjhBdx2221IJpP40pe+hMOHD5f0XhUKSYH5pdXbnFgAEolESuUaDodp6AvDMHjf+94nyOOwLIs//vGPRd0fGdip1WqMjo5CqVRiZmYGLpcLDQ0NaGpqWnMgRPTDDQ0NGBkZqRljB3CJ7EKhECYmJvLud3IJieQoE4cdV91RymtlWRZWqxVWqxXDw8MpeRC1AFKxt7W1YXBwEFKplGqlyXuSTStdyuuIJpJwBOKUTC/+RCnJOgIxrPqjWA1EEeO5MJRLJWjTKS8SLflJ/+/3fnTK3LVgNoUHF8TU8vLLL+Phhx/GX/7lX+KNN97AzTffjK1btxb8+pPJJMbHx/Ff//VfNGb2qaeeKuq+eJDXB1PXFTLRtrpcLoyMjKCtrQ1nzpwR/HEKQbq+mTuwGxsbA8uyGQMhACkVo0wmg9FoRCKRqLhtdy2Q9svi4iKGh4fpzsB8QSy8Op0uxWFHSNrpdKbYoLntjnxIuhLtiWLBtWNv27YtpQWQTStNCHp5eRmhUAgKhSJ1Q4tGC084wUuu9Oe923wR/qu9Fq0C7Tol9EoWfdo4rhhoQ19nE9obVGhvUKLjPZJt0iggFTC4C8gvEe93v/sdIpEIvnXu9usAACAASURBVPzlL5f0mGfOnMHo6CgNKTp48CBOnDghFCHnhdr5NgoIkpW7tLSEwcHBFG1rtTYaJBIJGI1GOBwOjI2Noa2tLesOu3RNMLlMc7vdOHv2LEKhENRqNVpbW6lhoVL611wgG1I6OjoEbZ1kI2lig3a5XNQGnd7uICQdi8XoiTCd7KoNbtspHzs2y7LwRxNwBBJwBOVwhHRwhBVwBBtg90ew4vVj1e+AKxSHN8ryXtZqlTJawY526LB/qOUSueqUlGxbdQokYlGcP38eKpUKY2PViztNJ2m73Y6vfvWrkEqleOCBB0q+f769eeXYLJ0LdUnILpcLALB///6q91NJxbiwsID+/v6iBnYSiYS6ovr6+tDT0wOGYeglLNG/cgdHZHhWiUvxUCiEmZkZSKXSvNxiQiCbDZovqwK4SMibN2/GyMiIoHbsUhEKhXD+/Hmo1Wps37kb3iiLtxZ9WOWpYLk/sWTmJbxcKqGkOtChx/R75NqilUMnTUKNGJRMBPJ4CBqlLEWGx9caY1kWCwsLWFpawvj4eM0sAWBZFr/61a/wgx/8AHfddRf++q//WpDveaX25uVCXRJyZ2dnRW2UfGtr0gd2xTjsANCQnfRdcTKZDM3NzSmyNjKV9vl8MJlMdOsGl6RL2dOWDpKr4Xa7MTY2VnXrajpJezweXLhwAXq9Hi0tLQgEAjh79ixisRgNcyfvS7nUCwmGgTsYv0SwwUu9WIvdDbsvgjCrgCscQCC6kvmaALTqFLSaHWprRhv57/eqWPLTpJbn/dkmk0mqYLBarQgEAgCQsoHEarXS1k61CxuClZUVfOUrX4FOp8N///d/F7UMIhsqtTcvF+qSkHOBBMEL9QUj5hAuIZNhm1qtxu7du6FUKgt22AUCgZS9fPlUnXwmBe5qJIPBgFAoBLlcnkLShQ6DWJaFzWbD/Pw8+vr6MDo6WlNDMW57Il2OBVwKFPL5fBkbsrmDw2yVNMuy8EYSaVVrlEq4uL1aVzDO2zLQyIFWrQLdzXoMpRErGYR1NCjRolNAXoZWFN8JnWEYejIn6g6n04lIJJKXVrqcYBgGzz77LH74wx/ie9/7Hj7+8Y8L/p3bu3cvZmdnYTKZ0NPTg+PHj+PJJ58U9DHWQl0Scq4PiuRZCEnIiUQCCoWCDuzITj0SEl8IEXONHYXu5eMD32qkWCxGSXplZQWhUIhugiY/2dbSezyelE3atbI+CchfPcENFOJuyHb7gpi3u3Fuzo4llxGOYAyhpAwhVg5/XApvlIEzdJGIE0wmzSplUkqqPc0a7OxtSlEZNKuk8Nmt0MkYXLZtsuayF3w+H86fP49NmzZh586d9IqOq5Umw2SNRlPwQLVYLC8v44477kBjYyNOnTqVYl0XEnK5HA899BCuueYaJJNJXH/99di2bVtZHisb6lL2xrIs7R2m480338TIyIhgyoS33noLfX19sNvtcLlc1GFHHEr5EnEymYTFYsHKygqGhoYqLsWKxWJUaubz+RAOh6FSqVIu6S0WC+LxOCYmJmpqKAZcPFHMzMygtbUVQ0ND9IQbTzJwBeMpKoPVLCqDEI+WSyoBmjVyNKkkaJCz0MsZtOkU6GzUYHNrA/raG9HdokOHXgW9ir9lkB4ElC1ov1pIJBIpuwDXsrIT1Qtpefh8vpQ2kFBaaYZh8Mwzz+D+++/HP/3TP+Gv/uqvaup9KxB5PfENR8hnz55FT0+PIJZihmFw5swZxGIxDA8PU1daKcaO/v7+qqslCKLRKL2k9/l8UCgUNDym3L3XXGBYFt7wRZK1eUI4a1iAIxCDTNcMb5RNIVl3iD9MqlEtz66T5agMWrQKyKSXPkOuzTddE5ze7pBIJAgEAjh//jz0ej1GRkZqSmYHXLJkDwwM0LyJYsDNlRZCK22z2XD77bejtbUV999/f80MFEvAxiVkAFlzimdmZtDS0oKOjo6i75ubgyGVSjE4OIiOjo6CB3bk8l+v12N4eLimjB1kKGk0GrFp0yb09/dDIpGkkJHP50M0Gk25fG1qair6dQRj6X1ZfpWBM8jfMlDJpVQTy0euF4dhKrQ3KKCSCzek4pI0IaRIJIJEIgGWZTEwMIDOzs6sbaBqIBaL4fz58wCAiYmJsp1Y09+XcDgMhUKR0u7gDpoZhsFTTz2FBx98EHfffTeuvfbamnnPSsTGJuRYLMYrYzEajdBoNFTPWiiIsUCr1WJkZASLi4sIh8Po7u7OO8GL9JqTySTGxsZqytgBXBooqlQqjI6O5jxYuQMy8kP0wI2NjdDoGhCTqjMq10wnWBzheGbLQCaRoK1BkeIA08tZJANu9LTpsWW4B12NmvfcX7KaOHhJLkVXVxd0Ol0KSXMrxly9+nKBO5AtNUe5WHBnGH6/H6FQCKdOncK5c+dgsViwefNmPPLII1V5bmWESMh8r81isUAikaQIwPNBOBzGzMwM4vE4ddgxDEPzMciXC7jkqmtqakoxbBCZGOk1l2s4USzi8TgMBgP8fn/OgSLDsnCH4ikV62oaya76L6oOfFH+8KUmtRzt+tQKNtWY8N4gTHvJ/cVVT9RiHzsSieDChQuQSqUYHx/nPZGltzvC4XDKQFUIC3Q2hMNhvPvuu9BoNBgbG6uZ9gnDMDh27BiefvppbN++HZFIBOfOncOxY8ewc+fOaj89obCxCTkej/P63202G8LhMLVH5nM/xH5NHHa5BnZE30kOOmLYkEqldC370NBQzfSJgUvmlbl5KxraeyDRNMKRNgjjVrWuYK7AGFVq20CnRNt77YIGOQMVE4UsHkI46EcikUjZQJIt3J5rx67FoRh5fktLS/Q7Ugj4LusJSZOTeykkzTAMTQKcmJioul6ci8XFRdx6663o6enBfffdVxebo7NAJGQ+Ql5dXYXb7cb4+HjOvydfYqvVioGBgaIGdsBFY8fc3BwN6g4EAgXJzEoFCYzJlmOw7AnB7gvDGwNiycyPV8jAmHRwc5MJISWTSeoga2xspKHs6eqJWgFpYQn9/LiqF3JZT/KBC9GP+/1+vPvuu2hra6upQoBhGDzxxBP42c9+hnvvvRfXXHNNTZ1kywCRkPkI2ePxYGlpKWtgCHdg19nZiYGBgaIcdqQPq1QqMTo6mmEyiMVi8Hq99KCLRCIpxoRcw7Ekw8IVSh1yFRoY06yRQy9n0aSWoL+jGV1NmoxBWEeDEo0auWCBMfmAbCBxuVw0qFylUlFlB5FVVZuYSXsnGAxicnKyIu0TQtKkkiZhQnwDsmQyCaPRCI/HU3MBVFarFbfccgsGBwdx7733lqy1XyfY2IScHsFJEAgEYDAYeHtTRMuq1WoxOjoKhUJRMBFHo1HqiBsbG8v7y8YwDBzeIOZX3FhweLHkDsIVjCPEyBFi5fDFJfBGGDiCcbhDMfCIDKAjgTGcnzZdamBMi0YG3+oivO6LLZhakxPxtSe4CXjcXj23ktbr9RVbLkpkiqVKxYQAd0BG2h1E9dHe3o7BwcGa2I4NXPxsH3/8cTzyyCO477778MEPfrAmnleFIBIyHyFHo1G88847mJ6epreRcByu6oFE/JVq7IjEkxlVa/oAzPlevzZbYEyrRoZGlQQNcgaNCqBDr8KmZi16WvXo62hEZ6MGbToltMrsVSPZ2mE2m9Hb24uenp6auXwlIJtcyOV1riqYu7eO9OrJphK+gaoQ4AYBjY2N1ZRLEbhYtc/OziIcDmPz5s1UF0zs8tx2h5CZJvnAYrHglltuwcjICH7wgx9UZR9glbGxCTmZTPJu80gmk3jllVewf/9+etlJwnHWGtgRJJhU95dhcRWGJQcSci2iUhWcHPIN8KgM0gNj0kNiuD3axrTAGHJJT9odayk7gEvZGsScUGtEEo1GMTs7S12AxS49zTZQ5VbRxZA0wzAwm81YXV3FxMRETQ6eSJttcHAQmzZtyvjuxuPxjHYHCZ4i351yVNIMw+DRRx/Fz3/+c/zwhz/E1VdfvZGqYi5EQuYjZJZl8Yc//AE9PT1YXFzMOrDzRhL4f28tY8WXOQxzh/gDY/QqeUa7oBKBMXxEJJPJoNVqEQqFwLIstmzZUnNVSSXUE9wEPJ/PRxPwuNViLiJyuVyYmZmh5phau6ogUjuZTIbx8fGCTDnc4CkuSef73qyF+fl53HzzzZicnMTRo0drqo9dBWxsQmYYhm6YJiD9v7feegsjIyMYGBhI2eHFrYqffMWK7z8/mxIYcykkRoJk0INGBYsdY4Po62hEm04JtaI2FAAMw8BkMsFms6G5uRnJZLKiyo58QNoTpM9ZySFdNiJKl9+Rqn1ysvaCgLj5GGNjY4LFUHLfG7/fz3sC02q1OU9MRFf86KOP4v7778eVV165UatiLjbuCic+kGzchoYGaLVaDA4O5uwTq+QXv3C//v/2oaf5okIiHo/DbDZfNHbsrj1jB5DqEjtw4EAK0XGVHUtLSwUpO4QCtz2xY8eOotsTpYAvAY9c0pPA/2AwCK1Wi7a2NrqVpVyGjUIRDAZx/vx5NDQ0CL6Giu+9SSQSlKC5OdvcICGtVguZTAaTyYRbbrkF27Ztw0svvVRz5p1aR90TMndgt3XrVuh0Orz66qs4e/YsmpqaaM81/UDTqS6+NaFYAgzDYHFxEVarFf39/TWX/wtcPEhJfvKuXbt4s3yVSiU6Ojpojkd6LrDFYkmxPTc1NQmWf8s1J5D2RC1BoVBApVLB6XSipaUFe/bsofnAPp8Py8vLGa66Sl9lMAyD+fl52O12TE5OVkwuJpfLeUmau7HmnnvuwdmzZ+HxePDFL34Rn/3sZ2tqM8t6Qd22LOLxOM6dO0dNIK2trXRgR2RUpFokPVdCQk1NTVgNs/jfOScu75LDZZunl9a1YjclIE5Cn8+HsbGxkgdO6VuffT5fhlmjUB1wNdsT+YBodt1uNyYnJ1P2GaYjPVwpPemNkLTQIAYU8h7WUi/baDTSqvjjH/84zp07h9deew3//M//XPZWTzKZxJ49e9DT04PnnnuurI9VIjZ2DzmRSMBqtaKrqysvhx33ktXr9SIYDCIej0OlUmFgYADt7e01l8ZGeojl1sMWo+wALpIXuToZHx+vSntiLRAnZU9PD3p7ewt+D9PjOEkCHrcVVEpMKXEq+nw+bNmypaZaAMlkEo888gh+8Ytf4Mc//jHe//73V/w5/OhHP8Krr74Kn88nOCG/8MILaGhowBVXXCHE3W3sHjIZ1CQSCUrCa20SaWtrQ0NDA8LhMA0RIpetS0tLiMfjNA+YXM5Xo9pzu92YnZ2t2Cp7qVRK+4UEXGUHWbLK7SuGQiE4HA6Mjo7WXHsCuKROkEgk2L17d9GEKZFIoFaroVaraTpZtlYQiSlN34idDWSLd09PD8bGxmqqTTY3N4dbbrkF09PT+P3vf1+Vk63VasWvf/1rfOtb38KPfvQjwe53fn4eN954IxKJBAKBAD74wQ/i1ltvrcj3uG4J+fz58/jqV78Kr9eLyclJTE9PY+/evVm3IqcbO7Zs2UIPAO6aH5K9YLPZMDMzA5ZlaQ5wtn60UIhEIpiZmQHDMNi+fXtVK85sS1YXFxdhMBggl1/UT1ssFng8nppQdgCpm5TJdhehkW1FFNmI7XK5YDabEY/Hab+e/CgUCsRiMczMzCCRSGSdB1QLyWQSP/vZz3D8+HH85Cc/Eap6LAq333477r33XnrFJgSi0Si+//3vY8+ePfje976Hubk5/P3f/z1WV1crQsh127IgiMfjOHv2LF5++WW88soreOONNyCVSrF7925MTU1hamoKv//979HV1YWpqSn09fUV1J/jVoqk1UGWiJL8hVJJKJlMUmNCuUikVEQiEZrxPDExQU96pWR2CI1yBQEVC75+PQm27+joQHd3N5qammpmbjEzM4Nbb70Vl19+Ob73ve9VVQr43HPP4eTJk/jpT3+KU6dO4b777iupZREMBqFSqSCXy3H27FlqOQeAQ4cO4aqrrsIXv/jFUp7yxu4hZwMZ6L322ms4fvw4nn32WfT29qKtrQ1TU1OYnp7G5Zdfjq6urqJJNB6PUxLyer2UhEgVTSqhfJ6r3W6H0WikPc5aGuYAF/vLFosFy8vLeakn1gq0F1LZQcANApqYmKhJg0IkEsG7774LhUKB3t5eStQkAS89prSSJ5NEIoGf/vSn+OUvf4kHH3wQ73vf+yr22NnwjW98A0888QTkcjn9Pn3iE5/Av/3bvxV8X9/97nfx/PPPY2hoCN/4xjewY8cOAKDb6a+66ip85zvfwQc+8IFSnrJIyLkQjUbxd3/3d/jGN76B8fFx2Gw2nDlzhlbSdrsdo6OjmJ6exp49e7B79240NDQURdJcEiJETbKAuSlm3IPM7/fToKORkZGaGigSuFwuzM7OoqOjAwMDA0WTRDmUHeR+yQmtFoKA+EA2ZS8uLlI1UDrIUJUbrsQwTMnvTz44f/48br31VvzZn/0Z7rrrrppqnxCUUiG/8MILuO+++/DUU0/hX//1X/E///M/+NrXvkZbMVKpFNdeey0eeeQR9PT00OKoiJnDxh7qrQWVSoX/v71zj4rqPN/2NTAKCooYz6IRlaOCESSaajUm9YCyPK8u/cXPE6ZdS43QhCoJ8YCmiahBqcYo2CRLrXQ1NjVVEqO2WttEQaJIUCAK0giiougwDKc57O8Pu3cGYxoOw8wG3+s/DMx+hwz3fvf9Ps/9fPTRR8rXffr0YcaMGcyYMQN4eHfMz88nPT2dw4cPs27dOoxGI8HBwYpIBwYGNmgn9zhP0fqP7ObNm+j1ejQaDW5ublRVVWGxWAgICPifJViOwtrLDg4Obvajq/y+3dzclNFa1pUd8u8HfrqyQ6aqqor8/Hzat29PaGioKm9olZWV5Obm0qVLF8LCwn5UUK0PVfv27Qv88PdTWVmJJEk2S8AzmUzs3LmTTz75hPfee4+RI0c2+X2qGZ1Oh6urK927d+e1116jurqao0eP0q1bNwIDA6mpqaFTp07odDpiY2Nxd3dny5YtLTaD8IndITeFqqoqLl68SEZGBhkZGVy5coVOnTopAh0WFtZkW0F+9C8uLqZLly5YLBYlpUveRXt4eDh0h2JtTzjCy/6xzA75UFX267/77jvu3LmDr6+vqqZjyMit7ffu3fvJuufGvu6jvx+gXtuzu7v7T34+c3NzeeWVVxg3bhzr1q1T5a64qVy9epW//vWvjBw5knHjxpGbm8uuXbt4+eWXCQ4OpqysjNjYWMaNG8eCBQs4d+4ckyZNwsfHh/DwcDZu3NjUSwvLoqWRJIl79+6RkZFBeno6GRkZ3Lhxg/79+xMWFkZoaCihoaF06dLlfz4qy+VN3bt3/0HjhBxKLlsdNTU1dOjQQRGghvrRzUUO2ZFD+x19ICYjt/VWVFRw9+5ddDodLi4u9OjRw2aHqrZEbuHv2bOnXcKK5JuYLNRyTOmjAUJOTk6YTCaSkpL429/+xq5duwgLC2vRtdmbzZs3k5qaysSJE/n2228JCQlh5cqVvP766zz77LPMnDkTDw8PDh48yI4dOzh79iyFhYVERkayefPm5v4+hCA7AovFQmFhIenp6aSnp5OZmYnBYCAwMJARI0YwYsQIgoODcXFx4dq1a+j1eiWpqyE7EdmP1ul0ikjLfqu1H22rP3TZnpAkCV9fX9WF7MD3w0/r6urw9/fH2dlZNZUdMiaTiWvXrmEwGAgICHBoyeKjTxppaWkcOnSI2tpaAgMD2bBhA0FBQaq56dqCc+fO8cc//pHVq1fj5eXFwYMHKS8vZ8WKFRw+fJjPP/+cyZMnM3PmTADCw8PZt28fnp6etqpyEYKsFurq6sjOzlZEOisri4qKCrRaLa+++io/+9nPGDx4cJNF1NpP1Ol0yi7IWoAaG0juaHuiIViPtB84cKAyFOBx32fvyg5rysrKuHbtGv3796dPnz6q2a3DwwqU7du38/nnn7No0SIMBgOZmZnMnTuX6dOnt9h1a2pqGDt2LLW1tZhMJubMmUN8fLzNr2OxWHBycsJoNKLRaNBqteTk5DBhwgReeOEFgoODWb16NQcOHODjjz8mMDCQK1euoNVqSU1NteWNWwiyGiktLSU8PJzIyEh8fHz4+uuvycjIoKCggN69eyt+9IgRI5qVDyyHv8g7RYPBQPv27ev50T92MCFbKD169FBdboJMZWWlkng2ePDgRu9irCs7dDrdYwesNrdyoa6ujvz8fCRJws/Pr8UOgppKTk4OK1euZOLEicTFxdl1fXKTlbu7O0ajkTFjxpCUlMSoUaNsdg1ZjK3R6XTs2bMHLy8vRo8eTVRUFN27dyclJYWrV6/y8ccfo9VqWbVqlc3W8V+EIKsR6w/io/9eXFzMuXPnlEPD8vJyfH19FYF+5plnmjV6R85ckEW6tra23i6xffv2FBYWqtqeMJvNXL9+nfLycpseiEHTMzsexXrnPmjQIKWlWi0YjUYSExM5duwY77//PiEhIQ5dT1VVFWPGjOH999+3eTVHcXExe/fu5fnnnyckJITOnTtTXV2tfLYLCgqIjIzkk08+oWvXro8VcRshBLm1YzKZyM3NVWqjL168iCRJDBs2TBFpPz+/JntccjvvgwcPKCkpQa/X0759ezw9PRUBasipvL2Qg4D69OlDv3797PLo35DKDuubZHV1Nbm5uXTo0EEZlKsmvvnmG1auXEl4eDhvvPGGQ8sBzWYzoaGhXLt2jeXLl5OQkNCs15PPCkwmk2I5vPvuuyxevJisrCzKy8v5y1/+Uu9ntmzZwvXr19m1a1ezrt0AhCC3NeTdtWxzZGRkkJ+fj6enZ73Su8b4lLI90bNnT55++mmAetGk8qGjtR9t76B2+WARaPDhZ0tiXdmh0+moqqpShLe6uhofH59mdXq2BHV1dWzdupWTJ0+ye/dunnnmGUcvSeHBgwfMnDmTHTt2MHTo0Ca9xpo1azhz5gzHjx9XrJdt27Yxd+5cioqKiI6OZurUqaxduxaTycSlS5fYvHkz9+/fJzExscnXbQRCkJ8E5G40+cDw/PnzlJaW4u3trQQqDR8+nM6dO9cTiMZUT8gCJIt0VVWVkgEst4O3xE7LHkFAtkCv13P58mU6duxIx44d0ev1qqjskLl06RJRUVFEREQQGxuryiaZ+Ph43NzciImJafTP7tixg0OHDuHp6Um/fv3YsWMHJpOJX/ziFwC4urqyatUqXnjhBW7evEnv3r356KOPePDgAb/5zW9s/VZ+DCHITyoWi4WrV68qfvSFCxeoqalh6NChDBs2jLy8PPr27cvSpUubPIbKumpBp9NRV1en5C3IVQvNKReqqKggLy9PNUFAj0P2s+Vge+t4UkdXdsDDM4MtW7Zw6tQp9uzZQ3BwcItdq7GUlZXRrl07unTpQnV1NRMnTmT16tVEREQ0+rUMBgMlJSV4eHgwZcoU4uPjiYiIYN++fbzyyivodDrg4dNgdHQ0K1as4Nlnn7X3E4wQZMH31NbWsmfPHrZu3Ur//v0xGAx06NCBkJAQxY9uTkWFXLVgXf8rt/JaR5P+1OvL9bqVlZX4+/urMggIvm+UaYyfbY/KDpmsrCyioqKYMWMGq1atUp2XnZ2dzcKFCzGbzVgsFn75y1+ydu3aZr/uvn372LlzJ2fOnMHV1ZVp06bh4uKCt7c3x48fZ86cObz55ps2eAeNRggyPAwPiYqKwmw2s3TpUmJjYx29JIeRkpLC5MmT6devH5Ikcf/+fc6fP68cGhYVFeHl5aUIdGhoKF27dm3yTsK6lVeuj350VJbcRWcdBKTGel0Zo9HI1atXqampISAgoNmVKLaq7JCpra0lISGBf/3rX+zevVtJLntSsFgsLFy4EA8PD3bu3ElZWRmZmZnk5uYyefJkAgMDHbU0Icjy6KATJ07g5eVFWFgYqampjvyfomrkIZrWfrRer68X8N/cMKFHR2VVV1fTrl07ZeyRv7+/qsYUWXPnzh0KCgoYMGAAvXr1arEbRmMrO2QuXLhAdHQ0s2fPJiYmRnW7Yntx7949pk6dyqBBg4CHh3sqKD0Ugnz27FnWr1/PF198AcA777wDPMxSFTQMo9FITk6O4kdnZ2fj7OysBPyHhYXh4+PTpMdsi8VCUVERt27dolevXpjNZnQ6nWpGZcnU1taSl5entLg74lDsxyo7vvnmG8xmM7m5uWRnZ5OcnMyQIUPsvj41cfv2bUaNGkWPHj1ISUlRi3cuBPnQoUMcO3aMvXv3ArB//37S09PZuXOng1fWepEkCb1ez9dff61YHXIwknXp3U+VfcmTqOWwIutHcutRWbLXaj0qS04ta2lLw3qQrI+Pj+qqPOrq6ti/fz8HDx5Ep9Oh1Wrp378/27dvZ/DgwY5ensOIiYnB09OTuLg4Ry/FGpGH/LibjRp9ydaEnJExfvx4ZYKCJEncvHlTCfjfs2cPZWVl+Pj4KIl3ISEhuLm5cfv2bUpKStBoNAQFBT02ZEej0eDu7o67uzt9+vQB6j/GFxUVtcioLGsMBoPSmm2PQbKNpaamhrfffpv09HQ++OADAgIClG7PplbONIYbN26wYMECbt26hZOTE7/61a+Iiopq8es2hISEBFVW5TQEdX3KbIyXlxc3btxQvi4uLlb+wAW2Q6PR0LdvX2bOnKmkZZnNZvLy8pSA/zVr1nD37l1MJhORkZFMmjSpUY/+jxuqKvvRckh7U0dlWSP76Hfu3MHf3x8PD49G/bw9yMjI4LXXXmPevHmcOnVKuVloNBr69etnlzVotVreffddQkJC0Ov1hIaGMmHCBFWcz7RWMYY2blmYTCZ8fX35+9//Tt++fQkLC+PgwYMt7rEtWbKEo0eP0qNHD3Jyclr0Wq0BSZKYNWsWvXr1Ytq0aeTm5pKRkUFubi6dO3euZ3X07du3WaV3jR2VZY1c+9ytWzdVhipVV1fz1ltvceHCBfbs2YO/v7+jUBBl7AAAChBJREFUl6Qwffp0VqxYwYQJExy9FLUiPGSAzz77jOjoaMxmM0uWLLGLr3TmzBnc3d1ZsGCBEOT/Ul5e/oN5cZIkcffu3XoB/8XFxTz99NP1Su88PDyaVXr3qB8tB7RbWx2FhYVUVFQQEBCgyiqPc+fOERMTw/z584mKilLVLrCoqIixY8eSk5OjypFjKqFtCLLJZOLNN9+kU6dOzJkzBz8/P3svoUkUFRUREREhBLmRWCwWCgoKFIHOzMykqqqqXsB/UFBQs6IiZT9ap9NRVlZWb66abHe4uLio4ryhqqqKjRs3kpWVRUpKCr6+vo5eUj0qKysZN24ccXFxzJo1y9HLUTNt41CvqqqKp556ir1791JUVERKSooSkdeCUXkCB+Hk5ISPjw8+Pj7Mnz8feFhNcOnSJdLT00lJSSEnJwcXFxeGDx+uiPSgQYMa/FlwdnbGzc2NkpIStFoto0ePxsnJ6Qd+tCNGZVnz1Vdf8dvf/paFCxeSmJioql0xPPTwZ8+ezUsvvSTE2EaofocsSRIajYbly5czZMgQli1bpsTrPfo99+/fJz8/X+n0ciRih9xySJKETqfj/Pnzyk66sLCQPn36KLXRI0aMoFu3bj/Y5UqSxO3bt7l+/XqDpozYa1SWNQaDgQ0bNpCTk0NycjI+Pj42v0ZzkSSJhQsX0rVrV7Zv3+7o5bQG2oZlAQ8L88ePH8+OHTsIDQ0FIDU1FScnJ0aNGqXERl6/fp2kpCROnz6NRqMhLi6OOXPmOGLJQpDtjCRJfPfdd4pAnz9/nvv379cL+Pfw8CAtLY2pU6fi6+vbpAoMW4/KevQ9fPnll6xevZolS5awbNky1e2KZf7973/z85//nKCgIOWm9PbbbzNlyhQHr0y1tB1BLiwsZPr06Vy8eBGtVsvKlSsJCgriypUrZGZmsm3bNkaMGFHvZ/7xj3/wpz/9ibVr1+Ll5UVaWhoffvih4snJwt5SCEF2PCaTicuXL3P27Fn2799Pbm4uQ4cOrSfSfn5+zRK95ozKssZgMLB+/Xry8vJITk5W2n4FbYYGCXKrMGCzs7Pp2rWrMqBw9+7d3Lt3j8WLFxMXF8emTZuU7y0sLMRgMNCuXTu++uorOnXqREFBAWvXrmXJkiVERESwa9cu9Ho9JpOJdevWERUVxYEDB3jw4IFN1jtv3jyee+458vPz8fLy4g9/+INNXvenuHHjBuPHjycgIIAhQ4aQlJRkl+uqFa1Wy7Bhw3jqqacYPXo0N27c4LPPPmP+/PlUVFSQkJDA6NGjmTJlCmvWrOHTTz/l5s2bj20o+l/X8PT0ZMCAAQQHB/Pcc88xdOhQPDw8qKioUG4Ily5dUkZPmUwm5eclSeLMmTNKDe+JEyeEGD/BtIodclJSEjdv3iQhIYGkpCQOHjzIkiVLOHLkCIWFhXh7e5OWloZer2fv3r1s2rQJb29vfv/73xMSEkJcXBw9e/bk1VdfBcDb25tLly5RV1fHF198QWlpKenp6fTu3ZutW7eqMsC7IZSWllJaWlqvWP/w4cOqKNZXK7KnbB2odOvWLQYOHFgv4L9Tp07NsiKqq6vrRZOmp6fzz3/+E6PRyIMHDzhw4IDqKigENqX1WxZpaWls376dixcvMmvWLJKTk4mNjaWuro7ExETl+6qqquq14JaXlxMfH8/LL7+MVqslMTGRl156iXHjxlFYWMjrr7/O+vXrCQgI+P6NSRL+/v6cPHnSbt1OLY0o1m8aFouFb7/9tl7Af11dHUOHDlVEesiQIU2uupAkiZMnT7Jp0yYGDhxIu3btyMnJYdGiRaxYscLG7+aHiMYlh9D6y97Gjx+Ps7MzZ8+e5eTJk2zcuJFFixbxu9/9jg0bNhAeHo6Liwt+fn5K2yxA165dKS8v58svvyQkJARXV1clfq+goIBu3bopj6VyxUZ6ejpBQUEYDAaHvV9bUlRUxMWLF20+xfdJwMnJCX9/f/z9/Vm0aBHwMDsiKyuLc+fO8d5773H58mXc3NzqBfw/GpL0OPR6PWvWrKGoqIjU1FQGDBig/LdGbo6ajCz8CxYssMv1BA1H1YLcsWNHJk+ezOTJk4mPj1f+PTIykiNHjhATE8PIkSNJSEggJSUFJycnXnzxRQwGA+fPn2fZsmV06NCBuro6xYbIysqie/fuSnKXVqvlwoULbNu2jUGDBqmqHbWpVFZWMnv2bLZv3y46p2yEq6sro0aNYtSoUcBD8SwvL1cC/v/85z/zn//8h379+tXrMvT09FQC+E+fPs0bb7zB8uXL2b179w/E216NKGPHjqWoqMgu1xI0DlVbFo3h+PHjpKamUlBQgJOTE9HR0UybNo2amhrGjBnDqVOn8PDwICQkhLfeeovw8HA0Gg0pKSkcOHCAefPmsWjRIlxdXZW65taI0WgkIiKCSZMmKZ65wD7I+c6y1ZGZmYler8fX15c7d+7QoUMHkpOT6d+/v6OXKqqA7E/r95CbSm1trVJqZLFY2LRpEx9++CEDBgzA29ub5ORkTCYT69evp7i4mLi4OFUW3zcWUayvPoxGI9nZ2Rw5coS1a9eqprNUCLLdeXIF+XGUlJRw/fp1Ro4ciVarZe/evfz6179mwoQJLF68mCFDhrT6+WOOLNavqalh7Nix1NbWYjKZmDNnTj2bSaAuhCDbHSHIP0VeXh45OTkcPXoUZ2dnYmNj28RO2RHIUz7c3d0xGo2MGTOGpKQkxXMVqAshyHan9VdZtDTySbqj2qvbEvKUD3j4mG40GlutD9/WmTdvHqdPn+bu3bt4eXkRHx9PZGSko5cl4AnfIQtsi9lsJjQ0lGvXrrF8+XISEhIcvSSBQC20ndZpQevA2dmZrKwsiouLycjIEI/DAkEjEYIssDldunTh+eef59ixY45eikDQqhCCLLAJZWVlSjhTdXU1J0+ebBNNNgKBPRGCLLAJpaWljB8/nuDgYMLCwpgwYQIRERF2XYPZbGb48OF2v649OXbsGH5+fgwePLheyqGgbSAO9QRthsTERDIzM6moqODo0aOOXo7NMZvN+Pr6cuLECby8vAgLCyM1NVWk+bUOxKGe4MmhuLiYtLQ0li5d6uiltBgZGRkMHjyYgQMH0r59e+bOncunn37q6GUJbIgQZEGbIDo6ms2bN6umNbklKCkpqRcN6+XlRUlJiQNXJLA1bffTK3hikLN9W3osl6N5nL0omm/aFo31kAUC1aHRaN4B/h9gAlyBzsAnkiTNd+jCbIxGo3kOWC9J0qT/fv06gCRJ7zh0YQKbIQRZ0KbQaDTPAzGSJLW5UguNRqMFvgVeBEqA88D/SZJ02aELE9iMJzrLQiBoTUiSZNJoNCuALwBn4AMhxm0LsUMWCAQClSAO9QQCgUAlCEEWCAQClSAEWSAQCFSCEGSBQCBQCUKQBQKBQCUIQRYIBAKVIARZIBAIVML/ByK+lx/tYZLnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Spring( 1, 1000 )#, 4, 100, 3, 4, 1.0, 2.0 )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hey,\n",
    "\n",
    "I've been looking at a lot of folk's Runge-Kutta codes. The most common mistake is to do something wrong with the k's . . . to add the wrong one.\n",
    "\n",
    "Lot of people are doing this\n",
    "\n",
    "k1L = h * Ldot(vLvec[i]); -- which is fine\n",
    "\n",
    "k1vL = h * vLdot(stuff) - also fine \n",
    "\n",
    "k2L = h * vLdot(vLvec[i] + k1L/2.0) -- which is not fine . . . see you added k1L/2.0 and should have added k1vL/2.0 . . . \n",
    "\n",
    "Same thing with the other functions - make sure that you are adding the right k.\n",
    "\n",
    "Also - note that in Runge-Kutta, the k4's take the full k3\n",
    "\n",
    "k4L = h * Ldot(vLvec[i] + k3vL) - is correct . . . dividing by 2.0 on the k3vL is not correct.\n",
    "\n",
    "For the most part, I got people's code to work if I fixed all the mistakes of this type.\n",
    "\n",
    "tpk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Old Code\n",
    "\n",
    "def LengthVelocity( previousLengthVelocity ):\n",
    "    return Truncate( previousLengthVelocity )\n",
    "\n",
    "def AngleVelocity( previousAngleVelocity ):\n",
    "    return Truncate( previousAngleVelocity )\n",
    "\n",
    "def LengthAcceleration( length, angle, angleVelocity, springConstant, mass, gravity, initialLength ):\n",
    "    return Truncate( length ) * Truncate( angleVelocity ) ** 2 - Truncate( springConstant / mass ) * ( Truncate( length ) - Truncate( initialLength ) ) + ( Truncate( gravity ) * Truncate( np.cos( Truncate( angle ) ) ) )\n",
    "\n",
    "def AngleAcceleration( length, angle, lengthVelocity, angleVelocity, gravity ):\n",
    "    return ( ( -gravity ) / Truncate( length ) ) * Truncate( np.sin( Truncate( angle ) ) ) - ( 2.0 / Truncate( length ) ) * Truncate( lengthVelocity ) * Truncate( angleVelocity )\n",
    "\n",
    "def Spring( step, timeEnd, mass = .04, springConstant = 100, initialLength = .3, initialAngle = .4, initialLengthVelocity = 0.0, initialAngleVelocity = 0.0 ):\n",
    "    time = 0\n",
    "    lengthZero = 0.2\n",
    "    gravity = 9.8\n",
    "    cells = int( np.ceil( timeEnd / step ) )\n",
    "    length = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    angle = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    lengthVelocity = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    angleVelocity = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    time = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    x = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    y = np.zeros( cells + 1, dtype = np.longfloat )\n",
    "    length[ 0 ] = initialLength\n",
    "    angle[ 0 ] = initialAngle\n",
    "    lengthVelocity[ 0 ] = initialLengthVelocity\n",
    "    angleVelocity[ 0 ] = initialAngleVelocity\n",
    "    x[ 0 ] = length[ 0 ] * np.sin( angle[ 0 ] )\n",
    "    y[ 0 ] = 1.0 - length[ 0 ] * np.cos( angle[ 0 ] )\n",
    "    for i in range( cells ):\n",
    "        #print( '<k1>' )\n",
    "        k1Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] ) )\n",
    "        k1Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] ) % 360.0 )\n",
    "        k1LengthVelocity = Truncate( step * LengthAcceleration( length[ i ], angle[ i ], angleVelocity[ i ], springConstant, mass, gravity, initialLength ) )\n",
    "        k1AngleVelocity = Truncate( step * AngleAcceleration( length[ i ], angle[ i ], lengthVelocity[ i ], angleVelocity[ i ], gravity ) )\n",
    "        #print( '</k1>' )\n",
    "        #########################################\n",
    "        #print( '<k2>' )\n",
    "        k2Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] + k1Length / 2.0 ) )\n",
    "        k2Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] + k1AngleVelocity  / 2.0 ) % 360.0 )\n",
    "        k2LengthVelocity = Truncate( step * LengthAcceleration( length[ i ] + k1LengthVelocity / 2.0, angle[ i ] + k1AngleVelocity / 2.0, angleVelocity[ i ] + k1AngleVelocity / 2.0, springConstant, mass, gravity, initialLength ) )\n",
    "        k2AngleVelocity = Truncate( step * AngleAcceleration( length[ i ] + k1LengthVelocity / 2.0, angle[ i ] + k1AngleVelocity / 2.0, lengthVelocity[ i ] + k1LengthVelocity / 2.0, angleVelocity[ i ] + k1AngleVelocity / 2.0, gravity ) )\n",
    "        #print( '</k2>' )\n",
    "        #########################################\n",
    "        #print( '<k3>' )\n",
    "        k3Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] + k2Length / 2.0 ) )\n",
    "        k3Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] + k2AngleVelocity / 2.0 ) % 360.0 )\n",
    "        k3LengthVelocity = Truncate( step * LengthAcceleration( length[ i ] + k2LengthVelocity / 2.0, angle[ i ] + k2AngleVelocity / 2.0, angleVelocity[ i ] + k2AngleVelocity / 2.0, springConstant, mass, gravity, initialLength ) )\n",
    "        k3AngleVelocity = Truncate( step * AngleAcceleration( length[ i ] + k2LengthVelocity / 2.0, angle[ i ] + k2AngleVelocity / 2.0, lengthVelocity[ i ] + k2LengthVelocity / 2.0, angleVelocity[ i ] + k2AngleVelocity / 2.0, gravity ) )\n",
    "        #print( '</k3>' )\n",
    "        #########################################\n",
    "        #print( '<k4>' )\n",
    "        k4Length = Truncate( step * LengthVelocity( lengthVelocity[ i ] + k3Length ) )\n",
    "        k4Angle = Truncate( step * AngleVelocity( angleVelocity[ i ] + k3AngleVelocity ) % 360.0 )\n",
    "        k4LengthVelocity = Truncate( step * LengthAcceleration( length[ i ] + k3LengthVelocity, angle[ i ] + k3AngleVelocity, angleVelocity[ i ] + k1AngleVelocity, springConstant, mass, gravity, initialLength ) )\n",
    "        k4AngleVelocity = Truncate( step * AngleAcceleration( length[ i ] + k3LengthVelocity, angle[ i ] + k3AngleVelocity, lengthVelocity[ i ] + k3LengthVelocity, angleVelocity[ i ] + k3AngleVelocity, gravity ) )\n",
    "        #print( '</k4>' )\n",
    "        #########################################\n",
    "        #print( '<taylor>' )\n",
    "        length[ i + 1 ] = Truncate( ( length[ i ] + ( k1Length + 2.0 * ( k2Length + k3Length ) + k4Length ) / 6.0 ) )\n",
    "        angle[ i + 1 ] = Truncate( ( ( angle[ i ] + ( k1Angle + 2.0 * ( k2Angle + k3Angle ) + k4Angle ) / 6.0 ) % 360.0 ) )\n",
    "        lengthVelocity[ i + 1 ] = Truncate( lengthVelocity[ i ] + ( k1LengthVelocity + 2.0 * ( k2LengthVelocity + k3LengthVelocity ) + k4LengthVelocity ) / 6.0 )\n",
    "        angleVelocity[ i + 1 ] = Truncate( angleVelocity[ i ] + ( k1AngleVelocity + 2.0 * ( k2AngleVelocity + k3AngleVelocity ) + k4AngleVelocity ) / 6.0 )\n",
    "        #print( '</taylor>' )\n",
    "        #########################################\n",
    "        #print( '<xyt>' )\n",
    "        x[ i + 1 ] = Truncate( length[ i + 1 ] * np.sin( angle[ i + 1 ] ) )\n",
    "        y[ i + 1 ] = Truncate( 1.0 - length[ i + 1 ] * np.cos( angle[ i + 1 ] ) )\n",
    "        time[ i + 1 ] = step + time[ i ]\n",
    "        #print( '</xyt>' )\n",
    "    fig = plt.figure()\n",
    "    ax = fig.gca( projection = '3d' )\n",
    "    plt.plot( x, y, time )\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\ipykernel_launcher.py:13: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  del sys.path[0]\n",
      "C:\\Program Files (x86)\\Microsoft Visual Studio\\Shared\\Anaconda3_64\\lib\\site-packages\\ipykernel_launcher.py:13: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  del sys.path[0]\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "cannot convert float NaN to integer",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-30-9755af4b081c>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mSpring\u001b[0m\u001b[1;33m(\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1000\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-29-ed5a80d42c9a>\u001b[0m in \u001b[0;36mSpring\u001b[1;34m(step, timeEnd, mass, springConstant, initialLength, initialAngle, initialLengthVelocity, initialAngleVelocity)\u001b[0m\n\u001b[0;32m     36\u001b[0m         \u001b[0mk1Angle\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTruncate\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mstep\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mAngleVelocity\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mangleVelocity\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m \u001b[1;33m)\u001b[0m \u001b[1;33m%\u001b[0m \u001b[1;36m360.0\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     37\u001b[0m         \u001b[0mk1LengthVelocity\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTruncate\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mstep\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mLengthAcceleration\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mlength\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mangle\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mangleVelocity\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mspringConstant\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmass\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgravity\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minitialLength\u001b[0m \u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 38\u001b[1;33m         \u001b[0mk1AngleVelocity\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mTruncate\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mstep\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mAngleAcceleration\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mlength\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mangle\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlengthVelocity\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mangleVelocity\u001b[0m\u001b[1;33m[\u001b[0m \u001b[0mi\u001b[0m \u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgravity\u001b[0m \u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     39\u001b[0m         \u001b[1;31m#print( '</k1>' )\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     40\u001b[0m         \u001b[1;31m#########################################\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-26-f87bb772b748>\u001b[0m in \u001b[0;36mTruncate\u001b[1;34m(number, digits)\u001b[0m\n\u001b[0;32m     15\u001b[0m     \u001b[1;31m#    number = 1.0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m     \u001b[1;31m#print( number )\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 17\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfloat64\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mnumber\u001b[0m \u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: cannot convert float NaN to integer"
     ]
    }
   ],
   "source": [
    "Spring( 1, 1000 )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}