The nearest-anchor rule (Reach x distance to the nearest other root) moves to WellTree.SubzoneReach(anchors, root) so the generator spells it once; the instance method returns max(it, floor). PullAround no longer gives a strengthless orbit on a ring a period of its own (binaries-orbited-from-afar); its test goes with it. Builds. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
435 lines
20 KiB
C#
435 lines
20 KiB
C#
using System;
|
|
using System.Collections.Generic;
|
|
using GameCult.Caching;
|
|
using CultMath;
|
|
using static CultMath.math;
|
|
using float2 = CultMath.float2;
|
|
|
|
/// <summary>
|
|
/// An orbit as the well tree reads it: where it hangs on the tree and the mass of the well-bearing body on it (0 for none).
|
|
/// Extent is the outermost distance from the parent of anything riding this orbit's rail (a belt orbit: its outermost
|
|
/// asteroid; any other orbit: its Distance). Anchor is a parentless orbit's stationary position, ignored with a parent.
|
|
/// </summary>
|
|
public readonly struct WellNode
|
|
{
|
|
public readonly CultRecordKey Key, Parent;
|
|
public readonly float Distance, Phase, Mass, Extent;
|
|
public readonly float2 Anchor;
|
|
|
|
public WellNode(CultRecordKey key, CultRecordKey parent, float distance, float phase, float mass, float extent, float2 anchor)
|
|
{
|
|
Key = key;
|
|
Parent = parent;
|
|
Distance = distance;
|
|
Phase = phase;
|
|
Mass = mass;
|
|
Extent = extent;
|
|
Anchor = anchor;
|
|
}
|
|
}
|
|
|
|
/// <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 extent: inside it the node pulls exactly as its body alone. B_F, the support, bounds every descendant term:
|
|
/// beyond it the family pulls exactly as one body of mass M_F, out to max(B_F, Gamma Lam_F), and is gone at its reach
|
|
/// Lam_F. A zone's sole massive root has infinite reach, so its far field is exactly c(M_zone). A zone of several massive
|
|
/// roots (subzones) has one zone node, which is not an orbit: at the roots' centre of mass it owns the only global term,
|
|
/// c(M_zone) U_global less c(M_zone) U(eps, A, Beta A), nothing inside A and the whole zone as one body beyond Beta A;
|
|
/// each root is bounded to its own subzone (SubzoneReach).
|
|
/// A body with no massive descendants is one term c(m) U(eps(m), Gamma Lam, Lam). Beyond a node's reach its whole
|
|
/// subtree contributes exactly nothing, so Force and Potential skip it: the walk is exact, not an approximation.
|
|
/// </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, _extent;
|
|
private float2[] _anchor;
|
|
private int[] _massiveRoots;
|
|
private float _packRadius;
|
|
private int _zoneStart, _zoneCount;
|
|
private Term[] _terms;
|
|
|
|
/// <summary>Where the zone node's terms sit (the family-mass-weighted mean of the massive roots' anchors); meaningful when it has terms.</summary>
|
|
public float2 Centre { get; private set; }
|
|
|
|
/// <summary>Every term of the zone, grouped by node, the zone node's last (its Node is unset: it sits at Centre): what the lowering uploads and the shader sums.</summary>
|
|
public IReadOnlyList<Term> Terms => _terms;
|
|
|
|
private WellTree() { }
|
|
|
|
/// <summary>A node's reach in its parent's field (+inf for a zone's sole massive root; a subzone's radius for a root of several); 0 for an orbit that carries no mass.</summary>
|
|
public float Reach(CultRecordKey node) => _reach[_index[node]];
|
|
|
|
/// <summary>A node's support B: every term of its descendants lies inside it, so beyond it the family is one body.</summary>
|
|
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];
|
|
tree._extent = new float[n];
|
|
tree._anchor = new float2[n];
|
|
tree._packRadius = packRadius;
|
|
for (var i = 0; i < n; i++)
|
|
{
|
|
tree._extent[i] = nodes[i].Extent;
|
|
tree._anchor[i] = nodes[i].Anchor;
|
|
tree._index[nodes[i].Key] = i;
|
|
tree._key[i] = nodes[i].Key;
|
|
tree._distance[i] = nodes[i].Distance;
|
|
tree._phase[i] = nodes[i].Phase;
|
|
}
|
|
|
|
var parent = new int[n];
|
|
var kids = new List<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];
|
|
|
|
// Massive roots are the zone's subzones; with two or more, their centre of mass is the zone node's position.
|
|
var massive = new List<int>();
|
|
foreach (var r in roots)
|
|
if (family[r] > 0) massive.Add(r);
|
|
tree._massiveRoots = massive.ToArray();
|
|
var centre = float2.zero;
|
|
var total = 0f;
|
|
foreach (var r in massive)
|
|
{
|
|
centre += tree._anchor[r] * family[r];
|
|
total += family[r];
|
|
}
|
|
tree.Centre = total > 0f ? centre / total : float2.zero;
|
|
|
|
tree._kids = new int[n][];
|
|
var ring = tree._ring = new int[n][];
|
|
for (var i = 0; i < n; i++)
|
|
{
|
|
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 extent.
|
|
private static void BuildReaches(WellTree t, int[] parent, float[] body, float[] family, int[][] ring, int[] order,
|
|
out float[] lower, out RingGroup[] group)
|
|
{
|
|
var n = parent.Length;
|
|
var d = t._distance;
|
|
lower = new float[n];
|
|
group = new RingGroup[n];
|
|
t._support = new float[n];
|
|
t._reach = new float[n];
|
|
foreach (var i in order)
|
|
{
|
|
var low = 0f;
|
|
foreach (var k in t._kids[i]) low = Math.Max(low, t._extent[k]);
|
|
var support = GravityField.Beta * low;
|
|
foreach (var k in t._kids[i]) support = Math.Max(support, t._extent[k] + t._reach[k]);
|
|
if (family[i] > 0 && body[i] <= 0)
|
|
{
|
|
// Bodyless: its ring group stands in for a body (it has massive children, or it would carry no mass).
|
|
var inner = float.PositiveInfinity;
|
|
var mass = 0f;
|
|
foreach (var j in ring[i])
|
|
{
|
|
inner = Math.Min(inner, d[j]);
|
|
mass += family[j];
|
|
}
|
|
var ringSupport = GravityField.Beta * inner;
|
|
foreach (var j in ring[i]) ringSupport = Math.Max(ringSupport, d[j] + t._reach[j]);
|
|
low = Math.Max(low, ringSupport);
|
|
support = Math.Max(support, GravityField.Beta * low);
|
|
group[i] = new RingGroup(mass, inner, ringSupport);
|
|
}
|
|
lower[i] = low;
|
|
t._support[i] = support;
|
|
if (family[i] <= 0) continue;
|
|
var floor = GravityField.Beta * support;
|
|
var p = parent[i];
|
|
float reach;
|
|
if (p < 0) reach = t.SubzoneReach(i, floor);
|
|
else if (body[p] <= 0 && ring[p] != null && ring[p].Length >= 2 && Array.IndexOf(ring[p], i) >= 0)
|
|
reach = Math.Max(GravityField.RingReach * d[i], floor);
|
|
else reach = Math.Max(GravityField.Reach * d[i], floor);
|
|
t._reach[i] = reach > 0f ? reach : 0f;
|
|
}
|
|
}
|
|
|
|
// A root's reach: the floor (Beta x support) or the subzone reach, whichever is more. The subzone shape is decided in the
|
|
// public static SubzoneReach (multi-root-isolated-subzones): the generator calls the same function, so Reach x the nearest
|
|
// root's distance is spelled once and the roots' windows are fitted to it (root-window-floor).
|
|
private float SubzoneReach(int root, float floor)
|
|
{
|
|
var anchors = new List<float2>();
|
|
foreach (var other in _massiveRoots) anchors.Add(_anchor[other]);
|
|
return Math.Max(SubzoneReach(anchors, Array.IndexOf(_massiveRoots, root)), floor);
|
|
}
|
|
|
|
/// <summary>
|
|
/// A subzone root's reach among the zone's massive roots' anchors: a sole massive root is the whole zone (+inf, so its
|
|
/// global term is the far field); with several, each root is bounded to GravityField.Reach x the distance to the nearest
|
|
/// other root (multi-root-isolated-subzones). The superposed-inside option is not a change to this return: it needs every
|
|
/// root's support, and roots are processed one at a time.
|
|
/// </summary>
|
|
public static float SubzoneReach(IReadOnlyList<float2> anchors, int root)
|
|
{
|
|
if (anchors.Count < 2) return float.PositiveInfinity;
|
|
var nearest = float.PositiveInfinity;
|
|
for (var other = 0; other < anchors.Count; other++)
|
|
if (other != root) nearest = Math.Min(nearest, length(anchors[other] - anchors[root]));
|
|
return GravityField.Reach * nearest;
|
|
}
|
|
|
|
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];
|
|
}
|
|
_zoneStart = terms.Count;
|
|
if (_massiveRoots.Length >= 2)
|
|
{
|
|
var total = 0f;
|
|
var span = 0f;
|
|
foreach (var r in _massiveRoots)
|
|
{
|
|
total += family[r];
|
|
span = Math.Max(span, length(_anchor[r] - Centre) + _reach[r]);
|
|
}
|
|
var (c, eps) = Curves(settings, total);
|
|
terms.Add(new Term(default, c, eps, float.PositiveInfinity, float.PositiveInfinity,
|
|
c * GravityField.Potential(eps, float.PositiveInfinity, float.PositiveInfinity, packRadius)));
|
|
if (span > 0f) terms.Add(new Term(default, -c, eps, span, GravityField.Beta * span, 0f));
|
|
}
|
|
_zoneCount = terms.Count - _zoneStart;
|
|
_terms = terms.ToArray();
|
|
}
|
|
|
|
// One node's own terms at distance r from it, toward its centre.
|
|
private float PullOf(int node, float r) => PullOf(_termStart[node], _termCount[node], r);
|
|
|
|
private float PullOf(int start, int count, float r)
|
|
{
|
|
var pull = 0f;
|
|
for (var t = start; t < start + count; t++)
|
|
{
|
|
var term = _terms[t];
|
|
pull += term.Strength * GravityField.Pull(term.Softening, term.Inner, term.Outer, r);
|
|
}
|
|
return pull;
|
|
}
|
|
|
|
/// <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;
|
|
var toZone = Centre - p;
|
|
var rz = length(toZone);
|
|
if (_zoneCount > 0 && rz > 0f) force += toZone / rz * PullOf(_zoneStart, _zoneCount, rz);
|
|
foreach (var root in _roots) ForceFrom(root, p, positionAt, time, ref force);
|
|
return force;
|
|
}
|
|
|
|
private void ForceFrom(int node, float2 p, Func<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;
|
|
var rz = length(Centre - p);
|
|
for (var t = _zoneStart; t < _zoneStart + _zoneCount; t++)
|
|
sum += _terms[t].Strength * GravityField.Potential(_terms[t].Softening, _terms[t].Inner, _terms[t].Outer, rz) - _terms[t].Reference;
|
|
foreach (var root in _roots) sum += PotentialFrom(root, p, positionAt, time);
|
|
return sum;
|
|
}
|
|
|
|
private float PotentialFrom(int node, float2 p, Func<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. 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);
|
|
return pull;
|
|
}
|
|
|
|
// 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;
|
|
}
|
|
}
|