150 lines
6.8 KiB
C#
150 lines
6.8 KiB
C#
using System;
|
|
using CultMath;
|
|
using static CultMath.math;
|
|
|
|
/// <summary>
|
|
/// 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).
|
|
/// </summary>
|
|
public static class GravityField
|
|
{
|
|
/// <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;
|
|
|
|
/// <summary>The band is summed as a series for tau above this. float32 error per unit strength against the double integral is at most 1.8e-6 anywhere (worst 1.7e-6 just above the seam, where the series has not reached float32 resolution); the closed form is float32-clean only to about 0.45.</summary>
|
|
public const float SeriesFrom = .4f;
|
|
|
|
/// <summary>The series branch's length: its ratio is below 1/1.4, so the truncated tail is below 1.4e-6 at the seam and falls geometrically above it.</summary>
|
|
public const int SeriesTerms = 40;
|
|
|
|
/// <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 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 <= SeriesFrom ? BandClosed(tau, t) : BandSeries(tau, t);
|
|
|
|
private static float BandClosed(float tau, float t)
|
|
{
|
|
// (1 - S(t)) = (t + tau) Q(t) + R by synthetic division; I is Q's antiderivative.
|
|
var n = OneMinusS.Length - 1;
|
|
Span<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;
|
|
}
|
|
}
|