02 - Particles.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "# preamble script to check and install AMUSE components if necessary\n",
    "\n",
    "# required packages for this tutorial:\n",
    "PACKAGES=\"amuse-framework\"\n",
    "# skip in case a full development install is present\n",
    "pip show amuse-devel && exit 0\n",
    "for package in ${PACKAGES} \n",
    "do\n",
    "  pip show ${package} || pip install ${package}\n",
    "done"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "AMUSE tutorial on particle sets\n",
    "====================\n",
    "\n",
    "\n",
    "This tutorial introduces you the the most used data store in AMUSE, a *particle set*. AMUSE particle sets are a handy tool for storing data of stellar or planetary systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.units import units, constants\n",
    "from amuse.lab import Particles, Particle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A particle set store the properties of a system under study, for example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sun_and_earth = Particles(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "of course at this point we have not specified anything about the system, lets change that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# declare the sun\n",
    "sun = sun_and_earth[0]\n",
    "sun.mass = 1 | units.MSun\n",
    "sun.position = (0, 0, 0) | units.au\n",
    "sun.velocity = (0, 0, 0) | units.kms\n",
    "print(f\"Sun={sun}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now declare the Earth\n",
    "def relative_orbital_velocity(mass, distance):\n",
    "    return (constants.G * mass / distance).sqrt()\n",
    "\n",
    "earth = sun_and_earth[1]\n",
    "earth.mass = 1 | units.MEarth\n",
    "earth.position = (1, 0, 0) | units.au\n",
    "\n",
    "vorb = relative_orbital_velocity(\n",
    "    sun_and_earth.total_mass(), \n",
    "    earth.position.length()\n",
    ")\n",
    "earth.velocity = (0, 1, 0) * vorb\n",
    "print(f\"Earth={earth}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# normally it is a good idea to study a system in its center of mass:\n",
    "sun_and_earth.move_to_center()\n",
    "print(sun_and_earth)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's give the particles a specific name (or other attribute)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sun_and_earth.name = [\"sun\", \"earth\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How we have declared the particles and moved them to the center of mass.\n",
    "We can also search for a specific particle.\n",
    "For example, the one with the \"earth\" in the attribute \"name\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earth = sun_and_earth[sun_and_earth.name==\"earth\"][0] # note the index! do you understand the difference when removing it?\n",
    "print(f\"Earth={earth}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add a moon in orbit around the earth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moon = Particle()\n",
    "moon.name = \"moon\"\n",
    "moon.mass = 7.34767309e+22 | units.kg\n",
    "moon.position = (384400, 0, 0) | units.km\n",
    "vorb = relative_orbital_velocity(\n",
    "    earth.mass + moon.mass, \n",
    "    moon.position.length()\n",
    ")\n",
    "moon.velocity = (0, 1, 0) * vorb\n",
    "print(f\"moon={moon}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Moon, however, is now somewhere inside the Sun with zero velocity - which is not good.\n",
    "We will have to relocate the Moon to make it orbit around the Earth.\n",
    "We do that by simply adding the positions and velocity of Earth to the Moon's."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "moon.position += earth.position\n",
    "moon.velocity += earth.velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can add the moon to the Sun and Earth system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sun_and_earth.add_particle(moon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the use of the singular here, because we only add a single particle to the particle set sun_and_earth.\n",
    "It is probably better to rename the sun_and_earth now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar_system = sun_and_earth"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important now to recenter the entire system, because by adding the moon we shifted the center of mass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar_system.move_to_center()\n",
    "print(f\"Solar system:\\n{solar_system}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The key uniquely identifies a particle, but is otherwise never explicitly used:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(moon in solar_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now manipulate the planetary system, or query it.\n",
    "for example by querying the masses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"mass=\", solar_system.mass.in_(units.MEarth))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us a list of the masses of all objects, in units the the earth's mass. In fact, each of the particle's attributes is a simple numpy array: it can be assigned and manipulated as such.\n",
    "\n",
    "Another way to acquire the same information could be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"mass={solar_system.mass/solar_system[1].mass}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we get a reference to the attributes..lets make a copy of the solar system and play around with it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar_system2=solar_system.copy()\n",
    "x=solar_system2.x\n",
    "print(x is solar_system2.x)\n",
    "print(solar_system2.center_of_mass())\n",
    "x[0]+=0.1| units.AU\n",
    "print(solar_system2.center_of_mass())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this is not necessarily true if the particle set is stored in a code, and can not be relied upon!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In some (hopefully rare) cases you may want to use the particle set or its attributes as simple numpy arrays, without the units.\n",
    "This is easily achieved by stripping the unit from the array. This can be realied by explicitely querying the selected parameter with that specific unit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solar_system.position.value_in(units.parsec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which, in this case, gives you a 2-dimensional *numpy.array* (a copy!) of the positions of star, planet and moon in units of a parsec."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, you may want to query the particle set solar_system.\n",
    "For example by asking what attributes it has.\n",
    "This can be done as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dir(solar_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or get some general help on the underlying particle class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(solar_system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You have performed some rudimentary operations on a particle set.\n",
    "It is now time to experiment a little for yourself.\n",
    "\n",
    "\n",
    "Assignments and questions:\n",
    "---------------\n",
    "\n",
    "### Assignment 1:\n",
    "Add the planet Jupiter (see [Wikipedia](https://en.wikipedia.org/wiki/Jupiter)) to your small planetary system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 2:\n",
    "Your planetary system is notoriously planar, and initialy the earth and moon are positioned along the Cartesian x-axis with the velocity vector in the Cartesian y-direction. \n",
    "\n",
    "Make the Sun-Earth-Moon system more realistic by introducing a small inclination to the Earth's and Moon's orbits and by giving them a random mean anomaly.\n",
    "\n",
    "The [Orbital element module](https://github.com/amusecode/amuse/blob/main/src/amuse/ext/orbital_elements.py) of AMUSE could come in handy. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.ext.orbital_elements import new_binary_from_orbital_elements\n",
    "binary=new_binary_from_orbital_elements(1| units.MSun, 1| units.MEarth, 1| units.AU, eccentricity=0.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 3:\n",
    "Calculate the total  gravitational binding energy of solarsystem.\n",
    "\n",
    "Now displace the entire particle set by 100 parsec and give it a linear velocity of 100km/s in the z-direction.\n",
    "Then calculate the binding energy of the system again.\n",
    "\n",
    "Did the binding energy of the Solar system change by this translation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1:\n",
    "Particle sets have the attribute *get_binaries()*.\n",
    "If you use this function to check the binaries in your system you will find that (without Jupiter) you have 3 binaries. Explain why the Sun is in a binary with the Moon. You may want to take a look at the [source code](https://github.com/amusecode/amuse/blob/main/src/amuse/datamodel/particle_attributes.py)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment 4:\n",
    "Generate another particle set with a 2 solar-mass star and two planets of 10 and 100 Earth masses in circular orbits at 0.1 and 0.6 au.\n",
    "Place this second planetary system at apocenter around your Solar system (true anomaly of 180 degrees) at a semimajor axis of 60 au with an eccentricity of 0.6.\n",
    "Then move the entire system to the center of mass."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2:\n",
    "Which of the orbits of the binary star with planets from *Assignment 4* has the highest binding energy?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Appendix: the Particle attributes\n",
    "\n",
    "particle sets carry a predefined set of methods and calculated attributes - which may save you some work! Below, we listed some of the most useful, you are encouraged to experiment\n",
    "to see if you can understand what they do."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from amuse.ic.plummer import new_plummer_model\n",
    "p=new_plummer_model(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following provide particle quantities (some of these are aliases):\n",
    "`position`\n",
    "`velocity`\n",
    "`acceleration`\n",
    "`potential`\n",
    "`specific_kinetic_energy`\n",
    "`angular_momentum`\n",
    "\n",
    "The following provide information about the physical state of whole particle set:\n",
    "`center_of_mass`\n",
    "`center_of_mass_velocity`\n",
    "`kinetic_energy`\n",
    "`potential_energy`\n",
    "`virial_radius`\n",
    "`total_mass`\n",
    "`total_radius`\n",
    "`total_momentum`\n",
    "`total_angular_momentum`\n",
    "`moment_of_inertia`\n",
    "`dynamical_timescale`\n",
    "\n",
    "The following change the particle set:\n",
    "`move_to_center`\n",
    "`scale_to_standard`\n",
    "`rotate`\n",
    "`add_spin`\n",
    "\n",
    "The following are (mainly) used for the analysis of star clusters:\n",
    "`get_binaries`\n",
    "`densitycentre_coreradius_coredens`\n",
    "`cluster_core`\n",
    "`bound_subset`\n",
    "`LagrangianRadii`\n",
    "`mass_segregation_Gini_coefficient`\n",
    "\n",
    "The following provides some more miscelleaneous analysis tools:\n",
    "`find_closest_particle_to`\n",
    "`distances_squared`\n",
    "`nearest_neighbour`\n",
    "`Qparameter`\n",
    "`connected_components`\n",
    "`minimum_spanning_tree_length`\n",
    "`mass_segregation_ratio`\n",
    "`mass_segregation_ratio_from_nearest_neighbour`\n",
    "`correlation_dimension`\n",
    "`box_counting_dimension`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(p.virial_radius())\n",
    "cm=p.center_of_mass()\n",
    "print(p.LagrangianRadii(cm=cm))"
   ]
  },
  {
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}