#!/usr/bin/env python3
"""ABLAQ Studios: courtyard sunlight v1.0.0.

Python 3.10+, standard library only. No inputs are downloaded.
Usage:
  python3 -B courtyard-study.v1.py > results.csv
  python3 -B courtyard-study.v1.py --check

Fixed area 100 m². Ratios 1, 2, 4; wall heights 2.5, 5, 7.5 m;
solar elevations 15, 30, 45, 60, 75 degrees;
azimuths 0..90 in 15-degree steps, relative to the long plan axis.
315 synthetic cases. Four equally high, opaque vertical walls, open top,
flat floor, collimated direct sunlight. No diffuse or reflected light,
weather, material properties, heat transfer or temperature prediction.
Angles are inputs, not a time/location ephemeris.
"""
import csv
import itertools
import math
import sys

VERSION = "1.0.0"
AREA = 100.0
RATIOS = (1.0, 2.0, 4.0)
HEIGHTS = (2.5, 5.0, 7.5)
ALTITUDES = (15.0, 30.0, 45.0, 60.0, 75.0)
AZIMUTHS = (0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0)


def sunlight(area=100.0, ratio=1.0, height=5.0, altitude=45.0, azimuth=45.0):
    values = (area, ratio, height, altitude, azimuth)
    if not all(math.isfinite(v) for v in values):
        raise ValueError("Inputs must be finite.")
    if area <= 0 or ratio <= 0 or height < 0 or not 0 <= altitude <= 90:
        raise ValueError("Invalid geometry or solar elevation.")
    if altitude == 0:
        return 0.0
    w, l = math.sqrt(area / ratio), math.sqrt(area * ratio)
    d = 0.0 if altitude == 90 else height / math.tan(math.radians(altitude))
    a = math.radians(azimuth % 360)
    return min(1.0, max(0.0, w - abs(d * math.sin(a))) *
               max(0.0, l - abs(d * math.cos(a))) / area)


def ray_grid(area, ratio, height, altitude, azimuth, n=160):
    """Independent midpoint quadrature: trace each ray to its first wall.

    This deliberately does not calculate projected-shadow rectangle overlaps.
    The ray is clear only if its height at the first wall exceeds wall height.
    A ray grazing a wall top counts as blocked (zero-area boundary).
    """
    if altitude <= 0:
        return 0.0
    w, l = math.sqrt(area / ratio), math.sqrt(area * ratio)
    a, b = math.radians(azimuth), math.radians(altitude)
    sx, sy, sz = math.cos(b) * math.sin(a), math.cos(b) * math.cos(a), math.sin(b)
    count = 0
    for ix in range(n):
        x = w * (ix + 0.5) / n
        tx = math.inf if abs(sx) < 1e-14 else ((w - x) / sx if sx > 0 else -x / sx)
        for iy in range(n):
            y = l * (iy + 0.5) / n
            ty = math.inf if abs(sy) < 1e-14 else ((l - y) / sy if sy > 0 else -y / sy)
            if min(tx, ty) * sz > height:
                count += 1
    return count / (n * n)


def cases():
    for ratio, height, altitude, azimuth in itertools.product(RATIOS, HEIGHTS, ALTITUDES, AZIMUTHS):
        yield dict(area=AREA, ratio=ratio, height=height, altitude=altitude, azimuth=azimuth)


def verify():
    import json
    checks = 0
    def check(value, message):
        nonlocal checks
        if not value:
            raise AssertionError(message)
        checks += 1
    check(abs(sunlight(azimuth=0) - 0.5) < 1e-12, "square cardinal reference")
    check(abs(sunlight(ratio=4, azimuth=0) - 0.75) < 1e-12, "rectangle long-axis reference")
    check(sunlight(ratio=4, azimuth=90) < 1e-12, "rectangle short-axis reference")
    check(sunlight(height=0) == 1.0, "no-wall limit")
    check(sunlight(altitude=90) == 1.0, "zenith limit")
    check(sunlight(altitude=0) == 0.0, "horizon convention")
    all_cases = list(cases())
    for p in all_cases:
        f = sunlight(**p)
        check(0 <= f <= 1, "bounded fraction")
        check(abs(f - sunlight(**dict(p, azimuth=p["azimuth"] + 180))) < 1e-12, "opposite-direction symmetry")
        check(sunlight(**dict(p, height=p["height"] + 1)) <= f + 1e-12, "wall-height monotonicity")
        check(sunlight(**dict(p, altitude=p["altitude"] + 5)) >= f - 1e-12, "elevation monotonicity")
        check(abs(f - sunlight(**dict(p, area=p["area"] * 4, height=p["height"] * 2))) < 1e-12, "geometric scale invariance")
    # Off-grid configurations avoid only testing the published sweep.
    probes = [
        (100, 1, 5, 45, 45), (100, 2, 4, 37, 19),
        (100, 4, 5, 60, 130), (225, 0.6, 6, 52, 227),
        (100, 4, 7.5, 15, 71), (100, 1, 2.5, 75, 319),
    ]
    grid_results = []
    for n in (40, 80, 160, 320):
        errors = [abs(sunlight(*p) - ray_grid(*p, n=n)) for p in probes]
        # Two axis-aligned boundary strips give a conservative O(1/n) area bound.
        check(max(errors) <= 2 / n, "independent ray-grid error bound")
        grid_results.append({"grid": n, "max_absolute_fraction_error": max(errors)})
    print(json.dumps({"version": VERSION, "status": "passed", "cases": len(all_cases),
                      "checks": checks, "ray_grid": grid_results,
                      "scope": "Numerical and geometric verification only; not empirical validation."}, indent=2))


def write_csv():
    writer = csv.writer(sys.stdout, lineterminator="\n")
    writer.writerow(("area_m2", "length_width_ratio", "wall_height_m", "solar_elevation_deg",
                     "azimuth_from_length_axis_deg", "sunlit_fraction", "shaded_fraction"))
    for p in cases():
        f = sunlight(**p)
        writer.writerow((100, p["ratio"], p["height"], p["altitude"], p["azimuth"],
                         f"{f:.8f}", f"{1-f:.8f}"))


if __name__ == "__main__":
    if sys.argv[1:] == ["--check"]:
        verify()
    elif not sys.argv[1:]:
        write_csv()
    else:
        raise SystemExit(__doc__)