running_an_Nbody_code.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to run a simple N-body code\n",
    "====================\n",
    "\n",
    "Here we will generate initial conditions for an N-body code, run a small simulation and analyse the results. This analysis is performed on a 100 star cluster in a 1 pc virial-radius King model. Stellar masses are taken randomly from a Salpeter distribution.\n",
    "Stellar evolution is ignored in this simulation.\n",
    "\n",
    "For reference you can read chapter 2 of Portegies Zwart & McMillan 2018 (2018araa.book.....P). \n",
    "\n",
    "With this tutorial you will learn\n",
    "- how to generate initial conditions (mass function, King-model)\n",
    "- to initialize a gravitational N-body code\n",
    "- channels and intra-code data transfer\n",
    "- detecting binaries\n",
    "- plotting results\n",
    "- making cumulative distributions"
   ]
  },
  {
   "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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by setting-up a simulation by specifying a stellar mass distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_stars = 100\n",
    "alpha_IMF = -2.35\n",
    "\n",
    "from amuse.lab import new_powerlaw_mass_distribution\n",
    "m_stars = new_powerlaw_mass_distribution(n_stars, 0.1|units.MSun, \n",
    "                                        100|units.MSun, alpha_IMF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now check to see if this indeed gives one the expected mass function by plotting the cumulative distribution against a theoretical distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_cdf(m, alpha_IMF=-2.35):\n",
    "    m = sorted(m.value_in(units.MSun))\n",
    "    fm = numpy.array([0])\n",
    "    for mi in m:\n",
    "        fm = numpy.append(fm, fm[-1] + mi)\n",
    "    fm /= max(fm)\n",
    "    from amuse.plot import plot\n",
    "    pyplot.plot(m, fm[:len(m)], lw=4, ls=\"-\")\n",
    "    x = 10**numpy.arange(-1.0, 1.0, 1./len(m))\n",
    "    y = x**(alpha_IMF+2)\n",
    "    fy = numpy.array([0])\n",
    "    for yi in y:\n",
    "        fy = numpy.append(fy, fy[-1] + yi)    \n",
    "    fy /= max(fy)\n",
    "    pyplot.plot(x, fy[:len(x)], lw=2, ls=\":\")\n",
    "    pyplot.semilogx()\n",
    "    pyplot.show()\n",
    "plot_cdf(m_stars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question:\n",
    "Can you explain why the dotted curve does not overplot with the solid curve?\n",
    "\n",
    "Now we can initialize the cluster size and declare the converter.\n",
    "\n",
    "### Converter\n",
    "N-body codes operate using special dimensionless N-body units (where the gravitational constant is set to unity). In order to make the code understand SI units, we must pass it a unit converter. This contains two quantities in the relevant unit system, from which AMUSE derives the proper unit conversion. Best practice is to keep those quantities on the order of the system's scale (e.g. the total mass and the radius of the outermost planet in a planetary system) to ensure numerical stability.\n",
    "\n",
    "Also see the notebook tutorial '03-generic_units.ipynb' on the importance of converters in AMUSE. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_cluster = 1.0 | units.parsec\n",
    "from amuse.units import nbody_system\n",
    "converter=nbody_system.nbody_to_si(m_stars.sum(),r_cluster)\n",
    "from amuse.ic.kingmodel import new_king_model\n",
    "W0 = 3.0\n",
    "bodies=new_king_model(n_stars, W0, convert_nbody=converter)\n",
    "bodies.scale_to_standard(converter)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check the result by plotting the X-Y positions of the bodies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import scatter\n",
    "def plot_snapshot(bodies):\n",
    "    v = (bodies.vx**2 + bodies.vy**2 + bodies.vz**2).sqrt()\n",
    "    scatter(bodies.x, bodies.y, c=v.value_in(units.kms), alpha=0.5)\n",
    "    pyplot.colorbar()\n",
    "    pyplot.show()\n",
    "plot_snapshot(bodies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.ph4.interface import ph4\n",
    "from amuse.ext.LagrangianRadii import LagrangianRadii\n",
    "\n",
    "gravity = ph4(converter)\n",
    "gravity.particles.add_particles(bodies)\n",
    "# A channel is a 'permanent' connection to a code's particle\n",
    "# set. Multiple calls to a code's particle set need to set up\n",
    "# a new connection every time; with a channel, we can copy\n",
    "# information back without opening a new connection.\n",
    "# This does not automatically update bodies! See below\n",
    "channel = gravity.particles.new_channel_to(bodies)\n",
    "\n",
    "\n",
    "times = numpy.arange(0, 100, 0.1) | units.Myr\n",
    "RL25 = [] | units.parsec\n",
    "Rvir = [] | units.parsec\n",
    "for time in times:\n",
    "    gravity.evolve_model(time)\n",
    "    channel.copy() # Copy from gravity.particles to bodies\n",
    "    Rvir.append(bodies.virial_radius())\n",
    "    L = LagrangianRadii(bodies)\n",
    "    RL25.append(LagrangianRadii(bodies)[5])\n",
    "\n",
    "    if not time.value_in(units.Myr)%10.0:\n",
    "        print(\"cluster at Time=\", time.in_(units.Myr), \n",
    "              \"Mass=\", bodies.mass.sum().in_(units.MSun),\n",
    "              \"Rvir=\", Rvir[-1].in_(units.parsec))\n",
    "    b = bodies.get_binaries()\n",
    "    if(len(b)>0):\n",
    "        print(\"Number of binaries found:\", len(b))\n",
    "        print(\"first binary:\", b[0])\n",
    "\n",
    "pyplot.plot(times.value_in(units.Myr), RL25.value_in(units.parsec))\n",
    "pyplot.plot(times.value_in(units.Myr), Rvir.value_in(units.parsec))\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have now calculated the dynamical evolution of a small cluster of stars including a realistic mass function.\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Question 1:\n",
    "Why do we have to *scale_to_standard()* after assigning the masses and positions to all the particles?\n",
    "\n",
    "### Question 2:\n",
    "Run the notebook and describe the evolution of the virial radius.\n",
    "\n",
    "### Assignment 2:\n",
    "There is a check for binary detection in the script.\n",
    "If no binary formed in your run, that is okay.\n",
    "\n",
    "Run the script again until a binary forms.\n",
    "\n",
    "### Question 3:\n",
    "One of the parameters in the binary is the *hardness*.\n",
    "What is the hardness of the binary from your last run, and explain what is meant with this parameter.\n",
    "\n",
    "### Assignment 3:\n",
    "Rewrite the script in such a way that the event loop over time continues indefinitely, but that the code stops as soon as the first binary is detected.\n",
    "\n",
    "Run the script several times until you have enough first binaries formed to make a histogram of the first-binary formation time-scale.\n",
    "\n",
    "Now, make a small change to your initial conditions.\n",
    "A variation could include the dimensionless depth of the potential of the initial King model, the slope of the initial-mass function, its lower limit, or the number of stars in the simulation).\n",
    "\n",
    "Now redo the calculation in which you generate a histogram of moments of first-binary formation.\n",
    "\n",
    "Make one figure with both distributions: the original distribution for first-binary formation time and the second series of runs (with the altered initial conditions.\n",
    "Present both distributions as histograms and as cumulative distributions. \n",
    "\n",
    "### Question 4:\n",
    "Can you understand the difference in first-binary formation-time based on the changes you introduced in the initial conditions?\n",
    "\n",
    "Did you perform enough runs to make a statistically significant statement about the time of first-binary formation?\n",
    "\n",
    "### Assignment 4:\n",
    "Make a cumulative distribution of the masses of the two stars for each first binary and compare it to the initial mass-function.\n",
    "\n",
    "### Question 5:\n",
    "Can you understand the difference between the typical masses of the binaries that formed first and the initial mass-function?"
   ]
  },
  {
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}