"""Original Gaussian/reference phase model; no experimental data.

Run from mirrored _media. Python 3.14, NumPy 2.3.5, Matplotlib 3.10.7.
The reference wavefront is planar and common optical-path phase is removed.
"""

from pathlib import Path

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np


OUTPUT = "gouy-phase-curvature-bias.png"


def phase(zeta, radial_ratio):
    """Relative phase at fixed r/w0 and axial position z/zR."""
    return -np.arctan(zeta) + radial_ratio**2 * zeta / (1 + zeta**2)


def main():
    np.testing.assert_allclose(phase(1, 1), -np.pi / 4 + 0.5)
    np.testing.assert_allclose(phase(1, 1) - phase(-1, 1), -np.pi / 2 + 1)
    np.testing.assert_allclose(phase(0, np.array([0, 0.5, 1])), 0)
    zeta = np.linspace(-3, 3, 601)
    plt.rcParams.update({"font.size": 10, "axes.spines.top": False, "axes.spines.right": False})
    fig, (left, right) = plt.subplots(1, 2, figsize=(8, 4.6), dpi=100, facecolor="white")
    fig.suptitle("Off-axis curvature can bias a Gouy-phase comparison", fontsize=13)
    for radial, color, label in [(0, "#007c91", "On axis"), (0.5, "#6554a4", r"$r/w_0=0.5$"), (1, "#ba4d00", r"$r/w_0=1$")]:
        left.plot(zeta, phase(zeta, radial), color=color, lw=2, label=label)
    left.set(xlabel=r"Axial position $z/z_R$", ylabel="Relative phase (rad)", xlim=(-3, 3), ylim=(-1.5, 1.5))
    left.set_title("Fixed-radius detector points", fontsize=11)
    left.legend(loc="upper right", frameon=False, fontsize=9)

    bias = zeta / (1 + zeta**2)
    right.plot(zeta, bias, color="#ba4d00", lw=2)
    right.fill_between(zeta, 0, bias, color="#ba4d00", alpha=0.1)
    right.scatter([-1, 1], [-0.5, 0.5], color="#ba4d00", zorder=4)
    right.annotate("+0.5 rad", (1, 0.5), xytext=(8, 10), textcoords="offset points", fontsize=9)
    right.annotate("−0.5 rad", (-1, -0.5), xytext=(-8, -18), ha="right", textcoords="offset points", fontsize=9)
    right.set(xlabel=r"Axial position $z/z_R$", ylabel="Curvature contribution (rad)", xlim=(-3, 3), ylim=(-0.65, 0.65))
    right.set_title(r"Bias at $r=w_0$: $\zeta/(1+\zeta^2)$", fontsize=11)
    for ax in (left, right):
        ax.axhline(0, color="#999999", lw=0.7)
        ax.grid(alpha=0.18)
    fig.text(0.5, 0.03, "Original model: planar reference, fixed detector radius and stable path phase. Not measured data.", ha="center", fontsize=9)
    fig.subplots_adjust(left=0.09, right=0.975, bottom=0.18, top=0.79, wspace=0.36)
    fig.savefig(Path.cwd() / OUTPUT, facecolor="white", transparent=False)
    plt.close(fig)


if __name__ == "__main__":
    main()