{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integration Methods\n",
    "\n",
    "In this notebook we will investigate different methods of numerical integration and see their stability. The Particle class by default uses the velocity Verlet method, a symplectic method (meaning it nearly conserves energy) with error $O(\\Delta t^2)$. The algorithm is as follows:\n",
    "$$$$\n",
    "$$ x_{n+1} = x_n + v_n \\Delta t + \\frac{1}{2} a_n \\Delta t^2 $$\n",
    "$$$$\n",
    "$$ v_{n+1} = v_n + \\frac{a_n+ a_{n+1}}{2} \\Delta t $$\n",
    "$$$$\n",
    "We will subclass the Particle class and create an EulerParticle class which implements the Euler method, and see how it compares with the velocity Verlet method. First, we import the necessary libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/jquery-ui.custom.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/glow.2.1.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/glowcomm\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/pako.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/pako_deflate.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require.undef(\"nbextensions/pako_inflate.min\");"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "require([\"nbextensions/glowcomm\"], function(){console.log(\"glowcomm loaded\");})"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#NAME: Integration Methods\n",
    "#DESCRIPTION: Investigation of different numerical integration techniques and their stability in dynamics simulations.\n",
    "\n",
    "from pycav.mechanics import *\n",
    "from vpython import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will create the EulerParticle class, which implements the Euler method:\n",
    "$$$$\n",
    "$$ x_{n+1} = x_n + v_n \\Delta t $$\n",
    "$$$$\n",
    "$$ v_{n+1} = v_n + a \\Delta t $$\n",
    "$$$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class EulerParticle(Particle):\n",
    "    def update(self, dt):\n",
    "        self.pos += (self.v * dt)\n",
    "        self.v += (dt * self.total_force) * self.inv_mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also subclass the System class to create a MeasuredSystem class which will be able to tell us interesting properties of a system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class MeasuredSystem(System):\n",
    "    @property\n",
    "    def kinetic_energy(self):\n",
    "        _ke = 0.\n",
    "        for particle in self.particles:\n",
    "            _ke += (0.5 * np.inner(particle.v, particle.v))/(particle.inv_mass)\n",
    "        return _ke\n",
    "    \n",
    "    @property\n",
    "    def gravitational_potential_energy(self):\n",
    "        _gpe = 0.\n",
    "        for particle1 in self.particles:\n",
    "            for particle2 in self.particles:\n",
    "                dist = particle1.pos - particle2.pos\n",
    "                if particle1 is not particle2:\n",
    "                    r = np.sqrt(np.inner(dist, dist))\n",
    "                    _gpe = -1. / (r * particle1.inv_mass * particle2.inv_mass)\n",
    "        _gpe /= 2.\n",
    "        return _gpe\n",
    "    \n",
    "    @property\n",
    "    def total_energy(self):\n",
    "        return (self.gravitational_potential_energy + self.kinetic_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will then set up the same system using normal Particles and EulerParticles. The EulerPacticle system will be colored white so that we can see the difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "planet_verlet = Particle(pos = np.array([200.,0.,0.]),\n",
    "                        v = np.array([0., 31.62, 0.]),\n",
    "                        inv_mass = 1.,\n",
    "                        make_trail = True,\n",
    "                        radius = 10.)\n",
    "star_verlet = Particle(pos = np.array([0., 0., 0.]),\n",
    "                      v = np.array([0., 0., 0.]),\n",
    "                      inv_mass = 1./200000.,\n",
    "                      radius = 20.)\n",
    "\n",
    "planet_euler = EulerParticle(pos = np.array([200.,0.,0.]),\n",
    "                        v = np.array([0., 31.62, 0.]),\n",
    "                        inv_mass = 1.,\n",
    "                        make_trail = True,\n",
    "                        color = [1., 1., 1.],\n",
    "                        radius = 10.)\n",
    "star_euler = EulerParticle(pos = np.array([0., 0., 0.]),\n",
    "                      v = np.array([0., 0., 0.]),\n",
    "                      inv_mass = 1./200000.,\n",
    "                      radius = 20.)\n",
    "\n",
    "\n",
    "verlet = [planet_verlet, star_verlet]\n",
    "euler = [planet_euler, star_euler]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now simulate these two systems side by side, on the same canvas. We will also plot the number of steps taken and the Star-Planet distance for each system. You should see that for any time step, the Verlet system stays stable while the Euler method has the orbits spiralling out. Try looking at other properties of the system (total energy, kinetic energy, etc.) and see how they are preserved by these methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-5-b3f4078dfd5d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m     \u001b[0mrate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m150\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     14\u001b[0m     \u001b[0mverlet_sys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msimulate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m     \u001b[0meuler_sys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msimulate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Applications/anaconda/lib/python3.5/site-packages/vpython/vpython.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, maxRate)\u001b[0m\n\u001b[1;32m    132\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrval\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mmaxRate\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mmaxRate\u001b[0m \u001b[0;34m>=\u001b[0m \u001b[0;36m1.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    133\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmaxRate\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 134\u001b[0;31m         \u001b[0msuper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mRateKeeper2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmaxRate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    135\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    136\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0msys\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mversion\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;34m'3'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Applications/anaconda/lib/python3.5/site-packages/vpython/rate_control.py\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, maxRate)\u001b[0m\n\u001b[1;32m    176\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcalls\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0;31m# if there were some calls to rate after the last render\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    177\u001b[0m                 \u001b[0mdt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlastSleep\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcalls\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0muserTime\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcallTime\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdelay\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0m_clock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 178\u001b[0;31m                 \u001b[0m_sleep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    179\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuildStrategy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmaxRate\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    180\u001b[0m             \u001b[0mnr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwhenToRender\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Applications/anaconda/lib/python3.5/site-packages/vpython/rate_control.py\u001b[0m in \u001b[0;36m_sleep\u001b[0;34m(dt)\u001b[0m\n\u001b[1;32m     47\u001b[0m         \u001b[0mdtsleep\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnticks\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0m_tick\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     48\u001b[0m         \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_clock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 49\u001b[0;31m         \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msleep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdtsleep\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     50\u001b[0m         \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_clock\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     51\u001b[0m         \u001b[0mdt\u001b[0m \u001b[0;34m-=\u001b[0m \u001b[0mt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "scene1 =  canvas(title='Comparison of integrating methods')\n",
    "verlet_sys = MeasuredSystem(collides = False, interacts = True, visualize = True, particles = verlet, canvas = scene1)\n",
    "euler_sys = MeasuredSystem(collides = False, interacts = True, visualize = True, particles = euler, canvas = scene1)\n",
    "\n",
    "graph1 = graph(x=0, y=0, \n",
    "      xtitle='Steps', ytitle='Star-Planet distance', \n",
    "      foreground=color.black, background=color.black)\n",
    "f1 = gcurve(color=color.white)\n",
    "f2 = gcurve(color=color.red)\n",
    "dt = 0.1\n",
    "\n",
    "while True:\n",
    "    rate(150)\n",
    "    verlet_sys.simulate(dt)\n",
    "    euler_sys.simulate(dt)\n",
    "    dis_verlet = planet_verlet.pos - star_verlet.pos\n",
    "    dis_euler = planet_euler.pos - star_euler.pos\n",
    "    f2.plot(verlet_sys.steps, np.sqrt(np.inner(dis_verlet, dis_verlet)))\n",
    "    f1.plot(verlet_sys.steps, np.sqrt(np.inner(dis_euler, dis_euler)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}