bridged_Nbody_with_Galaxy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run a star cluster orbiting the Galaxy\n",
    "====================\n",
    "\n",
    "Contrary to common believes and folklore, star clusters are not isolated, but embedded in the potential of the Galaxy at large.\n",
    "This notion became even more apparent with the measurements of the Gaia satellite. Here we will simulate the evolution of a star cluster in a background (Galactic) potential.\n",
    "\n",
    "For reference you can read chapter 7 of Portegies Zwart & McMillan 2018 (2018araa.book.....P). \n",
    "\n",
    "With this tutorial you will learn\n",
    "- to generate inital conditions for star clusters\n",
    "- to initialize stellar evolution and gravitational N-body codes\n",
    "- channels and intra-code data transfer\n",
    "- making a class in Python\n",
    "- bridge two codes\n",
    "- plotting results\n",
    "- making cumulative distributions\n",
    "- Single-directional hierarchical code coupling strategy (i.e. classic bridge).\n",
    "- Bridge timesteps.\n",
    "- Incorporating an external potential to an N-body simulation\n",
    "- Appreciate the role of get_gravity_at_point function in bridge.\n",
    "- Appreciate the role of get_potential_at_point function in bridge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy\n",
    "from matplotlib import pyplot\n",
    "from amuse.units import units, constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by creating a class in which we describe the background (static) potential of the Milky way Galaxy.\n",
    "This is a very simple background potential model, which you can find in [Binney & Trainaine](https://ui.adsabs.harvard.edu/abs/2008gady.book.....B/abstract).\n",
    "It has components for the bulge, disk and halo.\n",
    "The coding is a bit arcane, but the routine is short and it has a simple single purpose.\n",
    "\n",
    "The class definition was quite simple. But the two main routines get_potential_at_point and get_gravity_at_point are the ones that do the job. These routines are used in the integrator the calculate the potential and local gravity at a particular point in space. The format of the argument list of these functions is fixed, because it is used elsewhere in AMUSE for integrating the cluster (which we still have to define)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class MilkyWay_galaxy(object):\n",
    "    def __init__(self, \n",
    "                 Mb=1.40592e10| units.MSun,\n",
    "                 Md=8.5608e10| units.MSun,\n",
    "                 Mh=1.07068e11 | units.MSun):\n",
    "        self.Mb= Mb\n",
    "        self.Md= Md\n",
    "        self.Mh= Mh\n",
    "\n",
    "    def get_potential_at_point(self,eps,x,y,z):\n",
    "        r=(x**2+y**2+z**2)**0.5\n",
    "        R= (x**2+y**2)**0.5\n",
    "        # buldge\n",
    "        b1= 0.3873 |units.kpc\n",
    "        pot_bulge= -constants.G*self.Mb/(r**2+b1**2)**0.5 \n",
    "        # disk\n",
    "        a2= 5.31 |units.kpc\n",
    "        b2= 0.25 |units.kpc\n",
    "        pot_disk = -constants.G*self.Md/(R**2 + (a2+ (z**2+ b2**2)**0.5 )**2 )**0.5\n",
    "        #halo\n",
    "        a3= 12.0 |units.kpc\n",
    "        cut_off=100 |units.kpc\n",
    "        d1= r/a3\n",
    "        c=1+ (cut_off/a3)**1.02\n",
    "        pot_halo= -constants.G*(self.Mh/a3)*d1**1.02/(1+ d1**1.02) \\\n",
    "                  - (constants.G*self.Mh/(1.02*a3))\\\n",
    "                      * (-1.02/c +numpy.log(c) + 1.02/(1+d1**1.02) \\\n",
    "                           - numpy.log(1.0 +d1**1.02) )\n",
    "        return 2*(pot_bulge+pot_disk+ pot_halo) \n",
    "                # multiply by 2 because it is a rigid potential\n",
    "    \n",
    "    def get_gravity_at_point(self, eps, x,y,z): \n",
    "        r= (x**2+y**2+z**2)**0.5\n",
    "        R= (x**2+y**2)**0.5\n",
    "        #bulge\n",
    "        b1= 0.3873 |units.kpc\n",
    "        force_bulge= -constants.G*self.Mb/(r**2+b1**2)**1.5 \n",
    "        #disk\n",
    "        a2= 5.31 |units.kpc\n",
    "        b2= 0.25 |units.kpc\n",
    "        d= a2+ (z**2+ b2**2)**0.5\n",
    "        force_disk=-constants.G*self.Md/(R**2+ d**2 )**1.5\n",
    "        #halo\n",
    "        a3= 12.0 |units.kpc\n",
    "        d1= r/a3\n",
    "        force_halo= -constants.G*self.Mh*d1**0.02/(a3**2*(1+d1**1.02))\n",
    "       \n",
    "        ax= force_bulge*x + force_disk*x  + force_halo*x/r\n",
    "        ay= force_bulge*y + force_disk*y  + force_halo*y/r\n",
    "        az= force_bulge*z + force_disk*d*z/(z**2 + b2**2)**0.5 + force_halo*z/r \n",
    "\n",
    "        return ax,ay,az"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.lab import Particles\n",
    "from amuse.lab import nbody_system\n",
    "from amuse.couple import bridge\n",
    "star = Particles(1)\n",
    "star.mass = 1 |units.MSun\n",
    "star.position = (1.0,0,0) * (8.5 | units.kpc)\n",
    "star.velocity = (0,-1.0,0) * (22 | units.kms)\n",
    "converter=nbody_system.nbody_to_si(star.mass.sum(), \n",
    "                                   star.position.length())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The single star has a mass, position and velocity. We can now initiate the N-body code, which is not different than what we have already seen."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.hermite.interface import Hermite\n",
    "gravity_code = Hermite(converter)\n",
    "gravity_code.particles.add_particles(star)\n",
    "ch_g2l = gravity_code.particles.new_channel_to(star)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We not initialize the Milky-way back-ground potential code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MWG = MilkyWay_galaxy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the bridge.\n",
    "The latter, in the next snippet, is like an integrator, except that it takes two or more codes to integrate with respect to each other.\n",
    "In this case we initalize the bridge and then add a system (the N-body system) to it. It is being perturbed by the Milky Way galaxy code, which we add as the second argument (in a list of length one; there can be multiple perturbing codes). Since there are no particles in the Milky way, they do not require any interaction (they do not need to be updated). For this reason we do not need to add it as a separate system.\n",
    "\n",
    "Eventually, we provide a timestep to the bridge. This is the time-scale in which the two systems are integrated. Regretfully, and this is one of the major disadvantages of *bridge* its timestep has to be set separately. Of course, this can be done automatically in the script, but sometimes it requires some fine-tuning and pre-knowledge in order to find the optimium between perserving energy while integrating the equations of motion and the speed of the calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gravity = bridge.Bridge(use_threading=False)\n",
    "gravity.add_system(gravity_code, (MWG,) )\n",
    "gravity.timestep = 20|units.Myr\n",
    "\n",
    "times = numpy.arange(0., 250, 20) | units.Myr\n",
    "x = [] | units.kpc\n",
    "y = [] | units.kpc\n",
    "for time in times:\n",
    "    gravity.evolve_model(time)\n",
    "    ch_g2l.copy()\n",
    "    x.append(star[0].x)\n",
    "    y.append(star[0].y)\n",
    "gravity.stop()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import plot, scatter\n",
    "scatter(star.x, star.y, s=10)\n",
    "plot(x, y, lw=1)\n",
    "pyplot.gca().set_aspect(\"equal\", adjustable=\"box\")\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like we can integrate the orbit of a single star around the Galactic center, and with an orbital velocity of 200 km/s the orbit turns out to be circular (or close enough to being circular).\n",
    "\n",
    "Now, we used a direct N-body code (Hermite) for integrating the single star, which is a bit silly. Of course, we only need to update the star's position and velocity in the potential and integrate those. So long as there are no other stars around, there is no need for a direct N-body code.\n",
    "\n",
    "On the other hand, the N-body code offers us the opportunity to integrate a cluster of stars, rather than just a single star. In this way, we can study the evolution of a cluster in the Galactic potential.\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Assignment 1:\n",
    "The output of the figure is not a nice smooth orbit, but a bit jagged. Change the script in such a way that the orbit becomes more smooth.\n",
    "\n",
    "### Assignment 2:\n",
    "The integration was for an (almost) circular orbit. Sometimes, however, we are interested in more eccentric orbits.\n",
    "Integrate the orbit of the same star (at a starting distance of 8.5kpc from the Galactic center), but with a velocity in the *y*-direction of *22km/s*.\n",
    "\n",
    "### Question 1:\n",
    "If your star escaped the potentical of the Galaxy, why is that the case. What change do you have to make to the script to assure that the orbit is still correct?\n",
    "\n",
    "### Question 2:\n",
    "Can you describe the orbit of the star in the eccentric orbit around the Galactic center?\n",
    "\n",
    "### Assignment 2:\n",
    "Replace the single star with a cluster of 10 stars with masses from the Salpeter mass function (between 1MSun and 100Msun) in a Plummer sphere with a characteristic radius of 10pc. Run the simulation and plot the cluster center-of-mass while it orbits the Galactic center.\n",
    "\n",
    "### Question 3: \n",
    "Describe how the cluster center-of-mass orbit differs compared to the orbit of the single star.\n",
    "\n",
    "### Question 4:\n",
    "Does the cluster stay bound?\n",
    "\n",
    "### Assignment 3:\n",
    "Now change the orbital velocity of the cluster center-of-mass to 22km/s, as you did before in assignment 2, and redo the calculation.\n",
    "Run the script again and describe how different the cluster orbit is now compared to the orbit of the single star.\n",
    "\n",
    "### Question 5:\n",
    "At what point in time does the cluster become unbound?\n",
    "\n",
    "### Assignment 4:\n",
    "Redo the calculation of the cluster with 10 stars in an eccentric orbit around the Galactic center, but with a slightly different orbital velocity for each calculation, ranging from 22km/s, 44km/s, 66km/s, up to 440km/s.\n",
    "\n",
    "While doing the calculation plot the shortest distance of the cluster center to the Galactic center and measure the moment the cluster becomes unbound.\n",
    "\n",
    " - Plot the distance of closest approach as a function of the initial orbital velocity of the cluster.\n",
    " \n",
    " - Plot the moment the cluster dissolves as a function of its initial orbital velocity.\n",
    "\n",
    "**Make sure that you use the same cluster for each calculation.**\n",
    "\n",
    "### Question 6:\n",
    "At what speed does the cluster escape the Galaxy?\n",
    "And show with theoretical arguments that you could have expected this to happen at that particular speed.\n",
    "\n",
    "### Question 7:\n",
    "Explain the curve from plotting the cluster dissolution time as function of its orbital velocity in the Galactic potential."
   ]
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}