Skip to content

Instantly share code, notes, and snippets.

@rieder

rieder/amuse_binary_cluster_tools.ipynb Secret

Created September 11, 2023 06:34
Show Gist options
  • Save rieder/768b27959ee64d0f7ebe5af346972de5 to your computer and use it in GitHub Desktop.
Save rieder/768b27959ee64d0f7ebe5af346972de5 to your computer and use it in GitHub Desktop.
Download ZIP
Jupyter example of four AMUSE tools - MASC, Multiples, Ekster and Fresco
Display the source blob
Display the rendered blob
Raw
amuse_binary_cluster_tools.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "a88c4196",
"metadata": {},
"source": [
"# Simulating dynamics of multiple stellar systems in clusters\n",
"**NOTE: this notebook will be updated during the week - please check for updates!**\n",
"\n",
"In this notebook, I show how to create a star cluster with binary stars using AMUSE with MASC.\n",
"\n",
"Then, I simulate an encounter of two binary stars using Multiples.\n",
"\n",
"Finally, I show how to make a visualisation of the cluster with Fresco.\n",
"\n",
"Please feel free to use and modify the code in this notebook if you find it useful - but please let me know if you do!\n",
"\n",
"~ Steven Rieder"
]
},
{
"cell_type": "markdown",
"id": "6de1b5eb-dab5-44f2-888b-2e46ef561e7e",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"## MASC\n",
"*make a star cluster - choose an IMF and spatial distribution, supports initial binaries.*\n",
"\n",
"MASC is a tool that integrates several functions in AMUSE to facilitate creating initial conditions.\n",
"It supports multiple spatial distributions (like Plummer, King and Fractal models), initial mass functions (Kroupa, Salpeter, flat) and can now also create initial binaries with various properties (semi-major axis, eccentricity, mass relation between the primary and secondary).\n",
"It can be run from the command line to create an AMUSE file, or directly from an AMUSE script (with more functionality)."
]
},
{
"cell_type": "markdown",
"id": "23978ec7-5eca-442e-94a4-5b2c2c98aa61",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"### Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8bee3a5e",
"metadata": {},
"outputs": [],
"source": [
"# The following Python packages and their dependencies are required for this notebook\n",
"# You can find installation instructions for AMUSE at https://www.amusecode.org\n",
"\n",
"# uncomment line below to install via pip - note that you must already have the dependencies!\n",
"#!pip install amuse-framework amuse-masc amuse-fractalcluster amuse-ph4 amuse-smalln amuse-kepler"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3faa4895",
"metadata": {},
"outputs": [],
"source": [
"# Imports used here - also needed for the examples below\n",
"import time\n",
"import numpy as np\n",
"import matplotlib as mpl\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_format='retina'\n",
"mpl.rcParams['figure.dpi'] = 144\n",
"\n",
"from amuse.datamodel import Particles\n",
"from amuse.units import units, nbody_system, constants\n",
"from amuse.ext.masc.cluster import new_star_cluster"
]
},
{
"cell_type": "markdown",
"id": "cb86b904",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"### Create a star cluster\n",
"The cell below sets up a star cluster, where the semi-major axis of the binaries kept constant.\n",
"The binaries can also be given different semi-major axes.\n",
"\n",
"More options are available (e.g. for mass distributions) and more can/will be added in the future."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e4e4eb28",
"metadata": {},
"outputs": [],
"source": [
"np.random.seed(1701)\n",
"number_of_binaries = 100\n",
"number_of_systems = 5000\n",
"semi_major_axis = np.ones(number_of_binaries) * (10 | units.au)\n",
"\n",
"single_stars, binary_stars, binary_systems = new_star_cluster(\n",
" star_distribution='fractal', # note that currently, fractal ignores pre-set random seeds... Known issue.\n",
" number_of_stars=number_of_systems, # note that this really means 'number of systems'\n",
" number_of_binaries=number_of_binaries,\n",
" return_binaries=True,\n",
" binary_semi_major_axis_distribution=semi_major_axis,\n",
" eccentricity_distribution=0,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "63023037",
"metadata": {},
"outputs": [],
"source": [
"print(\n",
" f\"Singles: {len(single_stars)}, \"\n",
" f\"binary systems: {len(binary_systems)}, \"\n",
" f\"total stars: {len(single_stars)+len(binary_stars)}\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "4b12e407-3c6c-4db3-ab4e-8a4542f046e5",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"### Helper functions for plotting"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6b8675ab",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"# This cell contains plot functions\n",
"\n",
"figure_width_px = 1600\n",
"px = 1 / plt.rcParams['figure.dpi'] | units.inch # pixel in inches\n",
"\n",
"# Note the syntax used by AMUSE to create quantities (a variable with a unit)\n",
"# x.value_in(unit) will convert a quantity x back to a number\n",
"\n",
"figure_width_inch = figure_width_px * px # .value_in(units.inch)\n",
"\n",
"def plot_binary(b):\n",
" aspect = 1\n",
" figsize = (\n",
" figure_width_inch.value_in(units.inch),\n",
" aspect * figure_width_inch.value_in(units.inch),\n",
" )\n",
" l = (b[0].position - b[1].position).length()\n",
" fig = plt.figure(figsize=figsize, dpi=dpi.value_in(units.inch**-1))\n",
" ax_zoom = fig.add_subplot(111)\n",
" ax_zoom.scatter(\n",
" b.x.value_in(units.pc),\n",
" b.y.value_in(units.pc),\n",
" edgecolors=None, s=10,\n",
" )\n",
" ax_zoom.set_xlim(\n",
" (b[0].x - l).value_in(units.pc),\n",
" (b[0].x + l).value_in(units.pc),\n",
" )\n",
" ax_zoom.set_ylim(\n",
" (b[0].y - l).value_in(units.pc),\n",
" (b[0].y + l).value_in(units.pc),\n",
" )\n",
" ax_zoom.quiver(\n",
" b.x.value_in(units.pc),\n",
" b.y.value_in(units.pc),\n",
" b.vx.value_in(units.kms),\n",
" b.vy.value_in(units.kms),\n",
" scale=50\n",
" )\n",
" ax_zoom.set_aspect(1)\n",
" return fig\n",
"\n",
"\n",
"def plot_interaction(\n",
" stars,\n",
" length_unit=units.au,\n",
" speed_unit=units.kms,\n",
" time=0 | units.yr,\n",
" interactions=0,\n",
" time_requested=0 | units.yr,\n",
" min_width=100 | units.au,\n",
" return_width=False,\n",
"):\n",
" \"plot a set of interacting stars, seen from the center-of-mass frame\"\n",
" aspect = 1\n",
" figsize = (\n",
" 0.5 * figure_width_inch.value_in(units.inch),\n",
" 0.5 * aspect * figure_width_inch.value_in(units.inch),\n",
" )\n",
" fig = plt.figure(figsize=figsize)\n",
" ax = fig.add_subplot(111, aspect=aspect)\n",
" \n",
" com = stars.center_of_mass()\n",
" comv = stars.center_of_mass_velocity()\n",
"\n",
" minx = stars.x.min()\n",
" maxx = stars.x.max()\n",
" centerx = minx + (maxx-minx)/2\n",
" miny = stars.y.min()\n",
" maxy = stars.y.max()\n",
" centery = miny + (maxy-miny)/2\n",
" minz = stars.z.min()\n",
" maxz = stars.z.max()\n",
" centerz = minz + (maxz-minz)/2\n",
" width = max(\n",
" (\n",
" (maxx-minx)**2\n",
" + (maxy-miny)**2\n",
" + (maxz-minz)**2\n",
" )**0.5,\n",
" min_width\n",
" ).value_in(length_unit)\n",
"\n",
" ax.set_title(\n",
" f\"{time.value_in(units.yr):07.1f} yr \"\n",
" f\"({time_requested.value_in(units.yr):07.1f} yr), \"\n",
" f\"{interactions} interactions\"\n",
" )\n",
" stars = stars.sorted_by_attribute('z')\n",
" plt.scatter(\n",
" (stars.x - centerx).value_in(length_unit),\n",
" (stars.y - centery).value_in(length_unit),\n",
" c=stars.color,\n",
" )\n",
" plt.quiver(\n",
" (stars.x - centerx).value_in(length_unit),\n",
" (stars.y - centery).value_in(length_unit),\n",
" (stars.vx - comv.x).value_in(speed_unit),\n",
" (stars.vy - comv.y).value_in(speed_unit),\n",
" )\n",
" ax.set_xlim(\n",
" centerx.value_in(length_unit)-0.8*width,\n",
" centerx.value_in(length_unit)+0.8*width\n",
" )\n",
" ax.set_ylim(\n",
" centery.value_in(length_unit)-0.8*width,\n",
" centery.value_in(length_unit)+0.8*width\n",
" )\n",
" ax.set_xlabel(f'x [{length_unit}]')\n",
" ax.set_ylabel(f'y [{length_unit}]')\n",
" if return_width:\n",
" return fig, width\n",
" return fig\n",
"\n",
" \n",
"def plot_stars_binaries(s, b, width=5 | units.pc, title=''):\n",
" length_unit = width.unit\n",
" aspect = 1.7\n",
" figsize = (\n",
" figure_width_inch.value_in(units.inch),\n",
" figure_width_inch.value_in(units.inch) / aspect,\n",
" )\n",
" fig = plt.figure(figsize=figsize)\n",
" ax = fig.add_subplot(121)\n",
" ax.scatter(\n",
" s.x.value_in(length_unit),\n",
" s.y.value_in(length_unit),\n",
" edgecolors=None, s=1,\n",
" )\n",
" ax.set_xlim(-width.value_in(length_unit), width.value_in(length_unit))\n",
" ax.set_ylim(-width.value_in(length_unit), width.value_in(length_unit))\n",
" ax.set_xlabel(f'x [{length_unit}]')\n",
" ax.set_ylabel(f'y [{length_unit}]')\n",
" ax.set_title(title)\n",
" ax.set_aspect(1)\n",
" \n",
" ax_zoom = fig.add_subplot(122)\n",
" \n",
" n = s.nearest_neighbour()\n",
" d = (n.position - s.position).lengths()\n",
"\n",
" zoom_in_mass = (b[0].mass + n[0].mass)\n",
" zoom_in_position = (\n",
" b[0].mass * b[0].position\n",
" + n[0].mass * n[0].position\n",
" ) / zoom_in_mass\n",
" zoom_in_velocity = (\n",
" b[0].mass * b[0].velocity\n",
" + n[0].mass * n[0].velocity\n",
" ) / zoom_in_mass\n",
" system_projected_size = (\n",
" (b[0].x - n[0].x)**2\n",
" + (b[0].y - n[0].y)**2\n",
" )**0.5\n",
" \n",
" ax_zoom.scatter(\n",
" zoom_in_position.x.value_in(length_unit),\n",
" zoom_in_position.y.value_in(length_unit),\n",
" edgecolors=None, s=30,\n",
" label='Center-of-mass',\n",
" )\n",
" ax_zoom.scatter(\n",
" b.x.value_in(length_unit),\n",
" b.y.value_in(length_unit),\n",
" edgecolors=None, s=30,\n",
" label='Component stars',\n",
" )\n",
" ax_zoom.set_xlim(\n",
" (zoom_in_position.x - system_projected_size).value_in(length_unit),\n",
" (zoom_in_position.x + system_projected_size).value_in(length_unit)\n",
" )\n",
" ax_zoom.set_ylim(\n",
" (zoom_in_position.y - system_projected_size).value_in(length_unit),\n",
" (zoom_in_position.y + system_projected_size).value_in(length_unit)\n",
" )\n",
" ax_zoom.quiver(\n",
" zoom_in_position.x.value_in(length_unit),\n",
" zoom_in_position.y.value_in(length_unit),\n",
" (zoom_in_velocity.x).value_in(units.kms),\n",
" (zoom_in_velocity.y).value_in(units.kms),\n",
" scale=10,\n",
" )\n",
" ax_zoom.quiver(\n",
" b.x.value_in(length_unit),\n",
" b.y.value_in(length_unit),\n",
" (b.vx - zoom_in_velocity.x).value_in(units.kms),\n",
" (b.vy - zoom_in_velocity.y).value_in(units.kms),\n",
" scale=10,\n",
" )\n",
" ax_zoom.set_aspect(1)\n",
" ax_zoom.set_xlabel(f'x [{length_unit}]')\n",
" ax_zoom.set_ylabel(f'y [{length_unit}]')\n",
" ax_zoom.set_title('Zoom in on a binary pair')\n",
" ax_zoom.legend()\n",
" fig.tight_layout()\n",
" return fig"
]
},
{
"cell_type": "markdown",
"id": "33bf48dc-66c7-4c36-919b-33920183d857",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"### Plot the resulting cluster"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a2435d3b",
"metadata": {},
"outputs": [],
"source": [
"starplot = plot_stars_binaries(\n",
" (binary_stars | single_stars), # union of binary_stars and single_stars, i.e. all stars with the binaries first\n",
" binary_stars,\n",
" title='Fractal cluster'\n",
")\n",
"plt.savefig('masc_example.png', dpi=144)"
]
},
{
"cell_type": "markdown",
"id": "0a6676b6-240b-4940-8aa3-883806be09bb",
"metadata": {},
"source": [
"## Multiples\n",
"*takes care of interactions between a binary pair and any system that could disrupt it.*\n",
"\n",
"Multiples is a solution to the problem encountered by gravitational dynamics codes when simulating large systems with hard binaries in them: the time step would get too short, slowing the code down to a crawl.\n",
"Multiples detects systems of two (or more) stars within a specified interaction radius, and replaces these for the gravity code by a single particle.\n",
"When this particle encounters another particle, Multiples calculates the orbital phase of the stars in the system with a Kepler integrator, and then handles the encounter to its conclusion, e.g. the formation of (a) new multiple system(s) or the dissolution into individual stars, using a gravity code optimized for a small number of particles.\n",
"In this way, gravitational dynamics are sped up by a large factor compared to using a pure N-body code.\n",
"\n",
"In the cells below I explain how to use Multiples, taking two binaries from the MASC cluster and putting them in a more convenient frame where they will collide within a short time."
]
},
{
"cell_type": "markdown",
"id": "097acdad-6e4f-4291-83f0-f2a1ae846a7f",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"### Imports and setup\n",
"The cell below imports the required AMUSE modules and creates a GravityWithMultiples class."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8cc95995",
"metadata": {},
"outputs": [],
"source": [
"# Now we simulate the system further using Multiples.\n",
"# This module uses Ph4 for the integration of the systems, \n",
"# but integrates the binaries using a Kepler integrator\n",
"# and interactions between stars and binary systems using a small-N integrator\n",
"# As a result, it can speed things up by a large factor when dealing with tight binaries.\n",
"# When this is not the case however, it can be much slower!\n",
"from amuse.community.ph4 import Ph4\n",
"from amuse.community.smalln import Smalln\n",
"from amuse.community.kepler import Kepler\n",
"import amuse.couple.multiples as multiples\n",
"\n",
"\n",
"# Multiples requires some setting up - you can use this class for that\n",
"class GravityWithMultiples:\n",
" def __init__(\n",
" self, stars,\n",
" gravity=Ph4, smalln=Smalln, kepler=Kepler,\n",
" converter=None, time_step=0.01 | units.Myr,\n",
" number_of_workers=1,\n",
" ):\n",
" self.smalln_code = smalln\n",
" if converter is None:\n",
" self.converter = nbody_system.nbody_to_si(stars.total_mass(), time_step)\n",
" else:\n",
" self.converter = converter\n",
" self.model_time = 0 | units.Myr\n",
"\n",
" self.smalln = None\n",
" self.gravity = gravity(convert_nbody=self.converter, number_of_workers=number_of_workers)\n",
" self.gravity.parameters.epsilon_squared = (self.converter.to_si(0.0 | nbody_system.length))**2\n",
" self.gravity.particles.add_particles(stars)\n",
"\n",
" stopping_condition = self.gravity.stopping_conditions.collision_detection\n",
" stopping_condition.enable()\n",
"\n",
" self.init_smalln()\n",
" self.kepler = kepler(unit_converter=self.converter)\n",
" self.kepler.initialize_code()\n",
" self.multiples_code = multiples.Multiples(\n",
" self.gravity, self.new_smalln, self.kepler, constants.G\n",
" )\n",
" self.multiples_code.neighbor_perturbation_limit = 0.05\n",
" self.multiples_code.global_debug = 0\n",
"\n",
" # global_debug = 0: no output from multiples\n",
" # 1: minimal output\n",
" # 2: debugging output\n",
" # 3: even more output\n",
"\n",
" print('Settings for Multiples:')\n",
" print(f'multiples_code.neighbor_veto = {self.multiples_code.neighbor_veto}')\n",
" print(f'multiples_code.neighbor_perturbation_limit = {self.multiples_code.neighbor_perturbation_limit}')\n",
" print(f'multiples_code.retain_binary_apocenter = {self.multiples_code.retain_binary_apocenter}')\n",
" print(f'multiples_code.wide_perturbation_limit = {self.multiples_code.wide_perturbation_limit}')\n",
"\n",
" self.energy_0 = self.print_diagnostics(self.multiples_code)\n",
"\n",
" @property\n",
" def stars(self):\n",
" return self.multiples_code.stars\n",
"\n",
" @property\n",
" def particles(self):\n",
" return self.multiples_code.particles\n",
"\n",
" def init_smalln(self):\n",
" self.smalln = self.smalln_code(convert_nbody=self.converter)\n",
"\n",
" def new_smalln(self):\n",
" self.smalln.reset()\n",
" return self.smalln\n",
"\n",
" def stop_smalln(self):\n",
" self.smalln.stop()\n",
"\n",
" def print_diagnostics(self, grav, energy_0=None):\n",
" # Simple diagnostics.\n",
" energy_kinetic = grav.kinetic_energy\n",
" energy_potential = grav.potential_energy\n",
" (\n",
" self.number_of_multiples,\n",
" self.number_of_binaries,\n",
" self.energy_in_multiples,\n",
" ) = grav.get_total_multiple_energy()\n",
" energy = energy_kinetic + energy_potential + self.energy_in_multiples\n",
" print(f'Time = {grav.get_time().in_(units.Myr)}')\n",
" print(f' top-level kinetic energy = {energy_kinetic}')\n",
" print(f' top-level potential energy = {energy_potential}')\n",
" print(f' total top-level energy = {energy_kinetic + energy_potential}')\n",
" print(f' {self.number_of_multiples} multiples, total energy = {self.energy_in_multiples}')\n",
" print(f' uncorrected total energy ={energy}')\n",
"\n",
" energy_tidal = (\n",
" grav.multiples_external_tidal_correction\n",
" + grav.multiples_internal_tidal_correction\n",
" ) # tidal error\n",
" energy_error = grav.multiples_integration_energy_error # integration error\n",
"\n",
" energy -= energy_tidal + energy_error\n",
" print(f' corrected total energy = {energy}')\n",
"\n",
" if energy_0 is not None: print(' relative energy error=', (energy-energy_0)/energy_0)\n",
"\n",
" return energy\n",
"\n",
" def evolve_model(self, t_end, return_number_of_encounters=False):\n",
" if return_number_of_encounters:\n",
" encounters_resolved, encounters_ignored = self.multiples_code.evolve_model(\n",
" t_end, return_number_of_encounters=return_number_of_encounters\n",
" )\n",
" else:\n",
" self.multiples_code.evolve_model(t_end)\n",
" # self.print_diagnostics(self.multiples_code, self.energy_0)\n",
" self.model_time = self.multiples_code.model_time\n",
" if return_number_of_encounters:\n",
" return encounters_resolved, encounters_ignored\n",
" \n",
"\n",
" def stop(self):\n",
" self.gravity.stop()\n",
" self.kepler.stop()\n",
" self.stop_smalln()"
]
},
{
"cell_type": "markdown",
"id": "d4885a6b-942c-4965-a294-3f4b1b2d7e8a",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"### Create a pair of binaries\n",
"These will collide after a short time."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7b541dcb-c7be-47e8-b19b-2ce21bf08edc",
"metadata": {},
"outputs": [],
"source": [
"# Use two binaries as a nice test case for Multiples.\n",
"# We create a copy of the binaries, so that we don't change the original particles.\n",
"colliding_binaries = Particles()\n",
"binary_pair_one = Particles()\n",
"binary_pair_one.add_particle(binary_stars[0].copy())\n",
"binary_pair_one.add_particle(binary_stars[number_of_binaries].copy())\n",
"binary_pair_two = Particles()\n",
"binary_pair_two.add_particle(binary_stars[1].copy())\n",
"binary_pair_two.add_particle(binary_stars[1+number_of_binaries].copy())\n",
"\n",
"# shift binaries to their respective center-of-mass coordinate frame\n",
"binary_pair_one.move_to_center()\n",
"binary_pair_two.move_to_center()\n",
"\n",
"# set the binaries 400 au apart - gravity will do the rest\n",
"binary_pair_one.x += 200 | units.au\n",
"binary_pair_two.x -= 200 | units.au\n",
"\n",
"# (setting a moderate vx helps the collision happen a bit sooner - it will also change the outcome though)\n",
"binary_pair_one.vx -= 1.1 | units.au / units.yr\n",
"binary_pair_two.vx += 1.1 | units.au / units.yr\n",
"binary_pair_one.vy -= 0.01 | units.au / units.yr\n",
"binary_pair_two.vy += 0.01 | units.au / units.yr\n",
"\n",
"# add them to a common particle set\n",
"colliding_binaries.add_particles(binary_pair_one)\n",
"colliding_binaries.add_particles(binary_pair_two)\n",
"\n",
"# give all stars their own color, so that we can uniquely identify them in plots\n",
"colliding_binaries.color = [\"red\", \"blue\", \"green\", \"orange\"]\n",
"\n",
"# The radius we set here is the *interaction* radius: Multiples will present systems smaller than this as a single star to the gravity code\n",
"# So don't set it too large or too small!\n",
"colliding_binaries.radius = 3 | units.au\n",
"# uncommenting the line below effectively disables Multiples\n",
"# colliding_binaries.radius = 2 | units.RSun"
]
},
{
"cell_type": "markdown",
"id": "91244e6b-5cfa-421e-add8-dd34c7e3f0c6",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"### Plot the initial setup"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9507a7d4",
"metadata": {},
"outputs": [],
"source": [
"colfig, width = plot_interaction(colliding_binaries, time=0 | units.yr, return_width=True)\n",
"plt.show()\n",
"colfig.clf()"
]
},
{
"cell_type": "markdown",
"id": "ee23b2cf-fdb5-4c3d-aae7-2bfd3e155f2d",
"metadata": {},
"source": [
"### Simulate binary-binary interaction with Multiples\n",
"We now let the binary systems collide. Multiples will at times take over from the gravity code, each time this happens the 'interactions' counter will increase by 1.\n",
"This will create a lot of plots if save_plots is set to True!"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "434416b5",
"metadata": {},
"outputs": [],
"source": [
"# A converter sets the units used internally in the code, using two units out of mass, length and time.\n",
"converter_multiples = nbody_system.nbody_to_si(\n",
" 0.1 | units.au,\n",
" 0.1 | units.yr,\n",
")\n",
"grav = GravityWithMultiples(\n",
" colliding_binaries,\n",
" converter=converter_multiples,\n",
" number_of_workers=1, # if we simulate a large cluster, best increase this - but here we just do 4 stars\n",
")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bbc90e19",
"metadata": {},
"outputs": [],
"source": [
"save_plots = True\n",
"i = 0\n",
"width = 0\n",
"time_step = 0.8 | units.yr\n",
"time_end = 1000 | units.yr\n",
"encounters_resolved = 0\n",
"encounters_ignored = 0\n",
"clock_start = time.time() | units.s\n",
"while grav.model_time < time_end:\n",
" t = i * time_step\n",
" res, ign = grav.evolve_model(t, return_number_of_encounters=True)\n",
" encounters_resolved += res\n",
" encounters_ignored += ign\n",
" channel = grav.stars.new_channel_to(colliding_binaries)\n",
" channel.copy_attributes([\"x\", \"y\", \"z\", \"vx\", \"vy\", \"vz\"])\n",
" if save_plots:\n",
" colfig, width = plot_interaction(\n",
" colliding_binaries, time=grav.model_time, interactions=encounters_resolved, time_requested=t,\n",
" return_width=True,\n",
" )\n",
" plt.savefig(f'interacting_binary_{i:06d}.png')\n",
" plt.close(colfig)\n",
" i += 1\n",
" print(f\"{encounters_resolved}, {encounters_ignored}, {width}, {grav.model_time.in_(units.yr)}\")\n",
"clock_finish = time.time() | units.s\n",
"print(f\"This took {(clock_finish-clock_start).in_(units.hour)}\")\n",
"print(f\"There are now {len(grav.particles)} particles in the gravity code\")\n",
"print(f\"There are however {len(grav.stars)} stars.\")\n",
"grav.stop()"
]
},
{
"cell_type": "markdown",
"id": "cc6870c8-627c-4176-b8ba-a69b98f86e96",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"## Fresco\n",
"*visualisation of stars with extinction from dust - create mock-observations from simulations.*\n",
"\n",
"Fresco (Rieder & Pelupessy, https://github.com/rieder/amuse-fresco) is a visualisation tool, that calculates the visual appearance of stars using their luminosity and temperature.\n",
"It calculates the stellar flux in several color bands, convolves this with a point-spread function (the PSF for Hubble's WFC3 is included) and combines the result to a color image (mimicking visual observation).\n",
"Additionally, Fresco supports calculating extinction and reflection by gas/dust."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d768bb4d-9fbb-414c-a122-5f3dec93297d",
"metadata": {},
"outputs": [],
"source": [
"# uncomment to install Fresco\n",
"#!pip install amuse-fresco\n",
"from amuse.plot.fresco import make_fresco_image\n",
"from amuse.plot.fresco.fresco import evolve_to_age\n",
"\n",
"stars_for_fresco = (binary_stars | single_stars).copy()\n",
"evolve_to_age(stars_for_fresco, 1 | units.Myr)\n",
"\n",
"make_fresco_image(stars=stars_for_fresco, gas=None)"
]
},
{
"cell_type": "markdown",
"id": "5724bff4-e7d4-46c0-9881-31fce341fd6e",
"metadata": {},
"source": [
"## Ekster\n",
"*star formation from molecular clouds at moderate mass resolution.*\n",
"\n",
"Ekster (Rieder et al. 2021) is a tool for simulating the formation of star clusters from molecular clouds.\n",
"It combines gravitational dynamics, hydrodynamics and stellar evolution with a recipe for star formation that forms stars on a proto-cluster basis, about 200 M<sub>☉</sub> at a time.\n",
"In Rieder et al. 2021, we used Ekster to simulate star formation in a spiral galaxy, while in Sills et al. 2022 we used it to simulate a possible future of the Orion nebula cluster based on observed gas- and star distributions.\n",
"\n",
"Ekster is still being developed, and currently requires some manual attention - please contact me if you want to try it!"
]
},
{
"cell_type": "markdown",
"id": "7e6a9b1f-2884-4df7-9455-da578b2f6329",
"metadata": {
"jp-MarkdownHeadingCollapsed": true
},
"source": [
"## Extra cells"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "20b2e8c2-6db9-4bf6-8202-7575d441f682",
"metadata": {},
"outputs": [],
"source": [
"# Extra cell - this is how to simulate a full star cluster with just a gravity code\n",
"\n",
"converter = nbody_system.nbody_to_si(\n",
" single_stars.total_mass() + binary_stars.total_mass(),\n",
" 0.01 | units.Myr\n",
")\n",
"\n",
"# We can simulate the cluster with the 'Ph4' code (Parallel Hermite 4th order),\n",
"# as long as the binaries aren't very tight.\n",
"# If they are, Ph4 will slow down to a crawl as the required timestep becomes too small.\n",
"# We then need to find another solution - like Multiples (below).\n",
"# The `number_of_workers` keyword indicates how many processes we will use for Ph4.\n",
"from amuse.community.ph4 import Ph4\n",
"grav = Ph4(converter, number_of_workers=6)\n",
"s_in_grav = grav.particles.add_particles(single_stars)\n",
"b_in_grav = grav.particles.add_particles(binary_stars)\n",
"\n",
"time_start = time.time() | units.s\n",
"grav.evolve_model(0.01 | units.Myr)\n",
"time_end = time.time() | units.s\n",
"print(f\"This took {time_end-time_start}\")\n",
"\n",
"# A channel is a way to copy updated properties from an in-code particle set back to our model.\n",
"\n",
"channel = b_in_grav.new_channel_to(binary_stars)\n",
"channel.copy()\n",
"channel = s_in_grav.new_channel_to(single_stars)\n",
"channel.copy()\n",
"grav.stop()\n",
"\n",
"# Note how the binary progressed in its orbit while the cluster as a whole is basically unchanged\n",
"# - this illustrates how much the timescales differ.\n",
"starplot = plot(binary_stars | single_stars, binary_stars)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "974ee08d",
"metadata": {},
"outputs": [],
"source": [
"# Extra cell - this is how to simulate a full star cluster with Multiples\n",
"\n",
"converter = nbody_system.nbody_to_si(\n",
" single_stars.total_mass() + binary_stars.total_mass(),\n",
" 0.01 | units.Myr\n",
")\n",
"\n",
"grav = GravityWithMultiples(\n",
" binary_stars | single_stars,\n",
" gravity=Ph4, smalln=Smalln, kepler=Kepler,\n",
" converter=converter,\n",
" number_of_workers=6,\n",
")\n",
"time_start = time.time() | units.s\n",
"grav.evolve_model(0.01 | units.Myr)\n",
"time_end = time.time() | units.s\n",
"print(f\"This took {time_end-time_start}.\")\n",
"channel = grav.particles.new_channel_to(binary_stars)\n",
"channel.copy()\n",
"channel = grav.particles.new_channel_to(single_stars)\n",
"channel.copy()\n",
"grav.stop()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ea1a3fc8",
"metadata": {},
"outputs": [],
"source": [
"# Let's plot the system and the binary again\n",
"starplot = plot(binary_stars | single_stars, binary_stars)"
]
},
{
"cell_type": "markdown",
"id": "fff11c43",
"metadata": {},
"source": [
"Below are some potentially useful methods/functions, without much explanation - mostly here to remember that they exist. Some are very slow so execute with caution!"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "99856d63",
"metadata": {},
"outputs": [],
"source": [
"# the 'get_binaries' method will search for binaries the hard way - use\n",
"# with caution and probably never use it on a large particle set!\n",
"# it will only return binaries of a specified (or default) 'hardness'\n",
"\n",
"stars_in_ph4 = grav.particles\n",
"q = stars_in_ph4.get_binaries() # control via 'hardness' keyword\n",
"print(len(q))"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7216706e",
"metadata": {},
"outputs": [],
"source": [
"# use this function to extract the orbital parameters from a binary.\n",
"# the module has other useful functions for dealing with binaries.\n",
"from amuse.ext.orbital_elements import get_orbital_elements_from_binary\n",
"\n",
"mass1, mass2, a, e, c_o_ta, c_o_inc, c_o_lotan, c_o_aop = get_orbital_elements_from_binary(\n",
" binary_stars[0] + binary_stars[number_of_binaries]\n",
")\n",
"print(mass1, mass2, a.in_(units.au), e)"
]
}
],
"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.10"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment