Orekit radiation pressure fixtures

tests/orekit_comparison.rs

@@ -160,3 +160,54 @@ fn test_the_sun_and_moon_agree_with_orekit() {
     );
     println!("third body: {position:.3e} m, {velocity:.3e} m/s");
 }
+
+/// Radiation pressure on a cannonball, over a point mass.
+///
+/// The matched knobs are the reference pressure (this crate's 1361 W/m^2 rather than
+/// Orekit's 1367), the cross-section and coefficient, the mass, and a *spherical*
+/// occulting Earth, since this crate's eclipse geometry treats the occulter as a sphere.
+/// What is left unmatched is the solar radius, 6.957e8 here against Orekit's 6.95508e8,
+/// which changes only how wide the penumbra is.
+#[test]
+fn test_radiation_pressure_agrees_with_orekit() {
+    let Some(almanac) = almanac() else {
+        eprintln!("skipped: kernels not fetched, run scripts/fetch-kernels.sh");
+        return;
+    };
+    let path = Path::new("tests/fixtures/orekit/radiation-low-earth.json");
+    let arc = ReferenceArc::load(path).expect("the fixture is committed");
+    assert!(arc
+        .setup
+        .forces
+        .iter()
+        .any(|f| f.contains("SolarRadiationPressure")));
+
+    let start = &arc.states[0];
+    let mut propagation = Propagation::new(
+        Epoch::from_gregorian_utc_hms(2026, 1, 1, 0, 0, 0),
+        Orbit
+            .with_position(Vector3::from_column_slice(&start.position))
+            .with_velocity(Vector3::from_column_slice(&start.velocity)),
+        &[TwoBody { body: Body::EARTH }, Sunshine],
+    )
+    .frame(Frame::EARTH_J2000)
+    .spacecraft(Spacecraft::cannonball(500.0, 5.0, 2.2, 1.3))
+    .almanac(almanac)
+    .duration(Duration::from_seconds(arc.duration()))
+    .accuracy(Accuracy::High);
+
+    let solution = propagation.run();
+    let (mut position, mut velocity) = (0.0_f64, 0.0_f64);
+    for state in &arc.states {
+        let ours = Orbit.state(solution.interpolate(state.elapsed));
+        position =
+            position.max((ours.position() - Vector3::from_column_slice(&state.position)).norm());
+        velocity =
+            velocity.max((ours.velocity() - Vector3::from_column_slice(&state.velocity)).norm());
+    }
+    println!("radiation: {position:.3e} m, {velocity:.3e} m/s");
+    assert!(
+        position < 0.01,
+        "worst position disagreement {position:.3e} m over the day"
+    );
+}