{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving many body problem using Vpython\n",
    "We will solve the motion of three or more particles interacting through gravitational interaction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Accelaration due to Gravitational force"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vpython import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#define the value of G\n",
    "G=6.673e-11\n",
    "\n",
    "## Acceleration of object a due to object b because of their gravitational interaction\n",
    "def acc(a, b):\n",
    "    rel_pos = b.pos - a.pos\n",
    "    return G*b.mass * norm(rel_pos)/rel_pos.mag2\n",
    "\n",
    "## Accelaration of a due to all the objects b intracting with it\n",
    "def totalacc (a, objlist):\n",
    "    sum_acc = vector (0,0,0)\n",
    "    for b in objlist:\n",
    "        if (a!=b):\n",
    "            sum_acc = sum_acc + acc(a, b)\n",
    "    return sum_acc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Constants defining the solar system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solar system on a computer\n",
    "\n",
    "#Constants that we will need\n",
    "#define the value of G\n",
    "G=6.673e-11\n",
    "myPi = 3.141592\n",
    "sun_mass = 2e30\n",
    "earth_mass = 6e24\n",
    "\n",
    "boost=1.0 #allow for boosting Jupiter's mass\n",
    "\n",
    "jupiter_mass=boost*1.9e27\n",
    "#for initial conditions\n",
    "\n",
    "AU = 149.6e9       #mean earth sun orbital distance\n",
    "\n",
    "earth_vel = 2* myPi * AU/(365.25 *24. *60.*60.) # average velocity = 2*Pi*R/T\n",
    "jupiter_vel=2* myPi *AU*5.2/(11.86*365.25*24.*60.*60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up the scene for animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#setting for animations\n",
    "scene.background = color.white\n",
    "scene.autoscale = 0\n",
    "scene.range = 10*AU"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Making the objects that will form our solar system "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#objects making up our solar system\n",
    "sun = sphere(pos= vector(0,0,0), velocity = vector(0,0,0),\n",
    "             mass=sun_mass, radius = 0.1*AU, color =color.yellow)\n",
    "earth = sphere(pos= vector(AU, 0, 0), velocity = vector(0,earth_vel,0),\n",
    "               mass=earth_mass, radius=0.05*AU, color =color.cyan)\n",
    "jupiter=sphere(pos=vector(5.2*AU,0,0),velocity=vector(0,jupiter_vel,0),\n",
    "             mass=jupiter_mass, radius=0.15*AU, color=color.red)\n",
    "\n",
    "#note the radius of sun, jupiter  and earth are NOT\n",
    "# their true radius these are the radius of the spherical object\n",
    "#that will be draw on the computer scree"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding a new property of acceleration for our objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Create a list of objects making up our solar system \n",
    "#and add attributes for their accelaration and orbits\n",
    "\n",
    "bodies = [sun, earth, jupiter]\n",
    "\n",
    "for b in bodies:\n",
    "    b.acc = vector(0,0,0)\n",
    "    b.track=curve (color = b.color)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Going to the centre of mass frame\n",
    "In the centre of mass fram the total momentum of the system is zero"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set total momentum of system to zero (centre of mass frame) \n",
    "sum=vector(0,0,0)\n",
    "for b in bodies:\n",
    "    if (b!=sun):\n",
    "        sum=sum+b.mass*b.velocity\n",
    "\n",
    "sun.velocity=-sum/sun.mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dt corresponds to 3000 mins here\n",
    "dt=30.*60.*100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the equations of motion using leap frog method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initialize leap-frog by finding the velocites at t=dt/2\n",
    "\n",
    "for b in bodies:\n",
    "    b.velocity = b.velocity + totalacc(b, bodies)*dt/2.0\n",
    "\n",
    "#start leap-frog\n",
    "while True:\n",
    "    rate(50)  #not more than 100 time steps in a second\n",
    "    for b in bodies:\n",
    "        #update the positions\n",
    "        b.pos = b.pos + b.velocity*dt\n",
    "        b.track.append(pos=b.pos)\n",
    "\n",
    "        #update the velocities\n",
    "        b.velocity = b.velocity + totalacc(b, bodies)*dt\n",
    "\n",
    "    scene.center = vector(0,0,0) #view centered on the origin of CM coord system\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "VPython",
   "language": "python",
   "name": "vpython"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}