Basic Usage

Here we'll walk through an example of integrating a 3-body system.

Units

A quick note on the units used throughout the code.

Distance: AU
Time: Days
Mass: Solar Masses
Angles: Radians

Initial Conditions

First, we define the orbital elements of the system. This can be done by creating Elements for each body in the system.

Start with a 1 solar mass star

using NbodyGradient

a = Elements(m = 1.0);

We then define the orbital elements for second body, say an Earth analogue.

b = Elements(
    m = 3e-6,     # Mass [Solar masses]
    t0 = 0.0,     # Initial time of transit [days]
    P = 365.256,  # Period [days]
    ecosω = 0.01, # Eccentricity * cos(Argument of Periastron)
    esinω = 0.0,  # Eccentricity * sin(Argument of Periastron)
    I = π/2,      # Inclination [Radians]
    Ω = 0.0       # Longitude of Ascending Node [Radians]
);

(ecosω can be typed as ecos\omega and then hitting tab)

Next we'll create a Jupiter analogue for the final body. Here the orbital elements are specified for the Keplerian ((a,b),c), or c orbiting the center of mass of a and b. (While this might not need to be stated explicitly here, this convention is useful for more complicated hierarchical systems).

c = Elements(
    m = 9.54e-4,
    P = 4332.59,
    ecosω = 0.05,
    I = π/2
);

Here, we can simply omit any orbital elements which are zero. Unspecified elements are set to zero by default.

Now we need to pass our Elements to ElementsIC.

t0 = 0.0 # Initial time
N = 3    # Number of bodies
ic = ElementsIC(t0,N,a,b,c)

ElementsIC{Float64}
Orbital Elements: 
3×7 Matrix{Float64}:
 1.0          0.0    0.0  0.0   0.0  0.0     0.0
 3.0e-6     365.256  0.0  0.01  0.0  1.5708  0.0
 0.000954  4332.59   0.0  0.05  0.0  1.5708  0.0

Finally, we can pass the initial conditions to State, which converts orbital elements to cartesian coordinates (and calculates the derivatives with respect to the initial conditions).

s = State(ic);

The positions and velocities can be accessed by s.x and s.v, respectively. Each matrix contains the vector component (rows) for a particular body (columns).

s.x

3×3 Matrix{Float64}:
 -4.13804e-19   2.91042e-16   4.32842e-16
  3.03083e-19  -6.09279e-17  -3.17506e-16
  0.00494973   -0.995028     -5.18526

Integration

Now that we have initial conditions, we can construct and run the integrator. First, define an Integrator, specifying the integration scheme, the time step, and final time. We'll use the ahl21! mapping.

h = b.P/30.0 # We want at most 1/20th of the smallest period for a time step
tmax = 5*c.P # Integrate for 5 orbital periods of the outer body
intr = Integrator(ahl21!,h,tmax);

A quick aside on constructors

The types in NbodyGradient.jl have a number of constructors that enable to user to write as terse or as verbose code as they would like. For example, the above Integrator construction could specify each field as Integrator(scheme, h, t0, tmax) (see Integrator for argument definitions), or only the time step and integration time as Integrator(h,tmax). In the latter case, the default value for scheme and t0 are ahl21! and 0, respectively. If you're coming from Agol, Hernandez, and Langford (2021) you'll notice some discrepancy in the examples due to this.

Finally, run the Integrator by passing it the State.

intr(s)
s.x # Show final positions

3×3 Matrix{Float64}:
 -2.75581e-6   0.911197      2.32881e-5
  3.02828e-19  2.41368e-17  -3.17506e-16
  0.00494556   0.394184     -5.18526

This integrates from s.t to s.t+tmax, steping by h. If you'd rather step a certain number of time steps:

N = 1000
intr(s,N)

Re-running simulations

If you want to run the integration from the initial condtions again, you must 'reset' the State. I.e. run s = State(ic). Otherwise, the integration will begin from what ever s.t is currently equal to, and with those coordinates.

Transit Timing

If we wish to compute transit times, we need only to pass a TransitTimingstructure to the Integrator.

s = State(ic) # Reset to the initial conditions
tt = TransitTiming(tmax, ic)
intr(s,tt)

To see the first 5 transit times of our second body about the first body, run:

tt.tt[2,1:5]

5-element Vector{Float64}:
  365.2505305537288
  730.5007334346752
 1095.7505677053489
 1461.001226871598
 1826.2531823012507