"""
ch07_improper_infinite.py

Improper integral with an infinite upper limit.

The program evaluates

             infinity
                /
               |  exp(-x) sin(x) dx = 1/2
               /
               0

by transforming [0,infinity) to [0,1):

    x = t/(1-t),        dx = dt/(1-t)^2.

Gauss-Legendre quadrature is then applied on [0,1]. Since the
Gauss-Legendre nodes do not include endpoints, t = 1 is never sampled.

Based on the Fortran implementation by Alexander Godunov.
Python version prepared for the companion website, 2026.

Matplotlib is used only for visualization.
"""

import math
import matplotlib.pyplot as plt


def f(x):
    """Original integrand on the semi-infinite interval."""
    return math.exp(-x) * math.sin(x)


def g(t):
    """Transformed integrand after x = t/(1-t)."""
    one_minus_t = 1.0 - t
    x = t / one_minus_t
    return f(x) / (one_minus_t * one_minus_t)


def legendre_pair(n, z):
    if n == 1:
        return z, 1.0

    p0 = 1.0
    p1 = z

    for j in range(2, n + 1):
        p2 = ((2.0 * j - 1.0) * z * p1 - (j - 1.0) * p0) / j
        p0 = p1
        p1 = p2

    return p1, p0


def gauss_legendre(fun, a, b, n):
    """n-point Gauss-Legendre quadrature on [a,b]."""
    tol = 1.0e-14
    max_iter = 100

    if n < 1:
        raise ValueError("gauss_legendre: n must be positive.")

    midpoint = 0.5 * (a + b)
    halfwidth = 0.5 * (b - a)
    integral = 0.0

    m = (n + 1) // 2

    for i in range(1, m + 1):
        z = math.cos(math.pi * (i - 0.25) / (n + 0.5))
        converged = False

        for _ in range(max_iter):
            pn, pnm1 = legendre_pair(n, z)
            dp = 1.0 if n == 1 else n * (z * pn - pnm1) / (z * z - 1.0)

            zold = z
            z = zold - pn / dp

            if abs(z - zold) <= tol:
                converged = True
                break

        if not converged:
            raise RuntimeError(
                "gauss_legendre: Newton iteration did not converge."
            )

        pn, pnm1 = legendre_pair(n, z)
        dp = 1.0 if n == 1 else n * (z * pn - pnm1) / (z * z - 1.0)
        weight = 2.0 / ((1.0 - z * z) * dp * dp)

        xleft = midpoint - halfwidth * z
        xright = midpoint + halfwidth * z

        if abs(z) <= tol:
            integral += weight * fun(midpoint)
        else:
            integral += weight * (fun(xleft) + fun(xright))

    return halfwidth * integral


def plot_integrands():
    # Original integrand over a finite window.
    xmax = 20.0
    nplot = 2000
    xplot = [i * xmax / nplot for i in range(nplot + 1)]
    yplot = [f(x) for x in xplot]

    plt.figure()
    plt.plot(xplot, yplot, linewidth=1.5)
    plt.xlabel("x")
    plt.ylabel("f(x)")
    plt.title("Original integrand: exp(-x) sin(x)")
    plt.grid(True)

    # Transformed integrand. The endpoint t = 1 is excluded.
    tmax = 0.98
    tplot = [i * tmax / nplot for i in range(nplot + 1)]
    gplot = [g(t) for t in tplot]

    plt.figure()
    plt.plot(tplot, gplot, linewidth=1.5)
    plt.xlabel("t")
    plt.ylabel("g(t)")
    plt.title("Transformed integrand after x = t/(1-t)")
    plt.grid(True)


def main():
    a = 0.0
    b = 1.0
    exact = 0.5
    nvalues = [8, 16, 32, 64]

    plot_integrands()

    print("Improper integral: infinite interval")
    print("Integral of exp(-x) sin(x) from 0 to infinity")
    print(f"Exact value          = {exact:15.7e}\n")

    print("        n       Gauss-Legendre")
    print("--------------------------------")

    for n in nvalues:
        integral = gauss_legendre(g, a, b, n)
        print(f"{n:9d}{integral:21.7e}")

    plt.show()


if __name__ == "__main__":
    main()