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] : 29784.7 
──────────────────────────

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.2568983840419

meanmotion(orb) # Mean motion (radians/yr)

6.283066640494721

periastron(orb) # Epoch of periastron passage (MJD)

0.0

semiamplitude(orb) # radial velocity semi-amplitude (m/s)

29784.691829676925

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

plot(orb)

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.2568983840419, 6.283066640494721, 1.0, 0.0, 29784.691829676925), 0.1, 0.1, 0.9950041652780258, 5.813244087623438)

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] : 29784.7 
──────────────────────────

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.2568983840419, 6.283066640494721, 1.0, 0.0, 29784.691829676925), 1.5653910042026546, 1.5653910042026546, 0.0054052962706005155, 59045.0)

We can query specifics at this solution:

trueanom(sol) # true anomaly (radians)

1.5653910042026546

eccanom(sol) # eccentric anomaly (radians)

1.5653910042026546

plot(sol) # defaults to kind=:radvel for RadialVelocityOrbit

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

ω 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

plot(orb, kind=(:x,:y,:z))

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.677033584602116, -2.257801479678382, -0.1761103009311325, 0.9843704393702332, -0.5364115213389841, 1.3805347403121708, 60676.0)

eccanom(sol) # eccentric anomaly (radians)

-2.257801479678382

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

PlanetOrbits.posx(sol)

1.1470772815211459

PlanetOrbits.posy(sol)

-0.7558070007079656

PlanetOrbits.posy(sol)

-0.7558070007079656

plot(sol)

plot(sol, kind=:x)

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

radvel(sol)

10588.51498487476

plot(sol, kind=:radvel)

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)

posangle(sol) # position angle offset from barycentre (milliarcseconds)

projectedseparation(sol) # separation from barycentre (milliarcseconds)

150.95795341303815

raoff(sol) # right ascension offset from barycentre (milliarcseconds)

141.8181317990455

decoff(sol) # declination offset from barycentre (milliarcseconds)

-51.72930689349668

raoff(sol) # right ascension offset from barycentre (milliarcseconds)

141.8181317990455

decoff(sol) # declination offset from barycentre (milliarcseconds)

-51.72930689349668

pmra(sol) # instantaneous right ascension velocity from barycentre (milliarcseconds/year)

-136.78742482708498

pmdec(sol) # instantaneous declination velocity from barycentre (milliarcseconds/year)

-242.62899418753864

accra(sol) # instantaneous right ascension acceleration from barycentre (milliarcseconds/year^2)

-1482.3927234208916

accdec(sol) # instantaneous declination acceleration from barycentre (milliarcseconds/year^2)

540.7146967299267

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"))

21878.311530501967

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)

10489.02128926145

# Primary radial velocity
radvel(sol, 0.1) # 0.1 solar mass secondary

-1048.9021289261452

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

PlanetOrbits.posx(sol), PlanetOrbits.posx(sol, 0.1)

(-0.03746496348200512, 0.003746496348200512)

radvel(sol), radvel(sol, 0.1)

(10489.02128926145, -1048.9021289261452)

raoff(sol), raoff(sol, 0.1)

(-3.746496348200512, 0.3746496348200512)

accra(sol), accra(sol, 0.1)

(147.90016062911553, -14.790016062911555)

projectedseparation(sol), projectedseparation(sol, 0.1)

(94.56050254092096, 9.456050254092096)

posangle(sol), posangle(sol, 0.1)

(-3.1019621829495057, 0.03963047064028754)