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
τ         = 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)

58849.0

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

29785.890525825977

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, 58849.0, 365.2422, 6.283185307179586, 1.0, 0.0, 29785.890525825977), 0.1, 0.1, 0.9950041652780258, 58489.57081015557)

When constructing an orbit, the location of the planet along its orbit can be specified by τ. This is a unitless value between 0 and 1 that represents the fraction of the orbit completed at reference epoch tref which by default is the MJD 58849.0, or 2020-01-01.

orb = orbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    τ=0.2,
);

RadialVelocityOrbit{Float64}
─────────────────────────
a   [au ] = 1.0
e         = 0.0
ω   [°  ] = 0.0
τ         = 0.2
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, 0.2, 1.0, 58849.0, 365.2422, 6.283185307179586, 1.0, 0.0, 29785.890525825977), 2.1151100154549987, 2.1151100154549987, -0.5178310849003872, 59045.0)

We can query specifics at this solution:

trueanom(sol) # true anomaly (radians)

2.1151100154549987

eccanom(sol) # eccentric anomaly (radians)

2.1151100154549987

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)
    τ=0.2,
    # 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)
    τ=0.2,
    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
τ         = 0.2
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, 0.2, 1.0), -2.402898315438049, -1.823941726564961, 0.09694916543950415, 0.9952893344754498, -0.4436089926163764, 1.1502702083729568, 60676.0)

eccanom(sol) # eccentric anomaly (radians)

-1.823941726564961

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

PlanetOrbits.posx(sol)

0.7923981058899824

PlanetOrbits.posy(sol)

-0.8314215952320527

PlanetOrbits.posy(sol)

-0.8314215952320527

plot(sol)

plot(sol, kind=:x)

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

radvel(sol)

10818.488767445728

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)
    τ=0.2,
    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
τ         = 0.2
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)

143.5680289082937

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

125.7543106934933

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

-69.26494255117096

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

125.7543106934933

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

-69.26494255117096

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

-261.77533787758915

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

-185.71163886586308

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

-1620.351401869363

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

892.4827001496805

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)
    τ=0.2,
    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)
) # Not type stable

Use:

orb = VisualOrbit(
    a=1.0, # semi major axis (AU)
    M=1.0, # primary mass (solar masses)
    τ=0.2,
    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

This will also suppress the log message.

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

29783.14689439822

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.

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)

14278.801241174964

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

-1427.8801241174965

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

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

(0.5888745445951494, -0.058887454459514946)

radvel(sol), radvel(sol, 0.1)

(14278.801241174964, -1427.8801241174965)

raoff(sol), raoff(sol, 0.1)

(58.88745445951494, -5.8887454459514945)

accra(sol), accra(sol, 0.1)

(-2324.7835188103113, 232.47835188103113)

projectedseparation(sol), projectedseparation(sol, 0.1)

(99.99788289642939, 9.99978828964294)

posangle(sol), posangle(sol, 0.1)

(2.511911598223864, -0.6296810553659294)