#!/usr/bin/env python3
"""ABLAQ-R05: Weather-weighted courtyard solar access, version 1.0.0.

Author and publisher: ABLAQ Studios. Scientific lead: Osama Khaled Ahmed Abdel-Qader.
Run: python3 -B weather-study.v1.py --inputs inputs --output results
Requires numpy, pandas and pvlib; exact versions accompany the release.
Original ABLAQ code: CC0-1.0. PVGIS inputs retain their own source attribution.
No network requests, field observations, thermal model or reflections.
"""
import argparse
import calendar
import csv
import gzip
import hashlib
import itertools
import json
import math
from pathlib import Path

import numpy as np
import pandas as pd
import pvlib

VERSION = "1.0.0"
LATITUDE, LONGITUDE = 32.3, 36.2
CONFIGS = tuple(itertools.product((1, 2, 4), (2.5, 5, 7.5), (0, 45, 90)))
SOURCES = ("sarah3", "era5")


def sha(path):
    return hashlib.sha256(path.read_bytes()).hexdigest()


def write_json(path, value):
    path.write_text(json.dumps(value, indent=2, allow_nan=False) + "\n", encoding="utf-8")


def write_csv(path, rows):
    with path.open("w", encoding="utf-8", newline="") as stream:
        writer = csv.DictWriter(stream, fieldnames=list(rows[0]), lineterminator="\n")
        writer.writeheader()
        for row in rows:
            writer.writerow({k: f"{v:.12g}" if isinstance(v, (float, np.floating)) else v
                             for k, v in row.items()})


def dimensions(ratio, height):
    return math.sqrt(100 / ratio), math.sqrt(100 * ratio), height


def sky_exact(width, length, height):
    x, y = width / height, length / height
    return 2 / (math.pi * x * y) * (
        .5 * math.log((1 + x*x) * (1 + y*y) / (1 + x*x + y*y))
        + x * math.sqrt(1 + y*y) * math.atan(x / math.sqrt(1 + y*y))
        + y * math.sqrt(1 + x*x) * math.atan(y / math.sqrt(1 + x*x))
        - x * math.atan(x) - y * math.atan(y))


def sky_moments(width, length, height, order=400):
    """Projected hemispherical integral; u=sin(elevation)^2, one quadrant."""
    nodes, weights = np.polynomial.legendre.leggauss(order)
    u, wu = (nodes + 1) / 2, weights / 2
    az, wa = (nodes + 1) * math.pi / 4, weights / 2
    distance = height * np.sqrt((1-u)/u)
    clear = (np.maximum(0, 1-distance[:, None]*np.cos(az)/width)
             * np.maximum(0, 1-distance[:, None]*np.sin(az)/length))
    return (float(wu @ clear @ wa), float((wu*np.sqrt(u)) @ clear @ wa))


def sky_kernel_moment(width, length, height, order=128):
    """Independent floor-to-opening integral of the additional zenith weight."""
    nodes, weights = np.polynomial.legendre.leggauss(order)
    x, y = (nodes+1)*width/2, (nodes+1)*length/2
    wx, wy = weights*width/2, weights*length/2
    radius2 = height*height + x[:, None]**2 + y[None, :]**2
    kernel = ((width-x[:, None])*(length-y[None, :])*height**3
              / (math.pi * radius2**2.5))
    return float(4/(width*length) * (wx @ kernel @ wy))


def vectors(times, altitude, method="spa"):
    if method == "spa":
        position = pvlib.solarposition.spa_python(
            times, LATITUDE, LONGITUDE, altitude=altitude,
            pressure=0, delta_t=None, how="numpy")
        # Geometric, topocentric elevation: no atmospheric refraction.
        el = np.radians(position["elevation"].to_numpy())
        az = np.radians(position["azimuth"].to_numpy())
        return np.column_stack((np.cos(el)*np.sin(az),
                                np.cos(el)*np.cos(az), np.sin(el)))
    # Separately expressed NOAA fractional-year approximation.
    days = np.asarray(times.dayofyear, dtype=float)
    hour = np.asarray(times.hour + times.minute/60 + times.second/3600)
    year_days = np.where(times.is_leap_year, 366., 365.)
    g = 2*np.pi/year_days*(days-1+(hour-12)/24)
    decl = (.006918-.399912*np.cos(g)+.070257*np.sin(g)
            -.006758*np.cos(2*g)+.000907*np.sin(2*g)
            -.002697*np.cos(3*g)+.00148*np.sin(3*g))
    eq = 229.18*(.000075+.001868*np.cos(g)-.032077*np.sin(g)
                 -.014615*np.cos(2*g)-.040849*np.sin(2*g))
    angle = np.radians((hour*60+eq+4*LONGITUDE)/4-180)
    p = math.radians(LATITUDE)
    return np.column_stack((-np.cos(decl)*np.sin(angle),
                            math.cos(p)*np.sin(decl)-math.sin(p)*np.cos(decl)*np.cos(angle),
                            math.sin(p)*np.sin(decl)+math.cos(p)*np.cos(decl)*np.cos(angle)))


def fraction(vec, ratio, height, bearing):
    width, length, _ = dimensions(ratio, height)
    a = math.radians(bearing)
    e, n, z = vec.T
    safe = np.where(z > 0, z, 1)
    dx = height*(e*math.cos(a)-n*math.sin(a))/safe
    dy = height*(n*math.cos(a)+e*math.sin(a))/safe
    return np.where(z > 0, np.maximum(0, 1-np.abs(dx)/width)
                    * np.maximum(0, 1-np.abs(dy)/length), 0.)


def fraction_polygon(vec, ratio, height, bearing):
    """Independent clipping/area implementation; does not use overlap lengths."""
    e, n, z = map(float, vec)
    if z <= 0:
        return 0.
    w, l, _ = dimensions(ratio, height)
    a = math.radians(bearing)
    dx = height*(e*math.cos(a)-n*math.sin(a))/z
    dy = height*(n*math.cos(a)+e*math.sin(a))/z
    polygon = [(0., 0.), (w, 0.), (w, l), (0., l)]
    for axis, bound, sign in [(0, -dx, 1), (0, w-dx, -1),
                              (1, -dy, 1), (1, l-dy, -1)]:
        clipped = []
        if not polygon:
            return 0.
        previous = polygon[-1]
        pv = sign*(previous[axis]-bound)
        for current in polygon:
            cv = sign*(current[axis]-bound)
            if (cv >= 0) != (pv >= 0):
                t = pv/(pv-cv)
                clipped.append(tuple(previous[k]+t*(current[k]-previous[k]) for k in (0, 1)))
            if cv >= 0:
                clipped.append(current)
            previous, pv = current, cv
        polygon = clipped
    return abs(sum(p[0]*q[1]-p[1]*q[0]
                   for p, q in zip(polygon, polygon[1:]+polygon[:1])))/(2*w*l)


def load_inputs(folder, source):
    path = folder / f"pvgis-{source}-32.3-36.2-2005-2023.json.gz"
    receipt = json.loads((folder/f"pvgis-{source}.receipt.json").read_text())
    if sha(path) != receipt["compressed_sha256"]:
        raise ValueError("Input digest differs: " + source)
    raw = gzip.decompress(path.read_bytes())
    if hashlib.sha256(raw).hexdigest() != receipt["raw_sha256"]:
        raise ValueError("Decompressed input digest differs")
    data = json.loads(raw)
    loc = data["inputs"]["location"]
    meteo = data["inputs"]["meteo_data"]
    assert loc["latitude"] == LATITUDE and loc["longitude"] == LONGITUDE
    assert meteo["radiation_db"].lower() == "pvgis-"+source
    assert not meteo["use_horizon"]
    assert data["inputs"]["mounting_system"]["fixed"]["slope"]["value"] == 0
    frame = pd.DataFrame(data["outputs"]["hourly"])
    times = pd.DatetimeIndex(pd.to_datetime(frame["time"], format="%Y%m%d:%H%M", utc=True))
    expected = pd.date_range("2005-01-01", "2023-12-31 23:00", freq="h", tz="UTC")
    assert times.floor("h").equals(expected), "Missing, duplicated or reordered hourly labels"
    for key in ("Gb(i)", "Gd(i)", "Gr(i)", "Int"):
        value = frame[key].to_numpy(dtype=float)
        assert np.all(np.isfinite(value)) and np.all(value >= 0), key
    assert np.all(frame["Gr(i)"].to_numpy() == 0), "Horizontal reflected component must be zero"
    assert set(frame["Int"].unique()) <= {0, 1}
    counts = {str(y): int(np.count_nonzero(times.year == y)) for y in range(2005, 2024)}
    assert all(v == (8784 if calendar.isleap(int(y)) else 8760) for y, v in counts.items())
    quality = dict(source=source, records=len(frame), year_counts=counts,
                   timestamp_minutes=sorted(map(int, np.unique(times.minute))),
                   flagged_records=int(np.count_nonzero(frame["Int"])),
                   raw_sha256=receipt["raw_sha256"], compressed_sha256=sha(path),
                   provider_inputs=data["inputs"])
    return frame, times, loc["elevation"], quality


def interval_fractions(times, altitude, configs, steps, return_weights=False):
    """Conserve hourly horizontal beam under constant within-hour DNI."""
    denom = np.zeros(len(times))
    numerators = [np.zeros(len(times)) for _ in configs]
    for offset in ((np.arange(steps)+.5)*60/steps - 30):
        vec = vectors(times + pd.to_timedelta(offset, unit="min"), altitude)
        weight = np.maximum(0, vec[:, 2])
        denom += weight
        for i, cfg in enumerate(configs):
            numerators[i] += weight*fraction(vec, *cfg)
    values = [np.divide(n, denom, out=np.zeros_like(n), where=denom > 0) for n in numerators]
    return (values, denom/steps) if return_weights else values


def periods(times):
    months, years = np.asarray(times.month), np.asarray(times.year)
    result = []
    for year in range(2005, 2024):
        result.append(("annual", year, years == year))
        result.append(("summer", year, (years == year) & np.isin(months, [6, 7, 8])))
        if year >= 2006:
            result.append(("winter", year, ((years == year) & np.isin(months, [1, 2]))
                           | ((years == year-1) & (months == 12))))
    return result


def main(inputs, output):
    output.mkdir(parents=True, exist_ok=False)
    sky = {}
    sky_checks = []
    for ratio, height in itertools.product((1, 2, 4), (2.5, 5, 7.5)):
        w, l, h = dimensions(ratio, height)
        f0 = sky_exact(w, l, h)
        angular, j1 = sky_moments(w, l, h, 400)
        angular_fine, j1_fine = sky_moments(w, l, h, 800)
        kernel = sky_kernel_moment(w, l, h)
        assert abs(angular_fine-f0) < 3e-5
        assert abs(j1_fine-kernel) < 3e-5
        assert abs(j1-j1_fine) < 3e-5
        sky[ratio, height] = (f0, j1_fine)
        sky_checks.append(dict(aspect_ratio=ratio, height_m=height, exact_f0=f0,
                               angular_f0=angular_fine, j1_400=j1, j1_800=j1_fine,
                               kernel_j1=kernel, f0_error=abs(f0-angular_fine),
                               j1_error=abs(kernel-j1_fine)))
    rng = np.random.default_rng(1302026)
    polygon_errors = []
    for _ in range(1000):
        az = rng.uniform(0, 2*np.pi)
        z = rng.uniform(.001, 1.)
        vector = np.array([math.sqrt(1-z*z)*math.sin(az), math.sqrt(1-z*z)*math.cos(az), z])
        cfg = (float(rng.uniform(.5, 6)), float(rng.uniform(.1, 12)), float(rng.uniform(0, 180)))
        polygon_errors.append(abs(float(fraction(vector[None, :], *cfg)[0])-fraction_polygon(vector, *cfg)))
    assert max(polygon_errors) < 2e-13
    seasonal, sensitivity, quality, solar_checks = [], [], [], []
    suspect_days = set()
    for source in SOURCES:
        frame, times, altitude, q = load_inputs(inputs, source)
        quality.append(q)
        print(json.dumps({"phase":"loaded","source":source,"records":len(times)}), flush=True)
        vec = vectors(times, altitude)
        beam = frame["Gb(i)"].to_numpy(float)
        diffuse = frame["Gd(i)"].to_numpy(float)
        flagged = frame["Int"].to_numpy(bool)
        above_ten = vec[:, 2] > math.sin(math.radians(10))
        suspect = above_ten & ((beam+diffuse) == 0)
        q["zero_irradiance_above_10deg_records"] = int(suspect.sum())
        q["zero_irradiance_above_10deg_flagged"] = int((suspect & flagged).sum())
        q["zero_irradiance_above_10deg_dates"] = sorted(set(times[suspect].strftime("%Y-%m-%d")))
        if source == "sarah3":
            suspect_days = set(q["zero_irradiance_above_10deg_dates"])
        screened = np.asarray([day not in suspect_days for day in times.strftime("%Y-%m-%d")])
        q["matched_day_screen_excluded_dates"] = len(suspect_days)
        q["matched_day_screen_excluded_records"] = int((~screened).sum())
        masks = periods(times)
        point = [fraction(vec, *cfg) for cfg in CONFIGS]
        if source == "sarah3":
            effective, open_cos = point, np.maximum(0, vec[:, 2])
        else:
            effective, open_cos = interval_fractions(times, altitude, CONFIGS, 12, return_weights=True)
        for cfg, f in zip(CONFIGS, effective):
            ratio, height, bearing = cfg
            f0, j1 = sky[ratio, height]
            sky_horizon = (f0-.8*j1)/(1-.8*2/3)
            sky_zenith = (f0+2*j1)/(1+2*2/3)
            for period, year, mask in masks:
                b, d = float(beam[mask].sum()), float(diffuse[mask].sum())
                bt = float(np.dot(beam[mask], f[mask]))
                total = b+d
                # Baseline separates irradiance weighting at identical primary geometry/time representation.
                cosine_weight = open_cos[mask]
                baseline = float(np.dot(cosine_weight, f[mask])/cosine_weight.sum())
                good = mask & ~flagged
                good_total = float((beam[good]+diffuse[good]).sum())
                no_flags = float((np.dot(beam[good], f[good])+diffuse[good].sum()*f0)/good_total) if good_total else None
                screened_mask = mask & screened
                screened_open = float((beam[screened_mask]+diffuse[screened_mask]).sum())
                screened_transmission = float((np.dot(beam[screened_mask], f[screened_mask])+diffuse[screened_mask].sum()*f0)/screened_open)
                seasonal.append(dict(source=source, period=period, year=year, aspect_ratio=ratio,
                    height_m=height, bearing_deg=bearing, records=int(mask.sum()),
                    beam_horizontal_wh_m2=b, diffuse_horizontal_wh_m2=d,
                    transmitted_beam_wh_m2=bt, sky_isotropic=f0,
                    direct_transmission=bt/b if b else 0, diffuse_fraction=d/total,
                    total_transmission=(bt+d*f0)/total,
                    horizon_weighted_transmission=(bt+d*sky_horizon)/total,
                    zenith_weighted_transmission=(bt+d*sky_zenith)/total,
                    constant_dni_transmission=baseline,
                    unflagged_transmission=no_flags,
                    matched_day_screen_transmission=screened_transmission,
                    matched_day_screen_open_wh_m2=screened_open,
                    matched_day_screen_records=int(screened_mask.sum())))
        # Targeted full-year resolution/model sensitivity: all 27 geometries, both products.
        select = times.year == 2023
        t = times[select]
        v = vec[select]
        nvec = vectors(t, altitude, "noaa")
        coarse = interval_fractions(t, altitude, CONFIGS, 12)
        fine = interval_fractions(t, altitude, CONFIGS, 24)
        b, d = beam[select], diffuse[select]
        for cfg, f12, f24 in zip(CONFIGS, coarse, fine):
            fp, fn = fraction(v, *cfg), fraction(nvec, *cfg)
            f0 = sky[cfg[0], cfg[1]][0]
            for period, months in (("annual",tuple(range(1,13))),("summer",(6,7,8)),("winter_blocks",(12,1,2))):
                m = np.isin(t.month, months);den = float((b[m]+d[m]).sum())
                calc = lambda f: float((np.dot(b[m], f[m])+d[m].sum()*f0)/den)
                point_t, noaa_t, fine_t, coarse_t = map(calc, (fp, fn, f24, f12))
                sensitivity.append(dict(source=source,year=2023,period=period,aspect_ratio=cfg[0],
                    height_m=cfg[1],bearing_deg=cfg[2],point_spa=point_t,point_noaa=noaa_t,
                    interval_5min=coarse_t,interval_2_5min=fine_t,
                    quadrature_difference=abs(fine_t-coarse_t),
                    point_vs_interval_difference=abs(point_t-fine_t),
                    solar_algorithm_difference=abs(point_t-noaa_t)))
        dots = np.clip(np.sum(v*nvec,axis=1), -1, 1)
        angular = np.degrees(np.arccos(dots))
        daylight = v[:,2] > 0
        provider_elevation = frame.loc[select,"H_sun"].to_numpy(float)
        spa_elevation = np.degrees(np.arcsin(v[:,2]))
        solar_checks.append(dict(source=source,comparison_year=2023,
            noaa_spa_max_daylight_vector_angle_deg=float(angular[daylight].max()),
            noaa_spa_p95_daylight_vector_angle_deg=float(np.percentile(angular[daylight],95)),
            provider_spa_max_positive_elevation_difference_deg=float(np.abs(provider_elevation[daylight & (provider_elevation > 0)]-spa_elevation[daylight & (provider_elevation > 0)]).max()),
            provider_zero_elevation_daylight_records=int((daylight & (provider_elevation == 0)).sum()),
            positive_beam_below_spa_horizon_wh_m2=float(b[v[:,2]<=0].sum()),
            scope="Algorithm and timestamp sensitivity; not astronomical or irradiance validation."))
        print(json.dumps({"phase":"integrated","source":source}),flush=True)
    aggregate = []
    for source, period, cfg in itertools.product(SOURCES, ("annual","winter","summer"), CONFIGS):
        selected = [r for r in seasonal if r["source"]==source and r["period"]==period
                    and (r["aspect_ratio"],r["height_m"],r["bearing_deg"])==cfg]
        b = sum(r["beam_horizontal_wh_m2"] for r in selected)
        d = sum(r["diffuse_horizontal_wh_m2"] for r in selected)
        bt = sum(r["transmitted_beam_wh_m2"] for r in selected)
        f0, j1 = sky[cfg[0],cfg[1]]
        aggregate.append(dict(source=source,period=period,aspect_ratio=cfg[0],height_m=cfg[1],
            bearing_deg=cfg[2],years=len(selected),direct_transmission=bt/b,diffuse_fraction=d/(b+d),
            total_transmission=(bt+d*f0)/(b+d),
            minimum_year_transmission=min(r["total_transmission"] for r in selected),
            maximum_year_transmission=max(r["total_transmission"] for r in selected),
            horizon_weighted_transmission=(bt+d*(f0-.8*j1)/(1-.8*2/3))/(b+d),
            zenith_weighted_transmission=(bt+d*(f0+2*j1)/(1+2*2/3))/(b+d),
            matched_day_screen_transmission=sum(r["matched_day_screen_transmission"]*r["matched_day_screen_open_wh_m2"] for r in selected)/sum(r["matched_day_screen_open_wh_m2"] for r in selected)))
    pairs = []
    for source, period, height in itertools.product(SOURCES,("annual","winter","summer"),(2.5,5,7.5)):
        sq = {(r["year"]):r for r in seasonal if r["source"]==source and r["period"]==period
              and r["aspect_ratio"]==1 and r["height_m"]==height and r["bearing_deg"]==0}
        for bearing in (0,90):
            elong = [r for r in seasonal if r["source"]==source and r["period"]==period
                    and r["aspect_ratio"]==4 and r["height_m"]==height and r["bearing_deg"]==bearing]
            differences = [r["total_transmission"]-sq[r["year"]]["total_transmission"] for r in elong]
            pairs.append(dict(source=source,period=period,height_m=height,elongated_bearing_deg=bearing,
                years=len(differences),elongated_higher_years=sum(x>1e-8 for x in differences),
                minimum_difference=min(differences),maximum_difference=max(differences),
                mean_year_difference=sum(differences)/len(differences),
                screened_elongated_higher_years=sum(r["matched_day_screen_transmission"]>sq[r["year"]]["matched_day_screen_transmission"]+1e-8 for r in elong)))
    maximum_resolution = max(r["quadrature_difference"] for r in sensitivity)
    assert maximum_resolution < .001, maximum_resolution
    assert all(0 <= r["total_transmission"] <= 1 for r in seasonal)
    summary = dict(version=VERSION,report="ABLAQ-R05",author="ABLAQ Studios",
        coordinates=dict(latitude_n=LATITUDE,longitude_e=LONGITUDE),
        input_period=[2005,2023],source_records=sum(r["records"] for r in quality),
        geometries=len(CONFIGS),complete_winters=18,annual_years=19,seasonal_rows=len(seasonal),
        aggregate=aggregate,pairwise_comparisons=pairs)
    verification = dict(status="passed",scope="Internal numerical verification and model sensitivity; no external peer review or field validation",
        input_quality_screen="Post-inspection sensitivity: remove the same UTC dates from both products when SARAH3 reports zero direct plus diffuse irradiance above 10 degrees geometric solar elevation. Primary results retain original values; no imputation or assertion of provider error.",
        inputs=quality,software_versions=dict(numpy=np.__version__,pandas=pd.__version__,pvlib=pvlib.__version__),
        polygon_comparison=dict(cases=len(polygon_errors),maximum_absolute_fraction_difference=max(polygon_errors),tolerance=2e-13),
        sky_comparisons=sky_checks,solar_comparisons=solar_checks,
        temporal_quadrature=dict(comparison_year=2023,geometries=27,sources=2,
                                maximum_total_fraction_difference_5_vs_2_5_minutes=maximum_resolution,tolerance=.001),
        maximum_point_vs_interval_total_fraction_difference=max(r["point_vs_interval_difference"] for r in sensitivity),
        maximum_noaa_vs_spa_total_fraction_difference=max(r["solar_algorithm_difference"] for r in sensitivity))
    write_csv(output/"weather-seasonal.v1.csv",seasonal)
    write_csv(output/"weather-aggregate.v1.csv",aggregate)
    write_csv(output/"weather-sensitivity.v1.csv",sensitivity)
    write_csv(output/"weather-pairwise.v1.csv",pairs)
    write_json(output/"weather-summary.v1.json",summary)
    write_json(output/"weather-verification.v1.json",verification)
    files = {p.name:dict(bytes=p.stat().st_size,sha256=sha(p)) for p in sorted(output.iterdir())}
    files[Path(__file__).name]=dict(bytes=Path(__file__).stat().st_size,sha256=sha(Path(__file__)))
    write_json(output/"weather-computation-manifest.v1.json",dict(version=VERSION,files=files))
    print(json.dumps(dict(status="passed",seasonal_rows=len(seasonal),aggregate_rows=len(aggregate),
                          maximum_resolution_difference=maximum_resolution)),flush=True)


if __name__ == "__main__":
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--inputs",type=Path,required=True)
    parser.add_argument("--output",type=Path,required=True)
    args=parser.parse_args()
    main(args.inputs,args.output)