Package {orbitr}

Add a physical body to the system

Description

This function introduces a new celestial or physical body into your 'orbit_system'. You must provide a unique identifier and its mass. By default, the body will be placed at the origin (0, 0, 0) with zero initial velocity unless specified.

Usage

add_body(system, id, mass, x = 0, y = 0, z = 0, vx = 0, vy = 0, vz = 0)

Arguments

Value

The updated 'orbit_system' object containing the newly added body.

Examples

my_universe <- create_system() |>
  add_body(id = "Earth", mass = 5.97e24) |>
  add_body(id = "Moon", mass = 7.34e22, x = 3.84e8, vy = 1022)

Add a body using Keplerian orbital elements

Description

A convenience wrapper around [add_body()] that lets you specify an orbit using classical Keplerian elements instead of raw Cartesian state vectors. The elements are converted to position and velocity in the reference frame of the parent body, which must already exist in the system.

Usage

add_body_keplerian(
  system,
  id,
  mass,
  a,
  e = 0,
  i = 0,
  lan = 0,
  arg_pe = 0,
  nu = 0,
  parent
)

Arguments

Value

The updated 'orbit_system' with the new body added.

Keplerian Elements

Six numbers fully describe a Keplerian orbit:

'a' (semi-major axis)

The size of the orbit — half the longest diameter of the ellipse, in meters. Negative for a hyperbolic orbit.

'e' (eccentricity)

The shape of the orbit. 0 is a perfect circle; values between 0 and 1 are ellipses; values above 1 are hyperbolas (unbound flybys).

'i' (inclination)

The tilt of the orbital plane relative to the reference plane, in degrees.

'lan' (longitude of ascending node)

The angle from the reference direction to where the orbit crosses the reference plane going "upward," in degrees. Sometimes written as \Omega.

'arg_pe' (argument of periapsis)

The angle within the orbital plane from the ascending node to the closest-approach point, in degrees. Sometimes written as \omega.

'nu' (true anomaly)

Where the body currently sits along its orbit, measured as an angle from periapsis in degrees. 0 = at periapsis (closest), 180 = at apoapsis (farthest).

Examples

# Earth orbiting the Sun with real orbital elements
system <- create_system() |>
  add_sun() |>
  add_body_keplerian(
    "Earth", mass = mass_earth,
    a = distance_earth_sun, e = 0.0167, i = 0.00005,
    parent = "Sun"
  )

# An interstellar visitor on a hyperbolic orbit (e > 1 needs a < 0):
# 'Oumuamua-like, perihelion 0.255 AU, approaching from 140 degrees
# before perihelion
q <- 0.2553 * distance_earth_sun
e <- 1.2
system <- system |>
  add_body_keplerian(
    "Visitor", mass = 1e10,
    a = -q / (e - 1), e = e, i = 122.7, nu = -140,
    parent = "Sun"
  )

# Mars with its notable eccentricity
system <- system |>
  add_body_keplerian(
    "Mars", mass = mass_mars,
    a = distance_mars_sun, e = 0.0934, i = 1.85,
    lan = 49.6, arg_pe = 286.5, nu = 0,
    parent = "Sun"
  )

Add a known solar system body by name

Description

A convenience wrapper around [add_body_keplerian()] that looks up real orbital elements for well-known solar system bodies. Instead of typing out mass, semi-major axis, eccentricity, and inclination by hand, just give the name and parent:

Usage

add_planet(
  system,
  name,
  parent,
  nu = 0,
  a = NULL,
  e = NULL,
  i = NULL,
  lan = NULL,
  arg_pe = NULL,
  mass = NULL
)

Arguments

Details

“' create_system() |> add_sun() |> add_planet("Earth", parent = "Sun") |> add_planet("Moon", parent = "Earth") “'

Any Keplerian element can be overridden to explore "what if" scenarios while keeping the rest of the real values:

“' # What if Mars had zero eccentricity? add_planet("Mars", parent = "Sun", e = 0) “'

Value

The updated 'orbit_system' with the body added.

Examples

# Build the inner solar system
create_system() |>
  add_sun() |>
  add_planet("Mercury", parent = "Sun") |>
  add_planet("Venus",   parent = "Sun") |>
  add_planet("Earth",   parent = "Sun") |>
  add_planet("Mars",    parent = "Sun") |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year) |>
  plot_orbits()

# Earth-Moon system
create_system() |>
  add_body("Earth", mass = mass_earth) |>
  add_planet("Moon", parent = "Earth") |>
  simulate_system(time_step = seconds_per_hour, duration = seconds_per_day * 28) |>
  plot_orbits()

# What if Jupiter were twice as massive?
create_system() |>
  add_sun() |>
  add_planet("Jupiter", parent = "Sun", mass = mass_jupiter * 2)

Add the Sun to the system

Description

A convenience function that adds the Sun as the central body of a simulation. By default it is placed at the origin with zero velocity, which is the natural choice for a heliocentric reference frame. Position and velocity can be overridden for advanced use cases such as barycentric coordinates.

Usage

add_sun(system, mass = mass_sun, x = 0, y = 0, z = 0, vx = 0, vy = 0, vz = 0)

Arguments

Details

This pairs naturally with [add_planet()]:

“' create_system() |> add_sun() |> add_planet("Earth", parent = "Sun") |> add_planet("Mars", parent = "Sun") “'

Value

The updated 'orbit_system' with the Sun added.

Examples

# Typical usage — Sun at the origin
create_system() |>
  add_sun()


# Full solar system in three lines
create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  add_planet("Mars",  parent = "Sun") |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year) |>
  plot_orbits()

Animate the System Over Time (Smart 2D/3D Dispatch)

Description

Plays the simulation forward as an animation. Bodies move through their orbits frame by frame, optionally leaving a fading wake behind them. This is the animated counterpart to [plot_system()] — a moving snapshot rather than a single frozen one.

Usage

animate_system(
  sim_data,
  fps = 20,
  duration = 10,
  trails = FALSE,
  three_d = NULL
)

Arguments

Details

If any body has non-zero motion in the Z dimension (or 'three_d = TRUE'), [animate_system_3d()] is used; otherwise a 2D 'gganimate' animation is returned.

The 2D path requires the 'gganimate' package, which is in 'Suggests'. Install it with 'install.packages("gganimate")'. Rendering a 2D animation is much slower than a static plot — expect tens of seconds for typical simulations, since every frame is drawn and encoded as a GIF (or MP4).

The 3D path uses ‘plotly'’s built-in 'frame' aesthetic, which produces an interactive HTML widget with a play button and time slider. No GIF encoding is involved, so 3D animations render essentially instantly.

Value

A rendered 'gganimate' animation (2D) or a 'plotly' HTML widget with built-in play/pause controls (3D). The 2D return value can be saved to disk with [gganimate::anim_save()].

Examples

sim <- create_system() |>
  add_sun() |>
  add_body("Earth", mass = mass_earth, x = distance_earth_sun, vy = speed_earth) |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year)

# 2D fading-wake animation (requires gganimate)
anim <- animate_system(sim, fps = 20, duration = 8)
anim

# Save to disk
gganimate::anim_save(file.path(tempdir(), "earth_orbit.gif"), anim)

Animate the System Over Time in Interactive 3D

Description

The 3D counterpart to [animate_system()]. Builds a 'plotly' 3D scene with the bodies as moving markers and an interactive Play / Pause control plus a time slider. Optionally shows the full orbit paths drawn faintly behind.

Usage

animate_system_3d(sim_data, fps = 20, duration = 10, trails = FALSE)

Arguments

Value

A 'plotly' HTML widget with a built-in play button and time slider.

Examples

create_system() |>
  add_body("Earth", mass = mass_earth) |>
  add_body("Moon",  mass = mass_moon,
           x = distance_earth_moon, vy = speed_moon, vz = 100) |>
  simulate_system(time_step = seconds_per_hour, duration = seconds_per_day * 30) |>
  animate_system_3d()

Conserved quantities and how well a simulation conserved them

Description

Joins [get_energy()], [get_momentum()], and [get_angular_momentum()] into one tibble and adds the relative error of each quantity against its value at the first time step. In exact Newtonian gravity all three are constant, so the errors measure the integrator, and their shape tells you what is wrong when something is.

Usage

conserved_quantities(sim_data, G = NULL, softening = NULL)

Arguments

Details

'energy_error' is (E - E_0)/|E_0|. 'momentum_error' is |\mathbf{P} - \mathbf{P}_0| divided by \sum_j m_j |\mathbf{v}_j| at the first step (total momentum is often exactly zero, so it cannot be its own scale). 'angular_momentum_error' is |\mathbf{L} - \mathbf{L}_0| / |\mathbf{L}_0|. An error is 'NA' when its scale is zero.

Value

A tibble with one row per time step and columns 'time', 'kinetic', 'potential', 'energy', 'px', 'py', 'pz', 'Lx', 'Ly', 'Lz', 'energy_error', 'momentum_error', and 'angular_momentum_error'.

What to expect

The Velocity Verlet and Euler-Cromer integrators are built from "kicks" and "drifts" that conserve linear and angular momentum exactly, so for 'method = "verlet"' or '"euler_cromer"' those two errors should stay at the level of floating-point rounding (around 1e-15) for any time step. Energy is conserved only approximately: with Verlet its error oscillates once per orbit within a band whose width scales as the square of the time step, and does not grow. A steady drift in energy means the time step is too large for the fastest or most eccentric orbit in the system, or that 'method = "euler"' was used. A sudden jump marks a close encounter the step could not resolve. A drift in momentum or angular momentum under Verlet points to a problem in the setup rather than the integrator.

Examples

sim <- create_system() |>
  add_body("Star", mass = 1e30) |>
  add_body("Planet", mass = 1e24, x = 1e11, vy = 30000) |>
  simulate_system(time_step = seconds_per_hour * 6,
                  duration = seconds_per_year * 2)

cq <- conserved_quantities(sim)

# Verlet: bounded energy error, momenta at rounding level
range(cq$energy_error)
max(cq$angular_momentum_error)

plot(cq$time / seconds_per_day, cq$energy_error, type = "l",
     xlab = "Day", ylab = "Relative energy error")

Continue a simulation from its last state

Description

Rebuilds the system from the final time step of a simulation, runs it forward, and appends the new rows with 'time' continuing from where the previous run ended. Use it to extend a run without starting over, to change the time step partway through (small steps through a close approach, large ones elsewhere), or to apply a change between segments, such as a velocity kick or an added body, by editing the tibble before continuing.

Usage

continue_simulation(sim_data, time_step, duration, G = NULL, ...)

Arguments

Details

The last row of each body in 'sim_data' is a complete state, so the restart is exact: the new segment begins from precisely where the old one stopped. The new segment's first step duplicates the old segment's last and is dropped, so 'time' is strictly increasing in the result.

Each segment is integrated independently with a fixed step. Changing the step between segments disturbs the integration by about one step's worth of error at the switch, which is usually far smaller than the error saved by using a small step only where it is needed. A comet's perihelion passage is the typical use: run with a step of days out to a few AU, a step of hours through perihelion, and days again on the way out.

Value

A tibble with the same columns as 'sim_data' and the new time steps appended.

Examples

# A comet on a highly eccentric orbit, started at aphelion
comet <- create_system() |>
  add_sun() |>
  add_body_keplerian("Comet", mass = 1e14, parent = "Sun",
                     a = 5 * distance_earth_sun, e = 0.9, nu = 180)

# Coarse steps on the way in, fine steps through perihelion, coarse again
sim <- comet |>
  simulate_system(time_step = seconds_per_day * 5,
                  duration = seconds_per_year * 5) |>
  continue_simulation(time_step = seconds_per_hour * 2,
                      duration = seconds_per_year * 1.2) |>
  continue_simulation(time_step = seconds_per_day * 5,
                      duration = seconds_per_year * 5)

range(sim$time) / seconds_per_year

Initialize an orbitr simulation system

Description

Sets up the foundational data structure for an orbital simulation. Universal N-body gravity is automatically integrated into the system using the specified gravitational constant.

Usage

create_system(G = gravitational_constant)

Arguments

Value

An empty 'orbit_system' object ready for bodies to be added.

Examples

# Creates a system with standard gravity
my_universe <- create_system()

# Creates a universe with 10x stronger gravity
heavy_universe <- create_system(G = gravitational_constant * 10)

# Creates a zero-gravity sandbox
floating_universe <- create_system(G = 0)

Export body states to CSV

Description

Writes the body table (id, mass, position, and velocity) from an 'orbit_system' to a CSV file. This is useful for sharing initial conditions with collaborators or loading them into other tools like Python or Excel.

Usage

export_bodies(system, path)

Arguments

Value

'system', invisibly.

Examples

sys <- create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  add_planet("Mars",  parent = "Sun")

export_bodies(sys, file.path(tempdir(), "bodies.csv"))

Extract the body table from a system

Description

Returns the bodies in an 'orbit_system' as a standalone tibble. Useful when you want to inspect, filter, or save the body states without dealing with the full system object.

Usage

get_bodies(system)

Arguments

Value

A tibble with columns 'id', 'mass', 'x', 'y', 'z', 'vx', 'vy', 'vz'.

Examples

sys <- create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  add_planet("Mars",  parent = "Sun")

# Get the tibble
get_bodies(sys)

# Use with dplyr

get_bodies(sys) |>
  dplyr::filter(mass > 1e24)

Total energy of a simulated system

Description

Computes the total kinetic energy, total gravitational potential energy, and their sum at every time step of a simulation. In exact Newtonian gravity the total is constant; with a numerical integrator it is not, and how far it wanders is a direct measure of integration error (see [conserved_quantities()] for the error itself).

Usage

get_energy(sim_data, G = NULL, softening = NULL)

Arguments

Details

At each time step, 'kinetic' is \sum_j \tfrac{1}{2} m_j v_j^2, 'potential' is -\sum_{j<k} G m_j m_k / \sqrt{r_{jk}^2 + \varepsilon^2} summed over every unordered pair of bodies, and 'energy' is their sum. Potential energy is a property of pairs, not of individual bodies, which is why the function reports system totals only.

Value

A tibble with one row per time step and columns 'time', 'kinetic', 'potential', and 'energy', all in joules.

Examples

sim <- create_system() |>
  add_body("Star", mass = 1e30) |>
  add_body("Planet", mass = 1e24, x = 1e11, vy = 30000) |>
  simulate_system(time_step = seconds_per_hour * 6,
                  duration = seconds_per_year)

energy <- get_energy(sim)
energy

# Kinetic and potential trade off along the orbit; the total barely moves
plot(energy$time / seconds_per_day, energy$kinetic, type = "l",
     xlab = "Day", ylab = "Kinetic energy (J)")

Total linear and angular momentum of a simulated system

Description

'get_momentum()' computes the total linear momentum \mathbf{P} = \sum_j m_j \mathbf{v}_j and 'get_angular_momentum()' the total angular momentum about the origin, \mathbf{L} = \sum_j m_j \mathbf{r}_j \times \mathbf{v}_j, at every time step of a simulation.

Usage

get_momentum(sim_data)

get_angular_momentum(sim_data)

Arguments

Details

Both are conserved exactly by Newtonian gravity, because every force comes in an equal and opposite pair directed along the line between the two bodies. The Velocity Verlet and Euler-Cromer integrators also conserve both to floating-point rounding, at any time step, so a drift in either points to a problem in the setup rather than the step size. A system whose total momentum is not zero has a center of mass that drifts at a constant velocity; see 'shift_reference_frame(sim_data, "barycenter")'.

Value

A tibble with one row per time step: 'time, px, py, pz' (kg m/s) for 'get_momentum()', and 'time, Lx, Ly, Lz' (kg m^2/s) for 'get_angular_momentum()'.

Examples

# A binary built with zero total momentum
sim <- create_system() |>
  add_body("A", mass = 2e30, x = 5e10, vy = 15000) |>
  add_body("B", mass = 1e30, x = -1e11, vy = -30000) |>
  simulate_system(time_step = seconds_per_hour, duration = seconds_per_year)

get_momentum(sim)          # px, py, pz all zero to rounding
get_angular_momentum(sim)  # Lz constant to rounding

Osculating orbital elements of a body relative to a parent

Description

The inverse of [add_body_keplerian()]: computes the classical Keplerian elements of ‘body'’s orbit about 'parent' at every time step of a simulation, from their relative position and velocity. In a two-body system the elements are constant (up to integration error). With more bodies present they drift as the orbit is perturbed, and the value at each instant is the *osculating* orbit: the ellipse the body would follow from that moment on if every other perturbation were switched off.

Usage

get_orbital_elements(sim_data, body, parent, mu = NULL, G = NULL)

Arguments

Details

The computation follows the standard textbook route. The specific angular momentum \mathbf{h} = \mathbf{r} \times \mathbf{v} gives the inclination (\cos i = h_z / h) and, through the node vector \hat{\mathbf{z}} \times \mathbf{h}, the longitude of the ascending node. The eccentricity vector \mathbf{e} = (\mathbf{v} \times \mathbf{h})/\mu - \mathbf{r}/r points to periapsis and gives the eccentricity, the argument of periapsis (the angle from the node to \mathbf{e}), and the true anomaly (the angle from \mathbf{e} to \mathbf{r}). The semi-major axis comes from the vis-viva equation, a = 1 / (2/r - v^2/\mu).

Two cases are degenerate and need a convention. For an orbit in the reference plane (i = 0) the ascending node is undefined: 'lan' is reported as 0 and 'arg_pe' is measured from the x axis. For a circular orbit (e = 0) periapsis is undefined: 'arg_pe' is reported as 0 and 'nu' is measured from the ascending node (or the x axis).

For an unbound orbit (e \ge 1) 'a' is negative and 'period' is 'NA'.

Value

A tibble with one row per time step and columns 'time', 'a' (meters), 'e', 'i', 'lan', 'arg_pe', 'nu' (degrees; angles other than 'i' are in [0, 360)), and 'period' (seconds).

Examples

# Round trip: the elements that went in come back out at t = 0
sim <- create_system() |>
  add_sun() |>
  add_body_keplerian("Mars", mass = mass_mars, parent = "Sun",
                     a = distance_mars_sun, e = 0.0934, i = 1.85,
                     lan = 49.6, arg_pe = 286.5, nu = 120) |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_day * 687)

elements <- get_orbital_elements(sim, "Mars", "Sun")
elements[1, ]

# Two bodies: a and e are constant to integration error
range(elements$e)

# Earth's eccentricity drifting under Jupiter's pull
perturbed <- create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  add_planet("Jupiter", parent = "Sun", nu = 90) |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year * 12)

earth <- get_orbital_elements(perturbed, "Earth", "Sun")
plot(earth$time / seconds_per_year, earth$e, type = "l",
     xlab = "Years", ylab = "Osculating eccentricity")

Load a pre-built solar system

Description

A convenience function that creates a complete solar system with the Sun and all eight planets (plus optionally the Moon and Pluto) using real orbital data. Bodies are placed using Keplerian orbital elements from the JPL DE440 ephemeris (J2000 epoch), giving realistic eccentricities, inclinations, and orbital orientations out of the box.

Usage

load_solar_system(moon = TRUE, pluto = TRUE)

Arguments

Details

This is a quick way to get a physically accurate starting point without typing out a dozen [add_body()] calls. The returned system is a normal 'orbit_system' that you can modify further — add bodies, change parameters, or pipe straight into [simulate_system()].

Value

An 'orbit_system' object containing the Sun and planets, ready for simulation.

Examples

# Simulate the full solar system for one year
solar <- load_solar_system() |>
  simulate_system(
    time_step = seconds_per_day,
    duration  = seconds_per_year
  )

plot_orbits(solar)

# Just the Sun and planets, no Moon or Pluto
load_solar_system(moon = FALSE, pluto = FALSE)

Load an orbit_system from disk

Description

Restores an 'orbit_system' previously saved with [save_system()].

Usage

load_system(path)

Arguments

Value

An 'orbit_system' object.

Examples

sys <- create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun")

path <- file.path(tempdir(), "my_system.rds")
save_system(sys, path)
restored <- load_system(path)

Physical Constants for Orbital Mechanics

Description

A curated set of real-world masses and orbital distances for use as convenient starting points in 'orbitr' simulations. All values are in SI units (kilograms and meters).

Usage

gravitational_constant

seconds_per_hour

seconds_per_day

seconds_per_year

mass_sun

mass_earth

mass_moon

mass_mars

mass_jupiter

mass_saturn

mass_venus

mass_mercury

mass_uranus

mass_neptune

mass_pluto

distance_earth_sun

distance_earth_moon

distance_mars_sun

distance_jupiter_sun

distance_venus_sun

distance_mercury_sun

distance_saturn_sun

distance_uranus_sun

distance_neptune_sun

distance_pluto_sun

speed_earth

speed_moon

speed_mars

speed_jupiter

speed_venus

speed_mercury

speed_saturn

speed_uranus

speed_neptune

speed_pluto

Format

Numeric scalar in kilograms.

Details

‘gravitational_constant': Newton’s gravitational constant (6.6743 x 10^-11 m^3 kg^-1 s^-2). Source: CODATA 2018 recommended value. Use this with 'create_system()' to scale gravity: 'create_system(G = gravitational_constant * 10)'.

'seconds_per_hour': 3,600 seconds. Convenient for setting 'time_step' in lunar or close-orbit simulations.

'seconds_per_day': 86,400 seconds. Convenient for setting 'time_step' in planetary-scale simulations.

'seconds_per_year': 31,557,600 seconds (365.25 days, the Julian year). Convenient for setting 'duration' in 'simulate_system()'.

'mass_sun': Mass of the Sun (1.989 x 10^30 kg). Source: IAU 2015 nominal solar mass.

'mass_earth': Mass of the Earth (5.972 x 10^24 kg). Source: IAU 2015 nominal Earth mass.

'mass_moon': Mass of the Moon (7.342 x 10^22 kg). Source: JPL DE440 ephemeris.

'mass_mars': Mass of Mars (6.417 x 10^23 kg). Source: JPL DE440 ephemeris.

'mass_jupiter': Mass of Jupiter (1.898 x 10^27 kg). Source: JPL DE440 ephemeris.

'mass_saturn': Mass of Saturn (5.683 x 10^26 kg). Source: JPL DE440 ephemeris.

'mass_venus': Mass of Venus (4.867 x 10^24 kg). Source: JPL DE440 ephemeris.

'mass_mercury': Mass of Mercury (3.301 x 10^23 kg). Source: JPL DE440 ephemeris.

'mass_uranus': Mass of Uranus (8.681 x 10^25 kg). Source: JPL DE440 ephemeris.

'mass_neptune': Mass of Neptune (1.024 x 10^26 kg). Source: JPL DE440 ephemeris.

'mass_pluto': Mass of Pluto (1.309 x 10^22 kg). Source: JPL DE440 ephemeris. Pluto is a dwarf planet but is included for convenience.

‘distance_earth_sun': Semi-major axis of Earth’s orbit around the Sun (1.496 x 10^11 m, ~149.6 million km). Earth's actual distance varies between ~147.1 million km (perihelion) and ~152.1 million km (aphelion).

‘distance_earth_moon': Semi-major axis of the Moon’s orbit around Earth (3.844 x 10^8 m, ~384,400 km). The Moon's actual distance varies between ~363,300 km (perigee) and ~405,500 km (apogee).

‘distance_mars_sun': Semi-major axis of Mars’s orbit around the Sun (2.279 x 10^11 m, ~227.9 million km). Mars has a notably eccentric orbit (e = 0.093), ranging from ~206.7 million km to ~249.2 million km.

‘distance_jupiter_sun': Semi-major axis of Jupiter’s orbit around the Sun (7.785 x 10^11 m, ~778.5 million km).

‘distance_venus_sun': Semi-major axis of Venus’s orbit around the Sun (1.082 x 10^11 m, ~108.2 million km). Venus has the most circular orbit of any planet (e = 0.007).

‘distance_mercury_sun': Semi-major axis of Mercury’s orbit around the Sun (5.791 x 10^10 m, ~57.9 million km). Mercury has the most eccentric planetary orbit (e = 0.206), ranging from ~46.0 million km to ~69.8 million km.

‘distance_saturn_sun': Semi-major axis of Saturn’s orbit around the Sun (1.434 x 10^12 m, ~1.434 billion km).

‘distance_uranus_sun': Semi-major axis of Uranus’s orbit around the Sun (2.871 x 10^12 m, ~2.871 billion km).

‘distance_neptune_sun': Semi-major axis of Neptune’s orbit around the Sun (4.495 x 10^12 m, ~4.495 billion km).

‘distance_pluto_sun': Semi-major axis of Pluto’s orbit around the Sun (5.906 x 10^12 m, ~5.906 billion km). Pluto has a highly eccentric orbit (e = 0.249), ranging from ~4.437 billion km to ~7.376 billion km.

'speed_earth': Mean orbital speed of Earth around the Sun (29,780 m/s).

'speed_moon': Mean orbital speed of the Moon around Earth (1,022 m/s).

'speed_mars': Mean orbital speed of Mars around the Sun (24,070 m/s).

'speed_jupiter': Mean orbital speed of Jupiter around the Sun (13,060 m/s).

'speed_venus': Mean orbital speed of Venus around the Sun (35,020 m/s).

'speed_mercury': Mean orbital speed of Mercury around the Sun (47,360 m/s).

'speed_saturn': Mean orbital speed of Saturn around the Sun (9,680 m/s).

'speed_uranus': Mean orbital speed of Uranus around the Sun (6,800 m/s).

'speed_neptune': Mean orbital speed of Neptune around the Sun (5,430 m/s).

'speed_pluto': Mean orbital speed of Pluto around the Sun (4,740 m/s).

A Note on "Distance" Constants

Orbital distances are not truly constant. Every orbit is an ellipse, so the separation between two bodies changes continuously. The distances provided here are **semi-major axes** — the average of the closest approach (periapsis) and farthest point (apoapsis). The semi-major axis is the single most characteristic length scale of an elliptical orbit: it determines the orbital period via Kepler's Third Law, and when paired with the circular velocity at that distance, it produces a near-circular orbit that closely approximates the real trajectory.

For example, the Earth-Sun distance varies from about 147.1 million km (perihelion in January) to 152.1 million km (aphelion in July). The semi-major axis of 149.598 million km sits right in the middle and gives the correct orbital period of one year.

Plot Orbital Trajectories (Smart 2D/3D Dispatch)

Description

Plot Orbital Trajectories (Smart 2D/3D Dispatch)

Usage

plot_orbits(sim_data, three_d = NULL)

Arguments

Value

A 'ggplot' object (2D) or a 'plotly' HTML widget (3D) showing the orbital trajectories of all bodies in the simulation.

Plot 3D Interactive Orbital Trajectories

Description

Generates an interactive 3D visualization of the orbital system using plotly. You can click, drag to rotate, and scroll to zoom in on the trajectories.

Usage

plot_orbits_3d(sim_data)

Arguments

Value

A plotly HTML widget displaying the 3D orbits.

Examples

create_system() |>
  add_body("Earth", mass = mass_earth) |>
  add_body("Moon", mass = mass_moon,
           x = distance_earth_moon, vy = speed_moon, vz = 150) |>
  simulate_system(time_step = seconds_per_hour, duration = seconds_per_day * 30) |>
  plot_orbits_3d()

Plot System Snapshot at a Single Time (Smart 2D/3D Dispatch)

Description

Plots the position of every body in the system at a single time step, optionally with the full orbital trajectories drawn faintly behind. This is the snapshot counterpart to [plot_orbits()], which draws full trajectories.

Usage

plot_system(sim_data, time = NULL, trails = FALSE, three_d = NULL)

Arguments

Details

If any body has non-zero motion in the Z dimension (or 'three_d = TRUE'), [plot_system_3d()] is used; otherwise a 2D 'ggplot2' plot is returned.

Value

A 'ggplot' object (2D) or a 'plotly' HTML widget (3D).

Examples

sim <- create_system() |>
  add_sun() |>
  add_body("Earth", mass = mass_earth, x = distance_earth_sun, vy = speed_earth) |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year)

# Final state with faint orbit trails
plot_system(sim)

# State at day 100, no trails
plot_system(sim, time = seconds_per_day * 100, trails = FALSE)

Plot 3D Interactive System Snapshot at a Single Time

Description

The 3D counterpart to [plot_system()]. Draws every body's position at a chosen time as a sphere in an interactive plotly scene, optionally with the full orbital trajectories shown faintly behind.

Usage

plot_system_3d(sim_data, time = NULL, trails = FALSE)

Arguments

Value

A 'plotly' HTML widget.

Examples

create_system() |>
  add_body("Earth", mass = mass_earth) |>
  add_body("Moon",  mass = mass_moon,
           x = distance_earth_moon, vy = speed_moon, vz = 100) |>
  simulate_system(time_step = seconds_per_hour, duration = seconds_per_day * 30) |>
  plot_system_3d()

Print an orbit_system

Description

Displays a compact, human-readable summary of an 'orbit_system' showing the gravitational constant and a tibble of body states.

Usage

## S3 method for class 'orbit_system'
print(x, ...)

Arguments

Value

'x', invisibly.

Examples

create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  add_planet("Mars",  parent = "Sun")

Remove one or more bodies from the system

Description

Drops bodies by name from an 'orbit_system'. This is the counterpart to [add_body()] — use it to strip out bodies you no longer need before simulating, or to prune a system built with [load_solar_system()].

Usage

remove_body(system, id)

Arguments

Details

“' load_solar_system() |> remove_body(c("Pluto", "Moon")) |> simulate_system(time_step = seconds_per_day, duration = seconds_per_year) “'

Value

The updated 'orbit_system' with the specified bodies removed.

Examples

# Remove a single body
create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  add_planet("Mars",  parent = "Sun") |>
  remove_body("Mars")


# Remove multiple bodies from the full solar system
load_solar_system() |>
  remove_body(c("Pluto", "Moon")) |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year) |>
  plot_orbits()

Save an orbit_system to disk

Description

Saves the full 'orbit_system' object (bodies, forces, and time) to an '.rds' file so it can be restored later with [load_system()]. This preserves everything — the gravitational constant, body states, and class — exactly as it was.

Usage

save_system(system, path)

Arguments

Value

'system', invisibly.

Examples

sys <- create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun")

save_system(sys, file.path(tempdir(), "my_system.rds"))

Shift the coordinate reference frame of the simulation

Description

Recalculates the positions and velocities of all bodies relative to a specific target body, or to the system's center of mass. This effectively "anchors the camera" to the chosen point, placing it at the origin (0, 0, 0) for all time steps.

Usage

shift_reference_frame(sim_data, center_id, keep_center = TRUE)

Arguments

Details

The shift is a Galilean transformation: at every time step the chosen point's position and velocity are subtracted from every body. No physics changes; the same forces and accelerations produced the data, and you are only choosing where to stand when you look at it.

The barycentric frame is the natural one for binary stars and any other system where no single body dominates. In it the total momentum is zero and the center of mass sits at the origin for the whole run, which removes the slow drift you get when a system is built with one body at rest but nonzero total momentum (for example, a planet given an orbital velocity around a star that was not given the balancing recoil).

If a body in the system is itself named '"barycenter"', that body is used as the center rather than the center of mass.

Value

A tidy 'tibble' with updated 'x', 'y', 'z', 'vx', 'vy', and 'vz' columns.

Examples

# Simulate Sun-Earth-Moon
orbit_data <- create_system() |>
  add_sun() |>
  add_body("Earth", mass = mass_earth, x = distance_earth_sun, vy = speed_earth) |>
  add_body("Moon", mass = mass_moon, x = distance_earth_sun + distance_earth_moon,
           vy = speed_earth + speed_moon) |>
  simulate_system(time_step = seconds_per_hour, duration = seconds_per_year)

# Shift view to Earth and plot
orbit_data |>
  shift_reference_frame(center_id = "Earth") |>
  plot_orbits()

# The Sun started at rest with Jupiter in orbit: the pair's center of mass
# drifts, because the total momentum is not zero. The barycentric frame
# removes the drift and shows the Sun's own small orbit.
sun_jupiter <- create_system() |>
  add_sun() |>
  add_planet("Jupiter", parent = "Sun") |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year * 12)

sun_jupiter |>
  shift_reference_frame("barycenter") |>
  plot_orbits(three_d = FALSE)

Simulate kinematics for an orbitr system

Description

Propagates the physical state of an 'orbit_system' through time using numerical integration. This engine supports multiple mathematical methods, defaulting to the energy-conserving Velocity Verlet algorithm to ensure highly stable orbital trajectories.

Usage

simulate_system(
  system,
  time_step = seconds_per_hour,
  duration = seconds_per_year,
  method = "verlet",
  softening = 0,
  use_cpp = TRUE
)

Arguments

Value

A tidy 'tibble' containing the physical state (time, id, mass, x, y, z, vx, vy, vz) of every body at every time step. The run's settings are recorded as attributes ('"G"', '"softening"', '"method"', '"time_step"'), which [get_energy()], [conserved_quantities()], and [continue_simulation()] use as defaults.

Examples

my_universe <- create_system() |>
  add_body("Earth", mass = mass_earth) |>
  add_body("Moon", mass = mass_moon, x = distance_earth_moon, vy = speed_moon) |>
  simulate_system(time_step = seconds_per_hour, duration = seconds_per_day * 28)

Rebuild an orbit_system from a simulation snapshot

Description

Takes the state of every body at one time step of a simulation and returns a new 'orbit_system' with those positions and velocities as its initial conditions. This is the bridge from a finished run back to a system you can modify (add a body, remove one, change a velocity) and simulate again.

Usage

system_from_simulation(sim_data, time = NULL, G = NULL)

Arguments

Value

An ‘orbit_system' whose bodies have the snapshot’s positions and velocities.

Examples

sim <- create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_day * 100)

# Where everything was on day 100, as a system ready to simulate again
later <- system_from_simulation(sim)
later