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

Well law: the band term, WellTree owns every term; the old law, per-body sums and Gauss ring law deleted

Browse source 
GravityField is now the band term U(eps, inner, outer) (C4 smoothstep fade, closed form and a series for tau > 1).
WellTree is built once per Zone from the orbit tree: node terms, reach, support, the pruned Force/Potential walk and
PullAround (rails). Zone.GetHeight/GetForce/PeriodAround delegate to it; WellReference, _wellOn, _strengths,
Strength() and DrawAhead.Well are gone; GravityLowering and the shader template take the flat term list;
HeightExaggeration halves to 2.0491803.

Does not pass yet: Gravity.cginc is not regenerated and PotentialOwnerTests/OrbitsFromFieldTests still spell the old law.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Meta Vens 2026-10-16 17:56:25 +02:00
commit eeeb6db290
11 changed files with 664 additions and 213 deletions
Download patch file Download diff file Expand all files Collapse all files

18
Assets/Scripts/ServerShared/DrawAhead.cs
View file

@ -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;

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

@ -1,32 +1,147 @@
using System;
using CultMath;
using static CultMath.math;
using float2 = CultMath.float2;
/// <summary>
/// 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 &lt;= a^2, 1 - S(t) on the band with t = (s - a^2) / (b^2 - a^2),
/// 0 for s &gt;= 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 &gt; 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).
/// </summary>
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;
/// <summary>The smoothstep's order: the potential is C^Order on every band edge.</summary>
public const int Order = 4;
/// <summary>A family's charge is fully on by its support, Beta x its outermost child.</summary>
public const float Beta = 1.5f;
/// <summary>A family is felt whole out to max(support, Gamma x its reach) and fades from there to its reach.</summary>
public const float Gamma = .5f;
/// <summary>A body's reach in its parent's field, in its distance (the nominal reach before the fade's room).</summary>
public const float Reach = .5f;
/// <summary>A ring member of a bodyless node reaches RingReach x its distance: whole out to the ring's diameter.</summary>
public const float RingReach = 4f;
// 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);
/// <summary>The series branch's length: its ratio is below 1/2, so 40 terms reach float32's floor.</summary>
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)
/// <summary>S(t) as polynomial coefficients in t, lowest first: t^Order sum_j C(Order-1+j, j) (1-t)^j.</summary>
public static readonly float[] S = BuildS();
/// <summary>1 - S(t) as polynomial coefficients in t, lowest first: the band's profile B.</summary>
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;
}
/// <summary>The pull toward the centre at distance r, per unit strength.</summary>
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);
}
/// <summary>The potential at distance r, per unit strength: 0 at and beyond outer, the global log when outer is +inf.</summary>
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<float> 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<float> 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<float> 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;
}
}

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

@ -1,10 +1,12 @@
using System.Globalization;
using System.Linq;
/// <summary>
/// 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.
/// </summary>
public static class GravityShader
@ -12,48 +14,110 @@ public static class GravityShader
public const string Path = "Assets/Shaders/Gravity.cginc";
/// <summary>The whole file, with LF line ends.</summary>
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<GravityWell> _GravityWells;
StructuredBuffer<GravityTerm> _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;
}

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

@ -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.

355
Assets/Scripts/ServerShared/WellTree.cs Normal file
View file

@ -0,0 +1,355 @@
using System;
using System.Collections.Generic;
using GameCult.Caching;
using CultMath;
using static CultMath.math;
/// <summary>An orbit as the well tree reads it: where it hangs on the tree and the mass of the well-bearing body on it (0 for none).</summary>
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;
}
}
/// <summary>
/// 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.
/// </summary>
public sealed class WellTree
{
/// <summary>Orbits whose distances differ by no more than this are on one ring (a binary, a ring of siblings).</summary>
public const float RingWindow = .1f;
/// <summary>A band term of the field: Strength x GravityField's U(Softening, Inner, Outer) at Node's position, less Reference.</summary>
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<CultRecordKey, int> _index = new Dictionary<CultRecordKey, int>();
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;
/// <summary>Every term of the zone, grouped by node: what the lowering uploads and the shader sums.</summary>
public IReadOnlyList<Term> Terms => _terms;
private WellTree() { }
/// <summary>A node's reach in its parent's field (+inf for a root); 0 for an orbit that carries no mass.</summary>
public float Reach(CultRecordKey node) => _reach[_index[node]];
/// <summary>A node's support B: every term of its descendants lies inside it, so beyond it the family is one body.</summary>
public float Support(CultRecordKey node) => _support[_index[node]];
/// <summary>
/// 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.
/// </summary>
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);
}
/// <summary>
/// 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).
/// </summary>
public static WellTree Build(IReadOnlyList<WellNode> 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<int>[n];
var depth = new int[n];
for (var i = 0; i < n; i++) kids[i] = new List<int>();
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<int>();
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<int> 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<int>();
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<Term>();
_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;
}
/// <summary>
/// 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.
/// </summary>
public float2 Force(float2 p, Func<CultRecordKey, double, float2> 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<CultRecordKey, double, float2> 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);
}
/// <summary>The field's potential at p: every node's terms less their references, walked as Force walks.</summary>
public float Potential(float2 p, Func<CultRecordKey, double, float2> 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<CultRecordKey, double, float2> 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;
}
/// <summary>
/// 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.
/// </summary>
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;
}
}

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

@ -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<CultRecordKey>();
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<CultRecordKey, float>();
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<WellNode>();
foreach (var orbit in Orbits)
wellNodes.Add(new WellNode(orbit.Key, orbit.Value.Data.Parent.Key, orbit.Value.Data.Distance, orbit.Value.Data.Phase,
wellMass.TryGetValue(orbit.Key, out var mass) ? mass : 0f));
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<CultRecordKey, Planet> _wellOn = new Dictionary<CultRecordKey, Planet>();
private readonly Dictionary<CultRecordKey, List<CultRecordKey>> _childOrbits = new Dictionary<CultRecordKey, List<CultRecordKey>>();
private readonly Dictionary<CultRecordKey, float> _strengths = new Dictionary<CultRecordKey, float>();
/// <summary>The zone's gravity: the terms of its orbit tree (WellTree owns which exist and the walk that sums them).</summary>
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<CultRecordKey, double, float2> _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);
}
}