"""
ch07_nested_simpson_2d.py

Two-dimensional integration using nested Simpson integration.

The program evaluates

        1   sin(x)
        /      /
       |      |       x^2
   I = |      |   --------- dy dx .
       |      |    y^2 + 2
       /      /
        0      0

The two-dimensional problem is reduced to repeated one-dimensional
integrations. The same Simpson refinement routine is used for the
inner and outer integrations.

The number of intervals is doubled until

    |I_(2n)-I_n|/15 < tolerance.

This is automatic global refinement, not local adaptive subdivision.

Based on code written by Alexander Godunov, October 2009.
Python version prepared for the companion website, 2026.

NumPy and Matplotlib are used only for visualization.
"""

import math
import numpy as np
import matplotlib.pyplot as plt


def simpson_refine(fun, a, b, eps):
    """
    Composite Simpson integration with automatic global step doubling.

    Returns
    -------
    integral
    error_estimate
    n_used
    """
    nmax = 1048576

    if eps <= 0.0:
        raise ValueError("simpson_refine: tolerance must be positive.")

    if abs(b - a) <= sys_float_min():
        return 0.0, 0.0, 0

    h = (b - a) / 2.0
    sn = h * (fun(a) + 4.0 * fun(a + h) + fun(b)) / 3.0

    n = 4

    while n <= nmax:
        h = (b - a) / n
        s2n = fun(a) + fun(b)

        for i in range(1, n):
            x = a + i * h
            s2n += (2.0 if i % 2 == 0 else 4.0) * fun(x)

        s2n = h * s2n / 3.0
        error_estimate = abs(s2n - sn) / 15.0

        if error_estimate <= eps:
            return s2n, error_estimate, n

        sn = s2n
        n *= 2

    raise RuntimeError(
        "simpson_refine: maximum number of intervals reached."
    )


def sys_float_min():
    """Smallest positive normalized float, used for a zero-width test."""
    return 2.2250738585072014e-308


def plot_integrand():
    xplot = np.linspace(0.0, 1.0, 140)
    yplot = np.linspace(0.0, math.sin(1.0), 140)
    xgrid, ygrid = np.meshgrid(xplot, yplot)

    zgrid = xgrid ** 2 / (ygrid ** 2 + 2.0)
    zgrid = np.where(ygrid <= np.sin(xgrid), zgrid, np.nan)

    fig = plt.figure()
    ax3d = fig.add_subplot(111, projection="3d")
    ax3d.plot_surface(xgrid, ygrid, zgrid, linewidth=0)
    ax3d.set_xlabel("x")
    ax3d.set_ylabel("y")
    ax3d.set_zlabel("f(x,y)")
    ax3d.set_title("Integrand over 0 <= y <= sin(x)")


def main():
    a = 0.0
    b = 1.0

    eps_outer = 1.0e-7
    eps_inner = 1.0e-8

    reference = 1.0344649764293148e-1
    nfun_2d = 0

    def integrand(x, y):
        nonlocal nfun_2d
        nfun_2d += 1
        return x * x / (y * y + 2.0)

    def y_lower(x):
        return 0.0

    def y_upper(x):
        return math.sin(x)

    def inner_integral(x):
        c = y_lower(x)
        d = y_upper(x)

        value, _, _ = simpson_refine(
            lambda y: integrand(x, y),
            c,
            d,
            eps_inner,
        )
        return value

    plot_integrand()

    integral, error_estimate, n_used = simpson_refine(
        inner_integral,
        a,
        b,
        eps_outer,
    )

    absolute_error = abs(integral - reference)

    print("Two-dimensional nested Simpson integration")
    print("Integral of x^2/(y^2+2), 0 <= x <= 1,")
    print("with 0 <= y <= sin(x)")
    print(f"Reference value      = {reference:15.7e}")
    print(f"Calculated value     = {integral:15.7e}")
    print(f"Estimated outer error= {error_estimate:12.4e}")
    print(f"Absolute error       = {absolute_error:12.4e}")
    print(f"Outer intervals      = {n_used:10d}")
    print(f"2D function calls    = {nfun_2d:10d}")

    plt.show()


if __name__ == "__main__":
    main()