From 59cb8889fa33e5e7fed42acd2b9d4d8e08e38ff0 Mon Sep 17 00:00:00 2001
From: Mattias Holm <lorrden@openorbit.org>
Date: Thu, 22 Oct 2026 18:52:30 +0200
Subject: [PATCH] Fix numerical propagation bug

---
 .../CelMek/NumericalSpacecraftPropagator.swift | 14 +++++++++++++-
 .../NumericalPropagationTests.swift            | 18 ++++++++++++++++++
 2 files changed, 31 insertions(+), 1 deletion(-)

diff --git a/Sources/CelMek/NumericalSpacecraftPropagator.swift b/Sources/CelMek/NumericalSpacecraftPropagator.swift
index 4c75c96..f0ec980 100644
--- a/Sources/CelMek/NumericalSpacecraftPropagator.swift
+++ b/Sources/CelMek/NumericalSpacecraftPropagator.swift
@@ -158,7 +158,19 @@ public struct NumericalSpacecraftPropagator: StateProvider {
             let remaining = try forward ? target.duration(since: current) : current.duration(since: target)
             let magnitude = min(nextStep,remaining.seconds)
             let final = magnitude >= remaining.seconds
-            let end = final ? target : try current.advanced(by: Duration(seconds: forward ? magnitude : -magnitude))
+            var end = final ? target : try current.advanced(by: Duration(seconds: forward ? magnitude : -magnitude))
+            // Split epochs canonicalize their fractional seconds in 0..<1. A
+            // normal step can land a few rounding units before the endpoint
+            // (e.g. 0.1 + 0.5 on a 0.6 s query), leaving an unresolvable RK
+            // remainder. Integrate that rounding sliver with this step instead
+            // of dropping elapsed time or fabricating a state at a later epoch.
+            // Do not merge genuinely tiny requested/attempted steps: their stage
+            // resolution checks below must still reject collapsed timestamps.
+            let endpointRoundoffSeconds = 8*Double.ulpOfOne
+            if !final, magnitude > 16*endpointRoundoffSeconds,
+               abs(try target.duration(since:end).seconds) <= endpointRoundoffSeconds {
+                end = target
+            }
             guard forward ? end > current : end < current else { throw NumericalPropagationError.epochResolutionExceeded }
             let h = try end.duration(since: current).seconds
             var stages: [[Double]] = [], candidate = y, depleted = false
diff --git a/Tests/CelMekTests/NumericalPropagationTests.swift b/Tests/CelMekTests/NumericalPropagationTests.swift
index 58420cd..ad260c2 100644
--- a/Tests/CelMekTests/NumericalPropagationTests.swift
+++ b/Tests/CelMekTests/NumericalPropagationTests.swift
@@ -137,4 +137,22 @@ import Testing
         }
     }
 
+    @Test func fractionalStepRemaindersReachExactEndpointWithoutCollapsedStages() throws {
+        let frame=Frame(name:"fractional endpoint",pseudoInertial:true)
+        let dynamics=try SpacecraftDynamics(models:[Constant(acceleration:SIMD3(2,0,0))],context:ForceEvaluationContext(frame:frame,origin:.earth))
+        for fraction in [0.5,0.2693127683618036,0.99] {
+            for sign in [1.0,-1.0] {
+                let start=Epoch(offsetFromJ2000:try Duration(wholeSeconds:800000000,fractionalSecond:fraction))
+                let seed=SpacecraftState(cartesian:try CartesianState(epoch:start,frame:frame,position:.zero,velocity:SIMD3(3,0,0)),mass:try SpacecraftMassProperties(dryMassKilograms:1))
+                let options=try NumericalPropagationOptions(initialStepSeconds:0.1,maximumStepSeconds:2)
+                let model=try NumericalSpacecraftPropagator(initialState:seed,dynamics:dynamics,options:options)
+                let target=try start.advanced(by:Duration(seconds:sign*0.6))
+                let result=try model.propagation(at:target)
+                #expect(result.state.epoch==target)
+                #expect(abs(result.state.cartesian.position.x-(3*sign*0.6+0.36))<1e-12)
+                #expect(abs(result.state.cartesian.velocity.x-(3+2*sign*0.6))<1e-12)
+            }
+        }
+    }
+
 }