inst/doc/checking-a-simulation.R

In orbitr: A Tidy Physics Engine for Building and Visualizing Orbital Simulations

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## -----------------------------------------------------------------------------
library(orbitr)
library(dplyr)
library(ggplot2)

## -----------------------------------------------------------------------------
star_planet <- create_system() |>
  add_body("Star", mass = 1e30) |>
  add_body("Planet", mass = 1e24, x = 1e11, vy = 30000)

sim <- simulate_system(star_planet, time_step = seconds_per_hour * 6,
                       duration = seconds_per_year * 3)

get_energy(sim)

## ----energy-tradeoff----------------------------------------------------------
energy <- get_energy(sim)

bind_rows(
  energy |> transmute(time, name = "kinetic",   value = kinetic),
  energy |> transmute(time, name = "potential", value = potential),
  energy |> transmute(time, name = "total",     value = energy)
) |>
  ggplot(aes(x = time / seconds_per_day, y = value, color = name)) +
  geom_line() +
  labs(x = "Day", y = "Energy (J)", color = NULL) +
  theme_minimal()

## -----------------------------------------------------------------------------
get_momentum(sim)
get_angular_momentum(sim)

## -----------------------------------------------------------------------------
cq <- conserved_quantities(sim)
cq

## ----conservation-compare-----------------------------------------------------
runs <- bind_rows(lapply(c("verlet", "euler_cromer", "euler"), function(m) {
  simulate_system(star_planet, time_step = seconds_per_hour * 6,
                  duration = seconds_per_year * 3, method = m) |>
    conserved_quantities() |>
    mutate(method = m)
}))

runs |>
  ggplot(aes(x = time / seconds_per_year, y = energy_error, color = method)) +
  geom_line() +
  facet_wrap(~ method, scales = "free_y", ncol = 1) +
  labs(x = "Years", y = "Relative energy error", color = NULL) +
  theme_minimal()

## -----------------------------------------------------------------------------
runs |>
  group_by(method) |>
  summarise(max_angular_momentum_error = max(angular_momentum_error),
            max_energy_error = max(abs(energy_error)))

## -----------------------------------------------------------------------------
mars <- create_system() |>
  add_sun() |>
  add_planet("Mars", parent = "Sun", nu = 45) |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_day * 687)

elements <- get_orbital_elements(mars, "Mars", "Sun")
elements

## -----------------------------------------------------------------------------
elements |>
  summarise(a_spread = max(a) - min(a),
            e_spread = max(e) - min(e),
            period_days = mean(period) / seconds_per_day)

## ----osculating---------------------------------------------------------------
earth_alone <- create_system() |>
  add_sun() |>
  add_planet("Earth", parent = "Sun") |>
  simulate_system(time_step = seconds_per_day, duration = seconds_per_year * 12)

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

bind_rows(
  get_orbital_elements(earth_alone, "Earth", "Sun") |> mutate(system = "Sun + Earth"),
  get_orbital_elements(earth_jupiter, "Earth", "Sun") |> mutate(system = "Sun + Earth + Jupiter")
) |>
  ggplot(aes(x = time / seconds_per_year, y = e, color = system)) +
  geom_line() +
  labs(x = "Years", y = "Osculating eccentricity", color = NULL) +
  theme_minimal()

## -----------------------------------------------------------------------------
comet <- create_system() |>
  add_sun() |>
  add_body_keplerian("Comet", mass = 1e14, parent = "Sun",
                     a = 5 * distance_earth_sun, e = 0.9, nu = 180)

# Coarse on the way in, fine through perihelion, coarse on the way out
segmented <- comet |>
  simulate_system(time_step = seconds_per_day * 5, duration = seconds_per_year * 4.9) |>
  continue_simulation(time_step = seconds_per_hour * 2, duration = seconds_per_year * 1.4) |>
  continue_simulation(time_step = seconds_per_day * 5, duration = seconds_per_year * 4.9)

n_distinct(segmented$time)

## ----continue-compare---------------------------------------------------------
uniform <- simulate_system(comet, time_step = seconds_per_day * 5,
                           duration = seconds_per_year * 11.2)

bind_rows(
  conserved_quantities(segmented) |> mutate(run = "5 day / 2 hour / 5 day"),
  conserved_quantities(uniform)   |> mutate(run = "5 day throughout")
) |>
  ggplot(aes(x = time / seconds_per_year, y = energy_error, color = run)) +
  geom_line() +
  labs(x = "Years", y = "Relative energy error", color = NULL) +
  theme_minimal()

## -----------------------------------------------------------------------------
later <- system_from_simulation(segmented, time = seconds_per_year * 5)
later