Solutions_3_IntroToSimulations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solutions: Making Orbital Dreams Come True\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Try it yourself: \n",
    "<font color='maroon'> Copy everything in the cell above. This time, set the eccentricity to 0, and all the other values to zero as well. Run it to make sure it works. Then, add a *SECOND* planet, and change just one of the values to see how the orbit changes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import rebound\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# so we can make plots\n",
    "%matplotlib inline \n",
    "%run tools.ipynb \n",
    "\n",
    "\n",
    "def start_new_sim(star_mass, star_name): # We have to give our def a unique name, and we can pass parameters to it\n",
    "    sim = rebound.Simulation()\n",
    "    sim.G = 4.*pi**2.\n",
    "    sim.units = ('yr', 'AU', 'Msun')\n",
    "    sim.add(m=star_mass, hash=star_name)  \n",
    "    return sim     # We always need to return from our definition; here, we want to pass back our simulation\n",
    "\n",
    "t_final = 5 # Length of simulation in years\n",
    "N_frames = 50 # Number of frames in the animation\n",
    "\n",
    "print('Changing the eccentricity:')\n",
    "sim = start_new_sim(1.,'Sun') \n",
    "sim.add(m=3e-6, a=1., e=0, hash='Planet')\n",
    "sim.add(m=3e-6, a=2., e=0, hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e=0.8, hash='Planet')\n",
    "sim.move_to_com()\n",
    "integrate_with_widget(sim, t_final, N_frames, 3)\n",
    "\n",
    "print('Changing the eccentricity even more!:')\n",
    "sim = start_new_sim(1.,'Sun') \n",
    "sim.add(m = 3e-6, a = 1, e = 0.0, hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.1,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.2,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.3,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.4,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.5,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.6,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.7,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.8,  hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.9,  hash='Planet')\n",
    "sim.move_to_com()\n",
    "integrate_with_widget(sim, t_final, N_frames, 3)\n",
    "\n",
    "print('Changing the eccentricity and adding inclination! Move this around a bunch:')\n",
    "sim = start_new_sim(1.,'Sun') \n",
    "sim.add(m = 3e-6, a = 1, e = 0.0, inc = radians(0), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.1, inc = radians(10), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.2, inc = radians(20), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.3, inc = radians(30), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.4, inc = radians(40), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.5, inc = radians(50), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.6, inc = radians(60), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.7, inc = radians(70), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.8, inc = radians(80), hash='Planet')\n",
    "sim.add(m = 3e-6, a = 1, e = 0.9, inc = radians(90), hash='Planet')\n",
    "sim.move_to_com()\n",
    "integrate_with_widget(sim, t_final, N_frames, 3)\n",
    "\n",
    "\n",
    "print('Changing the inclination:')\n",
    "sim = start_new_sim(1.,'Sun') \n",
    "sim.add(m=3e-6, a=1., hash='Planet')\n",
    "sim.add(m=3e-6, a=2., hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, inc=radians(20), hash='Planet')\n",
    "sim.move_to_com()\n",
    "integrate_with_widget(sim, t_final, N_frames, 3)\n",
    "\n",
    "\n",
    "print('Changing the argument of periapsis (rotates the ellipse):')\n",
    "sim = start_new_sim(1.,'Sun')\n",
    "sim.add(m=3e-6, a=1., hash='Planet')\n",
    "sim.add(m=3e-6, a=2., hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), omega=radians(45.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), omega=radians(25.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), omega=radians(35.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), omega=radians(15.), hash='Planet')\n",
    "sim.move_to_com()\n",
    "integrate_with_widget(sim, t_final, N_frames, 3)\n",
    "\n",
    "\n",
    "print('Changing the argument of the ascending node:')\n",
    "sim = start_new_sim(1.,'Sun') \n",
    "sim.add(m=3e-6, a=1., hash='Planet')\n",
    "sim.add(m=3e-6, a=2., hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), Omega=radians(65.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), Omega=radians(45.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), Omega=radians(25.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.2, inc=radians(20), Omega=radians(5.), hash='Planet')\n",
    "sim.move_to_com()\n",
    "integrate_with_widget(sim, t_final, N_frames, 3)\n",
    "\n",
    "\n",
    "print('Changing the true anomaly:')\n",
    "sim = start_new_sim(1.,'Sun') \n",
    "sim.add(m=3e-6, a=1., hash='Planet')\n",
    "sim.add(m=3e-6, a=2., hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.5, inc=radians(20), f=radians(270.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.5, inc=radians(20), f=radians(180.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.5, inc=radians(20), f=radians(90.), hash='Planet')\n",
    "sim.add(m=3e-6, a=2.5, e = 0.5, inc=radians(20), f=radians(0.), hash='Planet')\n",
    "sim.move_to_com()\n",
    "integrate_with_widget(sim, t_final, N_frames, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='maroon'>Challenge: </font>\n",
    "Try to make a simulation of the solar system! Only include the planets: Jupiter, Saturn, Uranus, Neptune, and the dwarf planet Pluto! Include the correct masses (already defined above), the correct distances (already defined above), the correct eccentricities (again, use data from the link provided above), and the correct names!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# code here: \n",
    "import rebound\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "# so we can make plots\n",
    "%matplotlib inline \n",
    "%run tools.ipynb # Some behind-the-scenes tools to simplify some complicated code\n",
    "\n",
    "# masses are all going to be in kg. In order to get masses in solar masses, just divide each one by the mass of the sun!\n",
    "\n",
    "mass_Sun     = 1.989e30\n",
    "\n",
    "mass_Mercury = 3.285e23\n",
    "mass_Venus   = 4.867e24 # similar to the mass of the earth!\n",
    "mass_Earth   = 5.972e24\n",
    "mass_Mars    = 6.390e23\n",
    "mass_Jupiter = 1.898e27\n",
    "mass_Saturn  = 5.683e26\n",
    "mass_Uranus  = 8.681e25\n",
    "mass_Neptune = 1.024e26\n",
    "mass_Pluto   = 1.309e22 # not a planet, but still fun to plot!\n",
    "\n",
    "def start_new_sim(star_mass, star_name): # We have to give our def a unique name, and we can pass parameters to it\n",
    "    sim = rebound.Simulation()\n",
    "    sim.G = 4.*pi**2.\n",
    "    sim.units = ('yr', 'AU', 'Msun')\n",
    "    sim.add(m=star_mass, hash=star_name)  # Here, we use the variables we passed to the definition to initialize our star\n",
    "    return sim # We always need to return from our definition; here, we want to pass back our simulation\n",
    "\n",
    "\n",
    "# code here: \n",
    "# code here! \n",
    "t_final = 50 # Length of simulation in years\n",
    "N_frames = 200 # Number of frames in the animation\n",
    "\n",
    "# Start new simulation - you only get 1 widget per simulation, so if you want a new widget, you need a new simulation\n",
    "sim = start_new_sim(1.,'Sun') # We call the definition, and we must pass the required elements\n",
    "sim.add(m = mass_Jupiter/mass_Sun, a = 5.205, e = 0.0485, inc = radians(1.30), hash = 'Jupiter')\n",
    "sim.add(m = mass_Saturn/mass_Sun, a = 9.582, e = 0.0556, inc = radians(2.49), f=np.pi*6/4, hash = 'Saturn')\n",
    "sim.add(m = mass_Uranus/mass_Sun, a = 19.20, e = 0.046, inc = radians(0.77), f=np.pi/4, hash = 'Uranus')\n",
    "sim.add(m = mass_Neptune/mass_Sun, a = 30.05, e = 0.011, inc = radians(1.77), f=np.pi, hash = 'Neptune')\n",
    "sim.add(m = mass_Pluto/mass_Sun, a = 39.48, e = 0.244, inc = radians(17.2), f=np.pi/2, hash = 'Pluto')\n",
    "sim.move_to_com()\n",
    "\n",
    "# And now integrate!\n",
    "integrate_with_widget(sim, t_final, N_frames, 26)\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<br><br><br><br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is Jupiter's Orbital Period? \n",
    "\n",
    "<font color='maroon'>Since the only information that is used in constructing this simulation is info about the orbit, the amount of time it takes for Jupiter to orbit the sun is only related to it's orbital parameters! In particular, the orbital period is closely-related to the distance from the sun. Can you use the animation below to determine Jupiter's orbital period in terms of Earth-years (normal years, to us).</font>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture \n",
    "# This cell produces a lot of unnecessary output, so we'll capture that output instead of printing it out.\n",
    "# The next code cell will actually display the animation.\n",
    "\n",
    "dt = 0.1 # Let's integrate a tenth of a year each step\n",
    "N_steps = 118 # Let's integrate it forward for 20 steps, so 20*0.1 = 2 years\n",
    "\n",
    "# Start a fresh simulation\n",
    "sim = start_new_sim(1.,'Sun') \n",
    "sim.add(m=3.E-6, a=1, e=0.0167,  hash='Earth')\n",
    "sim.add(m = mass_Jupiter/mass_Sun, a = 5.205, e = 0.0485, inc = radians(1.30), hash = 'Jupiter')\n",
    "sim.move_to_com()\n",
    "\n",
    "# We need to create a figure (fig) and axes (ax) for the animation\n",
    "fig,ax = plt.subplots(1,1,figsize=(5,5))\n",
    "\n",
    "# The animation is created from a list of frames\n",
    "frames = []\n",
    "\n",
    "for i in range(N_steps):\n",
    "    \n",
    "    new_time = sim.t+dt # The time we want to integrate to\n",
    "    sim.integrate(new_time) # And integrate!\n",
    "    \n",
    "    frames.append(make_rebound_frame(ax,sim)) # Add the new frame to the list for animation\n",
    "       \n",
    "# Create the animation\n",
    "ani = animation.ArtistAnimation(fig, frames)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ani"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The orbital period of Jupiter is ~11.8 Earth years!"
   ]
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}