import math
import numpy
from mpmath import mp

def ellipse_S_partial(alpha, beta, phi):
    a = mp.sin(alpha)
    b = mp.sin(beta)
    a2 = a*a
    b2 = b*b
    a2m = 1.0 - a2
    b2m = 1.0 - b2
    c = (b * a2m) / (a * mp.sqrt(b2m))
    m = (a2 - b2) / b2m
    n = m / a2

    t = mp.atan((mp.tan(alpha)/mp.tan(beta)) * mp.tan(phi))

    full = mp.pi/2 - c * mp.ellippi(n, m)
    partial = phi - c * mp.ellippi(n, t, m)
    return partial / full

N = 256
table = numpy.empty((N, N), dtype=numpy.float32)
# Perform precomputation in high precision
with mp.workdps(100):
    for yi in range(N):
        for xi in range(N):
            # Apply coordinate mapping
            # x=0 is implicitly zero, no need to tabulate.
            x = (xi + 1) / N
            y = yi / (N-1)

            # Convert to angles
            alpha = 0.5
            beta = max(1e-10, y*alpha)
            phi = x * math.pi/2

            # Tabulate
            table[xi][yi] = ellipse_S_partial(alpha, beta, phi)

# Print table in C syntax
print(f"static const float table_ellipse_CDF[{N}*{N}] = {{")
for row in table:
    print("  " + " ".join(f"{float(p):.8f}f," for p in row))
print("};")