Syzygy.multibodysystem Function

multibodysystem(masses::Vector{<:Quantity}; <kwargs>)

Compute and return a multibody system with given structures, orbital elements, and specified hierarchy.

Structure and element arguments can be given as vectors of length # of bodies or # of binaries respectively, or as numbers, in which case the same value will be used for each body and/or binary.

...

Binary elements arguments

• a = 1.0u"AU": semi-major axis.

• e = 0.1: eccentricity.

• ω = 0.0u"rad": argument of periapsis.

• i = 0.0u"rad": (mutual) inclination. The inclination of the first binary is with respect to the xy-plane.

• Ω = 0.0u"rad": longitude of the ascending node.

• ν = (π)u"rad": true anomaly, with ν = 0u"rad" corresponding to periapsis.

Particle structure arguments

• R = 1.0u"Rsun": radius of each particle.

• S = 0.0u"1/yr": spin magnitude of each particle. If given as a negative number, the spin will be calculated using Syzygy.stellar_spin.

• L = 1.0u"Lsun": luminosity of each particle.

• stellar_types = 1: type of each particle. See Syzygy.stellar_types for all the supported types.

• R_core = 0.0u"Rsun": stellar core radius. Only used if tidal potential is included.

• m_core = 0.0u"Msun": stellar core mass. Only used if tidal potential is included.

• R_env = 0.0u"Rsun": stellar envelope radius. Only used if tidal potential is included.

• m_env = 0.0u"Msun": stellar envelope mass. Only used if tidal potential is included.

Other arguments

• hierarchy: specification of the hierarchy structure. First element is total number of bodies, followed by number of binaries on each level. If given a number, system is assumed to be hierarchichal.

• time = 0.0u"s": time of the system.

• verbose = false

...

Examples

julia> binary = multibodysystem([1.0, 1.0]u"Msun");
julia> binary = multibodysystem([1.0, 1.3]u"Msun", R=[1.0, 0.13]u"Rsun", a=1.0u"AU", e=0.4);
julia> triple = multibodysystem([1.0, 1.0, 1.0]u"Msun", a=[0.1, 1.0]u"AU", e=[0.1, 0.7], i=[π/2, 0.0]u"rad");
julia> quadruple = multibodysystem([1.0, 1.0, 1.0, 1.0]u"Msun");
julia> quadruple = multibodysystem([1.0, 1.0, 1.0, 1.0]u"Msun", hierarchy=[4, 2, 1]);

source

multibodysystem(masses, positions, velocities; kwargs...)

Create a NonHierarchicalSystem from masses, positions, and velocities.

Examples

julia> masses = rand(3)u"kg"
julia> positions = [rand(3)u"m" for i = 1:3]
julia> velocities = [rand(3)u"m/s" for i = 1:3]
julia> multibodysystem(masses, positions, velocities);
julia> multibodysystem(masses, positions, velocities, radii=ones(3)u"m") # you can also supply any of the stellar structure arguments

source

multibodysystem(system::HierarchicalMultiple; new_params...)

Remake a given system with new parameters.

source

multibodysystem(sol::MultiBodySolution, time)

Remake a system at a specific point in time from a solution object.

Examples

julia> triple = multibodysystem([3, 2, 1]u"Msun")
julia> result = simulate(triple, t_sim=10u"yr", save_everystep=true)
julia> solution = to_solution(result)
julia> triple_new = multibodystem(solution, 5.0u"yr")

source

Syzygy.simulation Function

simulation(system::MultiBodyInitialConditions; <simulation kwargs>)

Setup a simulation with a given system and simulation arguments. Returns a Syzygy.MultiBodySimulation object.

...

Arguments

• alg = DPRKN8(): ODE solver to use for simulation. See DiffEq website for list of all solvers.

• t_sim = 1.0: simulation time. Can either be a $Quantity$ of time, or a $Number$, in which case $t_{s}im$ is a multiple of the outermost binary in the system.

• npoints = 0: number of datapoints to save. The temporal locations of the snapshots are then $range(tspan...,length=npoints)$.

• dt = 1/10: time step for the ODE solver in multiples of the innermost binary. Only relevant if solver is symplectic.

• potential = PureGravitationalPotential(): potential to use in simulation. Can either be a single object of abstract type $MultiBodyPotential$, or a $VectorMultiBodyPotential$, in which case the total acceleration will be the sum of all acceleration functions for each potential.

• callbacks::Vector = ["collision"]: callbacks to use in the integration. Can be used to define stopping conditions or other checks. Value should be an array containing a mix of strings, referencing a pre-defined callback, or a custom callback from the DifferentialEquations.jl ecosystem. See Event Handling and Callback Functions for more.

• showprogress = false: whether to display the progress of the simulation.

• max_cpu_time = Inf: maximum time (in seconds) allowed for the simulation to run. Will terminate if it runs for longer than max_cpu_time.

• stellar_evolution = false: whether the stars will be evolved, meaning their structural properties will change. If true, the structural parameters of the stars will be types MVector, otherwise SVector.

...

The function also accepts all keyword arguments supported by the CommonSolve.solve interface from DifferentialEquations.jl. See 'Common Solver Options' for more information.

Example

julia> triple = multibodysystem([1.0, 1.0, 1.0]u"Msun")
julia> sim = simulation(triple, t_sim=10) # simulate for 10 outer orbits
julia> sim = simulation(triple, t_sim=500u"kyr") # simulate for 500 kyr
julia> sim = simulation(triple, t_sim=1, npoints=1000) # save 1000 datapoints
julia> sim = simulation(triple, t_sim=1, alg=Syzygy.McAte5(), dt=1.0) # use the McAte5 symplectic integrator with a timestep of 1*innermost period.

source

Syzygy.simulate Function

simulate(simulation::MultiBodySimulation)

Simulate the given simulation setup.

source

simulate(system::MultiBodyInitialConditions; args...)

Allows call to simulate directly without setting up simulation first.

source

Syzygy.CollisionCB Type

CollisionCB(check_every=1, grav_rad_multiple=10_000)

A callback for checking if collision has occured in system. Collision depends on the stellar types of the objects, and is triggered when the distance between two objects is smaller than some value.

• If the two objects are stars or planets, the callback checks for overlapping radii

• If one of the objects is a compact object (CO) and the other is a star, the tidal disruption radius is used.

• If both objects are COs, the code checks if the distance is smaller than grav_rad_multiple times their mutual gravitational radius.

• If one object is a planet and one is a star, the roche limit is used.

source

Syzygy.EscapeCB Type

EscapeCB(check_every=1, max_a_factor=100)

A callback for checking if a body in a triple system has escaped (become unbound). The criterium for escape is based on three checks:

1. The distance between the body and the center of mass of the remaining binary is
larger than the semi-major axis of the binary multiplied by `max_a_factor`.
2. The body is moving away from the center of mass of the binary.
3. The body has a positive total energy.

If criterium 1 and 2 are satisfied, but not 3, then the code checks if the body is more than 1 parsec away from the binary, in which case it is classified as a Drifter.

Arguments

• check_every: How often the solver should check for escape. Default is 1 (every time step).

• max_a_factor: See above.

source

Syzygy.PureGravitationalPotential Type

PureGravitationalPotential()

Newtonian gravitational potential. Corresponds to the acceleration function Syzygy.pure_gravitational_acceleration!.

source

Syzygy.DynamicalTidalPotential Type

DynamicalTidalPotential(n_t, polytropic_index)

Set up the dynamical tidal potential for a system as defined by Samsing, Leigh & Trani 2018. Corresponds to the acceleration function Syzygy.dynamical_tidal_acceleration!.

Arguments

• n_t: tidal force power index. Can be either 4 or 10.

• polytropic_index: vector of polytropic indices of each body in the system

source

Syzygy.TimeDependentEquilibriumTidalPotential Type

TimeDependentEquilibriumTidalPotential()

Equilibrium tidal potential, with the assumption that the stars are evolving over time, which means the internal structure changes. The envelope structure is calculated throughout the simulation, rather than being fixed from the start.

Corresponds to the acceleration function Syzygy.pure_gravitational_acceleration!.

source

Syzygy.EquilibriumTidalPotential Type

EquilibriumTidalPotential(system; Z=0.02, lb_multiplier=1.1, ub_multiplier=1.1,
                                  supplied_apsidal_motion_constants=nothing, 
                                  supplied_rotational_angular_velocities=nothing)

Set up the equilibrium tidal potential corresponding to the equilibrium tidal force prescription from Hut 1981. The apsidal motion constant k for each stellar body is acquired by interpolating a grid of MESA models, using data from Claret 2023. If the structural properties of one or more of the stellar objects fall outside the grid, the interpolator can extrapolate to the given boundaries, which can be adapted using the keyword arguments lb_multiplier for the lower bounds, and ub_multiplier for the upper bounds.

Alternatively, the apsidal motion constants and rotational angular velocities can be supplied using the keyword arguments supplied_apsidal_motion_constants and supplied_rotational_angular_velocities.

Corresponds to the acceleration function Syzygy.equilibrium_tidal_acceleration!.

Arguments

• system: an instance of a HierarchicalMultiple or NonHierarchichalSystem type.

Keyword arguments

• Z: metallicity

• lb_multiplier: multiplier for the lower bounds of the apsidal motion constant.

• ub_multiplier: multiplier for the upper bounds of the apsidal motion constant.

• supplied_apsidal_motion_constants

• supplied_rotational_angular_velocities

source

Missing docstring.

Missing docstring for Syzygy.PN1Potential. Check Documenter's build log for details.

Missing docstring.

Missing docstring for Syzygy.PN2Potential. Check Documenter's build log for details.

Missing docstring.

Missing docstring for Syzygy.PN2p5Potential. Check Documenter's build log for details.

Missing docstring.

Missing docstring for Syzygy.PNPotential. Check Documenter's build log for details.

Syzygy.semi_major_axis Function

semi_major_axis(d, v², M)

Return the semi-major axis of a binary system with total mass 'M'. Here, 'd' denotes the relative separation, while 'v²' is the relative velocity magnitude squared.

$$\begin{gathered} a= \\ fracGMd2GM-dv² \end{gathered}$$source

Syzygy.orbital_period Function

orbital_period(a, M)

Get the orbital period of binary with semi-major axis a and total mass M.

source

Syzygy.eccentricity_vector Function

eccentricity_vector(r, v, d, M)

Eccentricity vector of binary orbit with relative position vector r, separation d, relative velocity v, and total mass M.

$$\begin{gathered} vece= \\ frac \\ vecv \\ times( \\ vecr \\ times \\ vecv) \\ mu- \\ frac \\ vecrd \end{gathered}$$source

Syzygy.eccentricity Function

eccentricity(r, v, d, M)

Eccentricity of binary orbit with relative position vector r, relative velocity v, separation d, and total mass M, where M ≡ m₁ + m₂. The eccentricity is calculated as the norm of the eccentricity vector

$$\begin{gathered} e=| \\ vece|= \\ frac \\ vecv \\ times( \\ vecr \\ times \\ vecv) \\ mu- \\ frac \\ vecrd \end{gathered}$$source

Syzygy.mutual_inclination Function

mutual_inclination(h₁, h₂)

Return the mutual inclination angle between two binaries with orbital angular momenta h₁ and h₂.

source

Syzygy.longitude_of_ascending_node Function

longitude_of_ascending_node(h)

Return the longitude of ascending node of a body in an orbit with angular momentum h.

source

Syzygy.argument_of_periapsis Function

argument_of_periapsis(r, v, h, m, G)

Return the argument of periapsis of a body in an orbit with relative position r, velocity v, angular momentum h and mass m.

source

Syzygy.true_anomaly Function

true_anomaly(r, v, h, M, G)

Return the true anomaly of a body in an orbit with relative position r, velocity r, angular momentum h and mass m.

source

Syzygy.binary_elements Function

binary_elements(positions, velocities, masses)

Calculate binary properties of bodies with given positions and velocities. Assumes the bodies are gravitationally bound. Returns an instance of OrbitalElements.

source

Syzygy.is_binary Function

is_binary(r1, r2, v1, v2, m1, m1)

Check if two bodies with the given positions (r1, r2), velocities (v1, v2), and masses (m1, m2) are gravitationally bound and form a binary.

source

Syzygy.centre_of_mass Function

centre_of_mass(positions, masses)

Return the centre of mass position given a collection of positions and masses.

positions can either be a Vector of Vectors, or a Matrix in which each column represents the position of each body.

Examples

julia> r = [rand(3) for i =1:3]
julia> m = rand(3)

julia> com = Syzygy.centre_of_mass(r, m)

source

centre_of_mass(sol::MultiBodySolution, 
               bodies=eachindex(sol.ic.particles); 
               tspan=extrema(sol.t))

Return the centre of mass (COM) position given a MultiBodySolution object.

By default, the COM position calculated for each body for each time step in sol is returned. Specific bodies can be specified with their indices, and a specific time span can be set with a tuple of (t_start, t_end).

Examples

julia> com = Syzygy.centre_of_mass(sol)
julia> com12 = Syzygy.centre_of_mass(sol, [1, 2])
julia> # just get the COM position for the first two time steps:
julia> com_tspan = Syzygy.centre_of_mass(sol, tspan=(sol.t[1], sol.t[2]))

source

Syzygy.centre_of_mass_velocity Function

centre_of_mass_velocity(velocities, masses)

Return the center of mass velocity.

source

Syzygy.potential_energy Function

potential_energy(positions, masses)

Return the total potential energy of bodies with positions and masses.

positions can either be a Vector of Vectors, or a Matrix in which each column represents the position of each body.

Examples

julia> r = [rand(3) for i =1:3]
julia> m = rand(3)

julia> com = Syzygy.potential_energy(r, m)

source

potential_energy(sol::MultiBodySolution, 
                 bodies=keys(sol.ic.particles); 
                 tspan=extrema(sol.t))

Return the total potential energy for the given bodies and tspan.

source

Syzygy.kinetic_energy Function

kinetic_energy(velocity::AbstractVector, mass::Number)

Return the kinetic energy of a body with given velocity vector and mass.

source

kinetic_energy(velocity::AbstractVector{<:AbstractVector}, mass::Number)

Total kinetic energy of bodies with velocity vectors and masses, where velocities[1] is a Vector with the velocity components of particle 1.

source

kinetic_energy(velocity::AbstractMatrix, mass::AbstractVector)

velocity can also be a Matrix where each column represent the velocity vector of each body.

source

kinetic_energy(sol::MultiBodySolution, 
               bodies=keys(sol.ic.particles); 
               tspan=extrema(sol.t))

Return the total kinetic energy for the given bodies and tspan from a simulation solution.

source

Syzygy.total_energy Function

total_energy(positions, velocities, masses)

Total energy (kinetic + potential) of bodies with given positions, velocities, and masses.

source

Syzygy.specific_orbital_energy Function

specific_orbital_energy(r, v², μ)

Return the specific orbital energy of two orbiting bodies with relative position r, relative velocity squared v², and reduced mass μ

The specific orbital energy is sum of the mutual potential and kinetic energy divided by the reduced mass.

source

Syzygy.reduced_mass Function

reduced_mass(m₁, m₂)

Return the reduced mass given the two masses m₁ and m₂.

The reduced mass is defined as

$$\frac{m_1{m}_2}{m_1+m_2}$$source

Syzygy.gravitational_radius Function

gravitational_radius(mass::Unitful.Mass)

Return the gravitational radius of an object with the given mass.

The gravitational radius is defined as GM/c².

source

Syzygy.schwarzschild_radius Function

schwarzschild_radius(mass::Unitful.Mass)

Return the Schwarzschild radius of an object with the given mass.

The Schwarzschild radius is defined as 2GM/c².

source

Syzygy.roche_radius Function

roche_radius(a, m₁, m₂)

Return the volume-equivalent Eggleton Roche lobe radius of two bodies with masses m₁ and m₂, and separation a, where m₁ > m₂.

The volume-equivalent Eggleton Roche radius is defined as

$$R_L=a\frac{0.49q^{2/3}}{0.6q^{2/3}+\ln (1+q^{1/3})}$$source

roche_radius(a, q)

The mass ratio q can also be given directly, where q = m₁/m₂.

source

Syzygy.stellar_rotational_frequency Function

stellar_rotational_frequency(m::Unitful.Mass, R::Unitful.Length)

Return the stellar rotation frequency of a star with mass 'm [M⊙]' and radius 'R [R⊙]', as described by Hurley, Pols, & Tout 2000, eq 107-108.

source

stellar_rotational_frequency(m::Real, R::Real)

If m and R are given without units, they must be [M⊙] and [R⊙] respectively.

source

Syzygy.stability_criterion_ma01 Function

stability_criterion_ma01(m₂, m₂, m₃, i, eout)

Return the critical semi-major axis ratio (a_out/a_in)_crit as defined in Mardling & Aarset 1999, given a triple with masses m₁, m₂, m₃, mutual inclination i, and outer eccentricity eout.

source

stability_criterion_ma01(triple::HierarchicalMultiple)

The stability criterion can also be calculated for a triple system directly.

source

Syzygy.is_unstable Function

is_unstable(triple::HierarchicalMultiple)

Check if a given triple is unstable, i.e. whether aout/ain < (aout/ain)_crit.

source

Syzygy.octupole_parameter Function

octupole_parameter(triple::HierarchicalMultiple)

Return the octupole parameter ϵₒ for a given triple.

source