Define the initial conditions

Posidonius is composed by a N-body simulator written in Rust and a Python package. Rust is a great language but its learning curve is a bit steep at first, thus it is fantastic to have a more user-friendly scripting interface to define the kind of simulation that the user wants to execute. The user will only need to learn Rust for advance cases where the integrator or new forces have to be modified/added.

Included in Posidonius there are already several defined cases in the cases/ directory (including examples of planetary and binary systems), it is recommended to inspect them to better understand how the user can create its own scripts. Let's see an example of a brown dwarf host star with a earth-like planet around:

  • The script will create a JSON file with all the initial conditions, which will be read by Posidonius to run the simulation.
  • In the script, the first thing we need to do is define some general variable to specify:
    • The age at which the planetary system is when the simulation is going to start (used for stars/planets models)
    • The time step for the simulation
    • The time limit after which the simulation will stop (i.e., how many days we want to simulate)
    • How frequently we want to record the position of all the particles
    • How frequently we want to save a global snapshot that will allow us to resume our simulation if it was interruption for some random reason (e.g., power failure, unexpected system reboot)
    • The effects that we want to globally enable/disable 
      initial_time = 4.5e6*365.25 # time [days] where simulation starts
time_step = 0.08 # days
time_limit   = 365.25 * 1.0e2 # days
historic_snapshot_period = 100.*365.25 # days
recovery_snapshot_period = 100.*historic_snapshot_period # days
consider_effects = posidonius.ConsiderEffects({
"tides": True,
"rotational_flattening": True,
"general_relativity": True,
"disk": True,
"wind": True,
"evolution": True,
})
  • Now it's time to import Posidonius and create a universe: 
    import posidonius
universe = posidonius.Universe(initial_time, time_limit, time_step, 
                recovery_snapshot_period, historic_snapshot_period, 
                consider_effects)
  • We are going to define the host star of the planetary system, we will need:
    • Mass (in solar masses)
    • Radius factor (in solar radius)
    • Radius of gyration: computed in terms of the mass moment of inertia, which depends on the shape of the body and determines the torque needed for a desired angular acceleration.
    • Position and velocity: which should be zero for the three axes x, y, and z given that this is the central body.
    • Rotation period (hours): it will be determine the initial spin of the body. 
      star_mass = 0.08 # Solar masses
star_radius_factor = 0.845649342247916
star_radius = star_radius_factor * posidonius.constants.R_SUN
star_radius_of_gyration = 4.47e-01 # Brown dwarf
star_position = posidonius.Axes(0., 0., 0.)
star_velocity = posidonius.Axes(0., 0., 0.)
# Initialization of stellar spin
star_rotation_period = 70.0 # hours
star_angular_frequency = posidonius.constants.TWO_PI/(star_rotation_period/24.) # days^-1
star_spin = posidonius.Axes(0., 0., star_angular_frequency)
star = posidonius.Particle(star_mass, star_radius, star_radius_of_gyration, star_position, star_velocity, star_spin)
  • We can now define the effects that we will enable for the host star:
    • Tidal forces
      • The tidal model only allow interactions between one body (usually the most massive one, the host star in this example) and the rest of bodies. Therefore, to enable this force the simulation needs to have one body with this effect set to CentralBody and the rest with OrbitingBody (although not necessarily all of them). We can use Disabled when we do not want that this effect is taken into account. In this example, we are defining the effect for the host star, so we will use CentralBody and these two parameters:
        • Dissipation factor: fraction of the orbital energy that is dissipated as heat in either the surface or interior of a planetary body.
        • Love number of degree 2 (i.e., k2): dimensionless parameters that measure the rigidity of a planetary body and the susceptibility of its shape to change in response to a tidal potential. 
          star_tides_parameters = {
"dissipation_factor_scale": 1.0,
"dissipation_factor": 2.006*3.845764e4,
"love_number": 0.307,
}
star_tides = posidonius.effects.tides.CentralBody(star_tides_parameters)
#star_tides = posidonius.effects.tides.OrbitingBody(star_tides_parameters)
#star_tides = posidonius.effects.tides.Disabled()
    • Rotational flattening effects
      • This effects is defined very similarly to the tidal forces (i.e., CentralBody, OrbitingBody or Disabled), thus the same rules apply. In this example we use CentralBody for the host start, and we need to define the "fluid" love number, which is like the tidal one but for rotational flattening effects: 
        star_fluid_love_number = 0.307
star_rotational_flattening_parameters = {"love_number": star_fluid_love_number }
star_rotational_flattening = posidonius.effects.rotational_flattening.CentralBody(star_rotational_flattening_parameters)
#star_rotational_flattening = posidonius.effects.rotational_flattening.OrbitingBody(star_rotational_flattening_parameters)
#star_rotational_flattening = posidonius.effects.rotational_flattening.Disabled()
    • General relativity corrections
      • Again, this effects is defined very similarly to the tidal forces (i.e., CentralBody, OrbitingBody or Disabled), thus the same rules apply. There are three prescriptions available:
        • Kidder1995 as implemented in Mercury-T (it assumes the CentralBody is the most massive of the system)
        • Anderson1975, as implemented in REBOUND (it assumes the CentralBody is the most massive of the system)
        • Newhall1983 takes into consideration the effects from all bodies (no assumptions regarding the CentralBody), making Posidonius ready for binary stars 
          star_general_relativity = posidonius.effects.general_relativity.CentralBody("Kidder1995")
#star_general_relativity = posidonius.effects.general_relativity.CentralBody("Anderson1975")
#star_general_relativity = posidonius.effects.general_relativity.CentralBody("Newhall1983")
#star_general_relativity = posidonius.effects.general_relativity.OrbitingBody()
#star_general_relativity = posidonius.effects.general_relativity.Disabled()
    • Wind:
      • If the body produces winds, its effects on the spin can be taken into account by enabling it. It should always be disabled for the host star, only enable it for the planets if desired. In any case, this simulation we will not include wind effects: 
        #star_wind = posidonius.effects.wind.Interaction({
## Solar wind parametrisation (Bouvier 1997)
#"k_factor": 4.0e-18, # K_wind = 1.6d47 cgs, which is in Msun.AU2.day
#"rotation_saturation": 1.7592918860102842, # 14. * TWO_PI/25.0, in units of the spin of the Sun today
#})
star_wind = posidonius.effects.wind.Disabled()
    • Disk
      • A disk can be defined around one of the bodies in the simulation (CentralBody), so that the rest of bodies can feel the disk force (OrbitingBody) although some can be disabled too (Disabled). In this case, we will not enable a disk: 
        #disk_surface_density_normalization_gcm = 1000. # g.cm^-2
#disk_surface_density_normalization_SI = disk_surface_density_normalization_gcm * 1.0e-3 * 1.0e4 # kg.m^-2
#disk_properties = {
#'inner_edge_distance': 0.01,  # AU
#'outer_edge_distance': 100.0, # AU
#'lifetime': 1.0e5 * 365.25e0, # days
#'alpha': 1.0e-2,
#'surface_density_normalization': disk_surface_density_normalization_SI * (1.0/posidonius.constants.M_SUN) * posidonius.constants.AU**2, # Msun.AU^-2
#'mean_molecular_weight': 2.4,
#}
#star_disk = posidonius.effects.disk.CentralBody(disk_properties)
#star_disk = posidonius.effects.disk.OrbitingBody()
star_disk = posidonius.effects.disk.Disabled()
    • Body evolution
      • The radius and other parameters (depending on the model) can vary with time for each body in a simulation. Posidonius includes several models for different masses and from different authors, although we will not enable any for this body: 
        #star_evolution = posidonius.GalletBolmont2017(star_mass) # mass = 0.30 .. 1.40
#star_evolution = posidonius.BolmontMathis2016(star_mass) # mass = 0.40 .. 1.40
#star_evolution = posidonius.Baraffe2015(star_mass) # mass = 0.01 .. 1.40
#star_evolution = posidonius.Leconte2011(star_mass) # mass = 0.01 .. 0.08
#star_evolution = posidonius.Baraffe1998(star_mass) # Sun (mass = 1.0) or M-Dwarf (mass = 0.1)
#star_evolution = posidonius.LeconteChabrier2013(False) # Jupiter without dissipation of dynamical tides
#star_evolution = posidonius.LeconteChabrier2013(True) # Jupiter with dissipation of dynamical tides
star_evolution = posidonius.NonEvolving()
  • Once the effects are defined, we can assign them to the previously created star: 
    star.set_tides(star_tides)
star.set_rotational_flattening(star_rotational_flattening)
star.set_general_relativity(star_general_relativity)
star.set_wind(star_wind)
star.set_disk(star_disk)
star.set_evolution(star_evolution)
  • Finally, we are ready to add the star to the universe: 
    universe.add_particle(star)
  • Now we are going to define an earth-like planet, where we need the same parameters as the ones defined for the star, plus some differences (most of the following definitions/figures come from Wikipedia):
    • Position and velocities (stable orbits) will be computed from:
      • Semi-major axis (AU): one half of the longest diameter.
      • Eccentricity (between 0 and 1): amount by which its orbit around another body deviates from a perfect circle.
      • Inclination (degrees): vertical tilt of the ellipse with respect to the reference plane, measured at the ascending node (where the orbit passes upward through the reference plane, the green angle i in the diagram)
      • Argument of pericentre/periapsis (degrees): it defines the orientation of the ellipse in the orbital plane, as an angle measured from the ascending node to the periapsis (the closest point the satellite object comes to the primary object around which it orbits, the blue angle ω in the diagram).
      • Longitude of the ascending node (degrees): horizontally orients the ascending node of the ellipse (where the orbit passes upward through the reference plane) with respect to the reference frame's vernal point (the green angle omega in the diagram).
      • Mean anomaly (degrees): angular distance from the pericenter which a fictitious body would have if it moved in a circular orbit, with constant speed, in the same orbital period as the actual body in its elliptical orbit. This is a convenient mathematical angle which varies linearly with time, but which does not correspond to a real geometric angle. It can be converted into the true anomaly ν, which does represent the real geometric angle in the plane of the ellipse, between periapsis (closest approach to the central body) and the position of the orbiting object at any given time. Semi-major axis Eccentricity Orbital elements (source: Wikipedia)
    • Apart from the rotation period (hours), we can consider now also the obliquity (degrees), both will be determine the initial spin of the body. Obliquity 
      planet_mass_factor = 1.0
planet_mass = planet_mass_factor * posidonius.constants.M_EARTH # Solar masses (3.0e-6 solar masses = 1 earth mass)
#--------------------------------------------------------------------------------
# Earth-like => mass-radius relationship from Fortney 2007
planet_radius_factor = posidonius.tools.mass_radius_relation(planet_mass_factor, 
                                            planet_mass_type='factor', 
                                            planet_percent_rock=0.70)
planet_radius = planet_radius_factor * posidonius.constants.R_EARTH
planet_radius_of_gyration = 5.75e-01 # Earth type planet
#--------------------------------------------------------------------------------
#////////// Specify initial position and velocity for a stable orbit
#////// Keplerian orbital elements, in the `asteroidal' format of Mercury code
a = 0.018;                             # semi-major axis (in AU)
e = 0.1;                               # eccentricity
i = 5. * posidonius.constants.DEG2RAD; # inclination (degrees)
p = 0. * posidonius.constants.DEG2RAD; # argument of pericentre (degrees)
n = 0. * posidonius.constants.DEG2RAD; # longitude of the ascending node (degrees)
l = 0. * posidonius.constants.DEG2RAD; # mean anomaly (degrees)
p = (p + n);                           # Convert to longitude of perihelion !!
q = a * (1.0 - e);                     # perihelion distance
planet_position, planet_velocity = posidonius.calculate_cartesian_coordinates(
                                        planet_mass, q, e, i, p, n, l, 
                                        masses=[star_mass], 
                                        positions=[star_position], 
                                        velocities=[star_velocity])
#--------------------------------------------------------------------------------
#////// Initialization of planetary spin
planet_obliquity = 11.459156 * posidonius.constants.DEG2RAD # 0.2 rad
planet_rotation_period = 24. # hours
planet_angular_frequency = posidonius.constants.TWO_PI/(planet_rotation_period/24.) # days^-1
planet_keplerian_orbital_elements = posidonius.calculate_keplerian_orbital_elements(
                                            planet_mass, planet_position, planet_velocity, 
                                            masses=[star_mass], 
                                            positions=[star_position], 
                                            velocities=[star_velocity])
planet_inclination = planet_keplerian_orbital_elements[3]
planet_spin = posidonius.calculate_spin(planet_angular_frequency, 
                                    planet_inclination, planet_obliquity)
planet = posidonius.Particle(planet_mass, planet_radius, planet_radius_of_gyration, planet_position, planet_velocity, planet_spin)
  • Alternatively, in case we want the planet to have the pseudo-syncrhonization period corresponding to its orbit since the beginning of the simulation, we would compute the planet angular frequency by following these steps (an replacing the adequate code presented before): 
    # Pseudo-synchronization period
planet_keplerian_orbital_elements = posidonius.calculate_keplerian_orbital_elements(
                                        planet_mass, planet_position, planet_velocity, 
                                        masses=[star_mass], 
                                        positions=[star_position], 
                                        velocities=[star_velocity])
planet_semi_major_axis = planet_keplerian_orbital_elements[0]
planet_eccentricity = planet_keplerian_orbital_elements[2]
planet_semi_major_axis = a
planet_eccentricity = e
planet_pseudo_synchronization_period = posidonius.calculate_pseudo_synchronization_period(
                                                        planet_semi_major_axis, 
                                                        planet_eccentricity, 
                                                        star_mass, planet_mass)
planet_angular_frequency = posidonius.constants.TWO_PI/(planet_pseudo_synchronization_period/24.) # days^-1
  • We can now define the effects that we will enable for the planet:
    • Tidal forces
      • As mentioned above, the tidal model only allow interactions between one body (usually the most massive one, the host star in this example) and the rest of bodies. Therefore, to enable this force the simulation needs to have one body with this effect set to CentralBody and the rest with OrbitingBody (although not necessarily all of them). We can use Disabled when we do not want that this effect is taken into account. In this example, we already defined the effect for the host star, and we will enable it for the planet too so we will use OrbitingBody and these two parameters:
        • Dissipation factor: fraction of the orbital energy that is dissipated as heat in either the surface or interior of a planetary body.
        • Love number of degree 2 (i.e., k2): dimensionless parameters that measure the rigidity of a planetary body and the susceptibility of its shape to change in response to a tidal potential. 
          k2pdelta = 2.465278e-3 # Terrestrial planets (no gas)
planet_tides_parameters = {
"dissipation_factor_scale": 1.0,
"dissipation_factor": 2. * posidonius.constants.K2 * k2pdelta/(3. * np.power(planet_radius, 5)),
"love_number": 0.305,
}
#planet_tides = posidonius.effects.tides.CentralBody(planet_tides_parameters)
planet_tides = posidonius.effects.tides.OrbitingBody(planet_tides_parameters)
#planet_tides = posidonius.effects.tides.Disabled()
    • Rotational flattening effects
      • This effects is defined very similarly to the tidal forces (i.e., CentralBody, OrbitingBody or Disabled), thus the same rules apply. In this example we enabled it for the host start (CentralBody), and we will also enable it for this planet (OrbitingBody), thus we need to define the "fluid" love number, which is like the tidal one but for rotational flattening effects: 
        planet_fluid_love_number = 0.307
planet_rotational_flattening_parameters = {"love_number": planet_fluid_love_number}
#planet_rotational_flattening = posidonius.effects.rotational_flattening.CentralBody(planet_rotational_flattening_parameters)
planet_rotational_flattening = posidonius.effects.rotational_flattening.OrbitingBody(planet_rotational_flattening_parameters)
#planet_rotational_flattening = posidonius.effects.rotational_flattening.Disabled()
    • General relativity corrections
      • Again, this effects is defined very similarly to the tidal forces (i.e., CentralBody, OrbitingBody or Disabled), thus the same rules apply. There are three prescriptions available and we already selected Kidder1995 when we enabled the effect for the host star (CentralBody) so now it is only needed to enable it (OrbitingBody) for the planet: 
        #planet_general_relativity = posidonius.effects.general_relativity.CentralBody("Kidder1995")
#planet_general_relativity = posidonius.effects.general_relativity.CentralBody("Anderson1975")
#planet_general_relativity = posidonius.effects.general_relativity.CentralBody("Newhall1983")
planet_general_relativity = posidonius.effects.general_relativity.OrbitingBody()
#planet_general_relativity = posidonius.effects.general_relativity.Disabled()
    • Wind:
      • This simulation we will not include wind effects (Disabled): 
        #planet_wind = posidonius.effects.wind.Interaction({
## Solar wind parametrisation (Bouvier 1997)
#"k_factor": 4.0e-18, # K_wind = 1.6d47 cgs, which is in Msun.AU2.day
#"rotation_saturation": 1.7592918860102842, # 14. * TWO_PI/25.0, in units of the spin of the Sun today
#})
planet_wind = posidonius.effects.wind.Disabled()
    • Disk
      • A disk can be defined around one of the bodies in the simulation (CentralBody), so that the rest of bodies can feel the disk force (OrbitingBody) although some can be disabled too (Disabled). In this case, we will not enable a disk: 
        #disk_surface_density_normalization_gcm = 1000. # g.cm^-2
#disk_surface_density_normalization_SI = disk_surface_density_normalization_gcm * 1.0e-3 * 1.0e4 # kg.m^-2
#disk_properties = {
#'inner_edge_distance': 0.01,  # AU
#'outer_edge_distance': 100.0, # AU
#'lifetime': 1.0e5 * 365.25e0, # days
#'alpha': 1.0e-2,
#'surface_density_normalization': disk_surface_density_normalization_SI * (1.0/posidonius.constants.M_SUN) * posidonius.constants.AU**2, # Msun.AU^-2
#'mean_molecular_weight': 2.4,
#}
#planet_disk = posidonius.effects.disk.CentralBody(disk_properties)
#planet_disk = posidonius.effects.disk.OrbitingBody()
planet_disk = posidonius.effects.disk.Disabled()
    • Body evolution
      • The radius and other parameters (depending on the model) can vary with time for each body in a simulation. Posidonius includes several models for different masses and from different authors, although we will not enable any for this body: 
        #planet_evolution = posidonius.GalletBolmont2017(planet_mass) # mass = 0.30 .. 1.40
#planet_evolution = posidonius.BolmontMathis2016(planet_mass) # mass = 0.40 .. 1.40
#planet_evolution = posidonius.Baraffe2015(planet_mass) # mass = 0.01 .. 1.40
#planet_evolution = posidonius.Leconte2011(planet_mass) # mass = 0.01 .. 0.08
#planet_evolution = posidonius.Baraffe1998(planet_mass) # Sun (mass = 1.0) or M-Dwarf (mass = 0.1)
#planet_evolution = posidonius.LeconteChabrier2013(False) # Jupiter without dissipation of dynamical tides
#planet_evolution = posidonius.LeconteChabrier2013(True) # Jupiter with dissipation of dynamical tides
planet_evolution = posidonius.NonEvolving()
  • Once the effects are defined, we can assign them to the previously created planet: 
    planet.set_tides(planet_tides)
planet.set_rotational_flattening(planet_rotational_flattening)
planet.set_general_relativity(planet_general_relativity)
planet.set_wind(planet_wind)
planet.set_disk(planet_disk)
planet.set_evolution(planet_evolution)
  • Finally, we are ready to add the star to the universe: 
    universe.add_particle(planet)
  • Finally, the universe with the host star and the earth-like planet should be saved in a JSON file and in the process we will select the integrator we want to use:
    • WHFAST integrator (Rein & Tamayo 2015), which can use three kind of coordinates (i.e., DemocraticHeliocentric, WHDS, WHDS). This is a symplectic integrator that assumes a massive body that gravitationally dominates the system, thus it cannot be used when there are multiple stars. It can be used for the example here defined since it consists on one star and one planet.
    • IAS15 integrator (Rein & Spiegel, 2015). This integrator can be used for any case, no assumption is done about one body being the most massive. It could be used with the current example, but it is slower than WHFAST.
    • LeapFrog. This is a toy integrator for testing and educational purposes only. 
      output_filename = "target/custom.json"
whfast_alternative_coordinates="DemocraticHeliocentric"
#whfast_alternative_coordinates="WHDS"
#whfast_alternative_coordinates="Jacobi"
universe.write(output_filename, integrator="WHFast", whfast_alternative_coordinates=whfast_alternative_coordinates)
#universe.write(output_filename, integrator="IAS15")
#universe.write(output_filename, integrator="LeapFrog")

Once the script is saved (e.g., as custom_case.py), it can be easily run by executing:

mkdir -p target/
python custom_case.py

And it will create the JSON file in target/custom.json. The examples included in Posidonius are slightly more complete, since the output filename can be specified as an argument (it falls out of the scope of this document to explain that part of the python code):

python cases/Bolmont_et_al_2015/case3.py target/case3.json
python cases/Bolmont_et_al_2015/case4.py target/case4.json
python cases/Bolmont_et_al_2015/case7.py target/case7.json
python cases/example.py target/example.json

An important aspect of N-body simulations (and more in general, for any scientific study) is ensuring that your computations can be reproduced and reverified (Rein & Tamayo 2017). Thanks to the Python script and the JSON output, any study can be completely reproduced by indicating what Posidonius version was used and the original script or JSON file.