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..16a985e8d 100644 --- a/Assets/Scripts/Editor/GravityFieldParityVerify.cs +++ b/Assets/Scripts/Editor/GravityFieldParityVerify.cs @@ -120,10 +120,11 @@ namespace Aetheria.EditorTools 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. + // 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) { var cache = AetheriaStores.Open(Path.Combine(work, "catalog.cc"), Path.Combine(work, "run.cc"), catalogWritable: true); @@ -132,11 +133,20 @@ namespace Aetheria.EditorTools var items = new ItemManager(cache, new ProvenanceLedger(), settings.GameplaySettings, _ => { }); var positions = new[] { new Vector2(-130f, 85f), new Vector2(610f, -240f), new Vector2(670f, -270f), new Vector2(1700f, 300f) }; var masses = new[] { 14152f, 273.9f, 29.2f, 0f }; + var parents = new[] { -1, 0, 1, -1 }; var keys = new CultRecordKey[positions.Length]; var pack = new ZonePack { Radius = radius }; for (var i = 0; i < positions.Length; i++) { - var orbit = keys[i] = cache.Upsert(new OrbitData { FixedPosition = new float2(positions[i].x, positions[i].y) }).Key; + 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)); 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..9b7854b8e 100644 --- a/Assets/Scripts/ServerShared/GravityField.cs +++ b/Assets/Scripts/ServerShared/GravityField.cs @@ -1,32 +1,147 @@ +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 series branch's length: its ratio is below 1/2, so 40 terms reach float32's floor. + public const int SeriesTerms = 40; - // The acceleration at `offsetFromBody`, toward the body: zero at the body itself. - public static float2 Acceleration(float c, float coreRadius, float2 offsetFromBody) + /// 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 <= 1f ? 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..6172fa161 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,110 @@ 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("@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 <= 1, 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 > 1, 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 <= 1 ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0)); + float t = (s - a2) / (b2 - a2); + return tau <= 1 ? 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/Settings.cs b/Assets/Scripts/ServerShared/Settings.cs index cda5e1623..a9e018909 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..c0dd3d825 --- /dev/null +++ b/Assets/Scripts/ServerShared/WellTree.cs @@ -0,0 +1,355 @@ +using System; +using System.Collections.Generic; +using GameCult.Caching; +using CultMath; +using static CultMath.math; + +/// 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). +public readonly struct WellNode +{ + public readonly CultRecordKey Key, Parent; + public readonly float Distance, Phase, Mass; + + public WellNode(CultRecordKey key, CultRecordKey parent, float distance, float phase, float mass) + { + Key = key; + Parent = parent; + Distance = distance; + Phase = phase; + Mass = mass; + } +} + +/// +/// 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 distance: inside it the node pulls exactly as its body alone. B_F, the support, bounds every descendant term: +/// beyond it the family pulls exactly as one body of mass M_F, out to max(B_F, Gamma Lam_F), and is gone at its reach +/// Lam_F. A root's reach is infinite, so the zone's far field is exactly c(M_zone) however the generator split it. +/// 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; + private Term[] _terms; + + /// Every term of the zone, grouped by node: what the lowering uploads and the shader sums. + public IReadOnlyList Terms => _terms; + + private WellTree() { } + + /// A node's reach in its parent's field (+inf for a root); 0 for an orbit that carries no mass. + 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]; + for (var i = 0; i < n; i++) + { + 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]; + + 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 distance. + 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, d[k]); + var support = GravityField.Beta * low; + foreach (var k in t._kids[i]) support = Math.Max(support, d[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 = float.PositiveInfinity; + 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; + } + } + + 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]; + } + _terms = terms.ToArray(); + } + + // One node's own terms at distance r from it, toward its centre. + private float PullOf(int node, float r) + { + var pull = 0f; + for (var t = _termStart[node]; t < _termStart[node] + _termCount[node]; 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; + 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; + 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/Zone.cs b/Assets/Scripts/ServerShared/Zone.cs index 87dd238ad..9ec7b9eb4 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,16 @@ 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); + var wellNodes = new List(); + foreach (var orbit in Orbits) + wellNodes.Add(new WellNode(orbit.Key, orbit.Value.Data.Parent.Key, orbit.Value.Data.Distance, orbit.Value.Data.Phase, + wellMass.TryGetValue(orbit.Key, out var mass) ? mass : 0f)); + 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 +560,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 +685,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 +716,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 +737,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 +792,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..e9da4178a 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 = DrawAhead.Orbit(zone, term.Node, lead); + _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/tests/Aetheria.Shared.Tests/DrawnPoseTests.cs b/tests/Aetheria.Shared.Tests/DrawnPoseTests.cs index be622e298..355c57432 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 { @@ -126,7 +127,7 @@ public sealed class DrawnPoseTests : IDisposable private static float2 Around(float2 center, double time, float distance, float phase, float strength) { - var period = 2 * PI * distance / sqrt(strength); + var period = 2 * PI * sqrt(distance / (strength * distance / (distance * distance + 4f))); var wrapped = (float) frac(time / period); return center + OrbitData.Evaluate(frac(wrapped + phase)) * distance; } @@ -135,13 +136,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; } @@ -342,34 +344,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/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")); } }