Introduction

This package is structured around a representation of an orbit (PlanetOrbits.AbstractOrbit, and a representation of a "solved" orbit (PlanetOrbits.AbstractOrbitSolution).

You start by creating an orbit with known information, e.g. the semi-major axis and eccentricity. You can then query information from this orbit, like its orbital period, mean motion, or periastron (closest approach). Then, you can "solve" the orbit one more times for a given time, eccentric anomaly, or true anomaly.

Let's see how this works.

using PlanetOrbits, Plots
orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
)
RadialVelocityOrbit{Float64}
─────────────────────────
a   [au ] = 1.0
e         = 0.0
ω   [°  ] = 0.0
tp         = 0.0
M   [M⊙ ] = 1.0 
──────────────────────────
period        [yrs ] : 1.0 
mean motion   [°/yr] : 360.0 
semiamplitude₂ [m/s] : 29785.9 
──────────────────────────

The orbit function accepts many combinations of orbital parameters and returns a subtype of PlanetOrbits.AbstractOrbit.

We can now query some basic properties about the orbit:

period(orb) # orbital period (days)
365.2422
meanmotion(orb) # Mean motion (radians/yr)
6.283185307179586
periastron(orb) # Epoch of periastron passage (MJD)
0.0
semiamplitude(orb) # radial velocity semi-amplitude (m/s)
29785.890525825977

We can plot the orbit (more on this in Plotting):

plot(orb)
Example block output

And we can solve the orbit for a given orbital location

sol = orbitsolve_ν(orb, 0.1)        # true anomaly (radians)
sol = orbitsolve_eccanom(orb, 0.1)  # eccentric anomaly (radians)
sol = orbitsolve_meananom(orb, 0.1) # mean anomaly (radians)
OrbitSolutionRadialVelocity{Float64, RadialVelocityOrbit{Float64}}(RadialVelocityOrbit{Float64}(1.0, 0.0, 0.0, 0.0, 1.0, 365.2422, 6.283185307179586, 1.0, 0.0, 29785.890525825977), 0.1, 0.1, 0.9950041652780258, 5.813010155575866)

When constructing an orbit, the location of the planet along its orbit can be specified by tp. This is the (or a) time the planet made its closest approach to the star.

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    tp=mjd("2020-04-15"),
);
RadialVelocityOrbit{Float64}
─────────────────────────
a   [au ] = 1.0
e         = 0.0
ω   [°  ] = 0.0
tp         = 59000.0
M   [M⊙ ] = 1.0 
──────────────────────────
period        [yrs ] : 1.0 
mean motion   [°/yr] : 360.0 
semiamplitude₂ [m/s] : 29785.9 
──────────────────────────

We can now meaningfully solve the location of the planet at a specific time:

t = mjd("2020-07-15") # date as in modified julian days.
sol = orbitsolve(orb, t) # can optionally pass `tref=...`
OrbitSolutionRadialVelocity{Float64, RadialVelocityOrbit{Float64}}(RadialVelocityOrbit{Float64}(1.0, 0.0, 0.0, 58954.0, 1.0, 365.2422, 6.283185307179586, 1.0, 0.0, 29785.890525825977), 1.5654539999850572, 1.5654539999850572, 0.005342301397802048, 59045.0)

We can query specifics at this solution:

trueanom(sol) # true anomaly (radians)
1.5654539999850572
eccanom(sol) # eccentric anomaly (radians)
1.5654539999850572
plot(sol) # defaults to kind=:radvel for RadialVelocityOrbit
Example block output

Notice that we now see a marker at the location found by orbitsolve.

We can create an orbit with some eccentricity. If not specified, eccentricity and the argument or periapsis default to 0 for any orbit type.

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    tp=mjd("2020-04-15"),
    # New:
    e=0.6, # eccentricity
    ω=2.5, # argument of periapsis (radians)
)
plot(orb) # defaults to kind=:radvel for RadialVelocityOrbit
Example block output
ω convention

The convention used in this package is that ω, the argument of periapsis, refers to the secondary body. This is in contrast to the typical standard adopted in the radial velocity literature where ω refers to the primary. You can convert by adding or subtracting 180°.

Since we only provided very minimal information to the orbit function, we've been receiving a RadialVelocityOrbit. This object contains sufficient information to calculate the above radial velocity plots, orbital period, etc., but not the 3D position in space.

Let's create a new orbit with a specified inclination and longitude of ascending node.

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    tp=mjd("2020-04-15"),
    e=0.6, # eccentricity
    ω=2.5, # argument of periapsis (radians)
    # New:
    i=0.6, # inclination (radians)
    Ω=2.3, # inclination (radians)
)
KepOrbit{Float64}
─────────────────────────
a   [au ] = 1.0
e         = 0.6
i   [°  ] = 34.4
ω   [°  ] = 143.0
Ω   [°  ] = 132.0
tp  [day] = 59000.0
M   [M⊙ ] = 1.0 
period      [yrs ] : 1.0 
mean motion [°/yr] : 360.0 
──────────────────────────

This time we received a full KepOrbit. This has the necessary information to solve the orbit in 2/3D.

plot(orb) # defaults to kind=(:x,:y) for KepOrbit
Example block output
plot(orb, kind=(:x,:y,:z))
Example block output
Cartesian convention

The convention used in this package is that x increases to the left (just like right-ascension), and the z increases away from the observer.

We can solve for a time or location as usual.

sol = orbitsolve(orb, mjd("2025-01"))
PlanetOrbits.OrbitSolutionKep{Float64, KepOrbit{Float64}}(KepOrbit(1.0, 0.6, 0.6, 2.5, 2.3, 59000.0, 1.0), -2.6765330609472615, -2.256937867171756, -0.1756175782017388, 0.9844584634338598, -0.536276904791952, 1.380133977827578, 60676.0)
eccanom(sol) # eccentric anomaly (radians)
-2.256937867171756

We can also query the cartesian position of the planet in AU:

PlanetOrbits.posx(sol)
1.1464609368326173
PlanetOrbits.posy(sol)
-0.7560870613717289
PlanetOrbits.posy(sol)
-0.7560870613717289
plot(sol)
Example block output
plot(sol, kind=:x)
Example block output

We can still of course calculate the radial velocity as well.

radvel(sol)
10590.791651219175
plot(sol, kind=:radvel)
Example block output

Finally, we'll specify the parallax distance to the system. This will allow us to plot orbits with angular units as they would appear in the sky from the Earth.

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    tp=mjd("2020-04-15"),
    e=0.6, # eccentricity
    ω=2.5, # argument of periapsis (radians)
    i=0.6, # inclination (radians)
    Ω=2.3, # inclination (radians)
    # New:
    plx=100.0 # parallax distance (milliarcseconds)
)
KepOrbit{Float64}
─────────────────────────
a   [au ] = 1.0
e         = 0.6
i   [°  ] = 34.4
ω   [°  ] = 143.0
Ω   [°  ] = 132.0
tp  [day] = 59000.0
M   [M⊙ ] = 1.0 
period      [yrs ] : 1.0 
mean motion [°/yr] : 360.0 
──────────────────────────
plx [mas] = 100.0 
distance    [pc  ] : 10.0 
──────────────────────────
sol = orbitsolve(orb, 0.0)
plot(sol)
Example block output
posangle(sol) # position angle offset from barycentre (milliarcseconds)
projectedseparation(sol) # separation from barycentre (milliarcseconds)
151.229956039581
raoff(sol) # right ascension offset from barycentre (milliarcseconds)
142.67541135317006
decoff(sol) # declination offset from barycentre (milliarcseconds)
-50.14206416709735
raoff(sol) # right ascension offset from barycentre (milliarcseconds)
142.67541135317006
decoff(sol) # declination offset from barycentre (milliarcseconds)
-50.14206416709735
pmra(sol) # instantaneous right ascension velocity from barycentre (milliarcseconds/year)
-127.18425300129124
pmdec(sol) # instantaneous declination velocity from barycentre (milliarcseconds/year)
-246.07331820336822
accra(sol) # instantaneous right ascension acceleration from barycentre (milliarcseconds/year^2)
-1475.3673389701755
accdec(sol) # instantaneous declination acceleration from barycentre (milliarcseconds/year^2)
518.505347761442

Performance

The orbit function is a convenience only for interactive use. It is inefficient since it is not type-stable. Instead, one should use one of the orbit constructors directly. For example, instead of

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    tp=mjd("2020-04-15"),
    e=0.6, # eccentricity
    ω=2.5, # argument of periapsis (radians)
    i=0.6, # inclination (radians)
    Ω=2.3, # inclination (radians)
) # Not type stable
KepOrbit{Float64}
─────────────────────────
a   [au ] = 1.0
e         = 0.6
i   [°  ] = 34.4
ω   [°  ] = 143.0
Ω   [°  ] = 132.0
tp  [day] = 59000.0
M   [M⊙ ] = 1.0 
period      [yrs ] : 1.0 
mean motion [°/yr] : 360.0 
──────────────────────────

Use:

orb = KepOrbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    tp=mjd("2020-04-15"),
    e=0.6, # eccentricity
    ω=2.5, # argument of periapsis (radians)
    i=0.6, # inclination (radians)
    Ω=2.3, # inclination (radians)
    plx=100.0 # parallax distance (milliarcseconds)
) # Type stable
KepOrbit{Float64}
─────────────────────────
a   [au ] = 1.0
e         = 0.6
i   [°  ] = 34.4
ω   [°  ] = 143.0
Ω   [°  ] = 132.0
tp  [day] = 59000.0
M   [M⊙ ] = 1.0 
period      [yrs ] : 1.0 
mean motion [°/yr] : 360.0 
──────────────────────────

This will prevent ‪unnecessary‬ allocations in some cases.

Convenience

All functions described above that apply to orbit solutions can be called directly on an orbit along with a time in days:

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses),
)
radvel(orb, mjd("2025-01-01"))
21011.191599677164

If you need to calculate many different properties, e.g. both x and y position at a given time/location, it is more efficient to calculate the orbit solution a single time and query the result as needed.

Host calculations

The above calculations treat the planet as a test particle and calculate their displacement/velocity/etc. compared to the two-particle system's barycentre. If you wish to calculate the same properties for the host object, you can additionally supply the mass of the planet.

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    i=0.5,
    Ω=2.5,
    plx=100.0
)
sol = orbitsolve(orb, mjd("2025-01-01"))
# Secondary radial velocity
radvel(sol)
10073.30184939133
# Primary radial velocity
radvel(sol, 0.1) # 0.1 solar mass secondary
-1007.330184939133

The following show pairs of results for the secondary and the primary:

PlanetOrbits.posx(sol), PlanetOrbits.posx(sol, 0.1)
(-0.07617039392126124, 0.007617039392126124)
radvel(sol), radvel(sol, 0.1)
(10073.30184939133, -1007.330184939133)
raoff(sol), raoff(sol, 0.1)
(-7.617039392126124, 0.7617039392126124)
accra(sol), accra(sol, 0.1)
(300.7086620311962, -30.07086620311962)
projectedseparation(sol), projectedseparation(sol, 0.1)
(94.04913029444482, 9.404913029444483)
posangle(sol), posangle(sol, 0.1)
(-3.06051384700244, 0.08107880658735325)