r"""
The Airy disk and the diffraction limit
============================================

George Biddell Airy worked out the exact Fraunhofer diffraction pattern of
a circular aperture, finding a bright central disk surrounded by faint
concentric rings rather than the sharp geometric image predicted by ray
optics. Lord Rayleigh turned Airy's pattern into a practical
resolving-power criterion: two point sources are just resolvable when the
central peak of one Airy pattern falls on the first dark ring of the
other, giving the diffraction limit :math:`\theta_{\min} \approx
1.22\,\lambda/D` that bounds every microscope, telescope, and camera lens.
:func:`~physicskit.optics.wave.circular_aperture` and
:func:`~physicskit.optics.wave.fraunhofer_diffraction` reproduce the Airy
pattern directly from a numerical Fourier transform -- no closed-form
Bessel function is needed -- and its first dark ring lands exactly at the
Rayleigh angle computed independently.
"""

import matplotlib.pyplot as plt
import numpy as np

from physicskit.optics.wave import circular_aperture, fraunhofer_diffraction, intensity

# %%
# A circular aperture, propagated to its far-field diffraction pattern
# -------------------------------------------------------------------------

wavelength = 0.5e-3  # mm
D = 0.2  # aperture diameter, mm
dx = 0.002
N = 512
z = 1000.0  # mm

aperture = circular_aperture((N, N), dx=dx, radius=D / 2.0)
U = fraunhofer_diffraction(aperture, wavelength=wavelength, z=z, dx=dx)
I = intensity(U)
I /= I.max()

x = (np.arange(N) - N // 2) * (wavelength * z / (N * dx))
x_first_zero = 1.22 * wavelength / D * z  # Rayleigh criterion, mapped to the screen

# %%
# The central disk's radius, found directly from the computed pattern
# (its first minimum), against the Rayleigh angle computed independently
# from :math:`D` and :math:`\lambda` alone.

profile = I[N // 2]
right_half = profile[N // 2 :]
first_min_index = np.argmin(right_half[: np.argmax(right_half[1:] > right_half[:-1]) + 2])
x_measured = x[N // 2 :][first_min_index]

fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(x, profile)
ax.axvline(x_first_zero, color="r", ls="--", label=r"$1.22\,\lambda z/D$")
ax.axvline(-x_first_zero, color="r", ls="--")
ax.set_xlabel("screen position x (mm)")
ax.set_ylabel("normalized intensity")
ax.legend()
ax.set_title("Airy pattern: central disk + first dark ring (Rayleigh criterion)")
fig.tight_layout()

print(f"aperture diameter D = {D} mm, wavelength = {wavelength} mm")
print(f"Rayleigh angular limit: 1.22 lambda/D = {1.22 * wavelength / D:.6f} rad")
print(f"predicted first dark ring at x = {x_first_zero:.4f} mm")
print(f"first minimum found directly in the computed pattern at x = {x_measured:.4f} mm")

# %%
# The full 2D Airy pattern: central disk and concentric diffraction rings
# ----------------------------------------------------------------------------
# The 1D cross-section above only cuts through the pattern's center; the
# actual diffraction-limited image of a point source is the full 2D Airy
# disk, with concentric bright rings surrounding the central maximum --
# exactly what limits the *angular* (not just linear) resolution of a real
# telescope or microscope aperture. Plotted on a log scale, both the bright
# central disk and the faint outer rings (each roughly 1.75 magnitudes
# fainter than the last) are visible together.

theta = np.linspace(0, 2 * np.pi, 200)
fig2, ax2 = plt.subplots(figsize=(5, 4.5))
im = ax2.imshow(
    np.log10(I + 1e-8),
    extent=(x[0], x[-1], x[0], x[-1]),
    origin="lower",
    cmap="inferno",
    vmin=-6,
    vmax=0,
)
fig2.colorbar(im, ax=ax2, label="log10 normalized intensity")
ax2.plot(x_first_zero * np.cos(theta), x_first_zero * np.sin(theta), "c--", lw=1, label="first dark ring (Rayleigh)")
ax2.set_xlabel("x (mm)")
ax2.set_ylabel("y (mm)")
ax2.legend(fontsize=8, loc="upper right")
ax2.set_title("Full 2D Airy pattern: central disk + concentric rings")
fig2.tight_layout()