Orbital Perturbations & Event Detection¶

Real orbits deviate from Kepler's ideal due to perturbations. Here we model the J2 oblateness perturbation (Earth's equatorial bulge) and use CVODE's event detection to find periapsis/apoapsis passages and eclipse boundaries.

In [36]:
load("numerics")$
load("numerics-sundials")$
load("ax-plots")$

1. J2 Perturbation¶

Earth's oblateness is characterised by $J_2 \approx 1.08263 \times 10^{-3}$. The perturbing acceleration in Cartesian coordinates is:

$$a_x = -\frac{\mu x}{r^3}\left[1 - J_2\frac{3R_E^2}{2r^2}\left(\frac{5z^2}{r^2} - 1\right)\right]$$

and similarly for $a_y$, $a_z$ (with $5z^2/r^2 - 3$ for the $z$ component).

We use Earth-standard units: $\mu = 398600.4418$ km$^3$/s$^2$, $R_E = 6378.137$ km.

In [37]:
/* Earth constants */
mu_e : 398600.4418$    /* km^3/s^2 */
R_E : 6378.137$        /* km */
J2 : 1.08263e-3$

/* State: [x, y, z, vx, vy, vz] in km and km/s */
r2_expr : x_j^2 + y_j^2 + z_j^2$
r_expr : sqrt(r2_expr)$
r5 : r2_expr^(5/2)$

/* J2 perturbation factor */
fac_j2 : 3*J2*R_E^2 / (2*r2_expr)$
z2_over_r2 : z_j^2 / r2_expr$

ax_j2 : -mu_e * x_j / r_expr^3 * (1 - fac_j2 * (5*z2_over_r2 - 1))$
ay_j2 : -mu_e * y_j / r_expr^3 * (1 - fac_j2 * (5*z2_over_r2 - 1))$
az_j2 : -mu_e * z_j / r_expr^3 * (1 - fac_j2 * (5*z2_over_r2 - 3))$

rhs_j2 : [vx_j, vy_j, vz_j, ax_j2, ay_j2, az_j2]$
vars_j2 : [t, x_j, y_j, z_j, vx_j, vy_j, vz_j]$
In [38]:
/* LEO orbit: 400 km altitude, 45 deg inclination */
alt : 400$
r0_leo : R_E + alt$
v_circ : sqrt(mu_e / r0_leo)$

/* Inclined velocity */
inc : float(%pi/4)$   /* 45 degrees */
y0_j2 : [r0_leo, 0.0, 0.0, 0.0, v_circ*cos(inc), v_circ*sin(inc)]$

T_orb : 2 * float(%pi) * sqrt(r0_leo^3 / mu_e)$
print("Orbital period:", T_orb, "seconds =", T_orb/60, "minutes")$

/* Integrate for 10 orbits */
n_orbits : 10$
tspan_j2 : np_linspace(0.0, n_orbits * T_orb, 5000)$

orbit_j2 : np_cvode(rhs_j2, vars_j2, y0_j2, tspan_j2, 1e-10, 1e-10, bdf)$
print("J2 orbit shape:", np_shape(orbit_j2))$

Orbital period: 5553.624271252228 seconds = 92.56040452087046 minutes
J2 orbit shape: [5000,7]
In [39]:
/* Also integrate pure Keplerian (no J2) for comparison */
rhs_kep : [vx_j, vy_j, vz_j,
           -mu_e*x_j/(x_j^2+y_j^2+z_j^2)^(3/2),
           -mu_e*y_j/(x_j^2+y_j^2+z_j^2)^(3/2),
           -mu_e*z_j/(x_j^2+y_j^2+z_j^2)^(3/2)]$

orbit_kep : np_cvode(rhs_kep, vars_j2, y0_j2, tspan_j2, 1e-10, 1e-10, bdf)$

/* Plot x-y projection: J2 vs Keplerian */
ax_draw2d(
  aspect_ratio = true,
  line_width = 1,
  color = "#cccccc",
  name = "Keplerian",
  lines(np_slice(orbit_kep, all, 1), np_slice(orbit_kep, all, 2)),
  color = royalblue,
  name = "With J2",
  lines(np_slice(orbit_j2, all, 1), np_slice(orbit_j2, all, 2)),
  color = seagreen,
  marker_size = 6,
  name = "",
  points([0], [0]),
  title = "J2 Perturbation: 10 Orbits (LEO, 45 deg inclination)"
)$
No description has been provided for this image

RAAN Precession¶

The J2 perturbation causes the right ascension of ascending node (RAAN) to precess. For a circular orbit:

$$\dot{\Omega} = -\frac{3}{2} J_2 \left(\frac{R_E}{r_0}\right)^2 n \cos i$$

where $n = \sqrt{\mu/r_0^3}$ is the mean motion. Let's measure the RAAN drift by tracking the ascending node.

In [40]:
/* Compute angular momentum direction over time */
t_j2 : np_slice(orbit_j2, all, 0)$
n_pts : first(np_shape(orbit_j2))$

/* Predicted RAAN precession rate for circular orbit */
n_motion : sqrt(mu_e / r0_leo^3)$
raan_rate : -1.5 * J2 * (R_E/r0_leo)^2 * n_motion * cos(inc)$
print("Predicted RAAN rate:", raan_rate, "rad/s =", raan_rate * 180/float(%pi) * 86400, "deg/day")$

/* Compute actual RAAN from angular momentum vector direction */
/* h = r x v; RAAN = atan2(-hx, hy) for inclined orbit */
/* Phase-unwrap to remove atan2 discontinuities at +/-pi */
raan_list : []$
prev_omega : false$
for i : 0 thru n_pts - 1 step 50 do (
  rx : np_ref(orbit_j2, i, 1),
  ry : np_ref(orbit_j2, i, 2),
  rz : np_ref(orbit_j2, i, 3),
  vx : np_ref(orbit_j2, i, 4),
  vy : np_ref(orbit_j2, i, 5),
  vz : np_ref(orbit_j2, i, 6),
  r_v : ndarray([rx, ry, rz], [3]),
  v_v : ndarray([vx, vy, vz], [3]),
  h_v : np_cross(r_v, v_v),
  hx : np_ref(h_v, 0),
  hy : np_ref(h_v, 1),
  omega : atan2(-hx, hy),
  if prev_omega # false then (
    while omega - prev_omega > float(%pi) do omega : omega - 2*float(%pi),
    while prev_omega - omega > float(%pi) do omega : omega + 2*float(%pi)
  ),
  prev_omega : omega,
  raan_list : append(raan_list, [[np_ref(orbit_j2, i, 0), omega]])
)$

raan_t : map(first, raan_list)$
raan_v : map(second, raan_list)$

/* Subtract initial RAAN so drift starts at zero (matches prediction) */
raan0 : first(raan_v)$
raan_v : map(lambda([r], r - raan0), raan_v)$

/* Predicted linear drift */
raan_pred : map(lambda([ti], raan_rate * ti), raan_t)$

ax_draw2d(
  line_width = 2,
  color = royalblue,
  name = "Measured RAAN",
  points(raan_t, raan_v),
  color = tomato,
  dash = "dash",
  name = "Predicted (linear)",
  lines(raan_t, raan_pred),
  title = "RAAN Precession due to J2",
  xlabel = "Time (s)",
  ylabel = "\u0394RAAN (rad)"
)$

Predicted RAAN rate: -1.150341252712621e-6 rad/s = -5.694597974611462 deg/day
No description has been provided for this image

2. Periapsis / Apoapsis Detection¶

Use CVODE's rootfinding to detect periapsis and apoapsis passages. The event function is the radial velocity $\dot{r} = \vec{r} \cdot \vec{v} / |\vec{r}|$, which crosses zero at periapsis and apoapsis.

In [41]:
/* 2D elliptical orbit for clarity */
mu_n : 1.0$
r3_n : (xn^2 + yn^2)^(3/2)$
rhs_n : [vxn, vyn, -mu_n*xn/r3_n, -mu_n*yn/r3_n]$
vars_n : [t, xn, yn, vxn, vyn]$

/* Elliptical orbit: e ≈ 0.5 */
v0_e : 0.7$
y0_e : [1.0, 0.0, 0.0, v0_e]$

/* Event: radial velocity r_dot = (x*vx + y*vy) */
rdot_event : xn*vxn + yn*vyn$

/* Integrate for 3 orbits */
eps_e : v0_e^2/2 - mu_n$
a_e : -mu_n/(2*eps_e)$
T_e : 2*float(%pi)*sqrt(a_e^3/mu_n)$

tspan_ev : np_linspace(0.0, 3*T_e, 1000)$
sol_ev : np_cvode(rhs_n, vars_n, y0_e, tspan_ev,
                  1e-12, 1e-12, adams, [rdot_event])$

traj_ev : first(sol_ev)$
events_ev : second(sol_ev)$

print("Number of events:", length(events_ev))$
for i : 1 thru length(events_ev) do (
  ev : events_ev[i],
  t_ev : first(ev),
  y_ev : second(ev),
  r_ev : sqrt(y_ev[1]^2 + y_ev[2]^2),
  type_str : if r_ev < a_e then "periapsis" else "apoapsis",
  print(sconcat("  ", type_str, " at t = ", t_ev, ", r = ", r_ev))
)$

Number of events: 6
  periapsis at t = 1.6931071811827705, r = 0.3245033111734518
  apoapsis at t = 3.3862143621390297, r = 0.9999999999477329
  periapsis at t = 5.079321543078222, r = 0.32450331112246633
  apoapsis at t = 6.772428723853829, r = 0.9999999998963812
  periapsis at t = 8.465535904640339, r = 0.32450331109196534
  apoapsis at t = 10.15864308530667, r = 0.9999999998731922
In [42]:
/* Plot orbit with periapsis/apoapsis markers */
xt : np_slice(traj_ev, all, 1)$
yt : np_slice(traj_ev, all, 2)$

/* Separate periapsis and apoapsis events */
peri_x : []$  peri_y : []$
apo_x : []$   apo_y : []$
for i : 1 thru length(events_ev) do (
  ev : events_ev[i],
  y_ev : second(ev),
  r_ev : sqrt(y_ev[1]^2 + y_ev[2]^2),
  if r_ev < a_e then (
    peri_x : append(peri_x, [y_ev[1]]),
    peri_y : append(peri_y, [y_ev[2]])
  ) else (
    apo_x : append(apo_x, [y_ev[1]]),
    apo_y : append(apo_y, [y_ev[2]])
  )
)$

ax_draw2d(
  aspect_ratio = true,
  line_width = 2,
  color = "#888888",
  name = "Orbit",
  lines(xt, yt),
  color = red,
  marker_size = 6,
  name = "Periapsis",
  points(peri_x, peri_y),
  color = blue,
  marker_size = 6,
  name = "Apoapsis",
  points(apo_x, apo_y),
  color = gold,
  marker_size = 8,
  name = "",
  points([0], [0]),
  title = "Periapsis/Apoapsis Detection (3 orbits)"
)$
No description has been provided for this image

3. Eclipse Detection¶

Detect when a satellite enters Earth's shadow (eclipse). Using a cylindrical shadow model: the satellite is eclipsed when $x < 0$ and $\sqrt{y^2 + z^2} < R_E$.

We define two event functions:

  • $g_1 = x$ (crossing the terminator plane)
  • $g_2 = \sqrt{y^2 + z^2} - R_E$ (entering/leaving the shadow cylinder)
In [43]:
/* Circular LEO orbit, 400 km altitude, 0 deg inclination (equatorial) */
alt_ecl : 400$
r0_ecl : R_E + alt_ecl$
vc_ecl : sqrt(mu_e / r0_ecl)$

/* 3D Keplerian ODE */
rhs_ecl : [vx_j, vy_j, vz_j,
           -mu_e*x_j/(x_j^2+y_j^2+z_j^2)^(3/2),
           -mu_e*y_j/(x_j^2+y_j^2+z_j^2)^(3/2),
           -mu_e*z_j/(x_j^2+y_j^2+z_j^2)^(3/2)]$

/* Equatorial orbit */
y0_ecl : [r0_ecl, 0.0, 0.0, 0.0, vc_ecl, 0.0]$

/* Eclipse events:
   g1 = x  (terminator crossing)
   We'll track x-crossings to find eclipse entry/exit */
T_ecl : 2*float(%pi)*sqrt(r0_ecl^3/mu_e)$
tspan_ecl : np_linspace(0.0, 2*T_ecl, 2000)$

sol_ecl : np_cvode(rhs_ecl, vars_j2, y0_ecl, tspan_ecl,
                   1e-10, 1e-10, adams, [x_j])$

traj_ecl : first(sol_ecl)$
events_ecl : second(sol_ecl)$

print("Eclipse-related events (x=0 crossings):", length(events_ecl))$
for i : 1 thru length(events_ecl) do (
  ev : events_ecl[i],
  t_ev : first(ev),
  y_ev : second(ev),
  vx_ev : y_ev[4],
  direction : if vx_ev < 0 then "entering shadow" else "leaving shadow",
  print(sconcat("  t = ", t_ev/60, " min: ", direction))
)$

Eclipse-related events (x=0 crossings): 4
  t = 23.140101112070635 min: entering shadow
  t = 69.42030287472667 min: leaving shadow
  t = 115.7005041771527 min: entering shadow
  t = 161.9807060194475 min: leaving shadow
In [44]:
/* Plot orbit with eclipse arcs highlighted */
xe_ecl : np_slice(traj_ecl, all, 1)$
ye_ecl : np_slice(traj_ecl, all, 2)$
n_ecl : first(np_shape(traj_ecl))$

/* Split orbit into sunlit (x>=0) and eclipsed (x<0) segments.
   We colour points by whether they are in shadow. */
sunlit_x : []$  sunlit_y : []$
shadow_x : []$  shadow_y : []$
for i : 0 thru n_ecl - 1 do (
  xi : np_ref(traj_ecl, i, 1),
  yi : np_ref(traj_ecl, i, 2),
  if xi < 0 and abs(yi) < R_E then (
    shadow_x : endcons(xi, shadow_x),
    shadow_y : endcons(yi, shadow_y)
  ) else (
    sunlit_x : endcons(xi, sunlit_x),
    sunlit_y : endcons(yi, sunlit_y)
  )
)$

ax_draw2d(
  aspect_ratio = true,
  /* Sunlit arc */
  color = royalblue,
  name = "Sunlit",
  points(sunlit_x, sunlit_y),
  /* Eclipsed arc */
  color = "#333333",
  name = "Eclipsed",
  points(shadow_x, shadow_y),
  /* Earth */
  color = seagreen,
  line_width = 2,
  name = "Earth",
  parametric(R_E*cos(s), R_E*sin(s), s, 0, 2*%pi),
  /* Shadow cylinder boundary (two horizontal lines extending left) */
  color = "#aaaaaa",
  dash = "dash",
  line_width = 1,
  name = "Shadow boundary",
  lines([0, -15000], [R_E, R_E]),
  name = "",
  lines([0, -15000], [-R_E, -R_E]),
  title = "Eclipse Detection (Equatorial LEO)",
  xlabel = "x (km)",
  ylabel = "y (km)"
)$
No description has been provided for this image

Summary¶

Concept Technique
J2 oblateness Modified gravity in acceleration terms, compare with Keplerian
RAAN precession np_cross for angular momentum, predicted vs measured drift
Periapsis/apoapsis Event detection on radial velocity $\dot{r} = \vec{r} \cdot \vec{v}$
Eclipse detection Event detection on $x$ coordinate (terminator crossing)
Conservation check Compare perturbed vs unperturbed energy/momentum drift