"""Draw the cubic against the chord, and the arc against the chord.

Two figures. The first is the worked example of the page: four samples at
1, 2, 4, 8, the straight line through the middle pair, and the Catmull-Rom
cubic through all four, with both readings at the halfway position marked.
The second shows why a straight line between two directions is the wrong
path: eleven equal steps along the chord, pushed back out onto the arc,
against eleven equal steps of angle.

Both figures come from the page's own equations rather than from a library,
and the script prints the two numbers the page quotes so that a wrong
coefficient cannot reach the page unnoticed.

The colours become custom properties on the way out, so the drawings follow
the page between light and dark. See the note on `save`.

Writes cubic.svg and slerp.svg beside this file.
"""
import numpy as np
import matplotlib
matplotlib.use("Agg")
matplotlib.rcParams["svg.fonttype"] = "path"
import matplotlib.pyplot as plt

TEAL, INDIGO, ROSE, AMBER = "#00706B", "#3A49A6", "#A93356", "#B36A1E"
INK, MUTED, GRID = "#0D1620", "#5A6874", "#D3DBE3"

PROPERTY = {TEAL: "--teal", INDIGO: "--indigo", ROSE: "--rose",
            AMBER: "--amber", INK: "--ink", MUTED: "--muted", GRID: "--line"}


def hermite(ym1, y0, y1, y2):
    """Equation (4): the four coefficients, in the order c0 to c3."""
    return (y0,
            0.5 * (y1 - ym1),
            ym1 - 2.5 * y0 + 2.0 * y1 - 0.5 * y2,
            0.5 * (y2 - ym1) + 1.5 * (y0 - y1))


def horner(c, mu):
    """Equation (5), evaluated from the inside out."""
    c0, c1, c2, c3 = c
    return ((c3 * mu + c2) * mu + c1) * mu + c0


def slerp(a, b, t):
    """Equation (7), for unit vectors a and b."""
    omega = np.arccos(np.clip(np.dot(a, b), -1.0, 1.0))
    return (np.sin((1.0 - t) * omega) / np.sin(omega) * a
            + np.sin(t * omega) / np.sin(omega) * b)


def save(fig, name):
    """Write the figure, then make it the page's own drawing.

    Matplotlib writes every colour into the SVG as a literal, so a figure
    saved the ordinary way carries fixed values and cannot follow the page
    between light and dark. Each literal becomes var(--prop, literal) here,
    which keeps the file readable on its own and lets the page recolour it.

    Ids are prefixed with the file's stem because both figures are inlined
    into one page and Matplotlib numbers its own from one in each.
    """
    import re

    fig.tight_layout()
    # transparent, or Matplotlib paints an opaque white rectangle behind
    # everything and the figure reads as a light card on a dark page.
    fig.savefig(name, transparent=True, bbox_inches="tight", pad_inches=0.04)
    text = open(name, encoding="utf-8").read()

    stem = re.sub(r"[^A-Za-z0-9]+", "-", name.rsplit(".", 1)[0]).strip("-")
    text = re.sub(r'(\sid=")([^"]+)"', lambda m: f'{m.group(1)}{stem}-{m.group(2)}"', text)
    text = re.sub(r'((?:xlink:)?href="#)([^"]+)"',
                  lambda m: f'{m.group(1)}{stem}-{m.group(2)}"', text)
    text = re.sub(r'(url\(#)([^)]+)\)', lambda m: f'{m.group(1)}{stem}-{m.group(2)})', text)

    for literal, prop in PROPERTY.items():
        text = re.sub(re.escape(literal), f"var({prop}, {literal.lower()})",
                      text, flags=re.IGNORECASE)

    with open(name, "w", encoding="utf-8") as f:
        f.write(text)
    print(f"{name} written")


# The cubic against the chord, on the page's own four samples.
SAMPLES = [1.0, 2.0, 4.0, 8.0]
coef = hermite(*SAMPLES)
half = horner(coef, 0.5)
print(f"coefficients c0..c3 = {coef}")
print(f"cubic at mu = 0.5: {half}, straight line: {0.5 * (SAMPLES[1] + SAMPLES[2])}")

fig, (axL, axR) = plt.subplots(1, 2, figsize=(7.2, 3.3),
                               gridspec_kw={"width_ratios": [1.15, 1.0]})

mu = np.linspace(0.0, 1.0, 400)
POS = [-1, 0, 1, 2]

for ax in (axL, axR):
    ax.plot([0.0, 1.0], [SAMPLES[1], SAMPLES[2]], lw=1.6, color=MUTED, zorder=2)
    ax.plot(mu, horner(coef, mu), lw=2.0, color=TEAL, zorder=3)
    ax.plot([0.5], [3.0], "o", ms=5, color=MUTED, zorder=5)
    ax.plot([0.5], [half], "o", ms=5, color=ROSE, zorder=5)
    ax.plot([0.5, 0.5], [half, 3.0], lw=1.0, color=GRID, zorder=1)
    ax.tick_params(which="both", colors=MUTED, labelsize=9)
    for side in ("top", "right"):
        ax.spines[side].set_visible(False)
    for side in ("left", "bottom"):
        ax.spines[side].set_color(GRID)

# Left: all four samples, which is what the cubic is built from.
axL.plot(POS, SAMPLES, "o", ms=6, color=INDIGO, zorder=4)
for x, y in zip(POS, SAMPLES):
    axL.annotate(f"{y:.0f}", (x, y + 0.3), ha="center", va="bottom",
                 fontsize=10, color=INDIGO)
axL.plot([0.0, 0.0], [0.0, SAMPLES[1]], lw=0.8, color=GRID, zorder=1)
axL.plot([1.0, 1.0], [0.0, SAMPLES[2]], lw=0.8, color=GRID, zorder=1)
axL.set_xlim(-1.35, 2.35)
axL.set_ylim(0.0, 9.2)
axL.set_xticks(POS)
axL.set_xlabel("sample position", fontsize=10, color=INK)
axL.set_ylabel("value", fontsize=10, color=INK)
axL.annotate("the four samples", (-1.2, 8.6), ha="left", va="top",
             fontsize=9.5, color=MUTED)

# Right: the interval the value is read from, at a readable scale.
axR.plot([0, 1], [SAMPLES[1], SAMPLES[2]], "o", ms=6, color=INDIGO, zorder=4)
axR.annotate("3", (0.47, 3.02), ha="right", va="bottom", fontsize=10, color=MUTED)
axR.annotate(f"{half:.4f}", (0.53, half - 0.02), ha="left", va="top",
             fontsize=10, color=ROSE)
axR.annotate("straight line through two", (0.02, 3.86), ha="left", va="top",
             fontsize=9.5, color=MUTED)
axR.annotate("cubic through four", (0.98, 2.16), ha="right", va="bottom",
             fontsize=9.5, color=TEAL)
axR.set_xlim(-0.03, 1.03)
axR.set_ylim(1.85, 4.15)
axR.set_xticks([0.0, 0.5, 1.0])
axR.set_xlabel("position between the middle pair", fontsize=10, color=INK)

save(fig, "cubic.svg")


# The chord against the arc. A hundred and twenty degrees, because at sixty
# the two point sets differ by 1.1 degrees and the drawing shows nothing.
OMEGA = np.radians(120.0)
# Centred on straight up, so the span reads as an arc rather than a fan.
MID = np.pi / 2
A0, A1 = MID - OMEGA / 2, MID + OMEGA / 2
a = np.array([np.cos(A0), np.sin(A0)])
b = np.array([np.cos(A1), np.sin(A1)])
steps = np.linspace(0.0, 1.0, 11)

fig, ax = plt.subplots(figsize=(6.0, 3.9))

arc = np.array([[np.cos(u), np.sin(u)] for u in np.linspace(A0, A1, 300)])
ax.plot(arc[:, 0], arc[:, 1], lw=1.6, color=GRID, zorder=1)
ax.plot([a[0], b[0]], [a[1], b[1]], lw=1.4, color=MUTED, zorder=2)

for v in (a, b):
    ax.annotate("", xy=tuple(v), xytext=(0, 0),
                arrowprops=dict(arrowstyle="-|>", lw=1.8, color=INDIGO,
                                shrinkA=0, shrinkB=0), zorder=3)
ax.annotate("b", (b[0] - 0.05, b[1] - 0.03), ha="right", va="top",
            fontsize=12, color=INDIGO)
ax.annotate("a", (a[0] + 0.05, a[1] - 0.03), ha="left", va="top",
            fontsize=12, color=INDIGO)

# Equal steps along the chord, pushed back out onto the circle.
chord = np.array([(1 - t) * a + t * b for t in steps])
pushed = chord / np.linalg.norm(chord, axis=1)[:, None]
ax.plot(chord[:, 0], chord[:, 1], "o", ms=4, color=MUTED, zorder=4)
ax.plot(pushed[:, 0], pushed[:, 1], "o", ms=5.5, color=ROSE, zorder=5)

# Equal steps of angle, which is what equation (7) produces. Drawn on a
# wider ring so the two sets can be told apart at a glance.
RING = 1.14
even = np.array([slerp(a, b, t) for t in steps]) * RING
ax.plot(even[:, 0], even[:, 1], "o", ms=5.5, color=TEAL, zorder=5)
outer_arc = np.array([[np.cos(u), np.sin(u)]
                      for u in np.linspace(A0, A1, 300)]) * RING
ax.plot(outer_arc[:, 0], outer_arc[:, 1], lw=1.0, color=GRID, zorder=1)

ax.annotate("chord, equal steps", (0.0, 0.44), ha="center", va="top",
            fontsize=9.5, color=MUTED)
ax.annotate("pushed onto the arc", (0.0, 0.92), ha="center", va="top",
            fontsize=9.5, color=ROSE)
ax.annotate("equal steps of angle", (0.0, RING + 0.03), ha="center",
            va="bottom", fontsize=9.5, color=TEAL)

# How far the two disagree, in degrees, at the worst step. The outer ring is
# drawn at a larger radius to keep the two sets apart, so the angle is read
# from the construction rather than from the plotted point.
angle_pushed = np.degrees(np.arctan2(pushed[:, 1], pushed[:, 0]))
angle_even = np.degrees(A0 + steps * OMEGA)
worst = np.max(np.abs(angle_pushed - angle_even))
print(f"worst disagreement: {worst:.2f} degrees over {np.degrees(OMEGA):.0f}")
print("angle per step, pushed:", np.round(np.diff(angle_pushed), 2))
print("angle per step, even:  ", round(np.degrees(OMEGA) / (len(steps) - 1), 2))

ax.set_aspect("equal")
ax.set_xlim(-1.22, 1.22)
ax.set_ylim(-0.06, 1.34)
ax.axis("off")

save(fig, "slerp.svg")
