diff --git a/Assets/Resources/Settings.asset b/Assets/Resources/Settings.asset index 4b17fe993..4536e99f1 100644 --- a/Assets/Resources/Settings.asset +++ b/Assets/Resources/Settings.asset @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:248373994a9d0caf17b8dcddbedb7b8c1d36001c31ae0f1e695385e6589734e6 +oid sha256:2c67fcecfea500365c54e5c80e59b6f26d0526e877c266d9fa75a43053df225d size 12647 diff --git a/Assets/Scripts/Editor/GravityFieldParityVerify.cs b/Assets/Scripts/Editor/GravityFieldParityVerify.cs index 1753b95ab..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 and windows is 0.0065 u, in the last frame's planetoid window (float32 on both sides, the GPU's log against the CPU's, the - // quad's interpolated world position). The tolerance is 4.6x that. A 0.1% error in the tail constant (.5 -> .501) - // moves the drawn height by 0.053 u where all three tails reach (the planetoid and zone windows), 1.8x the - // tolerance; a reference of zero moves it by about 150 u. + // 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.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,22 +73,29 @@ 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); - // Three frames through the one buffer, as ZoneRenderer.Update runs it: the bodies as built, then two moves - // (an orbit's Position is what Zone.GetHeight reads). A lowering that mirrors an axis or stops refreshing - // the upload agrees with the sim in the first frame only where it is symmetric, and not after a move. + // Three frames through the one buffer, as ZoneRenderer.Update runs it: the bodies as built, then two moves (an + // orbit's Position is what Zone.GetHeight reads). A move keeps the tree's geometry, since a rail never leaves + // its parent's distance and the walk skips a subtree outside its reach: the sun is the root, the planet is + // its Distance from the sun, the planetoid its Distance from the planet. A lowering that mirrors an axis or + // stops refreshing the upload agrees with the sim in the first frame only where it is symmetric. CheckFrame("built", zone, orbits); - orbits[1].Position = new float2(-420f, 560f); - orbits[2].Position = new float2(900f, -700f); + 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[2].Position = new float2(-1800f, 1300f); + 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); @@ -108,35 +118,69 @@ namespace Aetheria.EditorTools static float Core(Zone zone, float mass) => zone.Settings.RocheLimit.Evaluate(mass); + // Where an orbit lies when it stands `turns` round its parent at its Distance. + static float2 Around(Orbit parent, Orbit child, float turns) => parent.Position + child.Data.Distance * OrbitData.Evaluate(turns); + 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()); } - // A sun (14152 at (-130, 85)), a planet (273.9 at (610, -240)), a planetoid (29.2 at (670, -270)) and a body of - // mass 0, whose well has no strength (a legal body; its curves give c = 0 and R = 0), at (1700, 300): parentless, - // fixed, 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. - internal static Zone BuildZone(GameSettings settings, string work, out Orbit[] orbits, float radius = 4182f) + // A sun (14152 at (-130, 85)), its planet (273.9, placed at (610, -240) at the pack's time), the planet's planetoid + // (29.2, placed at (670, -270)) and a body of mass 0, whose well has no strength (a legal body; its curves give + // 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, 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 keys = new CultRecordKey[positions.Length]; - var pack = new ZonePack { Radius = radius }; - for (var i = 0; i < positions.Length; i++) + 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) { - var orbit = keys[i] = cache.Upsert(new OrbitData { FixedPosition = new float2(positions[i].x, positions[i].y) }).Key; + // 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.Count; i++) + { + var data = new OrbitData { FixedPosition = new float2(positions[i].x, positions[i].y) }; + if (parents[i] >= 0) + { + var offset = positions[i] - positions[parents[i]]; + data.Parent = new CultRecordRef(keys[parents[i]]); + data.Distance = offset.magnitude; + data.Phase = Mathf.Atan2(offset.y, offset.x) / (2f * Mathf.PI); + } + var orbit = keys[i] = cache.Upsert(data).Key; var body = cache.Upsert(new PlanetData { Name = "Body " + i, Orbit = new CultRecordRef(orbit), Mass = masses[i] }).Key; pack.Orbits.Add(new CultRecordRef(orbit)); @@ -195,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); @@ -215,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"); } @@ -239,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/DrawAhead.cs b/Assets/Scripts/ServerShared/DrawAhead.cs index 49ca0a3d0..c3cde0ef3 100644 --- a/Assets/Scripts/ServerShared/DrawAhead.cs +++ b/Assets/Scripts/ServerShared/DrawAhead.cs @@ -51,24 +51,6 @@ public static class DrawAhead return float3(at.x + bearing.x * distance, at.y, at.z + bearing.y * distance); } - // One well as the shader takes it: its body's rail at the drawn time, with the strength, core and reference the sim's - // GetHeight sums with. GravityLowering copies it into the shader's record and decides nothing. - public readonly struct DrawnWell - { - public readonly float2 Position; - public readonly float Strength, CoreRadius, Reference; - public DrawnWell(float2 position, float strength, float coreRadius, float reference) - { - Position = position; - Strength = strength; - CoreRadius = coreRadius; - Reference = reference; - } - } - - public static DrawnWell Well(Zone zone, Planet body, float lead) => - new DrawnWell(Orbit(zone, body.BodyData.Orbit.Key, lead), body.WellStrength, body.RocheLimit, zone.WellReference(body)); - // A free body (a vapour cloud) drawn where the frame falls lead sim seconds after the last step, by the velocity // that step moved it with: the launch velocity and the drift together. public static float2 Position(KinematicBody body, float lead) => body.Position + body.TotalVelocity * lead; diff --git a/Assets/Scripts/ServerShared/GravityField.cs b/Assets/Scripts/ServerShared/GravityField.cs index 01a070e9d..084aa1f45 100644 --- a/Assets/Scripts/ServerShared/GravityField.cs +++ b/Assets/Scripts/ServerShared/GravityField.cs @@ -1,32 +1,150 @@ +using System; using CultMath; using static CultMath.math; -using float2 = CultMath.float2; /// -/// The law of one gravity well: a 2D point mass with a softened core. Outside the core radius the potential is -/// logarithmic and the pull is c / r; inside it the pull falls linearly to zero at the centre and meets the tail in -/// value and slope at the core radius. c is the well's strength (u^2/s^2), R its core radius (u). No threshold -/// anywhere: a well pulls at every distance. Zone.GetForce and Zone.GetHeight sum this over a zone's bodies. +/// The law of one band term U(r; eps, a, b), the only primitive of a zone's gravity (WellTree decides which terms +/// exist). A Plummer-softened log source whose pull toward the centre is exact up to a and fades to nothing over [a, b]: +/// pull r / (r^2 + eps^2) * B(s), s = r^2, B = 1 for s <= a^2, 1 - S(t) on the band with t = (s - a^2) / (b^2 - a^2), +/// 0 for s >= b^2; S is the order-Order smoothstep. The potential is 0 beyond b. With a = b = +inf it is the global term +/// 1/2 ln(r^2 + eps^2). The potential is a function of r^2 (smooth at the centre), and C^Order on the circles r = a and +/// r = b (C4 at Order 4: the player feels higher-order seams), so the drawn height is exactly the sim's potential. +/// On the band the closed form is 1/2 [R ln((t + tau) / (1 + tau)) - (I(1) - I(t))], tau = (a^2 + eps^2) / (b^2 - a^2), +/// (1 - S(t)) = (t + tau) Q(t) + R, I the antiderivative of Q; for tau > 1 it cancels catastrophically in float32 +/// (7e-3 at tau 3, 12 at tau 10), so the band is summed as a series in (1 - t) / (1 + tau), whose ratio is below 1/2. +/// Pull and Potential are per unit strength; the caller scales by the term's strength (WellStrength at a mass). +/// Gravity.cginc is generated from this file's constants and coefficients (GravityShader). /// public static class GravityField { - // The tail's additive constant: it makes the tail meet the core in value at the core radius. Gravity.cginc is generated from it. - public const float TailOffset = .5f; + /// The smoothstep's order: the potential is C^Order on every band edge. + public const int Order = 4; + /// A family's charge is fully on by its support, Beta x its outermost child. + public const float Beta = 1.5f; + /// A family is felt whole out to max(support, Gamma x its reach) and fades from there to its reach. + public const float Gamma = .5f; + /// A body's reach in its parent's field, in its distance (the nominal reach before the fade's room). + public const float Reach = .5f; + /// 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 potential at distance r from the body, defined up to a constant (callers subtract their reference). - // A well of no strength (a body whose curves are unauthored, so c = 0 and R = 0) digs nothing. - public static float Potential(float c, float coreRadius, float r) => - c == 0f ? 0f : - r < coreRadius - ? c * r * r / (2 * coreRadius * coreRadius) - : c * (log(r / coreRadius) + TailOffset); + /// The band is summed as a series for tau above this. float32 error per unit strength against the double integral is at most 1.8e-6 anywhere (worst 1.7e-6 just above the seam, where the series has not reached float32 resolution); the closed form is float32-clean only to about 0.45. + public const float SeriesFrom = .4f; - // The acceleration at `offsetFromBody`, toward the body: zero at the body itself. - public static float2 Acceleration(float c, float coreRadius, float2 offsetFromBody) + /// The series branch's length: its ratio is below 1/1.4, so the truncated tail is below 1.4e-6 at the seam and falls geometrically above it. + 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. + public static readonly float[] S = BuildS(); + + /// 1 - S(t) as polynomial coefficients in t, lowest first: the band's profile B. + public static readonly float[] OneMinusS = BuildOneMinusS(); + + private static float[] BuildS() { - var r = length(offsetFromBody); - if (r == 0f) return float2.zero; - var pull = r < coreRadius ? c * r / (coreRadius * coreRadius) : c / r; - return -offsetFromBody / r * pull; + var c = new double[2 * Order]; + for (var j = 0; j < Order; j++) + { + var cj = Binomial(Order - 1 + j, j); + for (var i = 0; i <= j; i++) + c[Order + i] += cj * Binomial(j, i) * (i % 2 == 0 ? 1 : -1); + } + var s = new float[2 * Order]; + for (var i = 0; i < s.Length; i++) s[i] = (float) c[i]; + return s; + } + + private static float[] BuildOneMinusS() + { + var f = new float[2 * Order]; + for (var i = 0; i < f.Length; i++) f[i] = -S[i]; + f[0] += 1f; + return f; + } + + private static double Binomial(int n, int k) + { + var r = 1.0; + for (var i = 1; i <= k; i++) r = r * (n - k + i) / i; + return r; + } + + // B(s) for the band term (a, b); b = +inf is no fade at all. + private static float Profile(float inner, float outer, float s) + { + if (float.IsPositiveInfinity(outer)) return 1f; + var a2 = inner * inner; + var b2 = outer * outer; + if (s <= a2) return 1f; + if (s >= b2) return 0f; + var t = (s - a2) / (b2 - a2); + var acc = 0f; + for (var i = OneMinusS.Length - 1; i >= 0; i--) acc = acc * t + OneMinusS[i]; + return acc; + } + + /// The pull toward the centre at distance r, per unit strength. + public static float Pull(float softening, float inner, float outer, float r) + { + var s = r * r; + return r * Profile(inner, outer, s) / (s + softening * softening); + } + + /// The potential at distance r, per unit strength: 0 at and beyond outer, the global log when outer is +inf. + public static float Potential(float softening, float inner, float outer, float r) + { + var s = r * r; + var e2 = softening * softening; + if (float.IsPositiveInfinity(outer)) return .5f * log(s + e2); + var a2 = inner * inner; + var b2 = outer * outer; + if (s >= b2) return 0f; + var tau = (a2 + e2) / (b2 - a2); + if (s <= a2) return .5f * log((s + e2) / (a2 + e2)) + Band(tau, 0f); + return Band(tau, (s - a2) / (b2 - a2)); + } + + // The band's potential at t in [0, 1], 0 at t = 1. + private static float Band(float tau, float t) => tau <= SeriesFrom ? BandClosed(tau, t) : BandSeries(tau, t); + + private static float BandClosed(float tau, float t) + { + // (1 - S(t)) = (t + tau) Q(t) + R by synthetic division; I is Q's antiderivative. + var n = OneMinusS.Length - 1; + Span q = stackalloc float[n]; + q[n - 1] = OneMinusS[n]; + for (var i = n - 1; i > 0; i--) q[i - 1] = OneMinusS[i] - tau * q[i]; + var remainder = OneMinusS[0] - tau * q[0]; + return .5f * (remainder * log((t + tau) / (1f + tau)) - (Antiderivative(q, 1f) - Antiderivative(q, t))); + } + + private static float Antiderivative(ReadOnlySpan q, float t) + { + var acc = 0f; + for (var i = q.Length - 1; i >= 0; i--) acc = (acc + q[i] / (i + 1)) * t; + return acc; + } + + private static float BandSeries(float tau, float t) + { + // With y = 1 - t, 1 - S(1 - y) = S(y) and 1 / (t + tau) = sum_k y^k / (1 + tau)^(k+1): + // band = -1/2 sum_k (1 + tau)^-(k+1) sum_j S_j y^(j+k+1) / (j+k+1). + var y = 1f - t; + Span yj = stackalloc float[S.Length]; + yj[0] = 1f; + for (var j = 1; j < yj.Length; j++) yj[j] = yj[j - 1] * y; + var w = 1f / (1f + tau); + var wk = w; + var yk = y; + var acc = 0f; + for (var k = 0; k < SeriesTerms; k++) + { + var inner = 0f; + for (var j = Order; j < S.Length; j++) inner += S[j] * yj[j] / (j + k + 1); + acc += wk * yk * inner; + wk *= w; + yk *= y; + } + return -.5f * acc; } } diff --git a/Assets/Scripts/ServerShared/GravityShader.cs b/Assets/Scripts/ServerShared/GravityShader.cs index ad163e26d..e4736ec95 100644 --- a/Assets/Scripts/ServerShared/GravityShader.cs +++ b/Assets/Scripts/ServerShared/GravityShader.cs @@ -1,10 +1,12 @@ using System.Globalization; +using System.Linq; /// /// The single author of Assets/Shaders/Gravity.cginc, the GPU lowering of the zone's potential. The law is -/// GravityField's (its constants are read here, not respelled) and the height is Zone.GetHeight: each well's -/// potential minus its reference, summed, divided by the exaggeration. Unity does not run C# when it imports a -/// shader, so the committed cginc is this output; GravityShaderTests fails unless it equals Source() exactly. +/// GravityField's (its constants and smoothstep coefficients are read here, not respelled) and the height is +/// Zone.GetHeight: each term's potential minus its reference, summed over the terms whose edge the point is inside, +/// divided by the exaggeration. Unity does not run C# when it imports a shader, so the committed cginc is this output; +/// GravityShaderTests fails unless it equals Source() exactly. /// To change the shader, change this file and regenerate: AETHERIA_REGENERATE_GRAVITY_CGINC=1 dotnet test --filter GravityShaderTests. /// public static class GravityShader @@ -12,48 +14,111 @@ public static class GravityShader public const string Path = "Assets/Shaders/Gravity.cginc"; /// The whole file, with LF line ends. - public static string Source() => (Template - .Replace("@TAIL@", GravityField.TailOffset.ToString("R", CultureInfo.InvariantCulture))) + public static string Source() => Template + .Replace("@ORDER@", GravityField.Order.ToString(CultureInfo.InvariantCulture)) + .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)) .Replace("\r\n", "\n"); + private static string Coefficients(float[] c) => + string.Join(", ", c.Select(x => x.ToString("R", CultureInfo.InvariantCulture))); + private const string Template = @"#ifndef AETHERIA_GRAVITY_INCLUDED #define AETHERIA_GRAVITY_INCLUDED // GENERATED by ServerShared/GravityShader.cs from GravityField; do not edit (GravityShaderTests fails on any difference). // The lowering of the zone's potential (aetheria-gravity, cut grid-lowering). The law is ServerShared's GravityField // and the height is Zone.GetHeight; this file decides nothing. GravityLowering.Lower fills the buffer and the globals -// from the zone (positions from Orbit.Position, reference from Zone.WellReference) each frame. +// from the zone's WellTree (positions from the nodes' rails at the drawn time) each frame. // GravityFieldParityVerify renders it and compares every texel with -Zone.GetHeight. -struct GravityWell +struct GravityTerm { float2 position; float strength; - float coreRadius; + float softening; + float inner; + float outer; float reference; }; -StructuredBuffer _GravityWells; +StructuredBuffer _GravityWells; int _GravityWellCount; float _HeightExaggeration; -// GravityField.Potential: c is the well's strength, R its core radius, r the distance. -float gravityPotential(float c, float R, float r) +// GravityField.OneMinusS and GravityField.S: the band profile B(t) = 1 - S(t) and the order-@ORDER@ smoothstep, lowest first. +static const float GravityOneMinusS[@N@] = { @ONEMINUSS@ }; +static const float GravityS[@N@] = { @S@ }; + +// GravityField.BandClosed: the band's potential at t for tau <= @SERIESFROM@, per unit strength. +float gravityBandClosed(float tau, float t) { - return c == 0 ? 0 : - r < R - ? c * r * r / (2 * R * R) - : c * (log(r / R) + @TAIL@); + float q[@Q@]; + q[@QLAST@] = GravityOneMinusS[@Q@]; + [unroll] for (int i = @QLAST@; i > 0; i--) q[i - 1] = GravityOneMinusS[i] - tau * q[i]; + float remainder = GravityOneMinusS[0] - tau * q[0]; + float whole = 0; + float part = 0; + [unroll] for (int j = @QLAST@; j >= 0; j--) + { + whole = whole + q[j] / (j + 1); + part = (part + q[j] / (j + 1)) * t; + } + return 0.5 * (remainder * log((t + tau) / (1 + tau)) - (whole - part)); } -// Zone.GetHeight at world xz: each well's potential minus its own reference, summed, divided by the exaggeration. +// 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; + float yj[@N@]; + yj[0] = 1; + [unroll] for (int j = 1; j < @N@; j++) yj[j] = yj[j - 1] * y; + float w = 1 / (1 + tau); + float wk = w; + float yk = y; + float acc = 0; + for (int k = 0; k < @SERIES@; k++) + { + float inner = 0; + [unroll] for (int m = @ORDER@; m < @N@; m++) inner += GravityS[m] * yj[m] / (m + k + 1); + acc += wk * yk * inner; + wk *= w; + yk *= y; + } + return -0.5 * acc; +} + +// GravityField.Potential: one band term, per unit strength; the global term when b is infinite. +float gravityTermPotential(float softening, float a, float b, float r) +{ + float s = r * r; + float e2 = softening * softening; + if (isinf(b)) return 0.5 * log(s + e2); + float a2 = a * a; + 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 <= @SERIESFROM@ ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0)); + float t = (s - a2) / (b2 - a2); + 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. float gravityHeight(float2 p) { float result = 0; for (int i = 0; i < _GravityWellCount; i++) { - GravityWell well = _GravityWells[i]; - result += gravityPotential(well.strength, well.coreRadius, length(p - well.position)) - well.reference; + GravityTerm term = _GravityWells[i]; + float r = length(p - term.position); + if (r < term.outer) + result += term.strength * gravityTermPotential(term.softening, term.inner, term.outer, r) - term.reference; } return result / _HeightExaggeration; } diff --git a/Assets/Scripts/ServerShared/Scenarios/LockOn.cs b/Assets/Scripts/ServerShared/Scenarios/LockOn.cs index 87e5de829..7dc850497 100644 --- a/Assets/Scripts/ServerShared/Scenarios/LockOn.cs +++ b/Assets/Scripts/ServerShared/Scenarios/LockOn.cs @@ -12,9 +12,9 @@ public sealed class LockOn : Scenario // The hostile starts about 4000 m off, ahead and to port of the player, facing the player. At that range its passive // sensors hold no track of the player (the test asks the sim's own PassiveInfo against the detection threshold). The - // hostile closes, sees the player some fourteen seconds in and locks. The zone's pull bends the hostile's run and moves the - // player, and the warning must still point at where the hostile was staged to within 5 degrees at the lock: this side keeps - // the bend under one degree (on the far side it is more than seven). At Duel's 900 m (and at 600 m) it sees and locks + // hostile closes, sees the player some fourteen seconds in and locks. Its flanking approach and the zone's pull carry it + // off the staged line (about 5 degrees by the lock on the merged tree), and the warning must point at where the hostile + // is at the lock, which the test checks against the live bearing. At Duel's 900 m (and at 600 m) it sees and locks // within a second. public static readonly CultMath.float2 HostileAt = float2(-2000, 3450); diff --git a/Assets/Scripts/ServerShared/Settings.cs b/Assets/Scripts/ServerShared/Settings.cs index af22d2840..235092f90 100644 --- a/Assets/Scripts/ServerShared/Settings.cs +++ b/Assets/Scripts/ServerShared/Settings.cs @@ -16,11 +16,12 @@ public class PlanetSettings public ExponentialCurve LightRadius; public ExponentialCurve BodyRadius; - // The closest a child may orbit its parent, in u, by the parent's mass, and the softened core radius of the parent's - // well (the steepest drawn slope sits at the Roche limit). Authored like BodyRadius, never per body. + // The closest a child may orbit its parent, in u, by the parent's mass, and the Plummer softening of the parent's + // well (the steepest slope, 0.5 c / eps, sits at the Roche limit). Authored like BodyRadius, never per body. public ExponentialCurve RocheLimit = new ExponentialCurve(); - // A body's well strength by its mass (GravityField): c in u^2/s^2 (the pull is c / r beyond the core, which is the Roche limit). + // A body's well strength by its mass (GravityField): c in u^2/s^2. The pull is c d / (d^2 + eps^2) inside a family and + // c(family mass) beyond it (WellTree); eps is the Roche limit. public ExponentialCurve WellStrength = new ExponentialCurve(); // The one global exaggeration of gravity height: height = potential / K, in u/s^2 per unit of height. diff --git a/Assets/Scripts/ServerShared/WellTree.cs b/Assets/Scripts/ServerShared/WellTree.cs new file mode 100644 index 000000000..bacf474de --- /dev/null +++ b/Assets/Scripts/ServerShared/WellTree.cs @@ -0,0 +1,431 @@ +using System; +using System.Collections.Generic; +using GameCult.Caching; +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). +/// 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, Extent; + public readonly float2 Anchor; + + 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; + } +} + +/// +/// The one owner of which gravity terms a zone has: the orbit tree turned into band terms (GravityField is the law of one +/// term). Built once from the pack's orbits and the authored curves, never reading time; the caller passes positions. +/// +/// 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 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 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. +/// +public sealed class WellTree +{ + /// Orbits whose distances differ by no more than this are on one ring (a binary, a ring of siblings). + public const float RingWindow = .1f; + + /// A band term of the field: Strength x GravityField's U(Softening, Inner, Outer) at Node's position, less Reference. + public readonly struct Term + { + public readonly CultRecordKey Node; + public readonly float Strength, Softening, Inner, Outer, Reference; + + public Term(CultRecordKey node, float strength, float softening, float inner, float outer, float reference) + { + Node = node; + Strength = strength; + Softening = softening; + Inner = inner; + Outer = outer; + Reference = reference; + } + } + + private readonly Dictionary _index = new Dictionary(); + private CultRecordKey[] _key; + private int[][] _kids; + private int[] _roots; + private int[][] _ring; + private int[] _termStart, _termCount; + private float[] _distance, _phase, _reach, _support, _extent; + private float2[] _anchor; + private int[] _massiveRoots; + private float _packRadius; + private int _zoneStart, _zoneCount; + private Term[] _terms; + + /// 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 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. + public float Support(CultRecordKey node) => _support[_index[node]]; + + /// + /// A well's strength and softening at a mass, from the authored curves (the one place they are read and checked). + /// A mass of 0 is no well. + /// + public static (float Strength, float Softening) Curves(PlanetSettings settings, float mass) + { + var strength = settings.WellStrength.Evaluate(mass); + var softening = settings.RocheLimit.Evaluate(mass); + if (!(strength >= 0f) || float.IsInfinity(strength)) + throw new InvalidOperationException("A body's well strength must be finite and not negative (PlanetSettings.WellStrength)."); + if (strength > 0 && !(softening > 0f && !float.IsInfinity(softening))) + throw new InvalidOperationException("A body's Roche limit, which is its softening, must be positive and finite (PlanetSettings.RocheLimit)."); + return (strength, softening); + } + + /// + /// The tree of a zone's orbits. `packRadius` is where the global terms are zero: each carries its own potential there + /// as its Reference (a bounded term is already 0 beyond its outer edge, so it needs none). + /// + public static WellTree Build(IReadOnlyList nodes, PlanetSettings settings, float packRadius) + { + var tree = new WellTree(); + var n = nodes.Count; + 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; + tree._phase[i] = nodes[i].Phase; + } + + var parent = new int[n]; + var kids = new List[n]; + var depth = new int[n]; + for (var i = 0; i < n; i++) kids[i] = new List(); + for (var i = 0; i < n; i++) + { + parent[i] = nodes[i].Parent.IsSet() && tree._index.TryGetValue(nodes[i].Parent, out var p) ? p : -1; + if (parent[i] >= 0) kids[parent[i]].Add(i); + } + for (var i = 0; i < n; i++) + for (var a = parent[i]; a >= 0; a = parent[a]) depth[i]++; + var roots = new List(); + for (var i = 0; i < n; i++) + if (parent[i] < 0) roots.Add(i); + tree._roots = roots.ToArray(); + var deepestFirst = new int[n]; + for (var i = 0; i < n; i++) deepestFirst[i] = i; + Array.Sort(deepestFirst, (x, y) => depth[y] != depth[x] ? depth[y].CompareTo(depth[x]) : x.CompareTo(y)); + + var body = new float[n]; // the body's own mass + var family = new float[n]; // the family's mass M_F + for (var i = 0; i < n; i++) family[i] = body[i] = nodes[i].Mass; + 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++) + { + tree._kids[i] = kids[i].ToArray(); + ring[i] = RingOf(kids[i], family, tree._distance); + } + + BuildReaches(tree, parent, body, family, ring, deepestFirst, out var lower, out var group); + for (var i = 0; i < n; i++) + if (body[i] > 0) ring[i] = null; // a ring stands in for a body only on a bodyless node + tree.EmitTerms(parent, body, family, ring, lower, group, settings, packRadius); + return tree; + } + + // The massive children of a node at its innermost massive distance, within RingWindow of it; null when there are none. + private static int[] RingOf(List kids, float[] family, float[] distance) + { + var inner = float.PositiveInfinity; + foreach (var k in kids) + if (family[k] > 0) inner = Math.Min(inner, distance[k]); + if (float.IsPositiveInfinity(inner)) return null; + var members = new List(); + foreach (var k in kids) + if (family[k] > 0 && distance[k] - inner <= RingWindow) members.Add(k); + return members.ToArray(); + } + + // The massive children of a bodyless node that stand in for its body: their mass, the ring's distance and its support. + private readonly struct RingGroup + { + public readonly float Mass, Inner, Support; + public RingGroup(float mass, float inner, float support) + { + Mass = mass; + Inner = inner; + Support = support; + } + } + + // 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) + { + var n = parent.Length; + var d = t._distance; + lower = new float[n]; + group = new RingGroup[n]; + t._support = new float[n]; + t._reach = new float[n]; + foreach (var i in order) + { + var low = 0f; + 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, 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). + var inner = float.PositiveInfinity; + var mass = 0f; + foreach (var j in ring[i]) + { + inner = Math.Min(inner, d[j]); + mass += family[j]; + } + var ringSupport = GravityField.Beta * inner; + foreach (var j in ring[i]) ringSupport = Math.Max(ringSupport, d[j] + t._reach[j]); + low = Math.Max(low, ringSupport); + support = Math.Max(support, GravityField.Beta * low); + group[i] = new RingGroup(mass, inner, ringSupport); + } + lower[i] = low; + t._support[i] = support; + if (family[i] <= 0) continue; + var floor = GravityField.Beta * support; + var p = parent[i]; + float reach; + 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); + t._reach[i] = reach > 0f ? reach : 0f; + } + } + + // 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 + // (Beta x max(pack radius, max over massive roots of |anchor - Centre| + support)) is not a change to this return: it needs + // every massive root's support, and roots are processed one at a time here, so it would take a second pass over the roots + // once their supports exist. + 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) + { + var n = parent.Length; + var terms = new List(); + _termStart = new int[n]; + _termCount = new int[n]; + for (var i = 0; i < n; i++) + { + _termStart[i] = terms.Count; + if (family[i] <= 0 || !(_reach[i] > 0f)) continue; + var global = float.IsPositiveInfinity(_reach[i]); + var fadeFrom = global ? float.PositiveInfinity : Math.Max(_support[i], GravityField.Gamma * _reach[i]); + void Add(float strength, float softening, float inner, float outer) => terms.Add(new Term(_key[i], strength, softening, inner, + outer, float.IsPositiveInfinity(outer) ? strength * GravityField.Potential(softening, inner, outer, packRadius) : 0f)); + if (body[i] > 0 && family[i] <= body[i]) + { + var (c, eps) = Curves(settings, body[i]); + Add(c, eps, fadeFrom, _reach[i]); + _termCount[i] = terms.Count - _termStart[i]; + continue; + } + if (body[i] > 0) + { + var (c, eps) = Curves(settings, body[i]); + Add(c, eps, lower[i], _support[i]); + } + else + { + var (c, eps) = Curves(settings, group[i].Mass); + Add(c, eps, lower[i], _support[i]); + Add(-c, eps, group[i].Inner, group[i].Support); + } + var (cf, ef) = Curves(settings, family[i]); + Add(cf, ef, fadeFrom, _reach[i]); + 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) => PullOf(_termStart[node], _termCount[node], r); + + private float PullOf(int start, int count, float r) + { + var pull = 0f; + for (var t = start; t < start + count; t++) + { + var term = _terms[t]; + pull += term.Strength * GravityField.Pull(term.Softening, term.Inner, term.Outer, r); + } + return pull; + } + + /// + /// The field's force at p (toward the wells, u/s^2): every node's terms, walked from the roots; a subtree is skipped + /// the moment p is outside its reach, which is exact. `positionAt(orbit, time)` is the caller's, so this never owns time. + /// + 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; + } + + private void ForceFrom(int node, float2 p, Func positionAt, double time, ref float2 force) + { + var offset = positionAt(_key[node], time) - p; + var r = length(offset); + if (r > 0f) + force += offset / r * PullOf(node, r); + foreach (var kid in _kids[node]) + if (length(positionAt(_key[kid], time) - p) < _reach[kid]) ForceFrom(kid, p, positionAt, time, ref force); + } + + /// The field's potential at p: every node's terms less their references, walked as Force walks. + 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; + } + + private float PotentialFrom(int node, float2 p, Func positionAt, double time) + { + var r = length(positionAt(_key[node], time) - p); + var sum = 0f; + for (var t = _termStart[node]; t < _termStart[node] + _termCount[node]; t++) + { + var term = _terms[t]; + sum += term.Strength * GravityField.Potential(term.Softening, term.Inner, term.Outer, r) - term.Reference; + } + foreach (var kid in _kids[node]) + if (length(positionAt(_key[kid], time) - p) < _reach[kid]) sum += PotentialFrom(kid, p, positionAt, time); + return sum; + } + + /// + /// The pull (toward the node, u/s^2) a rail at `distance` round `parentOrbit` is given: the node's own terms at that + /// distance. A ring member of a bodyless node adds each other member's own terms at their chord (positions at t = 0), + /// projected toward the node; a strengthless orbit on that ring (a station, a turret) takes the mean of what the ring's + /// members feel, so it keeps their period. Infinite period for a pull of 0 is the caller's. + /// + public float PullAround(CultRecordKey parentOrbit, float distance, CultRecordKey self) + { + var node = _index[parentOrbit]; + var pull = PullOf(node, distance); + var members = _ring[node]; + if (members == null || members.Length < 2) return pull; + if (_index.TryGetValue(self, out var me) && Array.IndexOf(members, me) >= 0) return pull + PartnersPull(members, me); + var inner = float.PositiveInfinity; + foreach (var m in members) inner = Math.Min(inner, _distance[m]); + if (Math.Abs(distance - inner) > RingWindow) return pull; + var mean = 0f; + foreach (var m in members) mean += PartnersPull(members, m) / members.Length; + return pull + mean; + } + + // What the other members of a ring pull member m with, along the direction from m to the node. + private float PartnersPull(int[] members, int m) + { + var at = _distance[m] * OrbitData.Evaluate(_phase[m]); + var sum = 0f; + foreach (var q in members) + { + if (q == m) continue; + var chord = _distance[q] * OrbitData.Evaluate(_phase[q]) - at; + var r = length(chord); + if (r <= 0f || _distance[m] <= 0f) continue; + sum += PullOf(q, r) * dot(chord / r, -at / _distance[m]); + } + return sum; + } +} diff --git a/Assets/Scripts/ServerShared/WellTree.cs.meta b/Assets/Scripts/ServerShared/WellTree.cs.meta new file mode 100644 index 000000000..6834ea3c2 --- /dev/null +++ b/Assets/Scripts/ServerShared/WellTree.cs.meta @@ -0,0 +1,2 @@ +fileFormatVersion: 2 +guid: da20befcb0f606b4ab1165988b836096 \ No newline at end of file diff --git a/Assets/Scripts/ServerShared/Zone.cs b/Assets/Scripts/ServerShared/Zone.cs index 87dd238ad..02904e461 100644 --- a/Assets/Scripts/ServerShared/Zone.cs +++ b/Assets/Scripts/ServerShared/Zone.cs @@ -141,6 +141,7 @@ public class Zone Pack = pack; _itemManager = itemManager; Settings = settings; + _orbitPositionAt = OrbitPositionAt; CombatSeed = galaxyZone?.Name.StableHash() ?? 1337u; var cache = itemManager.ItemData; @@ -187,18 +188,29 @@ public class Zone } } - // Orbits, bodies and belts exist: every rail's period is derived once, from the field it orbits in. - foreach (var orbit in Orbits) - { - var parent = orbit.Value.Data.Parent; - if (!parent.IsSet()) continue; - if (!_childOrbits.TryGetValue(parent.Key, out var siblings)) - _childOrbits[parent.Key] = siblings = new List(); - siblings.Add(orbit.Key); - } + // Orbits, bodies and belts exist: the field's terms are built once, and every rail's period is derived from them. + var wellMass = new Dictionary(); foreach (var body in PlanetInstances) if (body.Value.WellStrength > 0) - _wellOn.TryAdd(Planets[body.Key].Orbit.Key, body.Value); + 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) + { + 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() ? PeriodAround(orbit.Value.Data.Parent.Key, orbit.Value.Data.Distance, orbit.Key) @@ -561,67 +573,19 @@ public class Zone return time == _time ? GetOrbitPosition(orbitID) : RailAt(orbitID, time); } - // orbit -> the well-bearing body on it, and orbit -> the orbits whose parent it is. Derived from the pack in the - // constructor, read only by PeriodAround and Strength. - private readonly Dictionary _wellOn = new Dictionary(); - private readonly Dictionary> _childOrbits = new Dictionary>(); - private readonly Dictionary _strengths = new Dictionary(); + /// The zone's gravity: the terms of its orbit tree (WellTree owns which exist and the walk that sums them). + public WellTree Wells { get; } - // The one owner of a rail's period (Orbit.Period, AsteroidBelt.Periods): 2 pi sqrt(d / a), a being the pull the - // field gives at distance d from the node being orbited. Round a body, that is the body's own well through - // GravityField. Round a bodyless node it is the ring law: the strength of the node's orbits inside d, plus the - // ring term for the orbits on d itself (within RingWindow of it), taken as a point at the node. A bodied self - // feels half the strength of its ring partners (as a binary does); a strengthless self (a station, a turret) - // takes the period its bodied ring partners have, the mean of what each of them feels. No pull gives an - // infinite period. - private const float RingWindow = 0.1f; + private readonly Func _orbitPositionAt; + // The one owner of a rail's period (Orbit.Period, AsteroidBelt.Periods): 2 pi sqrt(d / a), a being the pull the field + // gives at distance d from the node being orbited (WellTree.PullAround). No pull gives an infinite period. private float PeriodAround(CultRecordKey parentOrbit, float distance, CultRecordKey self) { - float a; - if (_wellOn.TryGetValue(parentOrbit, out var body)) - a = length(GravityField.Acceleration(body.WellStrength, body.RocheLimit, float2(distance, 0))); - else - { - var inside = 0f; - var ring = new List<(CultRecordKey Orbit, float Strength)>(); - if (_childOrbits.TryGetValue(parentOrbit, out var children)) - foreach (var child in children) - { - var childDistance = Orbits[child].Data.Distance; - var childStrength = Strength(child); - if (childDistance < distance - RingWindow) inside += childStrength; - else if (childDistance <= distance + RingWindow && childStrength > 0) ring.Add((child, childStrength)); - } - // Whose pull self feels: its own if it has strength, else each bodied member of the ring. - var ringSum = 0f; - foreach (var member in ring) ringSum += member.Strength; - var selfBodied = Strength(self) > 0; - var feelings = 0f; - var feelers = 0; - foreach (var member in ring) - if (!selfBodied || member.Orbit.Equals(self)) - { - feelings += (ringSum - member.Strength) / 2; - feelers++; - } - var strength = inside + (feelers > 0 ? feelings / feelers : 0f); - a = length(GravityField.Acceleration(strength, 0f, float2(distance, 0))); - } + var a = Wells.PullAround(parentOrbit, distance, self); return a > 0 && distance > 0 ? 2 * PI * sqrt(distance / a) : float.PositiveInfinity; } - // The well strength an orbit carries as a node: its own body's plus every orbit round it, recursively. - private float Strength(CultRecordKey orbit) - { - if (_strengths.TryGetValue(orbit, out var known)) return known; - var strength = _wellOn.TryGetValue(orbit, out var body) ? body.WellStrength : 0f; - if (_childOrbits.TryGetValue(orbit, out var children)) - foreach (var child in children) - strength += Strength(child); - return _strengths[orbit] = strength; - } - public float2 GetOrbitVelocity(CultRecordKey orbit) { if (Orbits.ContainsKey(orbit)) @@ -734,10 +698,6 @@ public class Zone return security; } - // The potential of a body's well at the zone's edge, subtracted so every well is zero there: the one owner of the - // reference. GetHeight sums with it and the renderer's lowering (GravityLowering) hands it to the shader unchanged. - public float WellReference(Planet body) => GravityField.Potential(body.WellStrength, body.RocheLimit, Pack.Radius); - // The extent of the populated area, from the zone's origin: owned here; presenters read it, nothing authors it. A static // bound on where anything can be, so it holds at every time and reads nothing the pack's Radius decides: for each body, // its orbit chain's reach plus its Roche limit; for each belt, the chain's reach to the orbit its asteroids circle (the @@ -769,19 +729,9 @@ public class Zone public float GetHeight(float2 position) => GetHeight(position, _time); - // The field at a time: each well sits where its body's rail puts it then. The sim reads it at _time, a presenter at + // The field at a time: each term sits where its node's rail puts it then. The sim reads it at _time, a presenter at // the drawn time, so a ship rides the same slope the drawn planet's well makes. - public float GetHeight(float2 position, double time) - { - var result = 0f; - foreach (var body in PlanetInstances.Values) - { - var offset = position - OrbitPositionAt(body.BodyData.Orbit.Key, time); - result += GravityField.Potential(body.WellStrength, body.RocheLimit, length(offset)) - WellReference(body); - } - - return result / Settings.HeightExaggeration; - } + public float GetHeight(float2 position, double time) => Wells.Potential(position, _orbitPositionAt, time) / Settings.HeightExaggeration; public float GetLight(float2 position) { @@ -800,16 +750,10 @@ public class Zone internal float SchematicCellSize => _itemManager.GameplaySettings.SchematicCellSize; internal float FieldCrossSection(AsteroidBeltData belt) => _itemManager.ItemData.Get(belt.Kind)?.CrossSection ?? 0f; - // The one source of gravitational acceleration (u/s^2): the sum of every body's well (GravityField). Ships, - // floating items, clouds and mines all read it, and so does Flow. It is the gradient of the potential GetHeight - // draws, so the two cannot disagree. - public float2 GetForce(float2 position) - { - var force = float2.zero; - foreach (var body in PlanetInstances.Values) - force += GravityField.Acceleration(body.WellStrength, body.RocheLimit, position - body.Orbit.Position); - return force; - } + // The one source of gravitational acceleration (u/s^2): the walk of the zone's terms (WellTree). Ships, floating items, + // clouds and mines all read it, and so does Flow. It is the gradient of the potential GetHeight draws, so the two + // cannot disagree. + public float2 GetForce(float2 position) => Wells.Force(position, _orbitPositionAt, _time); // The global flow velocity at `position` in zone time, in u/s (FlowField is the law). The slope term reads the // downhill direction, which is the force. @@ -861,12 +805,7 @@ public class Planet public void CalculateProperties() { BodyRadius = Settings.BodyRadius.Evaluate(BodyData.Mass) * BodyData.BodyRadiusMultiplier; - WellStrength = Settings.WellStrength.Evaluate(BodyData.Mass); - RocheLimit = Settings.RocheLimit.Evaluate(BodyData.Mass); - if (!(WellStrength >= 0f) || float.IsInfinity(WellStrength)) - throw new InvalidOperationException("A body's well strength must be finite and not negative (PlanetSettings.WellStrength)."); - if (WellStrength > 0 && !(RocheLimit > 0f && !float.IsInfinity(RocheLimit))) - throw new InvalidOperationException("A body's Roche limit, which is its core radius, must be positive and finite (PlanetSettings.RocheLimit)."); + (WellStrength, RocheLimit) = WellTree.Curves(Settings, BodyData.Mass); } } diff --git a/Assets/Scripts/Zone Display/GravityLowering.cs b/Assets/Scripts/Zone Display/GravityLowering.cs index b4ea1e34c..2c0a0c2b3 100644 --- a/Assets/Scripts/Zone Display/GravityLowering.cs +++ b/Assets/Scripts/Zone Display/GravityLowering.cs @@ -3,8 +3,8 @@ using UnityEngine; /// /// The renderer's lowering of the zone's potential (aetheria-gravity, cut grid-lowering). Zone.GetHeight owns the drawn -/// height: this hands its inputs (each well's position, strength, core radius and reference, and the exaggeration) to -/// Assets/Shaders/Gravity.cginc, which sums them per texel, and makes the one quad that carries the sum into every +/// height: this hands its inputs (each term's position, strength, softening, band edges and reference, and the +/// exaggeration) to Assets/Shaders/Gravity.cginc, which sums them per texel, and makes the one quad that carries the sum into every /// camera that culls layer 9. Nothing here decides a height. /// public static class GravityLowering @@ -14,30 +14,32 @@ public static class GravityLowering /// The quad's side in u: larger than any zone (10000) plus any view (4096); texels outside a camera cost nothing. public const float FieldExtent = 100000f; - // Mirrors the well record of Gravity.cginc: float2 position, strength, coreRadius, reference. + // Mirrors the term record of Gravity.cginc: float2 position, strength, softening, inner, outer, reference. [StructLayout(LayoutKind.Sequential)] - struct WellRecord + struct TermRecord { public Vector2 Position; public float Strength; - public float CoreRadius; + public float Softening; + public float Inner; + public float Outer; public float Reference; } - static readonly int Stride = Marshal.SizeOf(); + static readonly int Stride = Marshal.SizeOf(); static readonly int WellsId = Shader.PropertyToID("_GravityWells"); static readonly int CountId = Shader.PropertyToID("_GravityWellCount"); static readonly int ExaggerationId = Shader.PropertyToID("_HeightExaggeration"); - static WellRecord[] _records = new WellRecord[0]; + static TermRecord[] _records = new TermRecord[0]; /// - /// Fills one record per body of the zone (a null zone has none) and sets the shader globals. The buffer is the - /// caller's, (re)allocated to the body count (at least one, so the global is always bound). Wells are lowered at the - /// drawn time (DrawAhead.Time), the same T the planets and a ship's ride height use. + /// Fills one record per term of the zone's WellTree (a null zone has none) and sets the shader globals. The buffer is + /// the caller's, (re)allocated to the term count (at least one, so the global is always bound). Terms are lowered at + /// their node's rail at the drawn time (DrawAhead.Time), the same T the planets and a ship's ride height use. /// public static void Lower(Zone zone, float lead, ref ComputeBuffer wells) { - var count = zone?.PlanetInstances.Count ?? 0; + var count = zone?.Wells.Terms.Count ?? 0; var capacity = Mathf.Max(count, 1); if (wells == null || wells.count != capacity) { @@ -45,20 +47,21 @@ public static class GravityLowering wells = new ComputeBuffer(capacity, Stride); } - if (_records.Length != capacity) _records = new WellRecord[capacity]; - var i = 0; - if (zone != null) - foreach (var body in zone.PlanetInstances.Values) + if (_records.Length != capacity) _records = new TermRecord[capacity]; + for (var t = 0; t < count; t++) + { + var term = zone.Wells.Terms[t]; + 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 { - var well = DrawAhead.Well(zone, body, lead); - _records[i++] = new WellRecord - { - Position = new Vector2(well.Position.x, well.Position.y), - Strength = well.Strength, - CoreRadius = well.CoreRadius, - Reference = well.Reference - }; - } + Position = new Vector2(at.x, at.y), + Strength = term.Strength, + Softening = term.Softening, + Inner = term.Inner, + Outer = term.Outer, + Reference = term.Reference + }; + } wells.SetData(_records); Shader.SetGlobalBuffer(WellsId, wells); diff --git a/Assets/Shaders/Gravity.cginc b/Assets/Shaders/Gravity.cginc index 47a419518..5eacdec0f 100644 --- a/Assets/Shaders/Gravity.cginc +++ b/Assets/Shaders/Gravity.cginc @@ -4,38 +4,91 @@ // GENERATED by ServerShared/GravityShader.cs from GravityField; do not edit (GravityShaderTests fails on any difference). // The lowering of the zone's potential (aetheria-gravity, cut grid-lowering). The law is ServerShared's GravityField // and the height is Zone.GetHeight; this file decides nothing. GravityLowering.Lower fills the buffer and the globals -// from the zone (positions from Orbit.Position, reference from Zone.WellReference) each frame. +// from the zone's WellTree (positions from the nodes' rails at the drawn time) each frame. // GravityFieldParityVerify renders it and compares every texel with -Zone.GetHeight. -struct GravityWell +struct GravityTerm { float2 position; float strength; - float coreRadius; + float softening; + float inner; + float outer; float reference; }; -StructuredBuffer _GravityWells; +StructuredBuffer _GravityWells; int _GravityWellCount; float _HeightExaggeration; -// GravityField.Potential: c is the well's strength, R its core radius, r the distance. -float gravityPotential(float c, float R, float r) +// GravityField.OneMinusS and GravityField.S: the band profile B(t) = 1 - S(t) and the order-4 smoothstep, lowest first. +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 <= 0.4, per unit strength. +float gravityBandClosed(float tau, float t) { - return c == 0 ? 0 : - r < R - ? c * r * r / (2 * R * R) - : c * (log(r / R) + 0.5); + float q[7]; + q[6] = GravityOneMinusS[7]; + [unroll] for (int i = 6; i > 0; i--) q[i - 1] = GravityOneMinusS[i] - tau * q[i]; + float remainder = GravityOneMinusS[0] - tau * q[0]; + float whole = 0; + float part = 0; + [unroll] for (int j = 6; j >= 0; j--) + { + whole = whole + q[j] / (j + 1); + part = (part + q[j] / (j + 1)) * t; + } + return 0.5 * (remainder * log((t + tau) / (1 + tau)) - (whole - part)); } -// Zone.GetHeight at world xz: each well's potential minus its own reference, summed, divided by the exaggeration. +// 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; + float yj[8]; + yj[0] = 1; + [unroll] for (int j = 1; j < 8; j++) yj[j] = yj[j - 1] * y; + float w = 1 / (1 + tau); + float wk = w; + float yk = y; + float acc = 0; + for (int k = 0; k < 40; k++) + { + float inner = 0; + [unroll] for (int m = 4; m < 8; m++) inner += GravityS[m] * yj[m] / (m + k + 1); + acc += wk * yk * inner; + wk *= w; + yk *= y; + } + return -0.5 * acc; +} + +// GravityField.Potential: one band term, per unit strength; the global term when b is infinite. +float gravityTermPotential(float softening, float a, float b, float r) +{ + float s = r * r; + float e2 = softening * softening; + if (isinf(b)) return 0.5 * log(s + e2); + float a2 = a * a; + 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 <= 0.4 ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0)); + float t = (s - a2) / (b2 - a2); + 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. float gravityHeight(float2 p) { float result = 0; for (int i = 0; i < _GravityWellCount; i++) { - GravityWell well = _GravityWells[i]; - result += gravityPotential(well.strength, well.coreRadius, length(p - well.position)) - well.reference; + GravityTerm term = _GravityWells[i]; + float r = length(p - term.position); + if (r < term.outer) + result += term.strength * gravityTermPotential(term.softening, term.inner, term.outer, r) - term.reference; } return result / _HeightExaggeration; } diff --git a/tests/Aetheria.Shared.Tests/DrawnPoseTests.cs b/tests/Aetheria.Shared.Tests/DrawnPoseTests.cs index be622e298..acc7f8f72 100644 --- a/tests/Aetheria.Shared.Tests/DrawnPoseTests.cs +++ b/tests/Aetheria.Shared.Tests/DrawnPoseTests.cs @@ -16,7 +16,8 @@ using float3 = CultMath.float3; // Every drawn pose is evaluated at one drawn time, T = Zone.ExactTime + Clock.Lead (DrawAhead.Time). The operator's report: // coasting beside a planet and looking at it, the planet jumped back and forth every frame, because the planet was drawn at // the last step's orbit while the ship was drawn a lead ahead of its last step. Each rail below is spelled here from the -// field law outside the core (a = c / d, so T = 2 pi d / sqrt(c)), never by calling the production rail; the zone is built +// field law inside a node's outermost child (the node pulls as its body alone: a = c d / (d^2 + eps^2), T = 2 pi sqrt(d / a)), +// never by calling the production rail; the zone is built // by the production constructor and stepped by Zone.Update, and the clock is the production SimClock. public sealed class DrawnPoseTests : IDisposable { @@ -124,9 +125,12 @@ public sealed class DrawnPoseTests : IDisposable // ---- The rails, spelled from the field law ---- + // T = 2 pi sqrt(d / a), a the Plummer pull of a body alone (softening 2) at d. + private static float PeriodAt(float distance, float strength) => 2 * PI * sqrt(distance / (strength * distance / (distance * distance + 4f))); + private static float2 Around(float2 center, double time, float distance, float phase, float strength) { - var period = 2 * PI * distance / sqrt(strength); + var period = PeriodAt(distance, strength); var wrapped = (float) frac(time / period); return center + OrbitData.Evaluate(frac(wrapped + phase)) * distance; } @@ -135,13 +139,14 @@ public sealed class DrawnPoseTests : IDisposable private static float2 MoonAt(double time) => Around(PlanetAt(time), time, MoonDistance, MoonPhase, PlanetMass); private static float2 StationAt(double time) => Around(PlanetAt(time), time, StationDistance, StationPhase, PlanetMass); - // The field at a point at a time, with each well where its rail puts it then. + // The field at a point at a time: the zone's terms (WellTree owns which exist), each at the position this test spells + // for its node's rail at that time. private static float FieldAt(World w, float2 point, double time) { - var wells = new (float2 At, float Strength)[] { (SunCentre, SunMass), (PlanetAt(time), PlanetMass), (MoonAt(time), MoonMass) }; + float2 At(CultRecordKey node) => node.Equals(w.PlanetOrbit) ? PlanetAt(time) : node.Equals(w.MoonOrbit) ? MoonAt(time) : SunCentre; var sum = 0f; - foreach (var (at, strength) in wells) - sum += GravityField.Potential(strength, 2f, length(point - at)) - GravityField.Potential(strength, 2f, w.Zone.Pack.Radius); + foreach (var term in w.Zone.Wells.Terms) + sum += term.Strength * GravityField.Potential(term.Softening, term.Inner, term.Outer, length(point - At(term.Node))) - term.Reference; return sum / w.HeightExaggeration; } @@ -332,7 +337,7 @@ public sealed class DrawnPoseTests : IDisposable for (var i = 0; i < w.Asteroids.Length; i++) { var asteroid = w.Asteroids[i]; - var period = 2 * PI * asteroid.Distance / sqrt(PlanetMass); + var period = PeriodAt(asteroid.Distance, PlanetMass); var expected = OrbitData.Evaluate((float) frac(time / period + asteroid.Phase)) * asteroid.Distance + parent; AssertNear(expected, drawn[i].xy, 1e-3f, $"chunk {i} drawn"); Assert.Equal(w.Zone.ChunkPose(w.Belt, i, time), drawn[i]); @@ -342,34 +347,27 @@ public sealed class DrawnPoseTests : IDisposable } } - // The well the shader takes is its body's rail at the drawn time: the lowering copies DrawAhead.Well and decides - // nothing, so the drawn well and the drawn planet share one pose at a non-zero lead, and the strength, core and - // reference are the field law's (spelled here from the masses, the core 2 and the zone's edge), not the step's. + // The term the shader takes sits on its node's rail at the drawn time: the lowering copies DrawAhead.Orbit of the term's + // node and decides nothing, so a drawn term and the drawn planet share one pose at a non-zero lead. [Fact] - public void AWellIsLoweredAtTheDrawnTime() + public void ATermIsLoweredAtItsNodesDrawnRail() { var w = Build(); Step(w, 1f / 60f, 803); + Assert.Contains(w.Zone.Wells.Terms, t => t.Node.Equals(w.PlanetOrbit)); + Assert.Contains(w.Zone.Wells.Terms, t => t.Node.Equals(w.MoonOrbit)); foreach (var lead in new[] { 0f, .004f, .0166f }) { var time = w.Zone.ExactTime + lead; - foreach (var body in w.Zone.PlanetInstances.Values) + foreach (var term in w.Zone.Wells.Terms) { - var key = body.BodyData.Orbit.Key; - var (at, strength) = key.Equals(w.PlanetOrbit) ? (PlanetAt(time), PlanetMass) - : key.Equals(w.MoonOrbit) ? (MoonAt(time), MoonMass) - : (SunCentre, SunMass); - var well = DrawAhead.Well(w.Zone, body, lead); - AssertNear(at, well.Position, 1e-3f, $"well of the {strength} body at lead {lead}"); - Assert.Equal(strength, well.Strength, 3); - Assert.Equal(2f, well.CoreRadius, 3); - Assert.Equal(GravityField.Potential(strength, 2f, w.Zone.Pack.Radius), well.Reference, 3); + var at = term.Node.Equals(w.PlanetOrbit) ? PlanetAt(time) : term.Node.Equals(w.MoonOrbit) ? MoonAt(time) : SunCentre; + AssertNear(at, DrawAhead.Orbit(w.Zone, term.Node, lead), 1e-3f, $"term of a body at lead {lead}"); } } - var planet = w.Zone.PlanetInstances.Values.Single(b => b.BodyData.Orbit.Key.Equals(w.PlanetOrbit)); - Assert.True(length(DrawAhead.Well(w.Zone, planet, .0166f).Position - DrawAhead.Well(w.Zone, planet, 0f).Position) > .1f, "the well moves over a lead"); - AssertNear(w.Zone.Orbits[w.PlanetOrbit].Position, DrawAhead.Well(w.Zone, planet, 0f).Position, 0f, "at lead 0 the well is the step's memo"); + Assert.True(length(DrawAhead.Orbit(w.Zone, w.PlanetOrbit, .0166f) - DrawAhead.Orbit(w.Zone, w.PlanetOrbit, 0f)) > .1f, "the term moves over a lead"); + AssertNear(w.Zone.Orbits[w.PlanetOrbit].Position, DrawAhead.Orbit(w.Zone, w.PlanetOrbit, 0f), 0f, "at lead 0 the term is the step's memo"); } // A HUD marker round the ship the camera follows (an incoming lock, a contact's bearing, a ping's flash) is anchored on diff --git a/tests/Aetheria.Shared.Tests/FlowFieldTests.cs b/tests/Aetheria.Shared.Tests/FlowFieldTests.cs index 719ee1283..c5cb3b793 100644 --- a/tests/Aetheria.Shared.Tests/FlowFieldTests.cs +++ b/tests/Aetheria.Shared.Tests/FlowFieldTests.cs @@ -124,8 +124,8 @@ public sealed class FlowFieldTests : IDisposable } } - // Zone.Flow's slope term reads the owner's superposed force, not any one body's: two equal wells, probes off the - // axis between them where the nearest body's direction and the sum's differ. + // Zone.Flow's slope term reads the owner's force at the point: two equal wells 600 apart (two subzones, each felt within + // its reach of 300), probes inside each off the axis, where the direction is not along the axis. [Fact] public void TheSlopeTermFollowsTheSuperposedForce() { @@ -135,7 +135,7 @@ public sealed class FlowFieldTests : IDisposable WellFixture.Dig(_items.ItemData, pack, planets, 300f, 100f, new float2(300f, 0f)); var zone = new Zone(_items, planets, pack, new GalaxyZone { Name = "Two Wells", Owner = null }, null); var checkedPoints = 0; - foreach (var at in new[] { float2(100f, 400f), float2(-50f, -300f), float2(250f, 220f), float2(-420f, 90f) }) + foreach (var at in new[] { float2(-200f, 150f), float2(-350f, -120f), float2(250f, 220f), float2(420f, -90f) }) { var flow = zone.Flow(float3(at.x, 3f, at.y)); var expected = 1.5f * normalize(zone.GetForce(at)); diff --git a/tests/Aetheria.Shared.Tests/GravityShaderTests.cs b/tests/Aetheria.Shared.Tests/GravityShaderTests.cs index dcd0e5d22..963e6b1a7 100644 --- a/tests/Aetheria.Shared.Tests/GravityShaderTests.cs +++ b/tests/Aetheria.Shared.Tests/GravityShaderTests.cs @@ -4,7 +4,7 @@ using Xunit; // potential-owns-gravity (aetheria-gravity, cut grid-lowering): Assets/Shaders/Gravity.cginc is the output of // GravityShader.Source(), which reads GravityField's constants and writes the one lowering of the potential, the -// zero-strength guard, each well's reference and gravityHeight's sum. The committed file must equal it, character for +// band term's closed form and series, each term's reference and gravityHeight's sum. The committed file must equal it, character for // character (line ends aside): any added statement, respelling, flipped sign, swapped argument or macro fails here. // To change the shader, change GravityShader and run once with AETHERIA_REGENERATE_GRAVITY_CGINC=1. public sealed class GravityShaderTests @@ -17,7 +17,7 @@ public sealed class GravityShaderTests File.WriteAllText(path, GravityShader.Source()); var generated = GravityShader.Source(); // The law was generated, not an empty template agreeing with an empty file. - Assert.Contains("log(r / R) + 0.5", generated); + Assert.Contains("static const float GravityS[8] = { 0, 0, 0, 0, 35, -84, 70, -20 };", generated); Assert.Equal(generated, File.ReadAllText(path).Replace("\r\n", "\n")); } } diff --git a/tests/Aetheria.Shared.Tests/LockOnScenarioTests.cs b/tests/Aetheria.Shared.Tests/LockOnScenarioTests.cs index 5c0632787..0e326d7b0 100644 --- a/tests/Aetheria.Shared.Tests/LockOnScenarioTests.cs +++ b/tests/Aetheria.Shared.Tests/LockOnScenarioTests.cs @@ -35,6 +35,7 @@ public sealed partial class RunStartTests var lockedAt = -1; var bearingDot = 0f; + var stagedDrift = 0f; var strengths = new List(); for (var tick = 0; tick < 600 && (lockedAt < 0 || tick < lockedAt + 40); tick++) { @@ -44,13 +45,16 @@ public sealed partial class RunStartTests if (lockedAt < 0) { lockedAt = tick; - bearingDot = dot(normalize(paint.Bearing), toPlayer * -1f); + // The warning points at the hostile where it is at the lock; the hostile has left the staged line by then. + var live = normalize(hostile.Position.xz - player.Position.xz); + bearingDot = dot(normalize(paint.Bearing), live); + stagedDrift = degrees(acos(clamp(dot(live, -toPlayer), -1f, 1f))); } strengths.Add(paint.Strength); } } Assert.True(lockedAt >= 0, "the hostile never locked the player within 600 ticks of .1 s"); - Assert.True(bearingDot >= cos(radians(5f)), "the warning points at the hostile"); + Assert.True(bearingDot >= cos(radians(.5f)), $"the warning points at the hostile where it is (the hostile is {stagedDrift} deg off its staged line at the lock)"); Assert.True(strengths.First() < strengths.Last(), $"the lock tightens: Strength went {strengths.First()} to {strengths.Last()}"); } } diff --git a/tests/Aetheria.Shared.Tests/MiningCut3Tests.cs b/tests/Aetheria.Shared.Tests/MiningCut3Tests.cs index fdabaf6bd..fd522e357 100644 --- a/tests/Aetheria.Shared.Tests/MiningCut3Tests.cs +++ b/tests/Aetheria.Shared.Tests/MiningCut3Tests.cs @@ -134,7 +134,7 @@ public sealed partial class MiningCut3Tests : IDisposable var items = new ItemManager(cache, new ProvenanceLedger(), GameSettings(), _ => { }); var sunOrbit = cache.Upsert(new OrbitData { FixedPosition = sunAt }).Key; - var sun = cache.Upsert(new SunData { Orbit = new CultRecordRef(sunOrbit) }).Key; + var sun = cache.Upsert(new SunData { Orbit = new CultRecordRef(sunOrbit), Mass = 1f }).Key; var pack = new ZonePack { Orbits = { new CultRecordRef(sunOrbit) }, @@ -147,7 +147,7 @@ public sealed partial class MiningCut3Tests : IDisposable { var orbit = cache.Upsert(new OrbitData { FixedPosition = at }).Key; pack.Orbits.Add(new CultRecordRef(orbit)); - pack.Planets.Add(new CultRecordRef(cache.Upsert(new SunData { Orbit = new CultRecordRef(orbit) }).Key)); + pack.Planets.Add(new CultRecordRef(cache.Upsert(new SunData { Orbit = new CultRecordRef(orbit), Mass = 1f }).Key)); } var scene = new Scene { Cache = cache, Items = items, Hull = hull, Eye = eye, DimEye = dimEye, Pack = pack }; scene.Kinds.AddRange(kinds ?? new FieldKindData[0]); diff --git a/tests/Aetheria.Shared.Tests/OrbitsFromFieldTests.cs b/tests/Aetheria.Shared.Tests/OrbitsFromFieldTests.cs index 10de2c62e..fea0d84c6 100644 --- a/tests/Aetheria.Shared.Tests/OrbitsFromFieldTests.cs +++ b/tests/Aetheria.Shared.Tests/OrbitsFromFieldTests.cs @@ -9,76 +9,31 @@ using static CultMath.math; using float2 = CultMath.float2; using float3 = CultMath.float3; -// Aetheria gravity, orbits-from-field: every rail's period is derived from the field it orbits in. Round a body it is -// that body's own well; round a bodyless node it is the ring law (the strength inside the rail, plus half of what shares -// it). The zones here are hand-built packs read by the production constructor, with the shipped well curves -// (c = 45.08 m^.2, R = 7.32 m^.25). Each expected period is spelled here from the field law outside the core, T = 2 pi d / sqrt(c), -// with c and R read from the zone's own bodies and every distance outside the core. +// Aetheria gravity, orbits-from-field: every rail's period is derived from the field it orbits in (WellTree.PullAround). +// Round a body that is the body's own pull inside its outermost child; round a bodyless node a ring member also feels each +// other member's own pull at their chord, projected toward the node. The zones here are hand-built packs read by the +// production constructor, with the shipped well curves (c = 45.08 m^.2, eps = 7.32 m^.25). Each expected period is +// spelled here as T = 2 pi sqrt(d / a), a composed from GravityField.Pull and c and eps read from the zone's own bodies. public sealed class OrbitsFromFieldTests : IDisposable { - private readonly string _root = Path.Combine(Path.GetTempPath(), "aetheria-orbitsfield-" + Guid.NewGuid().ToString("N")); - private readonly List _caches = new List(); + private readonly ZoneWorlds _worlds = new ZoneWorlds(); - public OrbitsFromFieldTests() => Directory.CreateDirectory(_root); + public void Dispose() => _worlds.Dispose(); - public void Dispose() - { - foreach (var cache in _caches) cache.Dispose(); - Directory.Delete(_root, true); - } - - private static GameplaySettings GameSettings() => new GameplaySettings - { - DefaultEntitySettings = new EntitySettings(), - Tiers = new[] { new RarityTier { Name = "Common", Quality = .5f, Rarity = 0, Color = new float3(1, 1, 1) } }, - QualityPriceModifier = new ExponentialLerp() - }; - - private static PlanetSettings Curves() => new PlanetSettings - { - WellStrength = new ExponentialCurve { Exponent = .2f, Multiplier = 45.08f }, - RocheLimit = new ExponentialCurve { Exponent = .25f, Multiplier = 7.32f }, - BodyRadius = new ExponentialCurve { Exponent = .25f, Multiplier = 3f }, - LightRadius = new ExponentialCurve() - }; - - private sealed class World - { - public CultCache Cache; - public ItemManager Items; - public ZonePack Pack = new ZonePack { Radius = 20000f, Mass = 10000f }; - - public CultRecordKey Orbit(CultRecordKey parent, float distance, float phase = 0f) - { - var data = new OrbitData { Distance = distance, Phase = phase }; - if (parent.IsSet()) data.Parent = new CultRecordRef(parent); - var key = Cache.Upsert(data).Key; - Pack.Orbits.Add(new CultRecordRef(key)); - return key; - } - - public CultRecordKey Body(CultRecordKey orbit, float mass) - { - var key = Cache.Upsert(new PlanetData { Name = "Body", Orbit = new CultRecordRef(orbit), Mass = mass }).Key; - Pack.Planets.Add(new CultRecordRef(key)); - return key; - } - - public Zone Build() => new Zone(Items, Curves(), Pack, new GalaxyZone { Name = "OrbitsFromField", Owner = null }, null); - } - - private World NewWorld() - { - var dir = Path.Combine(_root, $"world{_caches.Count}"); - Directory.CreateDirectory(dir); - var cache = AetheriaStores.Open(Path.Combine(dir, "Aetheria.cc"), Path.Combine(dir, "run.cc"), catalogWritable: true); - _caches.Add(cache); - return new World { Cache = cache, Items = new ItemManager(cache, new ProvenanceLedger(), GameSettings(), _ => { }) }; - } + private ZoneWorlds.World NewWorld() => _worlds.New(); + private static PlanetSettings Curves() => ZoneWorlds.Curves(); private static float C(Zone zone, CultRecordKey body) => zone.PlanetInstances[body].WellStrength; private static float Core(Zone zone, CultRecordKey body) => zone.PlanetInstances[body].RocheLimit; - private static float Outside(float distance, float c) => 2 * PI * distance / sqrt(c); + + // The pull of a body's well as one body, at distance d: c times the band term's pull with no band (GravityField.Pull). + private static float Pull(float c, float roche, float d) => c * GravityField.Pull(roche, float.PositiveInfinity, float.PositiveInfinity, d); + private static float Pull(Zone zone, CultRecordKey body, float d) => Pull(C(zone, body), Core(zone, body), d); + + // A mass the zone has no body of (a family's total): its well from the authored curves, as the law reads them. + private static float FamilyPull(float mass, float d) => Pull(Curves().WellStrength.Evaluate(mass), Curves().RocheLimit.Evaluate(mass), d); + + private static float PeriodUnder(float pull, float distance) => 2 * PI * sqrt(distance / pull); private static void AssertPeriod(float expected, float actual, string what) => Assert.True(Math.Abs(actual / expected - 1f) <= 1e-4f, $"{what}: expected {expected}, got {actual}"); @@ -95,13 +50,12 @@ public sealed class OrbitsFromFieldTests : IDisposable w.Body(moonOrbit, 1f); var zone = w.Build(); - Assert.True(Core(zone, planet) < 150f && Core(zone, sun) < 3000f, "fixture: both rails run outside the core"); - AssertPeriod(Outside(150f, C(zone, planet)), zone.Orbits[moonOrbit].Period, "the moon, in its planet's field"); - AssertPeriod(Outside(3000f, C(zone, sun)), zone.Orbits[planetOrbit].Period, "the planet, in its sun's field"); + AssertPeriod(PeriodUnder(Pull(zone, planet, 150f), 150f), zone.Orbits[moonOrbit].Period, "the moon, in its planet's field"); + AssertPeriod(PeriodUnder(Pull(zone, sun, 3000f), 3000f), zone.Orbits[planetOrbit].Period, "the planet, in its sun's field"); } [Fact] - public void AChildInsideItsParentsCoreOrbitsAsASolidBody() + public void AChildInsideItsParentsRocheLimitOrbitsInThePlummerPull() { var w = NewWorld(); var sunOrbit = w.Orbit(default, 0f); @@ -113,9 +67,8 @@ public sealed class OrbitsFromFieldTests : IDisposable var far = w.Orbit(planetOrbit, coreRadius * .7f); var zone = w.Build(); - var solid = 2 * PI * coreRadius / sqrt(C(zone, planet)); // a = c d / R^2, so T = 2 pi R / sqrt(c) at every d inside - AssertPeriod(solid, zone.Orbits[near].Period, "the inner child"); - AssertPeriod(solid, zone.Orbits[far].Period, "the outer child"); + AssertPeriod(PeriodUnder(Pull(zone, planet, coreRadius * .3f), coreRadius * .3f), zone.Orbits[near].Period, "the inner child"); + AssertPeriod(PeriodUnder(Pull(zone, planet, coreRadius * .7f), coreRadius * .7f), zone.Orbits[far].Period, "the outer child"); } [Fact] @@ -130,11 +83,13 @@ public sealed class OrbitsFromFieldTests : IDisposable w.Body(outerOrbit, 100f); var zone = w.Build(); - AssertPeriod(Outside(3000f, C(zone, sun)), zone.Orbits[outerOrbit].Period, "the outer planet, with a heavy planet inside it"); + AssertPeriod(PeriodUnder(Pull(zone, sun, 3000f), 3000f), zone.Orbits[outerOrbit].Period, "the outer planet, with a heavy planet inside it"); } + // A binary on a bodyless node: a member feels its partner's own pull at their chord (twice the ring's radius), toward + // the node, and nothing of the node itself. [Fact] - public void ABinaryOrbitsUnderHalfItsPartnersPull() + public void ABinaryMemberFeelsItsPartnerAtTheChord() { var w = NewWorld(); var node = w.Orbit(default, 0f); @@ -145,10 +100,11 @@ public sealed class OrbitsFromFieldTests : IDisposable var zone = w.Build(); Assert.True(C(zone, bodyHeavy) > 1.5f * C(zone, bodyA), "fixture: unequal members"); - AssertPeriod(Outside(400f, C(zone, bodyHeavy) / 2), zone.Orbits[a].Period, "the light member"); - AssertPeriod(Outside(400f, C(zone, bodyA) / 2), zone.Orbits[heavy].Period, "the heavy member"); + AssertPeriod(PeriodUnder(Pull(zone, bodyHeavy, 800f), 400f), zone.Orbits[a].Period, "the light member"); + AssertPeriod(PeriodUnder(Pull(zone, bodyA, 800f), 400f), zone.Orbits[heavy].Period, "the heavy member"); } + // Conservation: the outer body feels one well of the whole mass inside it, c(M_pair + moon), not the sum of strengths. [Fact] public void AnOuterBodyFeelsTheWholePairAndItsMoons() { @@ -157,15 +113,14 @@ public sealed class OrbitsFromFieldTests : IDisposable var m1 = w.Orbit(node, 400f, 0f); var m2 = w.Orbit(node, 400f, .5f); var moon = w.Orbit(m1, 40f); - var body1 = w.Body(m1, 100f); - var body2 = w.Body(m2, 300f); - var moonBody = w.Body(moon, 5f); + w.Body(m1, 100f); + w.Body(m2, 300f); + w.Body(moon, 5f); var outer = w.Orbit(node, 3000f); w.Body(outer, 50f); var zone = w.Build(); - var whole = C(zone, body1) + C(zone, moonBody) + C(zone, body2); - AssertPeriod(Outside(3000f, whole), zone.Orbits[outer].Period, "the outer body, round the pair and the first one's moon"); + AssertPeriod(PeriodUnder(FamilyPull(405f, 3000f), 3000f), zone.Orbits[outer].Period, "the outer body, round the pair and the first one's moon"); } [Fact] @@ -191,7 +146,7 @@ public sealed class OrbitsFromFieldTests : IDisposable w.Body(planet2, 100f); var companion2 = w.Orbit(node, 500f, 1f / 6f); zone = w.Build(); - AssertPeriod(Outside(500f, C(zone, inner)), zone.Orbits[planet2].Period, "the planet"); + AssertPeriod(PeriodUnder(Pull(zone, inner, 500f), 500f), zone.Orbits[planet2].Period, "the planet"); Assert.Equal(zone.Orbits[planet2].Period, zone.Orbits[companion2].Period); } @@ -213,10 +168,10 @@ public sealed class OrbitsFromFieldTests : IDisposable var zone = w.Build(); Assert.True(C(zone, heavyBody) > 1.5f * C(zone, lightBody), "fixture: unequal members"); - AssertPeriod(Outside(400f + offset, (C(zone, lightBody) + C(zone, heavyBody)) / 4), zone.Orbits[companion].Period, + AssertPeriod(PeriodUnder((Pull(zone, lightBody, 800f) + Pull(zone, heavyBody, 800f)) / 2, 400f + offset), zone.Orbits[companion].Period, "the companion, between its two members"); - // An equal pair: the companion has exactly the members' period. + // An equal pair: the companion has the members' period. w = NewWorld(); node = w.Orbit(default, 0f); var a = w.Orbit(node, 400f, 0f); @@ -226,12 +181,13 @@ public sealed class OrbitsFromFieldTests : IDisposable companion = w.Orbit(node, 400f, .25f); zone = w.Build(); Assert.True(float.IsFinite(zone.Orbits[a].Period)); - Assert.Equal(zone.Orbits[a].Period, zone.Orbits[companion].Period); - Assert.Equal(zone.Orbits[b].Period, zone.Orbits[companion].Period); + AssertPeriod(zone.Orbits[a].Period, zone.Orbits[companion].Period, "the companion, against the first member"); + AssertPeriod(zone.Orbits[b].Period, zone.Orbits[companion].Period, "the companion, against the second member"); } - // The ring is the 0.1 u window about a rail's distance: a member just inside it is a ring partner (pulls half), one just - // outside is inside or beyond the rail (pulls whole, or not at all). + // The ring is the 0.1 u window about the innermost massive distance: a member just inside it is a ring partner (felt at + // the chord), one just outside it is not a partner and the rail has the node's own terms alone, which give a + // ring's own member next to no pull. [Theory] [InlineData(-.09f, true)] [InlineData(.09f, true)] @@ -246,8 +202,11 @@ public sealed class OrbitsFromFieldTests : IDisposable var partnerBody = w.Body(partner, 2000f); var zone = w.Build(); - var c = C(zone, partnerBody); - AssertPeriod(Outside(400f, onRing ? c / 2 : c), zone.Orbits[rail].Period, "the rail, with a partner at " + partnerOffset); + if (onRing) + AssertPeriod(PeriodUnder(Pull(zone, partnerBody, 800f + partnerOffset), 400f), zone.Orbits[rail].Period, + "the rail, with a partner at " + partnerOffset); + else + Assert.True(zone.Orbits[rail].Period > 1e6f, "a body outside the window is no partner: the rail has no pull of any size"); } [Fact] @@ -264,7 +223,8 @@ public sealed class OrbitsFromFieldTests : IDisposable Assert.True(float.IsPositiveInfinity(zone.Orbits[rail].Period), "a body beyond the rail does not pull it"); } - // A bodyless node's children read only that node's ring: the well of the body the node itself orbits does not enter. + // A bodyless node's children read only that node's own terms and their ring partners: the well of the body the node + // itself orbits does not enter, and a lone body on the node has no pull. [Fact] public void ABodylessNodesChildrenNeverReadItsParentsWell() { @@ -286,8 +246,8 @@ public sealed class OrbitsFromFieldTests : IDisposable var bodyA = w.Body(a, 100f); var bodyB = w.Body(b, 300f); zone = w.Build(); - AssertPeriod(Outside(200f, C(zone, bodyB) / 2), zone.Orbits[a].Period, "the first member, under its partner alone"); - AssertPeriod(Outside(200f, C(zone, bodyA) / 2), zone.Orbits[b].Period, "the second member, under its partner alone"); + AssertPeriod(PeriodUnder(Pull(zone, bodyB, 400f), 200f), zone.Orbits[a].Period, "the first member, under its partner alone"); + AssertPeriod(PeriodUnder(Pull(zone, bodyA, 400f), 200f), zone.Orbits[b].Period, "the second member, under its partner alone"); } [Fact] @@ -317,12 +277,12 @@ public sealed class OrbitsFromFieldTests : IDisposable var binaries = 0; for (uint seed = 1; seed <= 12; seed++) { - var dir = Path.Combine(_root, $"galaxy{seed}"); + var dir = Path.Combine(_worlds.Root, $"galaxy{seed}"); Directory.CreateDirectory(dir); var catalog = Path.Combine(dir, "Aetheria.cc"); File.Copy(TestCatalog.Repo, catalog); var cache = new CultCache(TestCatalog.Registry()); - _caches.Add(cache); + _worlds.Caches.Add(cache); cache.AddBackingStore(new GameCult.Caching.MessagePack.SingleFileMessagePackBackingStore(catalog), AetheriaStores.CatalogTypes); cache.AddBackingStore(new GameCult.Caching.MessagePack.SingleFileMessagePackBackingStore(Path.Combine(dir, "run.cc")), AetheriaStores.RunTypes); var authored = AuthoredSettings.Load(RunStartTests.FindRepoRoot()); diff --git a/tests/Aetheria.Shared.Tests/PotentialOwnerTests.cs b/tests/Aetheria.Shared.Tests/PotentialOwnerTests.cs index 3cd03747e..6805bc82e 100644 --- a/tests/Aetheria.Shared.Tests/PotentialOwnerTests.cs +++ b/tests/Aetheria.Shared.Tests/PotentialOwnerTests.cs @@ -73,8 +73,8 @@ public sealed class PotentialOwnerTests : IDisposable private static void Near(float expected, float actual, float relative) => Assert.InRange(actual, expected - Math.Abs(expected) * relative - 1e-5f, expected + Math.Abs(expected) * relative + 1e-5f); - // The law written out: potential c r^2 / 2R^2 inside the core, c (ln(r / R) + 1/2) beyond. - private static float Phi(float c, float core, float r) => r < core ? c * r * r / (2 * core * core) : c * (log(r / core) + .5f); + // The law of a lone body written out (its global term): the Plummer potential c ln(r^2 + eps^2) / 2, eps the Roche limit. + private static float Phi(float c, float eps, float r) => c * .5f * log(r * r + eps * eps); // A zone whose bodies have no authored gravity curves (most test zones) has no pull and no dip, and no NaN either. [Fact] @@ -103,12 +103,13 @@ public sealed class PotentialOwnerTests : IDisposable var zone = ZoneOf(planets, (mass, float2.zero)); var core = planets.RocheLimit.Evaluate(mass); foreach (var r in new[] { 3f * core, 10f * core, PackRadius / 2f }) - Near(planets.WellStrength.Evaluate(mass), length(zone.GetForce(new float2(r, 0))) * r, 1e-3f); + Near(planets.WellStrength.Evaluate(mass) * r * r / (r * r + core * core), length(zone.GetForce(new float2(r, 0))) * r, 1e-3f); } } + // The steepest pull of the Plummer well is at its softening, which is the Roche curve: it rises up to it and falls beyond. [Fact] - public void TheCoreRadiusIsTheRocheCurve() + public void TheSteepestPullIsAtTheRocheCurve() { var planets = Curves(); foreach (var mass in new[] { 29.2f, 14152.4f }) @@ -116,10 +117,8 @@ public sealed class PotentialOwnerTests : IDisposable var zone = ZoneOf(planets, (mass, float2.zero)); var core = planets.RocheLimit.Evaluate(mass); float Pull(float r) => length(zone.GetForce(new float2(0, r))); - var slope = Pull(.25f * core) / (.25f * core); - Near(slope, Pull(.5f * core) / (.5f * core), 1e-3f); - Near(slope, Pull(.9f * core) / (.9f * core), 1e-3f); - Assert.True(Pull(core) > Pull(.98f * core), "the pull still rises up to the core radius"); + Near(planets.WellStrength.Evaluate(mass) / (2f * core), Pull(core), 1e-3f); + Assert.True(Pull(core) > Pull(.98f * core), "the pull still rises up to the Roche limit"); Assert.True(Pull(core) > Pull(1.02f * core), "and falls beyond it"); } } @@ -146,10 +145,10 @@ public sealed class PotentialOwnerTests : IDisposable ZoneOf(Constants(50f), (100f, float2.zero)); } - // core-at-roche-and-exaggeration-up: growing the core to the Roche limit (3 m^.25 to 7.32 m^.25, 2.44x) and dividing the - // global K by that factor leaves the steepest drawn slope (|grad phi| / K, at the core) where it was, and draws every - // slope outside the old core 2.44x steeper. The "before" side is the previous authoring, stated here from the ruling: - // core 3 m^.25 and K 10. The "after" side is the shipped Settings.asset read through the production zone. + // core-at-roche-and-exaggeration-up: the Plummer well's steepest pull is 0.5 c / eps at the Roche limit (eps = 7.32 m^.25), + // half the harmonic core's c / R, so K halves with the law and the steepest drawn slope stays where the previous authoring + // had it. The "before" side is stated from the ruling: core 3 m^.25 and K 10. The "after" side is the shipped + // Settings.asset read through the production zone. Outside the old core the funnel draws steeper than it did. [Fact] public void ThePeakDrawnSlopeStaysAndTheFunnelDrawsSteeper() { @@ -165,10 +164,8 @@ public sealed class PotentialOwnerTests : IDisposable var peakBefore = c / oldCore / oldExaggeration; Assert.InRange(Drawn(planets.RocheLimit.Evaluate(mass)), peakBefore * .98f, peakBefore * 1.02f); - var rocheOverOldCore = planets.RocheLimit.Evaluate(mass) / oldCore; var farBefore = c / (3f * oldCore) / oldExaggeration; - Near(rocheOverOldCore, Drawn(3f * oldCore) / farBefore, 1e-2f); - Assert.InRange(Drawn(3f * oldCore) / farBefore, 2.44f * .99f, 2.44f * 1.01f); + Assert.InRange(Drawn(3f * oldCore) / farBefore, 2.9f, 2.98f); } } @@ -247,8 +244,9 @@ public sealed class PotentialOwnerTests : IDisposable } } + // A lone body is a global term: it pulls as the Plummer well at every range, with no rim and no gate. [Fact] - public void TheTailIsAPointMassWithNoRimAndNoGate() + public void TheGlobalTermIsPlummerAtEveryRange() { var planets = Curves(); var zone = ZoneOf(planets, (273.9f, float2.zero)); @@ -259,34 +257,11 @@ public sealed class PotentialOwnerTests : IDisposable { var at = new float2(r * cos(angle), r * sin(angle)); var force = zone.GetForce(at); - Near(c, length(force) * r, 1e-3f); + Near(c * r * r / (r * r + core * core), length(force) * r, 1e-3f); Assert.True(dot(normalize(force), normalize(-at)) > .9999f, "the pull points at the body"); } } - [Fact] - public void WellsSuperpose() - { - var planets = Curves(); - var a = (mass: 29.2f, at: new float2(-300f, 0f)); - var b = (mass: 14152.4f, at: new float2(300f, 0f)); - var both = ZoneOf(planets, a, b); - var onlyA = ZoneOf(planets, a); - var onlyB = ZoneOf(planets, b); - var probes = new[] - { - float2(-300f, 0f), float2(300f, 0f), float2(-310f, 4f), float2(305f, -9f), float2(0f, 0f), float2(0f, 250f), - float2(-150f, 100f), float2(150f, -100f), float2(900f, 700f), float2(-1200f, -300f), float2(40f, 3000f), float2(-300f, -40f) - }; - foreach (var p in probes) - { - var sum = onlyA.GetForce(p) + onlyB.GetForce(p); - Near(sum.x, both.GetForce(p).x, 1e-4f); - Near(sum.y, both.GetForce(p).y, 1e-4f); - Near(onlyA.GetHeight(p) + onlyB.GetHeight(p), both.GetHeight(p), 1e-4f); - } - } - [Fact] public void HeightIsThePotentialOverOneK() { @@ -297,7 +272,7 @@ public sealed class PotentialOwnerTests : IDisposable var c = planets.WellStrength.Evaluate(273.9f); var core = planets.RocheLimit.Evaluate(273.9f); Assert.Equal(0f, lone.GetHeight(new float2(PackRadius, 0)), 4); - Near(-c * (log(PackRadius / core) + .5f) / k, lone.GetHeight(float2.zero), 1e-4f); + Near(-c * .5f * log((PackRadius * PackRadius + core * core) / (core * core)) / k, lone.GetHeight(float2.zero), 1e-4f); var two = ZoneOf(planets, (29.2f, new float2(-300f, 0f)), (14152.4f, new float2(300f, 0f))); var c0 = planets.WellStrength.Evaluate(29.2f); @@ -305,65 +280,34 @@ public sealed class PotentialOwnerTests : IDisposable var c1 = planets.WellStrength.Evaluate(14152.4f); var core1 = planets.RocheLimit.Evaluate(14152.4f); const float step = .5f; - // Inside the two cores as well, out to their outer bands (0.5 R to R), where the potential is the quadratic bowl. + // Inside the two softenings as well, where the Plummer potential is smooth through the centre. var inCore = new[] { float2(-300f + .5f * core0, 0f), float2(-300f, .65f * core0), float2(-300f - .86f * core0, 0f), float2(-300f + .6f * core0, .6f * core0), float2(300f + .55f * core1, 0f), float2(300f, -.7f * core1), float2(300f - .86f * core1, 0f), float2(300f + .9f * core1, 0f) }; - foreach (var p in inCore.Concat(new[] { float2(-200f, 120f), float2(0f, 0f), float2(100f, -300f), float2(500f, 40f), float2(-290f, 12f), float2(310f, 20f), float2(700f, 700f), float2(-80f, -500f) })) + // Two suns 600 apart are two subzones: each is felt only inside its own reach (half the distance, 300), so the height + // there is that well's alone (plus the zone node's constant, which differences cancel) and the force is its alone. + var w0 = float2(-300f, 0f); + var w1 = float2(300f, 0f); + foreach (var p in inCore.Concat(new[] { float2(-200f, 120f), float2(100f, -300f), float2(500f, 40f), float2(-290f, 12f), float2(310f, 20f), float2(700f, 700f), float2(-80f, -500f) })) { - var expected = Phi(c0, core0, length(p - float2(-300f, 0f))) - Phi(c0, core0, PackRadius) - + Phi(c1, core1, length(p - float2(300f, 0f))) - Phi(c1, core1, PackRadius); - Near(expected, two.GetHeight(p) * k, 1e-3f); + var near0 = length(p - w0) < 140f; // inside Gamma x reach (150) the well is its Plummer pull whole + var near1 = length(p - w1) < 140f; + if (near0 || near1) + { + var (at, wc, wcore) = near0 ? (w0, c0, core0) : (w1, c1, core1); + var expected = Phi(wc, wcore, length(p - at)) - Phi(wc, wcore, 0f); + Near(expected, (two.GetHeight(p) - two.GetHeight(at)) * k, 1e-3f); + } // The downhill slope of the potential is the force: the central difference of height * K is -GetForce. var gradient = float2( (two.GetHeight(p + float2(step, 0)) - two.GetHeight(p - float2(step, 0))) * k / (2 * step), (two.GetHeight(p + float2(0, step)) - two.GetHeight(p - float2(0, step))) * k / (2 * step)); var force = two.GetForce(p); - Assert.InRange(length(gradient + force), 0f, length(force) * 1e-2f + 1e-4f); + Assert.InRange(length(gradient + force), 0f, length(force) * 1e-2f + 1e-4f + 3e-5f * k); // flat ground between the subzones reads the float32 resolution of the height, times K } } } - - // Inside the core the potential is the bowl c r^2 / 2R^2 all the way out to R: the drop from a lone body's centre to a point - // at r is exactly that, from the middle of the core to a hair inside its edge (the log form agrees with the bowl to second - // order at R, so only a tight tolerance sees a join that sits just inside it). - [Fact] - public void TheCoreIsTheBowlUpToItsEdge() - { - foreach (var k in new[] { 1f, 10f }) - { - var planets = Curves(k); - foreach (var mass in new[] { 29.2f, 14152.4f }) - { - var zone = ZoneOf(planets, (mass, float2.zero)); - var c = planets.WellStrength.Evaluate(mass); - var core = planets.RocheLimit.Evaluate(mass); - foreach (var fraction in new[] { .5f, .75f, .9f, .95f, .98f, .99f, .996f }) - foreach (var angle in new[] { 0f, 2.3f }) - { - var r = fraction * core; - var drop = (zone.GetHeight(new float2(r * cos(angle), r * sin(angle))) - zone.GetHeight(float2.zero)) * k; - Near(c * fraction * fraction / 2f, drop, 2e-5f); - } - } - } - } - - [Fact] - public void TheCoreMeetsTheTailAtItsRadius() - { - var planets = Curves(); - foreach (var mass in new[] { 29.2f, 14152.4f }) - { - var zone = ZoneOf(planets, (mass, float2.zero)); - var core = planets.RocheLimit.Evaluate(mass); - var inside = new float2(core * (1f - 1e-3f), 0); - var outside = new float2(core * (1f + 1e-3f), 0); - Assert.InRange(abs(zone.GetHeight(inside) - zone.GetHeight(outside)), 0f, abs(zone.GetHeight(outside)) * 1e-3f); - Assert.InRange(abs(length(zone.GetForce(inside)) - length(zone.GetForce(outside))), 0f, length(zone.GetForce(outside)) * 1e-3f); - } - } } diff --git a/tests/Aetheria.Shared.Tests/WellFixture.cs b/tests/Aetheria.Shared.Tests/WellFixture.cs index 08bed1ad4..60e2ea29e 100644 --- a/tests/Aetheria.Shared.Tests/WellFixture.cs +++ b/tests/Aetheria.Shared.Tests/WellFixture.cs @@ -2,7 +2,7 @@ using CultMath; using GameCult.Caching; // The one fixture that gives a test zone a slope: a body on a parentless orbit (the origin unless `at` says otherwise; call again for more bodies), with the -// curves a body needs. Gravity is the sum of body wells (GravityField), so a test that wants ground to slide on +// curves a body needs. Gravity is the zone's WellTree over its bodies' masses, so a body carries a mass (any: the constant curves ignore it), and a test that wants ground to slide on // digs it here and nowhere else. `strength` is the well's c (u^2/s^2) and `coreRadius` its softened core in u, which is // the Roche limit curve (the one curve of a body's core). internal static class WellFixture @@ -13,7 +13,7 @@ internal static class WellFixture settings.WellStrength = new ExponentialCurve { Constant = strength }; settings.RocheLimit = new ExponentialCurve { Constant = coreRadius }; var orbit = cache.Upsert(new OrbitData { FixedPosition = at }).Key; - var body = cache.Upsert(new PlanetData { Name = $"Well {at.x},{at.y}", Orbit = new CultRecordRef(orbit) }).Key; + var body = cache.Upsert(new PlanetData { Name = $"Well {at.x},{at.y}", Orbit = new CultRecordRef(orbit), Mass = 1f }).Key; pack.Orbits.Add(new CultRecordRef(orbit)); pack.Planets.Add(new CultRecordRef(body)); } diff --git a/tests/Aetheria.Shared.Tests/WellLawTests.cs b/tests/Aetheria.Shared.Tests/WellLawTests.cs new file mode 100644 index 000000000..44ac58041 --- /dev/null +++ b/tests/Aetheria.Shared.Tests/WellLawTests.cs @@ -0,0 +1,533 @@ +using System; +using System.Linq; +using CultMath; +using GameCult.Caching; +using Xunit; +using static CultMath.math; +using float2 = CultMath.float2; + +// aetheria-gravity, cut well-law: the conserving compact law. GravityField is the law of one band term; WellTree turns a zone's +// orbit tree into terms and walks them. The oracles here are written out from the law (the order-4 smoothstep, Simpson +// integration of the pull in double, merging a family into one body), never by calling the code under test twice. +public sealed class WellLawTests : IDisposable +{ + private readonly ZoneWorlds _worlds = new ZoneWorlds(); + + public void Dispose() => _worlds.Dispose(); + + // The Plummer pull of a term without a band: r / (r^2 + eps^2). + private static double Plummer(double eps, double r) => r / (r * r + eps * eps); + + // The order-4 smoothstep, spelled: S(t) = 35 t^4 - 84 t^5 + 70 t^6 - 20 t^7. + private static double S(double t) => t * t * t * t * (35 - 84 * t + 70 * t * t - 20 * t * t * t); + + private static double Simpson(Func f, double from, double to, int n = 4000) + { + var h = (to - from) / n; + var sum = f(from) + f(to); + for (var i = 1; i < n; i++) sum += f(from + i * h) * (i % 2 == 1 ? 4 : 2); + return sum * h / 3; + } + + public static readonly object[][] Terms = + { + new object[] { 79.8f, 1000f, 1500f }, + new object[] { 17f, 40f, 60f }, + new object[] { 79.8f, 0f, 1525f }, + new object[] { 1658f, 1000f, 1500f }, // tau 3: the series branch + new object[] { 59f, 11f, 22f }, // tau 10 + new object[] { 79.8f, float.PositiveInfinity, float.PositiveInfinity } + }; + + // The pull is the slope of the potential: the potential's drop between two radii is the integral of the pull between + // them, by Simpson in double, for terms on both branches, inside the band, across its edges and beyond it. + [Theory] + [MemberData(nameof(Terms))] + public void TheTermIsItsPotentialsSlope(float eps, float inner, float outer) + { + var span = float.IsPositiveInfinity(outer) ? 3000f : outer; + foreach (var (from, to) in new[] { (.2f, .55f), (.4f, .99f), (.9f, 1.3f), (.05f, 1.0f) }) + { + float r0 = from * span, r1 = Math.Min(to * span, float.IsPositiveInfinity(outer) ? 6000f : outer * 1.1f); + var drop = GravityField.Potential(eps, inner, outer, r1) - GravityField.Potential(eps, inner, outer, r0); + var integral = Simpson(r => GravityField.Pull(eps, inner, outer, (float) r), r0, r1); + Assert.InRange(drop, integral - 3e-5 - 1e-4 * Math.Abs(integral), integral + 3e-5 + 1e-4 * Math.Abs(integral)); + } + if (!float.IsPositiveInfinity(outer)) + { + Assert.Equal(0f, GravityField.Potential(eps, inner, outer, outer)); + Assert.Equal(0f, GravityField.Potential(eps, inner, outer, outer * 1.5f)); + Assert.Equal(0f, GravityField.Pull(eps, inner, outer, outer * 1.5f)); + } + } + + // The pull is the Plummer pull faded by 1 - S(t) across the band, S the order-4 smoothstep: a smoothstep of order 2 has a + // different profile, and the profile's first three derivatives vanish at both edges (the deficit goes as t^4). + [Fact] + public void ABandEdgeIsSmooth() + { + const float eps = 79.8f, a = 1000f, b = 1500f; + foreach (var t in new[] { 0f, .02f, .05f, .1f, .3f, .5f, .7f, .9f, .98f, 1f }) + { + var r = (float) Math.Sqrt(a * a + t * (b * b - a * a)); + var pull = GravityField.Pull(eps, a, b, r); + var expected = Plummer(eps, r) * (1 - S(t)); + Assert.InRange(pull, expected - 1e-5 * Plummer(eps, r), expected + 1e-5 * Plummer(eps, r)); + } + Assert.Equal((double) GravityField.Pull(eps, a, b, a), Plummer(eps, a), 6); + Assert.Equal(0f, GravityField.Pull(eps, a, b, b)); + var deficit = new Func(t => 1 - GravityField.Pull(eps, a, b, (float) Math.Sqrt(a * a + t * (b * b - a * a))) / Plummer(eps, Math.Sqrt(a * a + t * (b * b - a * a)))); + Assert.InRange(deficit(.1f) / deficit(.05f), 13.0, 15.0); // doubling t multiplies the deficit by about 14 (t^4), not by 4 + } + + // 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.8e-6 per unit strength of the double-precision Simpson integral of the pull. A closed form kept + // past tau 0.6 or a series cut short misses this by 1e-5 to 1e-3. Both branches are within the bound between tau 0.3 and 0.6, + // so the exact seam (SeriesFrom 0.4) inside that interval is pinned by the committed shader being the generated one, not here. + [Fact] + public void TheBandPotentialIsExactInSinglePrecision() + { + var checkedCases = 0; + foreach (var tau in new[] { .05, .3, .35, .39, .41, .45, .5, .55, .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.8e-6, oracle + 1.8e-6); + checkedCases++; + } + } + Assert.True(checkedCases > 400, "the sweep ran its cases"); + } + + private static float Plummer(PlanetSettings curves, float mass, float d) => + curves.WellStrength.Evaluate(mass) * d / (d * d + curves.RocheLimit.Evaluate(mass) * curves.RocheLimit.Evaluate(mass)); + + // A sun, a planet at 3000 and a moon at 150 round the planet, the moon's mass as given. + private Zone SunPlanetMoon(float moonMass, out CultRecordKey planetOrbit, out CultRecordKey moonOrbit, float planetMass = 100f) + { + var w = _worlds.New(); + var sunOrbit = w.Orbit(default, 0f); + w.Body(sunOrbit, 5000f); + planetOrbit = w.Orbit(sunOrbit, 3000f); + w.Body(planetOrbit, planetMass); + moonOrbit = w.Orbit(planetOrbit, 150f); + if (moonMass > 0f) w.Body(moonOrbit, moonMass); + return w.Build(); + } + + // Inside its outermost child a node pulls exactly as its body alone: the family terms cancel there, so a child's rail + // never feels its own mass through its parent. + [Fact] + public void InsideAFamilyItsParentPullsAsItsBodyAlone() + { + var curves = ZoneWorlds.Curves(); + foreach (var moonMass in new[] { 1f, 2f }) + { + var zone = SunPlanetMoon(moonMass, out var planetOrbit, out var moonOrbit); + foreach (var d in new[] { 20f, 60f, 150f }) + { + var expected = Plummer(curves, 100f, d); + Assert.InRange(zone.Wells.PullAround(planetOrbit, d, moonOrbit), expected * (1 - 1e-5f), expected * (1 + 1e-5f)); + } + } + } + + // Merging a family into one body of its mass in the same slot changes nothing beyond the family's support, and does + // inside it. + [Fact] + public void AFamilyMergedIntoOneBodyIsTheSameFieldOutsideIt() + { + var family = SunPlanetMoon(1f, out var planetOrbit, out _); + var merged = SunPlanetMoon(0f, out _, out _, planetMass: 101f); + var support = family.Wells.Support(planetOrbit); + Assert.True(support > 150f, "fixture: the support clears the moon"); + var planetAt = family.Orbits[planetOrbit].Position; + foreach (var (distance, angle) in new[] { (1.05f, .3f), (1.4f, 2.2f), (3f, 4.0f), (9f, 1.0f), (20f, 5.5f) }) + { + var p = planetAt + float2(support * distance * cos(angle), support * distance * sin(angle)); + var f = family.GetForce(p); + var m = merged.GetForce(p); + Assert.InRange(length(f - m), 0f, 1e-5f * length(f) + 1e-9f); + Assert.InRange(abs(family.GetHeight(p) - merged.GetHeight(p)), 0f, 1e-5f * (1f + abs(family.GetHeight(p)))); + } + var inside = planetAt + float2(100f, 0f); + 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"); + } + + // 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() + { + var curves = ZoneWorlds.Curves(); + var one = SunPlanetMoon(1f, out _, out _); + var w = _worlds.New(); + var sunOrbit = w.Orbit(default, 0f); + w.Body(sunOrbit, 5000f); + for (var i = 0; i < 5; i++) w.Body(w.Orbit(sunOrbit, 800f + 450f * i, i / 5f), 20.2f); + var five = w.Build(); + foreach (var zone in new[] { one, five }) + { + var far = 10f * zone.Pack.Radius; + foreach (var angle in new[] { 0f, 1.7f, 4.1f }) + Near(curves.WellStrength.Evaluate(5101f), length(zone.GetForce(float2(far * cos(angle), far * sin(angle)))) * far, 1e-3f); + } + } + + // 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, 3000f, .85f); + 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); + + // Each subzone reaches half the way to its nearest neighbour, or 1.5 of its support if more. + for (var i = 0; i < roots.Length; i++) + { + var nearest = Enumerable.Range(0, roots.Length).Where(j => j != i).Min(j => length(zone.Orbits[roots[j]].Position - zone.Orbits[roots[i]].Position)); + Assert.Equal(Math.Max(.5f * nearest, 1.5f * zone.Wells.Support(roots[i])), zone.Wells.Reach(roots[i]), 2); + } + + 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"); + + // The rim: from A to 1.5 A the zone's pull rises from nothing to its whole Plummer pull as the order-4 smoothstep of + // (d^2 - A^2) / ((1.5 A)^2 - A^2), A being the reach-derived radius above. + foreach (var (radius, angle) in new[] { (1.05f, .7f), (1.15f, 3.1f), (1.25f, 5.0f), (1.35f, 1.9f), (1.45f, 4.2f) }) + { + var p = centre + float2(span * radius * cos(angle), span * radius * sin(angle)); + var d = length(centre - p); + var t = ((double) d * d - (double) span * span) / (2.25 * span * span - (double) span * span); + var expected = (centre - p) / d * (float) (Plummer(curves, 6100f, d) * S(t)); + Assert.InRange(length(zone.GetForce(p) - expected), 0f, 2e-5f * Plummer(curves, 6100f, d)); + } + + 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, 87f, 43f); // the outermost asteroid is in the middle: neither the first nor the last listed + 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)); + var scale = Plummer(curves, 101f, r); // float rounding of the four cancelling terms is relative to the Plummer pull + Assert.InRange(zone.Wells.PullAround(planet, r, moon), expected - 1e-5 * scale, expected + 1e-5 * scale); + } + } + + // A wide family fades from Gamma x its reach, not its support: a planet at 2000 with a moon at 100 has support 150 and reach + // 1000 (half its distance), so it pulls as the whole family out to 500 and gives way from there as the order-4 smoothstep + // of (r^2 - 500^2) / (1000^2 - 500^2). + [Fact] + public void AWideFamilysPullFadesFromGammaTimesItsReach() + { + var w = _worlds.New(); + var sun = w.Orbit(default, 0f); + w.Body(sun, 5000f); + var planet = w.Orbit(sun, 2000f); + w.Body(planet, 100f); + var moon = w.Orbit(planet, 100f); + w.Body(moon, 1f); + var zone = w.Build(); + var curves = ZoneWorlds.Curves(); + + Assert.Equal(150f, zone.Wells.Support(planet), 3); + Assert.Equal(1000f, zone.Wells.Reach(planet), 3); + foreach (var r in new[] { 480f, 520f, 600f, 700f, 800f, 900f, 980f }) + { + var t = Math.Max(0.0, (r * r - 500.0 * 500) / (1000.0 * 1000 - 500.0 * 500)); + var expected = Plummer(curves, 101f, r) * (1 - S(t)); + var scale = Plummer(curves, 101f, r); + Assert.InRange(zone.Wells.PullAround(planet, r, moon), expected - 1e-5 * scale, expected + 1e-5 * scale); + } + } + + // 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); + var scale = Plummer(curves, 140f, r); // float rounding of the four cancelling terms is relative to the Plummer pull + Assert.InRange(zone.Wells.PullAround(node, r, default), expected - 1e-5 * scale, expected + 1e-5 * scale); + } + } + + // 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)); + + // A node's reach is half its distance, or one and a half times its support if that is more; the root's is infinite. + [Fact] + public void AReachIsHalfTheDistanceOrOneAndAHalfTheSupport() + { + var w = _worlds.New(); + var sunOrbit = w.Orbit(default, 0f); + w.Body(sunOrbit, 5000f); + var wide = w.Orbit(sunOrbit, 3000f); + w.Body(wide, 100f); + var tight = w.Orbit(sunOrbit, 100f); + w.Body(tight, 100f); + w.Body(w.Orbit(wide, 150f), 1f); + w.Body(w.Orbit(tight, 150f), 1f); + var zone = w.Build(); + + Assert.True(float.IsPositiveInfinity(zone.Wells.Reach(sunOrbit))); + Assert.Equal(1500f, zone.Wells.Reach(wide), 3); + Assert.Equal(225f, zone.Wells.Support(wide), 3); + Assert.Equal(337.5f, zone.Wells.Reach(tight), 3); + Assert.Equal(225f, zone.Wells.Support(tight), 3); + } + + // A lone body on a bodyless node is pulled by nothing of its own: the node's terms cancel where the ring is. + [Fact] + public void ALoneBodyOnABodylessNodeFeelsNoPullOfItsOwn() + { + var w = _worlds.New(); + var sunOrbit = w.Orbit(default, 0f); + w.Body(sunOrbit, 5000f); + var node = w.Orbit(sunOrbit, 3000f); + var alone = w.Orbit(node, 200f); + w.Body(alone, 100f); + var zone = w.Build(); + + Assert.Equal(0f, zone.Wells.PullAround(node, 200f, alone)); + Assert.True(float.IsPositiveInfinity(zone.Orbits[alone].Period)); + } + + // The walk skips a subtree outside its reach, which is exact: at points all over a tree with families, a ring and a belt of + // moons, the force and potential of the walk equal the flat sum over every term to float rounding of the terms. + [Fact] + public void TheWalkIsTheSum() + { + var w = _worlds.New(); + var sunOrbit = w.Orbit(default, 0f); + w.Body(sunOrbit, 5000f); + for (var i = 0; i < 5; i++) + { + var planet = w.Orbit(sunOrbit, 700f + 800f * i, frac(i * .3819f)); + w.Body(planet, 40f + 30f * i); + for (var j = 0; j < 3; j++) w.Body(w.Orbit(planet, 90f + 70f * j, frac(j * .6180339887f)), 1f + j); + } + 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() * 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 = (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; + var level = term.Strength * GravityField.Potential(term.Softening, term.Inner, term.Outer, r) - term.Reference; + potential += level; + scale += abs(pull) + abs(level); + } + Assert.InRange(length(zone.GetForce(p) - force), 0f, 2e-6f * scale + 1e-9f); + Assert.InRange(abs(zone.GetHeight(p) * zone.Settings.HeightExaggeration - potential), 0f, 2e-6f * scale + 1e-9f); + } + } + + // A ring of three on a bodyless node: each member feels both partners' own pulls at their chord (the side of the inscribed + // triangle), each projected toward the node (cos 30 degrees), so the projection and the chord both decide the period. + [Fact] + public void ARosetteMemberFeelsItsPartnersProjectedTowardTheNode() + { + var w = _worlds.New(); + var node = w.Orbit(default, 0f); + var members = new[] { w.Orbit(node, 400f, 0f), w.Orbit(node, 400f, 1f / 3f), w.Orbit(node, 400f, 2f / 3f) }; + var bodies = members.Select(m => w.Body(m, 120f)).ToArray(); + var zone = w.Build(); + + var chord = 400f * sqrt(3f); + var pull = 2f * zone.PlanetInstances[bodies[0]].WellStrength * GravityField.Pull(zone.PlanetInstances[bodies[0]].RocheLimit, + float.PositiveInfinity, float.PositiveInfinity, chord) * cos(PI / 6f); + foreach (var member in members) + Assert.InRange(zone.Orbits[member].Period, 2 * PI * sqrt(400f / pull) * (1 - 1e-4f), 2 * PI * sqrt(400f / pull) * (1 + 1e-4f)); + } +} diff --git a/tests/Aetheria.Shared.Tests/ZoneWorlds.cs b/tests/Aetheria.Shared.Tests/ZoneWorlds.cs new file mode 100644 index 000000000..d1eec2091 --- /dev/null +++ b/tests/Aetheria.Shared.Tests/ZoneWorlds.cs @@ -0,0 +1,85 @@ +using System; +using System.Collections.Generic; +using System.IO; +using System.Linq; +using CultMath; +using GameCult.Caching; +using Xunit; +using static CultMath.math; +using float2 = CultMath.float2; +using float3 = CultMath.float3; + +// The hand-built zones of the tree tests: packs of orbits and bodies read by the production Zone constructor with the shipped +// well curves (c = 45.08 m^.2, eps = 7.32 m^.25), each in its own store under one temporary root. +internal sealed class ZoneWorlds : IDisposable +{ + public readonly string Root = Path.Combine(Path.GetTempPath(), "aetheria-orbitsfield-" + Guid.NewGuid().ToString("N")); + public readonly List Caches = new List(); + + public ZoneWorlds() => Directory.CreateDirectory(Root); + + public void Dispose() + { + foreach (var cache in Caches) cache.Dispose(); + Directory.Delete(Root, true); + } + + private static GameplaySettings GameSettings() => new GameplaySettings + { + DefaultEntitySettings = new EntitySettings(), + Tiers = new[] { new RarityTier { Name = "Common", Quality = .5f, Rarity = 0, Color = new float3(1, 1, 1) } }, + QualityPriceModifier = new ExponentialLerp() + }; + + public static PlanetSettings Curves() => new PlanetSettings + { + WellStrength = new ExponentialCurve { Exponent = .2f, Multiplier = 45.08f }, + RocheLimit = new ExponentialCurve { Exponent = .25f, Multiplier = 7.32f }, + BodyRadius = new ExponentialCurve { Exponent = .25f, Multiplier = 3f }, + LightRadius = new ExponentialCurve(), + AsteroidSize = new ExponentialLerp { Minimum = 2f, Maximum = 9f, Exponent = 1.7f } + }; + + public sealed class World + { + public CultCache Cache; + public ItemManager Items; + public ZonePack Pack = new ZonePack { Radius = 20000f, Mass = 10000f }; + + public CultRecordKey Orbit(CultRecordKey parent, float distance, float phase = 0f) + { + var data = new OrbitData { Distance = distance, Phase = phase }; + if (parent.IsSet()) data.Parent = new CultRecordRef(parent); + var key = Cache.Upsert(data).Key; + Pack.Orbits.Add(new CultRecordRef(key)); + return key; + } + + public CultRecordKey Body(CultRecordKey orbit, float mass) + { + var key = Cache.Upsert(new PlanetData { Name = "Body", Orbit = new CultRecordRef(orbit), Mass = mass }).Key; + Pack.Planets.Add(new CultRecordRef(key)); + 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); + } + + public World New() + { + var dir = Path.Combine(Root, $"world{Caches.Count}"); + Directory.CreateDirectory(dir); + var cache = AetheriaStores.Open(Path.Combine(dir, "Aetheria.cc"), Path.Combine(dir, "run.cc"), catalogWritable: true); + Caches.Add(cache); + return new World { Cache = cache, Items = new ItemManager(cache, new ProvenanceLedger(), GameSettings(), _ => { }) }; + } +}