"""Independent focal-cell circle area by smooth Gauss quadrature, not the TS primitive.

Integrate in x = radius*sin(theta), split at intersections with cell edges.
The remaining integrands are smooth trig polynomials. The optical complex
field is still approximated as constant INSIDE each native focal cell.
"""
from functools import lru_cache
import math
import numpy as np


def cell_coverage(x, y, half, radius):
    if math.hypot(abs(x) + half, abs(y) + half) <= radius:
        return 1.0
    if math.hypot(max(0, abs(x) - half), max(0, abs(y) - half)) >= radius:
        return 0.0
    left, right = max(-radius, x - half), min(radius, x + half)
    lower, upper = y - half, y + half
    divisions = [left, right]
    for boundary in (lower, upper):
        if abs(boundary) < radius:
            crossing = math.sqrt(radius * radius - boundary * boundary)
            for point in (-crossing, crossing):
                if left < point < right:
                    divisions.append(point)
    divisions = sorted(set(divisions))
    nodes, weights = np.polynomial.legendre.leggauss(8)
    area = 0.0
    for first, second in zip(divisions, divisions[1:]):
        theta_first, theta_second = math.asin(first / radius), math.asin(second / radius)
        theta = (theta_first + theta_second) / 2 + nodes * (theta_second - theta_first) / 2
        height = radius * np.cos(theta)
        vertical = np.maximum(0, np.minimum(upper, height) - np.maximum(lower, -height))
        area += float(np.sum(weights * vertical * radius * np.cos(theta))) * (theta_second - theta_first) / 2
    return area / (4 * half * half)


@lru_cache(maxsize=100)
def circular_transmission(size, radius, sample):
    return np.asarray([[1 - cell_coverage((column - size / 2) * sample, (row - size / 2) * sample, sample / 2, radius) for column in range(size)] for row in range(size)])