Skip to content

Fix moon orbit generation #6330

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
sam52796 wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Fix moon orbit generation #6330

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
995ca33
Fix underflow in Hill Radius calculation
sam52796 Apr 20, 2026
53ae0b1
Ensure binary stars are in sensible orbits
sam52796 Apr 20, 2026
86dbed5
Don't create moons with 0kg mass
sam52796 Apr 21, 2026
c8b97ed
Move fixed-point PI into precision-specific template constants
sam52796 Apr 22, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
23 changes: 21 additions & 2 deletions src/fixed.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@ class fixedf {
public:
static const int FRAC = FRAC_BITS;
static const Uint64 MASK = (Uint64(1UL) << FRAC_BITS) - 1;
static const fixedf PI;

fixedf() :
v(0) {}
Expand Down Expand Up @@ -224,19 +225,37 @@ class fixedf {
return (fixedf(root));
}

static fixedf CubeRootOf(const fixedf &a)
static fixedf CubeRootOf(const fixedf &a_orig)

@sturnclaw sturnclaw Apr 21, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I do not have time at current to fully dig into the codepath and analyse the numerical inaccuracy at play here (and the effects of this change on its stability) - however, I'd strongly prefer if the rescaling version of this function be moved into a separate function, e.g. fixedf::CubeRootOfRescale. The automatic rescaling may be undesired behavior in codepaths that operate on values with known input bounds, and I'd prefer it to be opt-in.

That being said, I appreciate your tackling of this issue and I'm looking forward to reviewing it in-depth!

@sam52796 sam52796 Apr 21, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you aware that the cube root function is only used in one other place, and right next to that place there is a comment implying that underflow is a persistent problem for the function? In fact, that comment is what gave me the idea to use the bit shifting approach in the first place, since it explains how to use bit shifting to fix the problem.

I say this because I think it's better not to clutter up the codebase with duplicate library functions where the caller is expected to know which version to call. I realise that sometimes it's unavoidable when you really don't want to risk breaking things, but this doesn't seem to be one of those times.

@sturnclaw sturnclaw Apr 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are you aware that the cube root function is only used in one other place, and right next to that place there is a comment implying that underflow is a persistent problem for the function?

Yes - I wrote that comment.

I say this because I think it's better not to clutter up the codebase with duplicate library functions where the caller is expected to know which version to call.

I do agree that duplicating functionality is to be avoided. If this function has noticeable loss-of-precision relative to the original, then there is a material difference between the two functions and thus a reason to separate it as a separate entry point.

After having built a test harness for both versions of the function, the proposed implementation has a worst-case error magnitude increase of 20% when the while (a > fixedf(1024, 1)) branch is taken, due to discarding up to 24 bits of fractional precision.

The good news is that replacing that branch with if (a > fixedf(1024, 1)) x = a * fixedf(2, 3); as the initial estimate experimentally fixes a noticed issue where the x * x term can catastrophically overflow (with an input value of 0x4299E58D00000000L), while retaining the same error magnitude as the original version of the function.

When that is resolved, I agree with your intent that the implementation is a direct upgrade in terms of numerical accuracy and stability and retract my initial recommendation.

@sturnclaw sturnclaw Apr 22, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For clarity, the while (a > fixedf(1024, 1)) branch did not present any reduction in computation error (with regards to a double-precision floating-point cube root using std::pow) compared to computing the N-R method cube root of the original value directly. Additionally, this increase in error metric was observed in both 16.48 fixed-point as well as 32.32 fixed-point representations.

The only noticed outliers as mentioned above were specific values triggering an overflow that caused the first iteration of the method to wildly overshoot beyond the ability of 47 more iterations to converge. Changing the initial estimate of x to a rough estimate of the first iteration result allows large input values to nominally remain within the allowed limits of the fixed-point representation.

{
fixedf a = a_orig;
if (a == fixedf(0, 1)) return fixedf(0, 1);

@sturnclaw sturnclaw Apr 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sanity-check: should this be a <= fixedf(0, 1) instead? I believe computing the N-th root of a negative number is invalid by definition until you start involving complex numbers...


// In order to prevent over or underflow in the (x * x) calculation
// below, we first scale the parameter into a reasonable range
int k = 0;
while (a.v > (Sint64(1) << 60)) { a >>= 3; k++; }
while (a < fixedf(1,1024)) { a <<= 3; k--; }

/* NR method. XXX very bad initial estimate (we get there in
* the end... XXX */
fixedf x = a;
for (int i = 0; i < 48; i++)
x = fixedf(1, 3) * ((a / (x * x)) + 2 * x);
Comment thread
sturnclaw marked this conversation as resolved.
return x;

// We've got the cube root of the normalized parameter, now
// un-scale it back to get the real answer
return k >= 0 ? x << k : x >> -k;
}

Sint64 v;
};

typedef fixedf<32> fixed;

template<>
constexpr fixed fixed::PI = fixed(13493037705LL);

template<>
constexpr fixedf<48> fixedf<48>::PI = fixedf<48>(884279719003555LL);

#endif /* _FIXED_H */
23 changes: 23 additions & 0 deletions src/galaxy/CustomSystem.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1106,7 +1106,30 @@ void CustomSystemsDatabase::RunLuaSystemSanityChecks(CustomSystem *csys)
body->bodyData.m_volatileGas = rand.NormFixed(fixed(1050, 1000), fixed(8000, 1000)).Abs();
body->bodyData.m_atmosOxidizing = rand.NormFixed(fixed(0, 1), fixed(300, 1000)).Abs();
}
}

// This loop runs after the one above, to ensure that system bodies are all set up first.
// Now we ensure that any binary star systems have orbits set correctly.
for (CustomSystemBody *body : csys->bodies) {

if ((body->bodyData.m_type == SystemBody::TYPE_GRAVPOINT) && (body->children.size() == 2)) {
SystemBodyData *star_A = &body->children[0]->bodyData;
SystemBodyData *star_B = &body->children[1]->bodyData;

if ((star_A->m_type >= SystemBodyType::TYPE_STAR_MIN) && (star_A->m_type <= SystemBodyType::TYPE_STAR_MAX) &&
(star_B->m_type >= SystemBodyType::TYPE_STAR_MIN) && (star_B->m_type <= SystemBodyType::TYPE_STAR_MAX)) {

@sturnclaw sturnclaw Apr 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I believe this can be more succinctly expressed as star_A->GetSuperType() == SystemBodyType::SUPERTYPE_STAR etc.

In general the fixup-at-load approach is good. Ideally an extension to the editor for creating binary (and trinary) orbits from multiple bodies would be nice, though I suspect that would involve non-trivial engineering to track and express that a set of bodies are intended to be in a binary orbit.


// This is a binary star system. Ensure that the stars are opposite each other in their orbits

star_B->m_orbitalOffset = star_A->m_orbitalOffset;
star_B->m_inclination = star_A->m_inclination;
star_B->m_argOfPeriapsis = star_A->m_argOfPeriapsis;

star_B->m_orbitalPhaseAtStart = star_A->m_orbitalPhaseAtStart + fixed::PI;
if (star_B->m_orbitalPhaseAtStart > fixed(2, 1) * fixed::PI)
star_B->m_orbitalPhaseAtStart -= fixed(2, 1) * fixed::PI;
}
}
}
}

Expand Down
42 changes: 22 additions & 20 deletions src/galaxy/StarSystemGenerator.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@ static const fixed SAFE_DIST_FROM_BINARY = fixed(5, 1);
// very crudely
static const fixed AU_SOL_RADIUS = fixed(305, 65536);
static const fixed AU_EARTH_RADIUS = fixed(3, 65536); // XXX Duplication from StarSystem.cpp
static const fixed FIXED_PI = fixed(103993, 33102); // XXX Duplication from StarSystem.cpp

@sturnclaw sturnclaw Apr 22, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nitpick: if this is duplicated, the file it is duplicated from should be updated to also use the new fixed::PI values.

static const double CELSIUS = 273.15;
static const fixed ONEEUMASS = fixed::FromDouble(1);
static const fixed TWOHUNDREDEUMASSES = fixed::FromDouble(200.0);
Expand Down Expand Up @@ -385,8 +384,9 @@ fixedf<48> StarSystemLegacyGeneratorBase::CalcHillRadius(SystemBody *sbody) cons
fixedp a = sbody->GetSemiMajorAxisAsFixed();
fixedp e = sbody->GetEccentricityAsFixed();
fixedp pe = a * (fixedp(1, 1) - e); // periapsis in higher precision
fixedp mass_ratio = sbody->GetMassAsFixed() / (fixed(3, 1) * mprimary);

return pe * fixedp::CubeRootOf(sbody->GetMassAsFixed() / (fixed(3, 1) * mprimary));
return pe * fixedp::CubeRootOf(mass_ratio);

//fixed hr = semiMajorAxis*(fixed(1,1) - eccentricity) *
// fixedcuberoot(mass / (3*mprimary));
Expand Down Expand Up @@ -976,7 +976,6 @@ fixed StarSystemRandomGenerator::CalcBodySatelliteShellDensity(Random &rand, Sys
// gameplay purposes instead of fully representing reality
fixedf<48> hillSphereRad = CalcHillRadius(primary) * fixedf<48>(1, 4);
discMax = std::min(discMax, fixed(hillSphereRad));

return get_disc_density(primary, discMin, discMax, fixed(1, 500));
}
}
Expand Down Expand Up @@ -1129,6 +1128,9 @@ SystemBody *StarSystemRandomGenerator::MakeBodyInOrbitSlice(Random &rand, StarSy

mass *= discDensity;

if (mass == fixed(0, 1))
return nullptr; // We used our best math, but this body just turned out too small, sorry.

if (mass.v < 0) { // hack around overflow
Output("WARNING: planetary mass has overflowed! (child %d of %s)\n", primary->GetNumChildren(), primary->GetName().c_str());
mass = fixed(Sint64(0x7fFFffFFffFFffFFull));
Expand All @@ -1145,19 +1147,19 @@ SystemBody *StarSystemRandomGenerator::MakeBodyInOrbitSlice(Random &rand, StarSy
planet->m_rotationPeriod = fixed(rand.Int32(1, 200), 24);

// longitude of ascending node
planet->m_orbitalOffset = rand.Fixed() * 2 * FIXED_PI;
planet->m_orbitalOffset = rand.Fixed() * 2 * fixed::PI;
// inclination in the hemisphere above the equator, low probability of high-inclination orbits
fixed incl_scale = rand.Fixed() * fixed(666, 1000);
planet->m_inclination = rand.NormFixed().Abs() * incl_scale * FIXED_PI * fixed(1, 2);
planet->m_inclination = rand.NormFixed().Abs() * incl_scale * fixed::PI * fixed(1, 2);
// argument of periapsis, interval -PI .. PI
planet->m_argOfPeriapsis = rand.NormFixed() * FIXED_PI;
planet->m_argOfPeriapsis = rand.NormFixed() * fixed::PI;

// rare chance of reversed orbit
if (rand.Fixed() < fixed(1, 20))
planet->m_inclination = FIXED_PI - planet->m_inclination;
planet->m_inclination = fixed::PI - planet->m_inclination;

// true anomaly as rotation beyond periapsis
planet->m_orbitalPhaseAtStart = rand.Fixed() * 2 * FIXED_PI;
planet->m_orbitalPhaseAtStart = rand.Fixed() * 2 * fixed::PI;

planet->SetOrbitFromParameters();
planet->SetAtmFromParameters();
Expand Down Expand Up @@ -1269,10 +1271,10 @@ void StarSystemRandomGenerator::MakeBinaryPair(SystemBody *a, SystemBody *b, fix
b->m_orbit.SetPlane(matrix3x3d::RotateY(rotY - M_PI) * matrix3x3d::RotateX(rotX));

// store orbit parameters for later use to be accessible in other way than by rotMatrix
b->m_orbitalPhaseAtStart = b->m_orbitalPhaseAtStart + FIXED_PI;
b->m_orbitalPhaseAtStart = b->m_orbitalPhaseAtStart > 2 * FIXED_PI ? b->m_orbitalPhaseAtStart - 2 * FIXED_PI : b->m_orbitalPhaseAtStart;
a->m_orbitalPhaseAtStart = a->m_orbitalPhaseAtStart > 2 * FIXED_PI ? a->m_orbitalPhaseAtStart - 2 * FIXED_PI : a->m_orbitalPhaseAtStart;
a->m_orbitalPhaseAtStart = a->m_orbitalPhaseAtStart < 0 ? a->m_orbitalPhaseAtStart + 2 * FIXED_PI : a->m_orbitalPhaseAtStart;
b->m_orbitalPhaseAtStart = b->m_orbitalPhaseAtStart + fixed::PI;
b->m_orbitalPhaseAtStart = b->m_orbitalPhaseAtStart > 2 * fixed::PI ? b->m_orbitalPhaseAtStart - 2 * fixed::PI : b->m_orbitalPhaseAtStart;
a->m_orbitalPhaseAtStart = a->m_orbitalPhaseAtStart > 2 * fixed::PI ? a->m_orbitalPhaseAtStart - 2 * fixed::PI : a->m_orbitalPhaseAtStart;
a->m_orbitalPhaseAtStart = a->m_orbitalPhaseAtStart < 0 ? a->m_orbitalPhaseAtStart + 2 * fixed::PI : a->m_orbitalPhaseAtStart;
b->m_orbitalOffset = fixed(int(round(rotY * 10000)), 10000);
a->m_orbitalOffset = fixed(int(round(rotY * 10000)), 10000);

Expand Down Expand Up @@ -1448,8 +1450,8 @@ void PopulateStarSystemGenerator::PositionSettlementOnPlanet(SystemBody *sbody,

// store latitude and longitude to equivalent orbital parameters to
// be accessible easier
sbody->m_inclination = latitude * FIXED_PI * fixed(1, 2);
sbody->m_orbitalOffset = longitude * FIXED_PI * 2;
sbody->m_inclination = latitude * fixed::PI * fixed(1, 2);
sbody->m_orbitalOffset = longitude * fixed::PI * 2;

sbody->SetOrbitFromParameters();
}
Expand Down Expand Up @@ -1712,12 +1714,12 @@ void PopulateStarSystemGenerator::PopulateAddStations(SystemBody *sbody, StarSys
shells[1] = innerOrbit + ((orbMaxS - innerOrbit) * fixed(1, 2)); // med
shells[2] = orbMaxS; // high

shellIncl[0] = rand.NormFixed() * FIXED_PI;
shellIncl[1] = rand.NormFixed() * FIXED_PI;
shellIncl[2] = rand.NormFixed() * FIXED_PI;
shellIncl[0] = rand.NormFixed() * fixed::PI;
shellIncl[1] = rand.NormFixed() * fixed::PI;
shellIncl[2] = rand.NormFixed() * fixed::PI;
} else {
shells[0] = shells[1] = shells[2] = innerOrbit;
shellIncl[0] = shellIncl[1] = shellIncl[2] = rand.NormFixed() * FIXED_PI;
shellIncl[0] = shellIncl[1] = shellIncl[2] = rand.NormFixed() * fixed::PI;
}

Uint32 orbitIdx = 0;
Expand Down Expand Up @@ -1749,9 +1751,9 @@ void PopulateStarSystemGenerator::PopulateAddStations(SystemBody *sbody, StarSys
// slightly random min/max orbital distance
sp->m_eccentricity = rand.NormFixed().Abs() * fixed(1, 8);
// perturb the orbital plane to avoid all stations falling in line with each other
sp->m_inclination = currOrbitIncl + rand.NormFixed() * fixed(1, 4) * FIXED_PI;
sp->m_inclination = currOrbitIncl + rand.NormFixed() * fixed(1, 4) * fixed::PI;
// station spacing around the primary body
sp->m_argOfPeriapsis = rand.Fixed() * FIXED_PI * 2;
sp->m_argOfPeriapsis = rand.Fixed() * fixed::PI * 2;
// TODO: no axial tilt for stations / axial tilt in general is strangely modeled
sp->m_axialTilt = fixed();

Expand Down
Loading