#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's WellTree (positions from the nodes' rails at the drawn time) each frame.
// GravityFieldParityVerify renders it and compares every texel with -Zone.GetHeight.

struct GravityTerm
{
    float2 position;
    float strength;
    float softening;
    float inner;
    float outer;
    float reference;
};

StructuredBuffer<GravityTerm> _GravityWells;
int _GravityWellCount;
float _HeightExaggeration;

// GravityField.OneMinusS and GravityField.S: the band profile B(t) = 1 - S(t) and the order-4 smoothstep, lowest first.
static const float GravityOneMinusS[8] = { 1, -0, -0, -0, -35, 84, -70, 20 };
static const float GravityS[8] = { 0, 0, 0, 0, 35, -84, 70, -20 };

// GravityField.BandClosed: the band's potential at t for tau <= 0.4, per unit strength.
float gravityBandClosed(float tau, float t)
{
    float q[7];
    q[6] = GravityOneMinusS[7];
    [unroll] for (int i = 6; i > 0; i--) q[i - 1] = GravityOneMinusS[i] - tau * q[i];
    float remainder = GravityOneMinusS[0] - tau * q[0];
    float whole = 0;
    float part = 0;
    [unroll] for (int j = 6; j >= 0; j--)
    {
        whole = whole + q[j] / (j + 1);
        part = (part + q[j] / (j + 1)) * t;
    }
    return 0.5 * (remainder * log((t + tau) / (1 + tau)) - (whole - part));
}

// GravityField.BandSeries: the band's potential at t for tau > 0.4, a series in (1 - t) / (1 + tau).
float gravityBandSeries(float tau, float t)
{
    float y = 1 - t;
    float yj[8];
    yj[0] = 1;
    [unroll] for (int j = 1; j < 8; j++) yj[j] = yj[j - 1] * y;
    float w = 1 / (1 + tau);
    float wk = w;
    float yk = y;
    float acc = 0;
    for (int k = 0; k < 40; k++)
    {
        float inner = 0;
        [unroll] for (int m = 4; m < 8; m++) inner += GravityS[m] * yj[m] / (m + k + 1);
        acc += wk * yk * inner;
        wk *= w;
        yk *= y;
    }
    return -0.5 * acc;
}

// GravityField.Potential: one band term, per unit strength; the global term when b is infinite.
float gravityTermPotential(float softening, float a, float b, float r)
{
    float s = r * r;
    float e2 = softening * softening;
    if (isinf(b)) return 0.5 * log(s + e2);
    float a2 = a * a;
    float b2 = b * b;
    if (s >= b2) return 0;
    float tau = (a2 + e2) / (b2 - a2);
    if (s <= a2) return 0.5 * log((s + e2) / (a2 + e2)) + (tau <= 0.4 ? gravityBandClosed(tau, 0) : gravityBandSeries(tau, 0));
    float t = (s - a2) / (b2 - a2);
    return tau <= 0.4 ? gravityBandClosed(tau, t) : gravityBandSeries(tau, t);
}

// Zone.GetHeight at world xz: each term's potential minus its own reference, summed, divided by the exaggeration.
float gravityHeight(float2 p)
{
    float result = 0;
    for (int i = 0; i < _GravityWellCount; i++)
    {
        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;
}

#endif