Nemesis_and_the_Sun.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "AMUSE tutorial on modules and channels\n",
    "====================\n",
    "\n",
    "AMUSE is composed of domain-specific modules (i.e. the physics modules), which are often written in some native compiled language, and data-handling.\n",
    "The latter tends to be confusing for novel users.\n",
    "Data representation in AMUSE is often replicated.\n",
    "So there can be a parameter *mass* to indicate a property of a star in a stellar-evolution code, in a gravitational dynamics code and in your user script. These three parameters *mass* can mean the same, or they can have a different meaning. This makes running amuse somewhat confusing at times.\n",
    "Here we show how to formally separate these data streams.\n",
    "\n",
    "But first we will enter the realm of modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load in the amuse units module, the particle module and \n",
    "# some generator for producing some conditions.\n",
    "from amuse.units import units, constants\n",
    "from amuse.lab import Particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now want to generate the Solar system, as it was observed on April 5th 2063 using the [JPL emphemeris](https://ssd.jpl.nasa.gov/horizons.cgi#top).\n",
    "Let's start with the inner most planet, Mercury, in the Sun's barycenter, which gives the following output:\n",
    "\n",
    "*******************************************************************************\n",
    "$$SOE\n",
    "2474649.500000000 = A.D. 2063-Apr-05 00:00:00.0000 TDB \n",
    " X =-2.689701945882702E-01 Y = 1.947622508089924E-01 Z = 4.148911797144442E-02\n",
    " VX=-2.230751806031045E-02 VY=-2.157655548570704E-02 VZ= 2.791759037949523E-04\n",
    " LT= 1.932846451835678E-03 RG= 3.346619889187973E-01 RR= 5.406475142917982E-03\n",
    " \n",
    "We can use these to start an AMUSE particle set, as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'Particles' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_26547/1076903443.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mSun\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mParticles\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mSun\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmass\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1\u001b[0m \u001b[0;34m|\u001b[0m \u001b[0munits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMSun\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mSun\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mposition\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m|\u001b[0m \u001b[0munits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkm\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mSun\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvelocity\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m|\u001b[0m \u001b[0munits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mkm\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0munits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0mplanets\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mParticles\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'Particles' is not defined"
     ]
    }
   ],
   "source": [
    "Sun = Particles(1)\n",
    "Sun.mass = 1 | units.MSun\n",
    "Sun.position = (0,0,0) | units.km\n",
    "Sun.velocity = (0,0,0) | units.km/units.s\n",
    "planets = Particles(8)\n",
    "planets[0].mass = 3.302e+23 | units.kg # also according to JPL)\n",
    "planets[0].position = (-2.689701945882702E-01,\n",
    "                       1.947622508089924E-01,\n",
    "                       4.148911797144442E-02) | units.au\n",
    "planets[0].velocity = (-2.230751806031045E-02,\n",
    "                       -2.157655548570704E-02,\n",
    "                       2.791759037949523E-04) | units.au/units.yr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And continue doing this for the other 7 planets.\n",
    "Then add the Sun and we have the Solar system's particle set.\n",
    "It would be easier if there is a handy routine with the same effect, in particular because we will frequently be using the Solar system as some sort of template for a rather typical planetary system, or for specifically sutdying this planetary system. We therefore have a handy routine that allows us to initialize the Solar system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "|S11\n",
      "                 key         mass         name       radius           vx           vy           vz            x            y            z\n",
      "                   -         MSun         none         RSun       km / s       km / s       km / s           AU           AU           AU\n",
      "====================  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========\n",
      " 6651284024525347701    1.000e+00          SUN    1.000e+00   -1.874e-03   -1.560e-02    1.104e-04   -8.354e-03    1.923e-03    2.145e-04\n",
      "11361453784849012877    1.660e-07      MERCURY    5.654e-03    1.032e+01   -4.214e+01   -4.388e+00   -3.923e-01   -1.749e-01    2.101e-02\n",
      " 7973865882662013015    2.448e-06        VENUS    1.403e-02   -1.704e+01    3.049e+01    1.401e+00    6.251e-01    3.518e-01   -3.157e-02\n",
      " 9336316391706255048    3.040e-06     EARTHMOO    1.482e-02    2.844e+01    6.966e+00    1.301e-04    2.337e-01   -9.855e-01    2.100e-04\n",
      "12320338445353835314    3.227e-07         MARS    7.853e-03   -2.299e+01    5.985e+00    6.908e-01    2.435e-01    1.528e+00    2.599e-02\n",
      "15770409518373525325    9.548e-04      JUPITER    1.619e-01    2.872e+00    1.331e+01   -1.194e-01    4.833e+00   -1.158e+00   -1.034e-01\n",
      " 5855752668938755582    2.859e-04       SATURN    1.342e-01   -4.794e+00    8.639e+00    4.001e-02    8.335e+00    4.127e+00   -4.033e-01\n",
      "14241492025581500998    4.366e-05       URANUS    5.834e-02    5.131e+00    4.103e+00   -5.124e-02    1.289e+01   -1.511e+01   -2.231e-01\n",
      " 7611664477873830688    5.151e-05      NEPTUNE    5.572e-02    4.640e+00    2.804e+00   -1.647e-01    1.537e+01   -2.592e+01    1.795e-01\n",
      " 1465633314586103062    7.396e-09        PLUTO    3.413e-03    5.144e+00   -2.956e+00   -1.172e+00   -1.152e+01   -2.708e+01    6.229e+00\n",
      "====================  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========\n"
     ]
    }
   ],
   "source": [
    "from amuse.ext.solarsystem import new_solar_system\n",
    "sun_and_planets = new_solar_system()\n",
    "#sun_and_planets = sun_and_planets[0].as_set()\n",
    "print(sun_and_planets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We adopted the module `ext.new_solar_system` to generate a ready-made solar system. You can also make a solar system that includes it's moons. For this you will have to use the routine \n",
    "'new_lunar_system', which is a module from the 'amuse.ic.solar_system_moons' package.\n",
    "\n",
    "Now you can plot the positions of the sun and planets in the orbital plane."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In 1984  \n",
    "[Whitemire and Jackson](https://ui.adsabs.harvard.edu/abs/1984Natur.308..713W/abstract) and [Davies et al 1984](\n",
    "https://ui.adsabs.harvard.edu/abs/1984Natur.308..715D)([preprint](https://escholarship.org/uc/item/9gm5c682)) argued that the Solar system was orbited by a hypothetical star called nemesis, named after the goddess\n",
    "of Rhamnous who punished those who did not honor the Olympians.\n",
    "\n",
    "This companion to the Sun would have an orbital separation of about 95,000 au, and due to its interaction with the Oort cloud would initiate periodic comet shouwers in the inner Solar system; punishing life on Earth for being atheistic?\n",
    "\n",
    "Let's adopt the Solar system as we have just initialized and add Nemesis as an 0.2MSun main-sequence star in a 95,000 au orbit with an eccentricity of 0.7.\n",
    "\n",
    "Let's construct tha binary system of the Sun (without planets) and Nemesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 key         mass         name           vx           vy           vz            x            y            z\n",
      "                   -         MSun         none  3646245880.3572216 * m * s**-1  3646245880.3572216 * m * s**-1  3646245880.3572216 * m * s**-1  149597870691.0 * m  149597870691.0 * m  149597870691.0 * m\n",
      "====================  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========\n",
      "  388603917714193344    1.001e+00    Companion    0.000e+00   -1.151e-08   -0.000e+00   -4.745e+03   -0.000e+00   -0.000e+00\n",
      " 1364133191765162363    2.000e-01      Nemesis    0.000e+00    5.764e-08    0.000e+00    2.376e+04    0.000e+00    0.000e+00\n",
      "====================  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========\n"
     ]
    }
   ],
   "source": [
    "msun = sun_and_planets.mass.sum()\n",
    "mnemesis = 0.2 | units.MSun\n",
    "a = 95000 | units.au\n",
    "e = 0.7\n",
    "from amuse.ext.orbital_elements import new_binary_from_orbital_elements\n",
    "sun_and_nemesis = new_binary_from_orbital_elements(msun, mnemesis, \n",
    "                                          a, e, G=constants.G)\n",
    "setattr(sun_and_nemesis, \"name\", [\"Companion\", \"Nemesis\"])\n",
    "print(sun_and_nemesis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the sun in this binary-particle set and the sun in the Solar system particle-set do not have the same position."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.00835408821766, 0.0019229394648, 0.000214547107396]] au\n",
      "[[-4744.69448798, -0.0, -0.0]] au\n"
     ]
    }
   ],
   "source": [
    "sun = sun_and_planets[sun_and_planets.name==\"SUN\"]\n",
    "print(sun.position.in_(units.au))\n",
    "print(sun_and_nemesis[sun_and_nemesis.mass>0.9|units.MSun].position.in_(units.au))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And, equally they will have different velocity.\n",
    "\n",
    "However, we would like to adot these binary parameters for the Sun-Nemesis binary system. For this, we first move the entire Solar system to the Sun barycenter, move it to Nemesis' companion and add Nemesis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 key         mass         name       radius           vx           vy           vz            x            y            z\n",
      "                   -         MSun         none         RSun       km / s       km / s       km / s           AU           AU           AU\n",
      "====================  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========\n",
      " 6651284024525347701    1.000e+00          SUN    1.000e+00    0.000e+00   -4.198e-02    0.000e+00   -4.745e+03    0.000e+00    0.000e+00\n",
      "11361453784849012877    1.660e-07      MERCURY    5.654e-03    1.032e+01   -4.217e+01   -4.389e+00   -4.745e+03   -1.769e-01    2.080e-02\n",
      " 7973865882662013015    2.448e-06        VENUS    1.403e-02   -1.704e+01    3.046e+01    1.400e+00   -4.744e+03    3.499e-01   -3.179e-02\n",
      " 9336316391706255048    3.040e-06     EARTHMOO    1.482e-02    2.845e+01    6.939e+00    1.967e-05   -4.744e+03   -9.875e-01   -4.543e-06\n",
      "12320338445353835314    3.227e-07         MARS    7.853e-03   -2.298e+01    5.959e+00    6.907e-01   -4.744e+03    1.526e+00    2.578e-02\n",
      "15770409518373525325    9.548e-04      JUPITER    1.619e-01    2.874e+00    1.329e+01   -1.196e-01   -4.740e+03   -1.160e+00   -1.036e-01\n",
      " 5855752668938755582    2.859e-04       SATURN    1.342e-01   -4.792e+00    8.613e+00    3.990e-02   -4.736e+03    4.125e+00   -4.035e-01\n",
      "14241492025581500998    4.366e-05       URANUS    5.834e-02    5.133e+00    4.076e+00   -5.135e-02   -4.732e+03   -1.511e+01   -2.233e-01\n",
      " 7611664477873830688    5.151e-05      NEPTUNE    5.572e-02    4.641e+00    2.777e+00   -1.648e-01   -4.729e+03   -2.592e+01    1.793e-01\n",
      " 1465633314586103062    7.396e-09        PLUTO    3.413e-03    5.146e+00   -2.982e+00   -1.172e+00   -4.756e+03   -2.708e+01    6.229e+00\n",
      " 1364133191765162363    2.000e-01      Nemesis    0.000e+00    0.000e+00    2.102e-01    0.000e+00    2.376e+04    0.000e+00    0.000e+00\n",
      "====================  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========  ===========\n"
     ]
    }
   ],
   "source": [
    "sun_and_planets.position -= sun.position\n",
    "sun_and_planets.velocity -= sun.velocity\n",
    "companion = sun_and_nemesis[sun_and_nemesis.name==\"Companion\"]\n",
    "nemesis = sun_and_nemesis[sun_and_nemesis.name==\"Nemesis\"]\n",
    "sun_and_planets.position += companion.position\n",
    "sun_and_planets.velocity += companion.velocity\n",
    "sun_and_planets.add_particle(sun_and_nemesis[1])\n",
    "print(sun_and_planets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have created a hypothetical planetary system with the Sun, 8 planets and 1 dwarf planet (mimiking the Solar system), and a secondary star (Nemesis) in a wide and elliptic orbit.\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Assignment 1:\n",
    "To check the orbit of Nemesis around the Sun you can calculate the orbital elements from the Cartesian coordinates of the Sun and Nemesis. Perform this operation and check the orbital elements of the nemesis-Sun binary. \n",
    "\n",
    "### Question 1:\n",
    "When performing this assignment, we found that the orbit between the Sun and Nemesis is slightly off from what we put in. The semi-major axis is about 95,606 au with an eccentricity of about 0.702. \n",
    "\n",
    "Explain why the orbit is not what you anticipated it to be, and what can you do to mitigate this?\n",
    "\n",
    "### Assignment 2: \n",
    "Add 100 zero-mass Oort cloud object in circular orbits around the Sun with semi-major axes equally spaced between 10000au and 50000 au. Integrate those orbits for 10 full orbits of the Sun-Nemesis binary. Make a scatter plot of eccentricity as function of the semi-major axis for the initial and final Oort-cloud objects. Explain the figure.\n",
    "\n",
    "Be aware that 10 Nemesis periods is ~160 Myr, and running the system for that time will take a lot of time (about an hour). Start with shorter times and build up towards longer times when you know your code runs. You can also try different codes, and leaving out some planets. Is Mercury really that important? But leaving it out may let the code take greater timesteps.\n",
    "\n",
    "### Question 3:\n",
    "Do you think that Nemesis could be the cause of periodic cometary visitors?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}