running_stellarevolution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How to run a stellar evolution code\n",
    "====================\n",
    "\n",
    "Here we will generate initial conditions for a stellar evolution code. We will adopt a parametrized stellar evolution code.\n",
    "\n",
    "For reference you can read chapter 3 of Portegies Zwart & McMillan 2018 (2018araa.book.....P). \n",
    "\n",
    "With this tutorial you will learn\n",
    "- how to generate inital conditions (mass-function)\n",
    "- how to run a parametrized stellar evolution code\n",
    "- channels and data transfer\n",
    "- plotting results"
   ]
  },
  {
   "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 the initial mass function. What would you say is the biggest difference between these IMFs?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.datamodel import Particles\n",
    "n_stars = 1024\n",
    "mmin = 0.1 | units.MSun\n",
    "mmax = 100 | units.MSun\n",
    "\n",
    "from amuse.lab import new_kroupa_mass_distribution\n",
    "mkroupa = new_kroupa_mass_distribution(n_stars,\n",
    "                                           mass_min=mmin, \n",
    "                                           mass_max=mmax)\n",
    "k_stars = Particles(mass=mkroupa)\n",
    "from amuse.lab import new_salpeter_mass_distribution\n",
    "msalpeter = new_salpeter_mass_distribution(n_stars, \n",
    "                                           mass_min=mmin, \n",
    "                                           mass_max=mmax)\n",
    "s_stars = Particles(mass=msalpeter)\n",
    "print(\"mean-mass for Kroupa=\", numpy.mean(k_stars.mass))\n",
    "print(\"mean-mass for Salpeter=\", numpy.mean(s_stars.mass))\n",
    "\n",
    "pyplot.hist(numpy.log10(k_stars.mass.value_in(units.MSun)))\n",
    "pyplot.title('Kroupa IMF')\n",
    "pyplot.yscale('log')\n",
    "pyplot.xlim(-1, 2)\n",
    "pyplot.figure()\n",
    "pyplot.hist(numpy.log10(s_stars.mass.value_in(units.MSun)))\n",
    "pyplot.title('Salpeter IMF')\n",
    "pyplot.yscale('log')\n",
    "pyplot.xlim(-1, 2)\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the mass functions are determined, we can evolve the stars in time. For this we start two stellar evolution codes, both parametrized, and evolve the stars in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.community.seba.interface import SeBa\n",
    "def start_stellar_code(stars):\n",
    "    stellar = SeBa()\n",
    "    stellar.particles.add_particles(stars)\n",
    "    channels = {\"to_stars\": stellar.particles.new_channel_to(stars), \n",
    "                \"to_stellar\": stars.new_channel_to(stellar.particles)}\n",
    "    return stellar, channels\n",
    "kstellar, kchannels = start_stellar_code(k_stars)\n",
    "sstellar, schannels = start_stellar_code(s_stars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the stellar evolution code are declared. Now evolve in time and calculate the mean mass as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = 10**numpy.arange(0.0, 4.0, 0.1) | units.Myr\n",
    "mmean = []\n",
    "for time in times:\n",
    "    kstellar.evolve_model(time)\n",
    "    kchannels[\"to_stars\"].copy()\n",
    "    sstellar.evolve_model(time)\n",
    "    schannels[\"to_stars\"].copy()\n",
    "    mmean.append(numpy.mean(k_stars.mass)/numpy.mean(s_stars.mass))\n",
    "kstellar.stop()\n",
    "sstellar.stop()\n",
    "print(\"mean-mass for Kroupa=\", numpy.mean(k_stars.mass))\n",
    "print(\"mean-mass for Salpeter=\", numpy.mean(s_stars.mass))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import plot\n",
    "plot(times, mmean)\n",
    "pyplot.ylabel(\"relative mean mass\")\n",
    "pyplot.semilogx()\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.plot import scatter\n",
    "scatter(s_stars.temperature, s_stars.luminosity, c='r')\n",
    "scatter(k_stars.temperature, k_stars.luminosity, c='b', s=3)\n",
    "pyplot.xlim(2.e+4, 2000)\n",
    "pyplot.ylim(1.e-5, 1000)\n",
    "pyplot.loglog()\n",
    "pyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have calculated the evolution of a stellar population with solar metalicity.\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Assignment 1:\n",
    "Run the notebook and describe the evolution of the relative mean stellar mass. Why does it drop with time?\n",
    "\n",
    "### Assignment 2:\n",
    "Compare the mean mass of two stellar populations with different metalicity, 0.02 (solar) and 0.002 (one tenth times solar) but the same stellar masses.\n",
    "\n",
    "### Question 1:\n",
    "Is there a difference in the number of compact objects (white dwarfs, neutron stars, black holes; these types are defined [here](https://github.com/amusecode/amuse/blob/main/src/amuse/units/units.py), and e.g. 1.|stellar_type corresponds to a main sequene star) in the populations with different initial mass functions, and with different metallicities, and is this difference statistically significant?\n",
    "\n",
    "Is a similar difference also noticeable in the number of main sequence stars?\n",
    "\n",
    "### Assignment 3:\n",
    "Calculate, instead of the mean mass, the total luminosity and luminosity-weighted temperature ($T_{LW}=(\\Sigma_iL_iT_i)/(\\Sigma_iL_i)$, the effective color of the integrated stellar population) of the entire population of stars as a function of time.\n",
    "\n",
    "Make a plot of these parameters for both mass functions, and both metallicities. \n",
    "\n",
    "### Question 4:\n",
    "Can one, based on observational differences in luminosity and integrated temperature, make a distinction between the two mass functions?"
   ]
  },
  {
   "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
}