09 - Gravity_cascaded_bridges.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "AMUSE tutorial on multiple code in a single bridge\n",
    "====================\n",
    "\n",
    "A cascade of bridged codes to address the problem of running multiple planetary systems in, for example, a star cluster. This is just an example of how to initialize such a cascaded bridge without any stellar evolution, background potentials. The forces for one planetary system on the planets in the other systems are ignored to save computer time. This gives rise to some energy errors, and inconsistencies (for example when one star tries to capture planets from another system. The latter will not happen here.\n",
    "This can be addressed by intorducing some logic in checking what stars are nearby which planets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from amuse.units import (units, constants)\n",
    "from amuse.lab import Particles\n",
    "from amuse.units import nbody_system\n",
    "from matplotlib import pyplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## source https://en.wikipedia.org/wiki/TRAPPIST-1\n",
    "trappist= {\"b\": {\"m\": 1.374 | units.MEarth,\n",
    "                 \"a\": 0.01154 | units.au,\n",
    "                 \"e\": 0.00622,\n",
    "                 \"i\": 89.56},\n",
    "           \"c\": {\"m\": 1.308 | units.MEarth,\n",
    "                 \"a\": 0.01580 | units.au,\n",
    "                 \"e\": 0.00654,\n",
    "                 \"i\": 89.70},\n",
    "           \"d\": {\"m\": 0.388 | units.MEarth,\n",
    "                 \"a\": 0.02227 | units.au,\n",
    "                 \"e\": 0.00837,\n",
    "                 \"i\": 89.89},\n",
    "           \"e\": {\"m\": 0.692 | units.MEarth,\n",
    "                 \"a\": 0.02925 | units.au,\n",
    "                 \"e\": 0.00510,\n",
    "                 \"i\": 89.736},\n",
    "           \"f\": {\"m\": 1.039 | units.MEarth,\n",
    "                 \"a\": 0.03849 | units.au,\n",
    "                 \"e\": 0.01007,\n",
    "                 \"i\": 89.719},\n",
    "           \"g\": {\"m\": 1.321 | units.MEarth,\n",
    "                 \"a\": 0.04683 | units.au,\n",
    "                 \"e\": 0.00208,\n",
    "                 \"i\": 89.721},\n",
    "           \"h\": {\"m\": 0.326 | units.MEarth,\n",
    "                 \"a\": 0.06189 | units.au,\n",
    "                 \"e\": 0.00567,\n",
    "                 \"i\": 89.796}\n",
    "           }\n",
    "\n",
    "def trappist_system():\n",
    "    from amuse.ext.orbital_elements import new_binary_from_orbital_elements\n",
    "    from numpy.random import uniform\n",
    "    star = Particles(1)\n",
    "    setattr(star, \"name\", \"\")\n",
    "    setattr(star, \"type\", \"\")\n",
    "    star[0].mass = 0.898 | units.MSun\n",
    "    star[0].position = (0,0,0) | units.au\n",
    "    star[0].velocity = (0,0,0) | units.kms\n",
    "    star[0].name = \"trappist\"\n",
    "    star[0].type = \"star\"\n",
    "    bodies = Particles(len(trappist))\n",
    "    setattr(bodies, \"name\", \"\")\n",
    "    setattr(bodies, \"type\", \"\")\n",
    "\n",
    "    for bi, planet in zip(bodies, trappist):\n",
    "        true_anomaly = uniform(0, 360)\n",
    "        b = new_binary_from_orbital_elements(star.mass, \n",
    "                                             trappist[planet]['m'], \n",
    "                                             trappist[planet][\"a\"],\n",
    "                                             trappist[planet][\"e\"],\n",
    "                                             true_anomaly = true_anomaly,\n",
    "                                             inclination = trappist[planet][\"i\"],\n",
    "                                             G = constants.G)\n",
    "        bi.name = planet\n",
    "        bi.type = \"planet\"\n",
    "        bi.mass = b[1].mass\n",
    "        bi.position = b[1].position - b[0].position\n",
    "        bi.velocity = b[1].velocity - b[0].velocity\n",
    "    return star | bodies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.ph4 import Ph4\n",
    "from amuse.community.hermite import Hermite\n",
    "\n",
    "from amuse.ic.plummer import new_plummer_model\n",
    "import numpy.random\n",
    "numpy.random.seed(1624973942)\n",
    "\n",
    "converter = nbody_system.nbody_to_si(1 | units.MSun, 0.1 | units.parsec)\n",
    "\n",
    "t1 = trappist_system()\n",
    "t2 = trappist_system()\n",
    "t3 = trappist_system()\n",
    "p = new_plummer_model(3, convert_nbody=converter)\n",
    "t1.position += p[0].position\n",
    "t1.velocity += p[0].velocity\n",
    "t2.position += p[1].position\n",
    "t2.velocity += p[1].velocity\n",
    "t3.position += p[2].position\n",
    "t3.velocity += p[2].velocity\n",
    "\n",
    "converter=nbody_system.nbody_to_si(t1.mass.sum(), 0.1 | units.au)\n",
    "\n",
    "bodies = Particles(0)\n",
    "gravity1 = Ph4(converter)\n",
    "t = gravity1.particles.add_particles(t1)\n",
    "bodies.add_particles(t1)\n",
    "\n",
    "gravity2 = Ph4(converter)\n",
    "t = gravity2.particles.add_particles(t2)\n",
    "bodies.add_particles(t2)\n",
    "\n",
    "gravity3 = Hermite(converter)\n",
    "t = gravity3.particles.add_particles(t3)\n",
    "bodies.add_particles(t3)\n",
    "\n",
    "channel_from_g1 = gravity1.particles.new_channel_to(bodies)\n",
    "channel_from_g2 = gravity2.particles.new_channel_to(bodies)\n",
    "channel_from_g3 = gravity3.particles.new_channel_to(bodies)\n",
    "\n",
    "from amuse.plot import scatter\n",
    "from matplotlib import pyplot\n",
    "scatter(bodies.x-bodies[0].x, bodies.z-bodies[0].z)\n",
    "#pyplot.xlim(-0.1, 0.1)\n",
    "#pyplot.ylim(-0.1, 0.1)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot(bodies):\n",
    "    from amuse.plot import scatter\n",
    "    from matplotlib import pyplot\n",
    "    stars = bodies[bodies.type=='star']\n",
    "    planets = bodies-stars\n",
    "    pyplot.scatter((stars.x-bodies[0].x).value_in(units.au), \n",
    "                   (stars.z-bodies[0].z).value_in(units.au), c='r', s=100)\n",
    "    pyplot.scatter((planets.x-bodies[0].x).value_in(units.au), \n",
    "                   (planets.z-bodies[0].z).value_in(units.au), c='b', s=10)\n",
    "    pyplot.xlim(-0.1, 0.1)\n",
    "    pyplot.ylim(-0.1, 0.1)\n",
    "    pyplot.show()\n",
    "plot(bodies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.couple import bridge\n",
    "gravity = bridge.Bridge()\n",
    "gravity.add_system(gravity1, (gravity2,gravity3))\n",
    "gravity.add_system(gravity2, (gravity1,gravity3))\n",
    "gravity.add_system(gravity3, (gravity1,gravity2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.lab import zero\n",
    "Etot_init = gravity.kinetic_energy + gravity.potential_energy\n",
    "Etot_prev = Etot_init\n",
    "\n",
    "gravity.timestep = 100.0 | units.yr\n",
    "time = zero\n",
    "dt = 200.0 | units.yr\n",
    "t_end = 1000.0 | units.yr\n",
    "\n",
    "while time < t_end:\n",
    "    time += dt\n",
    "    gravity.evolve_model(time)\n",
    "\n",
    "    Etot_prev_se = gravity.kinetic_energy + gravity.potential_energy\n",
    "\n",
    "    channel_from_g1.copy()\n",
    "    channel_from_g2.copy()\n",
    "    channel_from_g3.copy()\n",
    "    plot(bodies)\n",
    "    print(bodies[1].position.in_(units.au))\n",
    "             \n",
    "    Ekin = gravity.kinetic_energy \n",
    "    Epot = gravity.potential_energy\n",
    "    Etot = Ekin + Epot\n",
    "    print(\"T=\", time.in_(units.yr), end=' ') \n",
    "    print(\"E= \", Etot/Etot_init, \"Q= \", Ekin/Epot, end=' ')\n",
    "    print(\"dE=\", (Etot_init-Etot)/Etot, \"ddE=\", (Etot_prev-Etot)/Etot) \n",
    "    Etot_prev = Etot\n",
    "gravity.stop()"
   ]
  },
  {
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}