GameCult/Aetheria
Watch
9
0
Fork
You've already forked Aetheria
 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Well law fix batch: zone node for several massive roots, belt extent, series from tau 0.4, Unity parity on the multi-root fixture

Browse source 
One zone node at the massive roots' centre of mass owns the only global term (SubzoneReach is the one place the subzone shape is decided: isolated-subzones). WellNode gains Extent and Anchor; a belt orbit's extent is its outermost asteroid. GravityField.SeriesFrom = 0.4 is emitted into the shader template and Gravity.cginc. Tests: the spec's WellLawTests; the old seam test is replaced by TheBandPotentialIsExactInSinglePrecision. GravityFieldParityVerify runs the multi-root fixture (second sun, moon inside the planet's Roche limit, far window): PASS, largest drawn tau 1.228, max error 0.01028 u. The dotnet suite has not been run at this commit.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-16 21:14:58 +02:00
commit 250b61be2c
9 changed files with 412 additions and 60 deletions
Download patch file Download diff file Expand all files Collapse all files

71
Assets/Scripts/Editor/GravityFieldParityVerify.cs
View file

@ -1,4 +1,5 @@
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using GameCult.Caching;
@ -27,14 +28,16 @@ namespace Aetheria.EditorTools
const string Tag = "[GravityFieldParityVerify]";
const int Size = 256;
// Measured on this tree (the "max error" line this run prints): the largest |texel + GetHeight| over the windows of the three
// frames is 0.0108 u, in the moved frame's planetoid window (float32 on both sides, the GPU's log and the closed form's
// frames is 0.0103 u, in the moved-again frame's sun window (float32 on both sides, the GPU's log and the closed form's
// single-precision floor of 1.5e-5 per unit strength against the CPU's, the quad's interpolated world position). The
// tolerance is 2.8x that. A band whose terms were drawn from the wrong edge (reference, inner or outer) moves the
// tolerance is 2.9x that. A band whose terms were drawn from the wrong edge (reference, inner or outer) moves the
// drawn height by whole units.
const float Tolerance = 0.03f;
// The worst angle between the texture's central difference and GetForce (B28: 1.0 degrees median zone, 2.9 worst).
const float GradientDegrees = 3f;
static float _maxError, _maxAngle;
// Below this force (u/s^2) the drawn height differs between neighbouring texels by less than float32 resolves.
const float FlatForce = 1e-3f;
static float _maxError, _maxAngle, _maxTau;
static bool _ok = true;
static void OnLog(string message, string stack, LogType type)
@ -70,7 +73,7 @@ namespace Aetheria.EditorTools
var material = AssetDatabase.LoadAssetAtPath<Material>("Assets/Materials/Brushes/Gravity Field.mat")
?? throw new InvalidOperationException("No Gravity Field.mat.");
var zone = BuildZone(settings, work, out var orbits);
var zone = BuildZone(settings, work, out var orbits, multiRoot: true);
GravityLowering.Lower(zone, 0f, ref wells);
field = GravityLowering.CreateField(material);
@ -82,13 +85,17 @@ namespace Aetheria.EditorTools
CheckFrame("built", zone, orbits);
orbits[1].Position = Around(orbits[0], orbits[1], .31f);
orbits[2].Position = Around(orbits[1], orbits[2], .77f);
orbits[5].Position = Around(orbits[1], orbits[5], .41f);
GravityLowering.Lower(zone, 0f, ref wells);
CheckFrame("moved", zone, orbits);
orbits[0].Position = new float2(210f, -95f);
orbits[1].Position = Around(orbits[0], orbits[1], .52f);
orbits[2].Position = Around(orbits[1], orbits[2], .05f);
orbits[5].Position = Around(orbits[1], orbits[5], .93f);
GravityLowering.Lower(zone, 0f, ref wells);
CheckFrame("moved-again", zone, orbits);
if (_maxTau <= 1f) Fail("series-branch", $"the largest tau over the drawn bounded terms is {_maxTau:G4}; no window draws the series branch");
Debug.Log($"{Tag} largest drawn tau {_maxTau:G4}");
Debug.Log($"{Tag} max error {_maxError:G4} u (tolerance {Tolerance:G4}); max gradient angle {_maxAngle:F2} degrees (limit {GradientDegrees:F0})");
CheckWiring(settings);
@ -117,12 +124,23 @@ namespace Aetheria.EditorTools
static Vector2 At(Orbit orbit) => new Vector2(orbit.Position.x, orbit.Position.y);
// Orbits 0-2 are the sun, the planet and the planetoid, each window centred where the body is now; orbit 3 is the
// body of no strength. The zone window holds all four.
// body of no strength, orbit 4 the second sun and orbit 5 the planet's moon inside its Roche limit (whose band has
// tau above 1: the series branch of the shader). The zone window holds all of them; the far window is centred on the
// zone node, so its band (nothing inside A, the zone's mass beyond Beta A) is drawn. A window draws a term when
// the term's disc reaches it; the largest tau over the drawn terms must exceed 1 or the series branch went unchecked.
static void CheckFrame(string frame, Zone zone, Orbit[] orbits)
{
CheckTexels(frame + "-planetoid", zone, At(orbits[2]), 4f * Core(zone, 29.2f));
CheckTexels(frame + "-planet", zone, At(orbits[1]), 4f * Core(zone, 273.9f));
CheckTexels(frame + "-sun", zone, At(orbits[0]), 4f * Core(zone, 14152f));
CheckTexels(frame + "-zone", zone, Vector2.zero, 1.05f * zone.Pack.Radius);
var centre = new Vector2(zone.Wells.Centre.x, zone.Wells.Centre.y);
var span = 0f;
foreach (var root in new[] { orbits[0], orbits[4] })
span = Math.Max(span, (At(root) - centre).magnitude + zone.Wells.Reach(zone.Orbits.First(o => o.Value == root).Key));
// The wells have their own windows: out here a texel is 40 u wide, so only the ring from 0.9 A outward, where the
// zone node's band rises and the zone becomes one body, is measured.
CheckTexels(frame + "-far", zone, centre, 1.6f * span, .9f * span);
CheckGradient(zone, orbits.Select(At).ToArray());
}
@ -131,18 +149,28 @@ namespace Aetheria.EditorTools
// c = 0 and R = 0), at (1700, 300): a family, so the band terms are drawn (the planet's and the sun's cusps and
// fades), off both axes and asymmetric, so the zone's wells sit where the verify says before any Update and a
// mirrored or swapped axis cannot agree with the sim. The sun and the zero-mass body are fixed roots.
internal static Zone BuildZone(GameSettings settings, string work, out Orbit[] orbits, float radius = 4182f)
internal static Zone BuildZone(GameSettings settings, string work, out Orbit[] orbits, float radius = 4182f, bool multiRoot = false)
{
var cache = AetheriaStores.Open(Path.Combine(work, "catalog.cc"), Path.Combine(work, "run.cc"), catalogWritable: true);
cache.Upsert(new VerseGrammar { Revision = 1 });
cache.FlushAsync().Wait();
var items = new ItemManager(cache, new ProvenanceLedger(), settings.GameplaySettings, _ => { });
var positions = new[] { new Vector2(-130f, 85f), new Vector2(610f, -240f), new Vector2(670f, -270f), new Vector2(1700f, 300f) };
var masses = new[] { 14152f, 273.9f, 29.2f, 0f };
var parents = new[] { -1, 0, 1, -1 };
var keys = new CultRecordKey[positions.Length];
var positions = new List<Vector2> { new Vector2(-130f, 85f), new Vector2(610f, -240f), new Vector2(670f, -270f), new Vector2(1700f, 300f) };
var masses = new List<float> { 14152f, 273.9f, 29.2f, 0f };
var parents = new List<int> { -1, 0, 1, -1 };
if (multiRoot)
{
// A second sun (a second subzone, so the zone node exists) and a moon of the planet 0.6 of its Roche limit away.
positions.Add(new Vector2(-2400f, -1900f));
masses.Add(3000f);
parents.Add(-1);
positions.Add(positions[1] + new Vector2(-.6f * settings.PlanetSettings.RocheLimit.Evaluate(273.9f), 0f));
masses.Add(1f);
parents.Add(1);
}
var keys = new CultRecordKey[positions.Count];
var pack = new ZonePack { Radius = radius };
for (var i = 0; i < positions.Length; i++)
for (var i = 0; i < positions.Count; i++)
{
var data = new OrbitData { FixedPosition = new float2(positions[i].x, positions[i].y) };
if (parents[i] >= 0)
@ -211,13 +239,29 @@ namespace Aetheria.EditorTools
}
}
static void CheckTexels(string window, Zone zone, Vector2 centre, float half)
// The largest tau = (a^2 + eps^2) / (b^2 - a^2) over the bounded terms whose disc reaches the window's square.
static void NoteTau(Zone zone, Vector2 centre, float half)
{
foreach (var term in zone.Wells.Terms)
{
if (float.IsPositiveInfinity(term.Outer)) continue;
var at = term.Node.IsSet() ? zone.Orbits[term.Node].Position : zone.Wells.Centre;
if (Math.Abs(at.x - centre.x) > half + term.Outer || Math.Abs(at.y - centre.y) > half + term.Outer) continue;
_maxTau = Math.Max(_maxTau, (term.Inner * term.Inner + term.Softening * term.Softening) / (term.Outer * term.Outer - term.Inner * term.Inner));
}
}
static void CheckTexels(string window, Zone zone, Vector2 centre, float half, float ringFrom = 0f)
{
Render(centre, half, out var depth, out var world);
NoteTau(zone, centre, half);
var worst = 0f;
var live = 0;
var measured = 0;
for (var i = 0; i < depth.Length; i++)
{
if ((world[i] - centre).magnitude < ringFrom) continue;
measured++;
var expected = -zone.GetHeight(new float2(world[i].x, world[i].y));
if (float.IsNaN(depth[i])) { Fail($"{window}-texels", $"texel {i}: NaN"); return; }
var error = Math.Abs(depth[i] - expected);
@ -231,7 +275,7 @@ namespace Aetheria.EditorTools
}
_maxError = Math.Max(_maxError, worst);
// A window that draws nothing agrees with a zone that is flat: the check must be live.
if (live < depth.Length / 2) Fail($"{window}-texels", $"only {live} of {depth.Length} texels are deeper than 1 u; the check is vacuous");
if (live < measured / 2) Fail($"{window}-texels", $"only {live} of {measured} texels are deeper than 1 u; the check is vacuous");
Debug.Log($"{Tag} {window}: {depth.Length} texels, max error {worst:G4} u");
}
@ -255,6 +299,7 @@ namespace Aetheria.EditorTools
var det = au.x * av.y - au.y * av.x;
var g = new Vector2((du * av.y - dv * au.y) / det, (dv * au.x - du * av.x) / det);
var f = zone.GetForce(new float2(here.x, here.y));
if (CultMath.math.length(f) < FlatForce) continue; // between the subzones the ground is flat: no direction to compare
var angle = Vector2.Angle(g, new Vector2(f.x, f.y));
if (float.IsNaN(angle) || angle > GradientDegrees)
{

7
Assets/Scripts/ServerShared/GravityField.cs
View file

@ -28,7 +28,10 @@ public static class GravityField
/// <summary>A ring member of a bodyless node reaches RingReach x its distance: whole out to the ring's diameter.</summary>
public const float RingReach = 4f;
/// <summary>The series branch's length: its ratio is below 1/2, so 40 terms reach float32's floor.</summary>
/// <summary>The band is summed as a series for tau above this. Measured (map B54): float32 within 1.0e-6 per unit strength everywhere with SeriesTerms terms; the closed form is float32-clean only to about 0.45.</summary>
public const float SeriesFrom = .4f;
/// <summary>The series branch's length: its ratio is below 1/1.4, so 40 terms reach float32's floor.</summary>
public const int SeriesTerms = 40;
/// <summary>S(t) as polynomial coefficients in t, lowest first: t^Order sum_j C(Order-1+j, j) (1-t)^j.</summary>
@ -102,7 +105,7 @@ public static class GravityField
}
// The band's potential at t in [0, 1], 0 at t = 1.
private static float Band(float tau, float t) => tau <= 1f ? BandClosed(tau, t) : BandSeries(tau, t);
private static float Band(float tau, float t) => tau <= SeriesFrom ? BandClosed(tau, t) : BandSeries(tau, t);
private static float BandClosed(float tau, float t)
{

9
Assets/Scripts/ServerShared/GravityShader.cs
View file

@ -19,6 +19,7 @@ public static class GravityShader
.Replace("@N@", GravityField.S.Length.ToString(CultureInfo.InvariantCulture))
.Replace("@Q@", (GravityField.S.Length - 1).ToString(CultureInfo.InvariantCulture))
.Replace("@QLAST@", (GravityField.S.Length - 2).ToString(CultureInfo.InvariantCulture))
.Replace("@SERIESFROM@", GravityField.SeriesFrom.ToString("R", CultureInfo.InvariantCulture))
.Replace("@SERIES@", GravityField.SeriesTerms.ToString(CultureInfo.InvariantCulture))
.Replace("@ONEMINUSS@", Coefficients(GravityField.OneMinusS))
.Replace("@S@", Coefficients(GravityField.S))
@ -54,7 +55,7 @@ float _HeightExaggeration;
static const float GravityOneMinusS[@N@] = { @ONEMINUSS@ };
static const float GravityS[@N@] = { @S@ };
// GravityField.BandClosed: the band's potential at t for tau <= 1, per unit strength.
// GravityField.BandClosed: the band's potential at t for tau <= @SERIESFROM@, per unit strength.
float gravityBandClosed(float tau, float t)
{
float q[@Q@];
@ -71,7 +72,7 @@ float gravityBandClosed(float tau, float t)
return 0.5 * (remainder * log((t + tau) / (1 + tau)) - (whole - part));
}
// GravityField.BandSeries: the band's potential at t for tau > 1, a series in (1 - t) / (1 + tau).
// GravityField.BandSeries: the band's potential at t for tau > @SERIESFROM@, a series in (1 - t) / (1 + tau).
float gravityBandSeries(float tau, float t)
{
float y = 1 - t;
@ -103,9 +104,9 @@ float gravityTermPotential(float softening, float a, float b, float r)
float b2 = b * b;
if (s >= b2) return 0;
float tau = (a2 + e2) / (b2 - a2);
if (s <= a2) return 0.5 * log((s + e2) / (a2 + e2)) + (tau <= 1 ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0));
if (s <= a2) return 0.5 * log((s + e2) / (a2 + e2)) + (tau <= @SERIESFROM@ ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0));
float t = (s - a2) / (b2 - a2);
return tau <= 1 ? gravityBandClosed(tau, t) : gravityBandSeries(tau, t);
return tau <= @SERIESFROM@ ? gravityBandClosed(tau, t) : gravityBandSeries(tau, t);
}
// Zone.GetHeight at world xz: each term's potential minus its own reference, summed, divided by the exaggeration.

101
Assets/Scripts/ServerShared/WellTree.cs
View file

@ -5,19 +5,26 @@ using CultMath;
using static CultMath.math;
using float2 = CultMath.float2;
/// <summary>An orbit as the well tree reads it: where it hangs on the tree and the mass of the well-bearing body on it (0 for none).</summary>
/// <summary>
/// An orbit as the well tree reads it: where it hangs on the tree and the mass of the well-bearing body on it (0 for none).
/// Extent is the outermost distance from the parent of anything riding this orbit's rail (a belt orbit: its outermost
/// asteroid; any other orbit: its Distance). Anchor is a parentless orbit's stationary position, ignored with a parent.
/// </summary>
public readonly struct WellNode
{
public readonly CultRecordKey Key, Parent;
public readonly float Distance, Phase, Mass;
public readonly float Distance, Phase, Mass, Extent;
public readonly float2 Anchor;
public WellNode(CultRecordKey key, CultRecordKey parent, float distance, float phase, float mass)
public WellNode(CultRecordKey key, CultRecordKey parent, float distance, float phase, float mass, float extent, float2 anchor)
{
Key = key;
Parent = parent;
Distance = distance;
Phase = phase;
Mass = mass;
Extent = extent;
Anchor = anchor;
}
}
@ -28,9 +35,12 @@ public readonly struct WellNode
/// Every node F (a body, or a bodyless node with its ring group standing in for the body) contributes at its own position
/// c(m_F) U(eps(m_F), L_F, B_F) + c(M_F) [ U(eps(M_F), max(B_F, Gamma Lam_F), Lam_F) - U(eps(M_F), L_F, B_F) ]
/// m_F is the body's mass, M_F the family's, c and eps the WellStrength and RocheLimit curves. L_F is the outermost
/// child's distance: inside it the node pulls exactly as its body alone. B_F, the support, bounds every descendant term:
/// child's extent: inside it the node pulls exactly as its body alone. B_F, the support, bounds every descendant term:
/// beyond it the family pulls exactly as one body of mass M_F, out to max(B_F, Gamma Lam_F), and is gone at its reach
/// Lam_F. A root's reach is infinite, so the zone's far field is exactly c(M_zone) however the generator split it.
/// Lam_F. A zone's sole massive root has infinite reach, so its far field is exactly c(M_zone). A zone of several massive
/// roots (subzones) has one zone node, which is not an orbit: at the roots' centre of mass it owns the only global term,
/// c(M_zone) U_global less c(M_zone) U(eps, A, Beta A), nothing inside A and the whole zone as one body beyond Beta A;
/// each root is bounded to its own subzone (SubzoneReach).
/// A body with no massive descendants is one term c(m) U(eps(m), Gamma Lam, Lam). Beyond a node's reach its whole
/// subtree contributes exactly nothing, so Force and Potential skip it: the walk is exact, not an approximation.
/// </summary>
@ -62,15 +72,22 @@ public sealed class WellTree
private int[] _roots;
private int[][] _ring;
private int[] _termStart, _termCount;
private float[] _distance, _phase, _reach, _support;
private float[] _distance, _phase, _reach, _support, _extent;
private float2[] _anchor;
private int[] _massiveRoots;
private float _packRadius;
private int _zoneStart, _zoneCount;
private Term[] _terms;
/// <summary>Every term of the zone, grouped by node: what the lowering uploads and the shader sums.</summary>
/// <summary>Where the zone node's terms sit (the family-mass-weighted mean of the massive roots' anchors); meaningful when it has terms.</summary>
public float2 Centre { get; private set; }
/// <summary>Every term of the zone, grouped by node, the zone node's last (its Node is unset: it sits at Centre): what the lowering uploads and the shader sums.</summary>
public IReadOnlyList<Term> Terms => _terms;
private WellTree() { }
/// <summary>A node's reach in its parent's field (+inf for a root); 0 for an orbit that carries no mass.</summary>
/// <summary>A node's reach in its parent's field (+inf for a zone's sole massive root; a subzone's radius for a root of several); 0 for an orbit that carries no mass.</summary>
public float Reach(CultRecordKey node) => _reach[_index[node]];
/// <summary>A node's support B: every term of its descendants lies inside it, so beyond it the family is one body.</summary>
@ -102,8 +119,13 @@ public sealed class WellTree
tree._key = new CultRecordKey[n];
tree._distance = new float[n];
tree._phase = new float[n];
tree._extent = new float[n];
tree._anchor = new float2[n];
tree._packRadius = packRadius;
for (var i = 0; i < n; i++)
{
tree._extent[i] = nodes[i].Extent;
tree._anchor[i] = nodes[i].Anchor;
tree._index[nodes[i].Key] = i;
tree._key[i] = nodes[i].Key;
tree._distance[i] = nodes[i].Distance;
@ -135,6 +157,20 @@ public sealed class WellTree
foreach (var i in deepestFirst)
if (parent[i] >= 0) family[parent[i]] += family[i];
// Massive roots are the zone's subzones; with two or more, their centre of mass is the zone node's position.
var massive = new List<int>();
foreach (var r in roots)
if (family[r] > 0) massive.Add(r);
tree._massiveRoots = massive.ToArray();
var centre = float2.zero;
var total = 0f;
foreach (var r in massive)
{
centre += tree._anchor[r] * family[r];
total += family[r];
}
tree.Centre = total > 0f ? centre / total : float2.zero;
tree._kids = new int[n][];
var ring = tree._ring = new int[n][];
for (var i = 0; i < n; i++)
@ -175,7 +211,7 @@ public sealed class WellTree
}
}
// Children before parents: a family's support needs its children's reaches. lower = L_F, the outermost child's distance.
// Children before parents: a family's support needs its children's reaches. lower = L_F, the outermost child's extent.
private static void BuildReaches(WellTree t, int[] parent, float[] body, float[] family, int[][] ring, int[] order,
out float[] lower, out RingGroup[] group)
{
@ -188,9 +224,9 @@ public sealed class WellTree
foreach (var i in order)
{
var low = 0f;
foreach (var k in t._kids[i]) low = Math.Max(low, d[k]);
foreach (var k in t._kids[i]) low = Math.Max(low, t._extent[k]);
var support = GravityField.Beta * low;
foreach (var k in t._kids[i]) support = Math.Max(support, d[k] + t._reach[k]);
foreach (var k in t._kids[i]) support = Math.Max(support, t._extent[k] + t._reach[k]);
if (family[i] > 0 && body[i] <= 0)
{
// Bodyless: its ring group stands in for a body (it has massive children, or it would carry no mass).
@ -213,7 +249,7 @@ public sealed class WellTree
var floor = GravityField.Beta * support;
var p = parent[i];
float reach;
if (p < 0) reach = float.PositiveInfinity;
if (p < 0) reach = t.SubzoneReach(i, floor);
else if (body[p] <= 0 && ring[p] != null && ring[p].Length >= 2 && Array.IndexOf(ring[p], i) >= 0)
reach = Math.Max(GravityField.RingReach * d[i], floor);
else reach = Math.Max(GravityField.Reach * d[i], floor);
@ -221,6 +257,19 @@ public sealed class WellTree
}
}
// A root's reach. The one place the subzone shape is decided (question multi-root-zone-shape): a sole massive root is the
// whole zone (+inf, so its global term is the far field); with several, isolated-subzones bounds each root to Reach x the
// distance to the nearest other massive root, or its floor (Beta x support) if more. The superposed-inside option is
// this one return: Beta x max(pack radius, max over massive roots of |anchor - Centre| + support).
private float SubzoneReach(int root, float floor)
{
if (_massiveRoots.Length < 2) return float.PositiveInfinity;
var nearest = float.PositiveInfinity;
foreach (var other in _massiveRoots)
if (other != root) nearest = Math.Min(nearest, length(_anchor[other] - _anchor[root]));
return Math.Max(GravityField.Reach * nearest, floor);
}
private void EmitTerms(int[] parent, float[] body, float[] family, int[][] ring, float[] lower, RingGroup[] group,
PlanetSettings settings, float packRadius)
{
@ -259,14 +308,32 @@ public sealed class WellTree
Add(-cf, ef, lower[i], _support[i]);
_termCount[i] = terms.Count - _termStart[i];
}
_zoneStart = terms.Count;
if (_massiveRoots.Length >= 2)
{
var total = 0f;
var span = 0f;
foreach (var r in _massiveRoots)
{
total += family[r];
span = Math.Max(span, length(_anchor[r] - Centre) + _reach[r]);
}
var (c, eps) = Curves(settings, total);
terms.Add(new Term(default, c, eps, float.PositiveInfinity, float.PositiveInfinity,
c * GravityField.Potential(eps, float.PositiveInfinity, float.PositiveInfinity, packRadius)));
if (span > 0f) terms.Add(new Term(default, -c, eps, span, GravityField.Beta * span, 0f));
}
_zoneCount = terms.Count - _zoneStart;
_terms = terms.ToArray();
}
// One node's own terms at distance r from it, toward its centre.
private float PullOf(int node, float r)
private float PullOf(int node, float r) => PullOf(_termStart[node], _termCount[node], r);
private float PullOf(int start, int count, float r)
{
var pull = 0f;
for (var t = _termStart[node]; t < _termStart[node] + _termCount[node]; t++)
for (var t = start; t < start + count; t++)
{
var term = _terms[t];
pull += term.Strength * GravityField.Pull(term.Softening, term.Inner, term.Outer, r);
@ -281,6 +348,9 @@ public sealed class WellTree
public float2 Force(float2 p, Func<CultRecordKey, double, float2> positionAt, double time)
{
var force = float2.zero;
var toZone = Centre - p;
var rz = length(toZone);
if (_zoneCount > 0 && rz > 0f) force += toZone / rz * PullOf(_zoneStart, _zoneCount, rz);
foreach (var root in _roots) ForceFrom(root, p, positionAt, time, ref force);
return force;
}
@ -299,6 +369,9 @@ public sealed class WellTree
public float Potential(float2 p, Func<CultRecordKey, double, float2> positionAt, double time)
{
var sum = 0f;
var rz = length(Centre - p);
for (var t = _zoneStart; t < _zoneStart + _zoneCount; t++)
sum += _terms[t].Strength * GravityField.Potential(_terms[t].Softening, _terms[t].Inner, _terms[t].Outer, rz) - _terms[t].Reference;
foreach (var root in _roots) sum += PotentialFrom(root, p, positionAt, time);
return sum;
}

17
Assets/Scripts/ServerShared/Zone.cs
View file

@ -193,10 +193,23 @@ public class Zone
foreach (var body in PlanetInstances)
if (body.Value.WellStrength > 0)
wellMass.TryAdd(Planets[body.Key].Orbit.Key, body.Value.BodyData.Mass);
// A belt's rail carries its asteroids, which spread beyond the rail's own distance: its extent is the outermost.
var beltExtent = new Dictionary<CultRecordKey, float>();
foreach (var belt in AsteroidBelts)
{
var extent = 0f;
foreach (var asteroid in belt.Value.Data.Asteroids) extent = Math.Max(extent, asteroid.Distance);
beltExtent[Planets[belt.Key].Orbit.Key] = extent;
}
var wellNodes = new List<WellNode>();
foreach (var orbit in Orbits)
wellNodes.Add(new WellNode(orbit.Key, orbit.Value.Data.Parent.Key, orbit.Value.Data.Distance, orbit.Value.Data.Phase,
wellMass.TryGetValue(orbit.Key, out var mass) ? mass : 0f));
{
var data = orbit.Value.Data;
wellNodes.Add(new WellNode(orbit.Key, data.Parent.Key, data.Distance, data.Phase,
wellMass.TryGetValue(orbit.Key, out var mass) ? mass : 0f,
beltExtent.TryGetValue(orbit.Key, out var reach) ? Math.Max(reach, data.Distance) : data.Distance,
data.FixedPosition + OrbitData.Evaluate(frac(data.Phase)) * data.Distance)); // a root's stationary place: Offset at an infinite period
}
Wells = WellTree.Build(wellNodes, settings, pack.Radius);
foreach (var orbit in Orbits)
orbit.Value.Period = orbit.Value.Data.Parent.IsSet()

2
Assets/Scripts/Zone Display/GravityLowering.cs
View file

@ -51,7 +51,7 @@ public static class GravityLowering
for (var t = 0; t < count; t++)
{
var term = zone.Wells.Terms[t];
var at = DrawAhead.Orbit(zone, term.Node, lead);
var at = term.Node.IsSet() ? DrawAhead.Orbit(zone, term.Node, lead) : zone.Wells.Centre; // an unset node is the zone node, stationary at Centre
_records[t] = new TermRecord
{
Position = new Vector2(at.x, at.y),

8
Assets/Shaders/Gravity.cginc
View file

@ -25,7 +25,7 @@ float _HeightExaggeration;
static const float GravityOneMinusS[8] = { 1, -0, -0, -0, -35, 84, -70, 20 };
static const float GravityS[8] = { 0, 0, 0, 0, 35, -84, 70, -20 };
// GravityField.BandClosed: the band's potential at t for tau <= 1, per unit strength.
// GravityField.BandClosed: the band's potential at t for tau <= 0.4, per unit strength.
float gravityBandClosed(float tau, float t)
{
float q[7];
@ -42,7 +42,7 @@ float gravityBandClosed(float tau, float t)
return 0.5 * (remainder * log((t + tau) / (1 + tau)) - (whole - part));
}
// GravityField.BandSeries: the band's potential at t for tau > 1, a series in (1 - t) / (1 + tau).
// GravityField.BandSeries: the band's potential at t for tau > 0.4, a series in (1 - t) / (1 + tau).
float gravityBandSeries(float tau, float t)
{
float y = 1 - t;
@ -74,9 +74,9 @@ float gravityTermPotential(float softening, float a, float b, float r)
float b2 = b * b;
if (s >= b2) return 0;
float tau = (a2 + e2) / (b2 - a2);
if (s <= a2) return 0.5 * log((s + e2) / (a2 + e2)) + (tau <= 1 ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0));
if (s <= a2) return 0.5 * log((s + e2) / (a2 + e2)) + (tau <= 0.4 ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0));
float t = (s - a2) / (b2 - a2);
return tau <= 1 ? gravityBandClosed(tau, t) : gravityBandSeries(tau, t);
return tau <= 0.4 ? gravityBandClosed(tau, t) : gravityBandSeries(tau, t);
}
// Zone.GetHeight at world xz: each term's potential minus its own reference, summed, divided by the exaggeration.

248
tests/Aetheria.Shared.Tests/WellLawTests.cs
View file

@ -80,25 +80,35 @@ public sealed class WellLawTests : IDisposable
Assert.InRange(deficit(.1f) / deficit(.05f), 13.0, 15.0); // doubling t multiplies the deficit by about 14 (t^4), not by 4
}
// tau = (a^2 + eps^2) / (b^2 - a^2): the closed form holds to tau 1 and the series beyond it, and the two meet to the closed
// form's single-precision floor (the research's bound of 1.5e-5 per unit strength; the double-precision gap is 7e-7); in single
// precision the tau 10 potential stays within 1e-5 of the double-precision integral of the pull.
// The band's potential is the integral of its pull, in single precision, on both sides of the series threshold: tau =
// (a^2 + eps^2) / (b^2 - a^2) is swept from 0.05 to 10 (every (a, b) where that tau has a real softening), and the
// potential stays within 1.2e-6 per unit strength of the double-precision Simpson integral of the pull. A closed form kept
// past tau 0.45 or a series cut short misses this by 1e-5 to 1e-3.
[Fact]
public void TheSeriesBranchMeetsTheClosedForm()
public void TheBandPotentialIsExactInSinglePrecision()
{
// tau 1 +- 1e-6: eps^2 = b^2 - 2 a^2 gives tau 1 exactly.
const float a = 0.8f, b = 1.28062487f;
var below = (float) Math.Sqrt(b * b - 2 * a * a) * (1 - 1e-6f);
var above = (float) Math.Sqrt(b * b - 2 * a * a) * (1 + 1e-6f);
foreach (var r in new[] { 0.5f, 0.85f, 1.0f, 1.2f })
Assert.InRange(GravityField.Potential(above, a, b, r) - GravityField.Potential(below, a, b, r), -2e-5f, 2e-5f);
const float eps = 59f, inner = 11f, outer = 22f;
foreach (var r in new[] { 12f, 15f, 18f, 21f })
{
var oracle = -Simpson(x => GravityField.Pull(eps, inner, outer, (float) x), r, outer);
Assert.InRange(GravityField.Potential(eps, inner, outer, r), oracle - 1e-5, oracle + 1e-5);
}
var checkedCases = 0;
foreach (var tau in new[] { .05, .3, .39, .41, .6, .7, 1, 1.5, 2.5, 10 })
foreach (var (a, b) in new[] { (1000.0, 1500.0), (40.0, 60.0), (11.0, 22.0), (0.0, 1525.0), (150.0, 4000.0) })
{
var e2 = tau * (b * b - a * a) - a * a;
if (e2 < 0) continue;
var eps = (float) Math.Sqrt(e2);
e2 = (double) eps * eps; // the oracle integrates the softening the code is given
for (var k = 0; k < 15; k++)
{
var r = b * (.02 + .979 * k / 14);
var oracle = -Simpson(x =>
{
var t = x <= a ? 0.0 : (x * x - a * a) / (b * b - a * a);
return x / (x * x + e2) * (1 - S(t));
}, r, b, 8000);
var actual = GravityField.Potential(eps, (float) a, (float) b, (float) r);
Assert.InRange(actual, oracle - 1.2e-6, oracle + 1.2e-6);
checkedCases++;
}
}
Assert.True(checkedCases > 400, "the sweep ran its cases");
}
private static float Plummer(PlanetSettings curves, float mass, float d) =>
@ -156,7 +166,7 @@ public sealed class WellLawTests : IDisposable
Assert.True(length(family.GetForce(inside) - merged.GetForce(inside)) > 1e-3f * length(merged.GetForce(inside)), "inside the support the moon is a body of its own");
}
// The root has no reach: far from everything the zone is one body of its total mass, however the generator split it.
// A sole root has no reach: far from everything the zone is one body of its total mass, however the generator split it.
[Fact]
public void TheZonesFarFieldIsItsTotalMass()
{
@ -175,6 +185,203 @@ public sealed class WellLawTests : IDisposable
}
}
// Three subzones: a sun of 5000 with a planet, a sun of 900 alone, and a bodyless pair; each is a root (a parentless orbit
// stands where its distance and phase put it). Masses are the test's own, the zone's total 6100.
private Zone ThreeSubzones(out CultRecordKey[] roots)
{
var w = _worlds.New();
var big = w.Orbit(default, 0f);
w.Body(big, 5000f);
w.Body(w.Orbit(big, 800f), 100f);
var small = w.Orbit(default, 9000f, .3f);
w.Body(small, 900f);
var pair = w.Orbit(default, 12000f, .7f);
w.Body(w.Orbit(pair, 300f, 0f), 60f);
w.Body(w.Orbit(pair, 300f, .5f), 40f);
roots = new[] { big, small, pair };
return w.Build();
}
private static readonly float[] ThreeSubzoneMass = { 5100f, 900f, 100f };
// The zone node sits at the family-mass-weighted mean of its roots' places and pulls the zone's total mass from there.
[Fact]
public void AZoneOfSeveralSunsPullsAsItsTotalMassFromTheirBarycentre()
{
var zone = ThreeSubzones(out var roots);
var centre = float2.zero;
for (var i = 0; i < roots.Length; i++) centre += zone.Orbits[roots[i]].Position * ThreeSubzoneMass[i] / 6100f;
Assert.InRange(length(zone.Wells.Centre - centre), 0f, 1e-3f);
var far = 10f * zone.Pack.Radius;
foreach (var angle in new[] { 0f, 2.1f, 4.4f })
Near(ZoneWorlds.Curves().WellStrength.Evaluate(6100f),
length(zone.GetForce(zone.Wells.Centre + float2(far * cos(angle), far * sin(angle)))) * far, 1e-3f);
}
// One sun with a planet, alone in its zone at the place a root of `distance` and `phase` stands.
private Zone OneSubzone(float sun, float distance, float phase, out CultRecordKey root)
{
var w = _worlds.New();
root = w.Orbit(default, distance, phase);
w.Body(root, sun);
w.Body(w.Orbit(root, 700f, .2f), 120f);
return w.Build();
}
// Isolated subzones: inside its own support a sun is felt exactly as if it were alone in the zone, the other suns and the
// zone node adding nothing; a root's reach is half the distance to the nearest other sun, or 1.5 of its support if more.
[Fact]
public void ASubzoneFeelsOnlyItsOwnSun()
{
var w = _worlds.New();
var a = w.Orbit(default, 0f);
w.Body(a, 4000f);
w.Body(w.Orbit(a, 700f, .2f), 120f);
var b = w.Orbit(default, 9000f, 0f);
w.Body(b, 2500f);
w.Body(w.Orbit(b, 700f, .2f), 120f);
var zone = w.Build();
var rng = new System.Random(3);
foreach (var (root, sun, distance) in new[] { (a, 4000f, 0f), (b, 2500f, 9000f) })
{
var alone = OneSubzone(sun, distance, 0f, out var aloneRoot);
Assert.True(float.IsPositiveInfinity(alone.Wells.Reach(aloneRoot)), "a sole massive root is the whole zone");
var support = zone.Wells.Support(root);
Assert.Equal(Math.Max(.5f * 9000f, 1.5f * support), zone.Wells.Reach(root), 2);
var at = zone.Orbits[root].Position;
for (var i = 0; i < 50; i++)
{
var angle = (float) (rng.NextDouble() * 2 * Math.PI);
var radius = support * (.2f + .8f * (float) rng.NextDouble());
var p = at + float2(radius * cos(angle), radius * sin(angle));
var expected = alone.GetForce(p);
Assert.InRange(length(zone.GetForce(p) - expected), 0f, 1e-5f * length(expected) + 1e-9f);
}
}
}
// Between the subzones the ground is flat; beyond the populated area the zone is one body of its total mass at its centre.
[Fact]
public void BeyondEverySubzoneTheZoneIsOneBody()
{
var zone = ThreeSubzones(out var roots);
var curves = ZoneWorlds.Curves();
var centre = zone.Wells.Centre;
var span = roots.Max(r => length(zone.Orbits[r].Position - centre) + zone.Wells.Reach(r)); // A, as the law defines it
var flat = 0;
for (var x = -14000f; x <= 14000f; x += 700f)
for (var y = -14000f; y <= 14000f; y += 700f)
{
var p = float2(x, y);
if (length(p - centre) >= span || roots.Any(r => length(p - zone.Orbits[r].Position) <= zone.Wells.Reach(r))) continue;
flat++;
Assert.InRange(length(zone.GetForce(p)), 0f, 1e-7f);
}
Assert.True(flat > 100, "fixture: there is open ground between the subzones");
foreach (var (radius, angle) in new[] { (1.6f, .4f), (3f, 2.5f), (8f, 5.1f) })
{
var p = centre + float2(span * radius * cos(angle), span * radius * sin(angle));
var expected = (centre - p) / length(centre - p) * Plummer(curves, 6100f, length(centre - p));
Assert.InRange(length(zone.GetForce(p) - expected), 0f, 1e-5f * length(expected));
}
}
// A belt is a child of its parent's family: the family's support covers the belt's outermost asteroid, so every asteroid
// circles the planet's body alone at its own distance (not at the belt orbit's, 43, where the fade used to start).
[Fact]
public void ABeltRidesInsideItsParentsBody()
{
var w = _worlds.New();
var sun = w.Orbit(default, 0f);
w.Body(sun, 5000f);
var planet = w.Orbit(sun, 3000f);
w.Body(planet, 172f);
var belt = w.Belt(w.Orbit(planet, 43f, .4f), 10f, 43f, 87f);
var zone = w.Build();
Assert.Equal(1.5f * 87f, zone.Wells.Support(planet), 3);
var curves = ZoneWorlds.Curves();
var periods = zone.AsteroidBelts[belt].Periods;
Assert.Equal(3, periods.Length);
for (var i = 0; i < 3; i++)
{
var d = zone.AsteroidBelts[belt].Data.Asteroids[i].Distance;
var expected = 2 * PI * sqrt(d / Plummer(curves, 172f, d));
Assert.InRange(periods[i], expected * (1 - 1e-5f), expected * (1 + 1e-5f));
}
}
// The fade of a tight family begins at its support: the planet (moon at 150, support 225, reach at its floor 337.5) pulls
// as the whole family out to 225 and then gives way as the order-4 smoothstep of (r^2 - 225^2) / (337.5^2 - 225^2).
[Fact]
public void ATightFamilysPullFadesFromItsSupport()
{
var w = _worlds.New();
var sun = w.Orbit(default, 0f);
w.Body(sun, 5000f);
var planet = w.Orbit(sun, 400f);
w.Body(planet, 100f);
var moon = w.Orbit(planet, 150f);
w.Body(moon, 1f);
var zone = w.Build();
var curves = ZoneWorlds.Curves();
Assert.Equal(225f, zone.Wells.Support(planet), 3);
Assert.Equal(337.5f, zone.Wells.Reach(planet), 3);
foreach (var r in new[] { 232f, 250f, 270f, 290f, 310f, 325f, 333f, 336f })
{
var t = (r * r - 225.0 * 225) / (337.5 * 337.5 - 225.0 * 225);
var expected = Plummer(curves, 101f, r) * (1 - S(t));
Assert.InRange(zone.Wells.PullAround(planet, r, moon), expected * (1 - 1e-5) - 1e-9, expected * (1 + 1e-5) + 1e-9);
}
}
// A ring group's pull rises from its ring: a pair at 300 round a bodyless node pulls from the node as the order-4 smoothstep
// of (r^2 - 300^2) / (1500^2 - 300^2) of the pair's Plummer, 1500 being the ring's distance and the members' RingReach.
[Fact]
public void ARingGroupsPullRisesFromItsRing()
{
var w = _worlds.New();
var sun = w.Orbit(default, 0f);
w.Body(sun, 5000f);
var node = w.Orbit(sun, 3000f);
w.Body(w.Orbit(node, 300f, 0f), 80f);
w.Body(w.Orbit(node, 300f, .5f), 60f);
var zone = w.Build();
var curves = ZoneWorlds.Curves();
foreach (var r in new[] { 310f, 400f, 600f, 900f, 1200f, 1400f, 1490f })
{
var t = (r * r - 300.0 * 300) / (1500.0 * 1500 - 300.0 * 300);
var expected = Plummer(curves, 140f, r) * S(t);
Assert.InRange(zone.Wells.PullAround(node, r, default), expected * (1 - 1e-5) - 1e-9, expected * (1 + 1e-5) + 1e-9);
}
}
// The drawn height is the force's potential everywhere, however far: its slope along the radius is the pull, for a zone of
// one sun and for one of three subzones, out to ten pack radii (a potential flattened at some radius reads zero here).
[Fact]
public void FarOutTheHeightIsTheForcesPotential()
{
var single = SunPlanetMoon(1f, out _, out _);
var several = ThreeSubzones(out _);
foreach (var (zone, centre) in new[] { (single, float2.zero), (several, several.Wells.Centre) })
foreach (var radius in new[] { 2f, 4f, 7f, 10f })
foreach (var angle in new[] { .3f, 2.4f, 4.6f })
{
var outward = float2(cos(angle), sin(angle));
var p = centre + outward * (radius * zone.Pack.Radius);
const float h = 10f;
var slope = (zone.GetHeight(p + outward * h) - zone.GetHeight(p - outward * h)) * zone.Settings.HeightExaggeration / (2 * h);
var pull = -dot(zone.GetForce(p), outward);
Assert.InRange(slope, pull * .98f, pull * 1.02f);
}
}
private static void Near(float expected, float actual, float relative) =>
Assert.InRange(actual, expected * (1 - relative), expected * (1 + relative));
@ -233,18 +440,19 @@ public sealed class WellLawTests : IDisposable
var node = w.Orbit(sunOrbit, 5200f, .13f);
w.Body(w.Orbit(node, 300f, 0f), 80f);
w.Body(w.Orbit(node, 300f, .5f), 60f);
w.Body(w.Orbit(default, 9000f, .3f), 900f); // a second massive root: the zone node joins the walk
var zone = w.Build();
var rng = new System.Random(7);
for (var n = 0; n < 200; n++)
{
var p = float2((float) (rng.NextDouble() * 12000 - 6000), (float) (rng.NextDouble() * 12000 - 6000));
var p = float2((float) (rng.NextDouble() * 24000 - 6000), (float) (rng.NextDouble() * 24000 - 12000));
var force = float2.zero;
var potential = 0f;
var scale = 0f;
foreach (var term in zone.Wells.Terms)
{
var offset = zone.Orbits[term.Node].Position - p;
var offset = (term.Node.IsSet() ? zone.Orbits[term.Node].Position : zone.Wells.Centre) - p; // an unset node is the zone node
var r = length(offset);
var pull = term.Strength * GravityField.Pull(term.Softening, term.Inner, term.Outer, r);
if (r > 0f) force += offset / r * pull;

9
tests/Aetheria.Shared.Tests/ZoneWorlds.cs
View file

@ -61,6 +61,15 @@ internal sealed class ZoneWorlds : IDisposable
return key;
}
// A belt of asteroids at the given distances riding `orbit`; returns the belt's body key.
public CultRecordKey Belt(CultRecordKey orbit, params float[] distances)
{
var asteroids = distances.Select((d, i) => new Asteroid { Distance = d, Phase = frac(i * .6180339887f), Size = .5f, RotationSpeed = .1f }).ToArray();
var key = Cache.Upsert(new AsteroidBeltData { Orbit = new CultRecordRef<OrbitData>(orbit), Asteroids = asteroids }).Key;
Pack.Planets.Add(new CultRecordRef<BodyData>(key));
return key;
}
public Zone Build() => new Zone(Items, Curves(), Pack, new GalaxyZone { Name = "OrbitsFromField", Owner = null }, null);
}