diff --git a/Assets/Scripts/Editor/GravityFieldParityVerify.cs b/Assets/Scripts/Editor/GravityFieldParityVerify.cs index fe1aef918..766a58a6f 100644 --- a/Assets/Scripts/Editor/GravityFieldParityVerify.cs +++ b/Assets/Scripts/Editor/GravityFieldParityVerify.cs @@ -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("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 { new Vector2(-130f, 85f), new Vector2(610f, -240f), new Vector2(670f, -270f), new Vector2(1700f, 300f) }; + var masses = new List { 14152f, 273.9f, 29.2f, 0f }; + var parents = new List { -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) { diff --git a/Assets/Scripts/ServerShared/GravityField.cs b/Assets/Scripts/ServerShared/GravityField.cs index 9b7854b8e..0eb3797de 100644 --- a/Assets/Scripts/ServerShared/GravityField.cs +++ b/Assets/Scripts/ServerShared/GravityField.cs @@ -28,7 +28,10 @@ public static class GravityField /// A ring member of a bodyless node reaches RingReach x its distance: whole out to the ring's diameter. public const float RingReach = 4f; - /// The series branch's length: its ratio is below 1/2, so 40 terms reach float32's floor. + /// 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. + public const float SeriesFrom = .4f; + + /// The series branch's length: its ratio is below 1/1.4, so 40 terms reach float32's floor. public const int SeriesTerms = 40; /// S(t) as polynomial coefficients in t, lowest first: t^Order sum_j C(Order-1+j, j) (1-t)^j. @@ -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) { diff --git a/Assets/Scripts/ServerShared/GravityShader.cs b/Assets/Scripts/ServerShared/GravityShader.cs index 6172fa161..e4736ec95 100644 --- a/Assets/Scripts/ServerShared/GravityShader.cs +++ b/Assets/Scripts/ServerShared/GravityShader.cs @@ -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. diff --git a/Assets/Scripts/ServerShared/WellTree.cs b/Assets/Scripts/ServerShared/WellTree.cs index 8f9cd50da..3448c8a9e 100644 --- a/Assets/Scripts/ServerShared/WellTree.cs +++ b/Assets/Scripts/ServerShared/WellTree.cs @@ -5,19 +5,26 @@ using CultMath; using static CultMath.math; using float2 = CultMath.float2; -/// 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). +/// +/// 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. +/// 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. /// @@ -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; - /// Every term of the zone, grouped by node: what the lowering uploads and the shader sums. + /// Where the zone node's terms sit (the family-mass-weighted mean of the massive roots' anchors); meaningful when it has terms. + public float2 Centre { get; private set; } + + /// 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. public IReadOnlyList Terms => _terms; private WellTree() { } - /// A node's reach in its parent's field (+inf for a root); 0 for an orbit that carries no mass. + /// 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. public float Reach(CultRecordKey node) => _reach[_index[node]]; /// A node's support B: every term of its descendants lies inside it, so beyond it the family is one body. @@ -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(); + 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 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 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; } diff --git a/Assets/Scripts/ServerShared/Zone.cs b/Assets/Scripts/ServerShared/Zone.cs index 9ec7b9eb4..02904e461 100644 --- a/Assets/Scripts/ServerShared/Zone.cs +++ b/Assets/Scripts/ServerShared/Zone.cs @@ -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(); + 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(); 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() diff --git a/Assets/Scripts/Zone Display/GravityLowering.cs b/Assets/Scripts/Zone Display/GravityLowering.cs index e9da4178a..2c0a0c2b3 100644 --- a/Assets/Scripts/Zone Display/GravityLowering.cs +++ b/Assets/Scripts/Zone Display/GravityLowering.cs @@ -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), diff --git a/Assets/Shaders/Gravity.cginc b/Assets/Shaders/Gravity.cginc index b6b4d3300..5eacdec0f 100644 --- a/Assets/Shaders/Gravity.cginc +++ b/Assets/Shaders/Gravity.cginc @@ -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. diff --git a/tests/Aetheria.Shared.Tests/WellLawTests.cs b/tests/Aetheria.Shared.Tests/WellLawTests.cs index 746712a60..096a57a78 100644 --- a/tests/Aetheria.Shared.Tests/WellLawTests.cs +++ b/tests/Aetheria.Shared.Tests/WellLawTests.cs @@ -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; diff --git a/tests/Aetheria.Shared.Tests/ZoneWorlds.cs b/tests/Aetheria.Shared.Tests/ZoneWorlds.cs index 98a43a6b5..aaa1a840a 100644 --- a/tests/Aetheria.Shared.Tests/ZoneWorlds.cs +++ b/tests/Aetheria.Shared.Tests/ZoneWorlds.cs @@ -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(orbit), Asteroids = asteroids }).Key; + Pack.Planets.Add(new CultRecordRef(key)); + return key; + } + public Zone Build() => new Zone(Items, Curves(), Pack, new GalaxyZone { Name = "OrbitsFromField", Owner = null }, null); }