I tracked down the root cause of this bug using a diagnostic Harmony mod to investigate the crash path in hope of perhaps saving the devs some time and effort (and of course mainly because i was also curious about how the calculation works
).
(Full LLM disclosure: Alongside my own brain i used Claude Code with Opus 4.6 to save some time and help with the initial code analysis of the decompiled source and to write the diagnostic mod, but all findings were verified empirically in-game and of course by my own brain too xD.)
Summary
The crash occurs when a porkchop plot cell produces a
near-parabolic flyby of Luna with an extremely small periapsis. The position calculation hits a
1/(1+cos(pi)) = 1/0 = Infinity singularity, and the subsequent quaternion rotation produces
NaN from
-Inf + Inf.
Methodology
1. Decompiled KSA.dll and traced the crash path from the exception through
CreateFromStateCci ->
CreateFromCalcData ->
CalculateEscapePatch ->
CalculatePatch ->
ComputeCompleteTrajectory ->
BuildFlightPlan
2. Built a Harmony mod that patches
CalculateEscapePatch,
CreateFromCalcData, and
CreateFromStateCci to log intermediate values when NaN is detected
3. Reproduced the crash multiple times with LEO-to-Luna Hohmann transfers (2-14 day range)
The Bug Chain
Step 1: Porkchop cell creates a near-parabolic Luna encounter
For certain departure/arrival time combinations, the Lambert solver produces a transfer orbit that encounters Luna with near-zero relative velocity. The resulting encounter orbit around Luna is classified as
parabolic (
e = 1.0):
Code:
PrimaryBody=Luna OrbitType=PARABOLIC
e = 1.000000000000406
p = 4.785e-06 m (periapsis = 2.4 micrometers!)
The periapsis of ~2.4 micrometers means this is essentially a radial trajectory through Luna's center - physically impossible but not caught before the escape calculation runs.
Step 2: SOI true anomaly approaches pi (the parabolic asymptote)
Code:
TrueAnomaly@SOI = 3.141592260 (pi = 3.141592654)
SOI TA - pi = -3.93e-7
Step 3: Newton solver (SolveForTheta) fails to converge
The Barker equation solver starts at
D=1 and needs to reach
D~1.7e10. With only 50 iterations and cubic convergence behavior, it does not converge:
Code:
SolveForTheta: D = 1.704e+10, converged = False
Step 4: True anomaly becomes indistinguishable from pi in double precision
Code:
ta = 2 * atan(1.7e10) = 3.141592653472
cos(ta) = -1.000000000000000 (exactly -1.0 in IEEE 754)
sin(ta) = 1.174e-10
1 + cos(ta) = 0.0 (exactly zero!)
Step 5: Position formula hits singularity (THE ROOT CAUSE)
In
Parabolic.GetPositionOrb(TrueAnomaly ta) (Orbit.cs):
C#:
double r = semiLatusRectum / (1.0 + Math.Cos(ta)); // p / 0 = Infinity!
double x = Math.Cos(ta) * r; // -1 * Inf = -Inf
double y = Math.Sin(ta) * r; // ~0 * Inf = Inf (not NaN yet)
// posOrb = (-Inf, +Inf, 0)
Step 6: Quaternion rotation converts Infinity to NaN
Code:
posOrb = (-Inf, +Inf, 0) // valid Infinity in orbital frame
posCci = (NaN, Inf, Inf) // after Transform(Orb2ParentCci)
The quaternion rotation mixes components:
posCci.X = -Inf * q_a + Inf * q_b = NaN (IEEE 754:
-Inf + Inf = NaN).
Step 7: Velocity remains valid
The velocity formula in
Parabolic.GetVelocityOrb uses
(1 + cos(ta)) as a
multiplier (producing ~0), not as a
divisor (which would produce Infinity):
Code:
velOrb = (-0.119, -0.0004, 0) // tiny but finite
velCci = valid // stays valid through rotation
This is why the error shows
PositionCci = NaN but
VelocityCci remains numeric.
Step 8: Crash
CreateFromStateCci receives NaN position, computes NaN eccentricity, throws "unclassifiable eccentricity: NaN".
Diagnostic Log Output (verbatim)
Code:
[TPD] === PARABOLIC CalculateEscapePatch (#167) ===
[TPD] PrimaryBody=Luna e=1,000000000000406E+000
[TPD] SOI TA=3,141592260187849 (pi=3,141592653589793)
[TPD] SOI TA - pi = -3,934019439100211E-007
[TPD] soiTime=395005,864118s elapsed=48436,683495s
[TPD] p=4,7852974777E-006 mu=4,9002710600E+012 t=48436,683495
[TPD] SolveForTheta: D=1,703861728369721E+010 converged=False at iter=-1
[TPD] ta=2*atan(D)=3,141592653472413
[TPD] cos(ta)=-1,000000000000000E+000 sin(ta)=1,173804501446189E-010
[TPD] 1+cos(ta)=0,000000000000000E+000
[TPD] r=p/(1+cos)=Infinity
[TPD] posOrb=(-Inf, +Inf, 0)
[TPD] posOrb hasInf=True hasNaN=False
[TPD] posCci after Transform=(NaN, Inf, Inf)
[TPD] posCci hasNaN=True
[TPD] ROOT CAUSE CONFIRMED: posOrb has Infinity components,
[TPD] quaternion Transform mixes -Inf + Inf = NaN!
[TPD] velOrb=(-1,1878219657E-001, -4,1052003222E-004, 0)
[TPD] velOrb hasInf=False hasNaN=False
[TPD] velCci=valid
Why It Doesn't Always Crash
In a successful (non-crashing) run, the diagnostic mod still detected
562 parabolic escape calculations in a single porkchop plot - but none of them hit the singularity. The difference is the periapsis size:
| Crash case | Non-crash cases |
|---|
| Periapsis | 2.4e-6 m (micrometers) | 4.9 - 92 m |
| SolveForTheta | Not converged, D=1.7e10 | Converged, D=2000-9000 |
| 1+cos(ta) | exactly 0.0 | 2.2e-8 to 4.2e-7 |
| Result | Infinity -> NaN -> crash | Large but finite -> OK |
This explains the ~10% crash rate: every LEO-to-Luna porkchop has hundreds of near-parabolic cells, but only occasionally does one produce a periapsis so small (micrometers) that the Newton solver can't converge and
cos(ta) rounds to exactly -1.0.
Contributing Factors
1. An orbit with periapsis of 2.4 micrometers around Luna (radius 1,737 km) is physically an impact, but
CalculateEscapePatch runs before any impact check catches it.
2.
Parabolic.GetRemainingTimeTo returns
GetTimeFromPeTo(ta) (time from periapsis to angle) instead of the actual remaining time from current position. This gives a wrong but non-NaN SOI encounter time, which can push the evaluation deeper toward the asymptote.
3.
GetPositionOrb doesn't check for the
1+cos(ta) -> 0 singularity.
4.
TransferTask.WorkerTask.ThreadPoolCallback has no try/catch, so a single NaN in one of the 85,264 porkchop cells crashes the entire worker thread.
Possible Fixes
The simplest fix would be in
CalculateEscapePatch: if the orbit is parabolic (or near-parabolic)
and the periapsis is below the parent body's radius (i.e., it's an impact trajectory), skip the escape calculation and mark it as
PatchTransition.Impact instead.
Alternatively,
GetPositionOrb could clamp
1/(1+cos(ta)) to avoid Infinity when
ta is within machine epsilon of pi.
Proof-of-Concept Fix
I wrote a Harmony prefix patch on
CalculateEscapePatch that skips the escape calculation when the orbit's periapsis is below the parent body's radius (i.e., it's an impact trajectory, not an escape):
C#:
[HarmonyPatch(typeof(FlightPlan), nameof(FlightPlan.CalculateEscapePatch))]
static class NaNFix
{
static bool Prefix(PatchedConic patch, ref PatchedConic? nextPatch, ref bool __result)
{
double ecc = patch.Orbit.Eccentricity;
if (ecc < 1.0)
return true; // run original
double periapsis = patch.Orbit.Periapsis;
if (!(patch.PrimaryBody is Astronomical body))
return true;
if (periapsis >= body.MeanRadius)
return true; // periapsis above surface, legitimate escape
// Periapsis below surface = impact trajectory, not escape.
// Skip to avoid r = p/(1+cos(pi)) = Inf -> NaN singularity.
patch.EndTransition = PatchTransition.Final;
nextPatch = null;
__result = false;
return false; // skip original
}
}
Verification
Tested by running multiple LEO-to-Luna Hohmann transfers.
- The fix intercepted ~10 Luna impact trajectories per porkchop plot (periapsis 70-82 km vs Luna radius 1,737 km)
- 4,358 total parabolic escape events per run, ~10 caught by the fix
The diagnostic mod (DiagnosticPatches.cs + NaNFix.cs) is attached to this comment.
