#!/usr/bin/env python3
"""Compare Azimuthal Map's R=6371 km spherical AEQD distances with WGS84.

Run: python -m pip install geographiclib==2.1
     python static/research/compare-sphere-wgs84.py --out static/research

Generated files have public, deterministic coordinates only. No site analytics are used.
"""
from __future__ import annotations

import argparse
import csv
import json
import math
import re
from pathlib import Path
from statistics import median

from geographiclib.geodesic import Geodesic

# Matches static/js/map-generator.js and static/js/map-generator-equal-area.js
EARTH_RADIUS_KM = 6371.0
LATITUDES = (-80, -60, -30, 0, 30, 60, 80)
LONGITUDES = (-120, 0, 120)
AZIMUTHS = (0, 45, 90, 135, 180, 225, 270, 315)
LENGTHS_KM = (1, 10, 100, 500, 1000, 2500, 5000, 10000, 15000, 19000)
FIELDS = ('group', 'center_lat', 'center_lon', 'target_lat', 'target_lon',
          'planned_wgs84_km', 'initial_azimuth', 'sphere_km', 'wgs84_km',
          'signed_error_km', 'absolute_error_km', 'percentage_error')


def spherical_distance(lat1: float, lon1: float, lat2: float, lon2: float) -> float:
    """Central angle of two points on the sphere * 6371 km.

    Mathematically equivalent to the D3 geoAzimuthalEquidistant radial
    coordinate (angular arc), with atan2 for floating-point stability.
    The map scales that angular arc using EARTH_RADIUS_KM. This does not
    claim to reproduce D3's rasterized screen or clipped antipodal output.
    """
    a, b = math.radians(lat1), math.radians(lat2)
    d = math.radians((lon2 - lon1 + 180) % 360 - 180)
    y = math.cos(b) * math.sin(d)
    x = math.cos(a) * math.sin(b) - math.sin(a) * math.cos(b) * math.cos(d)
    dot = math.sin(a) * math.sin(b) + math.cos(a) * math.cos(b) * math.cos(d)
    return EARTH_RADIUS_KM * math.atan2(math.hypot(x, y), dot)


def make_row(group, lat1, lon1, lat2, lon2, planned=None, azimuth=None):
    sphere = spherical_distance(lat1, lon1, lat2, lon2)
    reference = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2)['s12'] / 1000
    signed = sphere - reference
    return dict(zip(FIELDS, (
        group, lat1, lon1, lat2, lon2, planned, azimuth,
        sphere, reference, signed, abs(signed),
        (100 * signed / reference if reference > 0 else None))))


def cases():
    for lat in LATITUDES:
        for lon in LONGITUDES:
            for az in AZIMUTHS:
                for dist in LENGTHS_KM:
                    endpoint = Geodesic.WGS84.Direct(lat, lon, az, dist * 1000)
                    yield make_row('stratified', lat, lon, endpoint['lat2'],
                                   endpoint['lon2'], dist, az)
            yield make_row('coincident', lat, lon, lat, lon, 0, None)
            for delta in (0, 0.000001, 0.001, 0.1):
                antipode_lon = (lon + 180 + 180) % 360 - 180
                # Exact antipode, plus three specified signed perturbations.
                yield make_row('antipodal' if delta == 0 else 'near_antipodal',
                               lat, lon, -lat + delta, antipode_lon - delta, None, None)


def pct(values, quantile):
    if not values:
        return None
    vals = sorted(values)
    pos = (len(vals) - 1) * quantile
    lo = int(math.floor(pos))
    hi = int(math.ceil(pos))
    return vals[lo] + (vals[hi] - vals[lo]) * (pos - lo)


def describe(rows):
    errors = [r['signed_error_km'] for r in rows]
    absolute = list(map(abs, errors))
    percentages = [r['percentage_error'] for r in rows if r['percentage_error'] is not None]
    worst = max(rows, key=lambda r: r['absolute_error_km'])
    return {'n': len(rows), 'signed_error_min_km': min(errors),
            'signed_error_max_km': max(errors), 'signed_error_median_km': median(errors),
            'absolute_error_median_km': median(absolute),
            'absolute_error_p95_km': pct(absolute, .95),
            'absolute_error_max_km': max(absolute),
            'percentage_error_min': min(percentages) if percentages else None,
            'percentage_error_max': max(percentages) if percentages else None,
            'absolute_percentage_error_median': median(map(abs,percentages)) if percentages else None,
            'absolute_percentage_error_p95': pct(list(map(abs, percentages)), .95) if percentages else None,
            'worst_route': {k: worst[k] for k in FIELDS}}


def write_csv(path, rows):
    with path.open('w', newline='', encoding='utf-8') as handle:
        writer = csv.DictWriter(handle, FIELDS)
        writer.writeheader()
        for r in rows:
            writer.writerow({k: ('%.9f' % v if isinstance(v, float) else v) if v is not None else ''
                             for k, v in r.items()})


def validate_source(repo_root):
    source = (repo_root / 'static/js/map-generator.js')
    if not source.exists():
        return 'Source not present; compare this script with static/js/map-generator.js separately.'
    content = source.read_text(encoding='utf-8')
    m = re.search(r'const EARTH_RADIUS_KM\s*=\s*([\d.]+)', content)
    if not m or float(m.group(1)) != EARTH_RADIUS_KM:
        raise RuntimeError('Generator Earth radius changed; rerun review before using results')
    if 'd3.geoAzimuthalEquidistant()' not in content or 'maxDistance / EARTH_RADIUS_KM' not in content:
        raise RuntimeError('Generator distance or projection implementation changed')
    minified = repo_root / 'static/js/map-generator.min.js'
    if minified.exists():
        compact = minified.read_text(encoding='utf-8')
        if not re.search(r'EARTH_RADIUS_KM\s*=\s*6371\b', compact):
            raise RuntimeError('Minified generator Earth radius differs from source')
        if 'd3.geoAzimuthalEquidistant()' not in compact:
            raise RuntimeError('Minified generator projection type differs from source')
        return 'Matched source and minified generator radius, D3 projection type, and scaling.'
    return 'Matched source generator radius, D3 projection type, and scaling; minified build absent.'


def write_figures(rows, out):
    """Create accessible, original SVG figures without extra plot dependencies."""
    base = [r for r in rows if r['group'] == 'stratified']
    W, H, L, T, RW, BH = 960, 530, 94, 96, 820, 338

    def svg_start(title, desc):
        return [f'<svg xmlns="http://www.w3.org/2000/svg" width="{W}" height="{H}" viewBox="0 0 {W} {H}" role="img" aria-labelledby="title desc">',
                f'<title id="title">{title}</title><desc id="desc">{desc}</desc>',
                '<rect width="960" height="530" fill="#ffffff"/>',
                f'<text x="{L}" y="42" font-family="Arial,sans-serif" font-size="25" fill="#193c36" font-weight="700">{title}</text>',
                '<text x="94" y="66" font-family="Arial,sans-serif" font-size="15" fill="#4b5b57">Spherical distance minus WGS84 distance. Positive values are longer.</text>']

    def line(x1,y1,x2,y2, color, width=1, extra=''):
        return f'<line x1="{x1:.2f}" y1="{y1:.2f}" x2="{x2:.2f}" y2="{y2:.2f}" stroke="{color}" stroke-width="{width}" {extra}/>'

    def text_label(x,y,value,anchor='middle',size=14):
        return f'<text x="{x:.2f}" y="{y:.2f}" fill="#34453f" text-anchor="{anchor}" font-family="Arial,sans-serif" font-size="{size}">{value}</text>'

    # Figure 1: logarithmic route-length axis and signed-km error.
    by_length = {k:[r['signed_error_km'] for r in base if r['planned_wgs84_km']==k]
                 for k in LENGTHS_KM}
    xv = lambda dist: L + RW*math.log10(dist)/math.log10(max(LENGTHS_KM))
    yv = lambda err: T + BH*(40-err)/80
    outsvg = svg_start('Distance Error by WGS84 Route Length',
        'Ranges of signed error over the tested centers and bearings, with median by distance. Horizontal axis is logarithmic.')
    for e in (-40,-20,0,20,40):
        y=yv(e)
        outsvg.append(line(L,y,L+RW,y,'#dbe3df',1 if e else 2))
        outsvg.append(text_label(L-14,y+5,str(e),anchor='end'))
    for dist in (1,10,100,1000,5000,10000,19000):
        x=xv(dist)
        outsvg.append(line(x,T,x,T+BH,'#e4eae6'))
        outsvg.append(text_label(x,T+BH+28,f'{dist:,}'))
    upper=' '.join(f'{xv(k):.1f},{yv(max(v)):.1f}' for k,v in by_length.items())
    lower=' '.join(f'{xv(k):.1f},{yv(min(v)):.1f}' for k,v in reversed(tuple(by_length.items())))
    outsvg.append(f'<polygon points="{upper} {lower}" fill="#bfd9d0" fill-opacity="0.55"/>')
    medpoints=' '.join(f'{xv(k):.1f},{yv(median(v)):.1f}' for k,v in by_length.items())
    outsvg.append(f'<polyline points="{medpoints}" fill="none" stroke="#14654e" stroke-width="3"/>')
    for k,v in by_length.items():
        outsvg.append(f'<circle cx="{xv(k):.2f}" cy="{yv(median(v)):.2f}" r="4.5" fill="#14654e"/>')
    outsvg.extend([text_label(L+RW/2,H-52,'WGS84 route length (km, logarithmic axis)'),
                   text_label(L-62,T+BH/2,'Error (km)',size=14).replace('text-anchor="middle"','transform="rotate(-90,32,'+str(int(T+BH/2))+')" text-anchor="middle"'),
                   '<rect x="650" y="80" width="22" height="12" fill="#bfd9d0"/>',
                   text_label(680,91,'Tested min–max range',anchor='start',size=13),
                   line(650,108,673,108,'#14654e',3),
                   text_label(680,113,'Median',anchor='start',size=13),
                   text_label(L,H-13,'560 distinct geometries at three longitudes. Not a global error bound.',anchor='start',size=13),
                   '</svg>'])
    (out/'spherical-error-by-distance.svg').write_text('\n'.join(outsvg)+'\n',encoding='utf-8')

    # Figure 2: center latitude on a linear axis and percentage error.
    by_lat = {lat:[r['percentage_error'] for r in base if r['center_lat']==lat] for lat in LATITUDES}
    xv2=lambda lat: L+RW*(lat+80)/160
    yv2=lambda err: T+BH*(.7-err)/1.3
    outsvg=svg_start('Percentage Distance Error by Center Latitude',
        'Minimum, median and maximum signed percentage errors at each tested center latitude.')
    for e in (-.6,-.3,0,.3,.6):
        y=yv2(e)
        outsvg.append(line(L,y,L+RW,y,'#dbe3df',1 if e else 2))
        outsvg.append(text_label(L-14,y+5,f'{e:+.1f}%',anchor='end'))
    for lat, v in by_lat.items():
        x=xv2(lat)
        outsvg.append(line(x,T,x,T+BH,'#e4eae6'))
        outsvg.append(text_label(x,T+BH+28,str(lat)))
        outsvg.append(line(x,yv2(min(v)),x,yv2(max(v)),'#528979',4))
        for bound in (min(v),max(v)):
            outsvg.append(line(x-9,yv2(bound),x+9,yv2(bound),'#528979',3))
        outsvg.append(f'<circle cx="{x:.2f}" cy="{yv2(median(v)):.2f}" r="6" fill="#14654e" stroke="white" stroke-width="1.3"/>')
    outsvg.extend([text_label(L+RW/2,H-52,'Center latitude (degrees)'),
                   text_label(L,H-13,'80 distinct geometries per latitude, checked at three longitudes.',anchor='start',size=13),
                   '<circle cx="674" cy="87" r="6" fill="#14654e"/>',
                   text_label(692,92,'Median',anchor='start',size=13),
                   line(674,103,674,123,'#528979',4),
                   text_label(692,116,'Tested min–max range',anchor='start',size=13),
                   '</svg>'])
    (out/'spherical-error-by-center-latitude.svg').write_text('\n'.join(outsvg)+'\n',encoding='utf-8')

def crosscheck_pyproj(rows):
    """Optional reproducible cross-check of 100 fixed principal rows using PROJ."""
    try:
        from pyproj import Geod
    except ImportError as exc:
        raise RuntimeError('Optional check requires: python -m pip install pyproj==3.7.2') from exc
    base = [r for r in rows if r['group'] == 'stratified']
    indices = [round(i * (len(base) - 1) / 99) for i in range(100)]
    geod = Geod(ellps='WGS84')
    diffs = []
    for i in indices:
        r = base[i]
        _, _, meters = geod.inv(r['center_lon'], r['center_lat'],
                                r['target_lon'], r['target_lat'])
        diffs.append(abs(meters / 1000 - r['wgs84_km']))
    return {'comparator': 'pyproj.Geod WGS84', 'sample_n': len(indices),
            'max_absolute_difference_km': max(diffs)}


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--out', type=Path, default=Path('static/research'))
    parser.add_argument('--repo-root', type=Path, default=Path(__file__).resolve().parents[2])
    parser.add_argument('--no-figures', action='store_true')
    parser.add_argument('--verify-pyproj', action='store_true',
                        help='Optionally compare 100 deterministic pairs with pyproj/PROJ')
    args = parser.parse_args()
    args.out.mkdir(parents=True, exist_ok=True)
    validation = validate_source(args.repo_root)
    rows = list(cases())
    assert len(rows) == 21 * (8 * 10 + 1 + 4)
    base = [r for r in rows if r['group'] == 'stratified']
    for row in base:
        assert abs(row['wgs84_km'] - row['planned_wgs84_km']) < 1e-6
    for row in rows:
        assert math.isfinite(row['sphere_km'])
        assert math.isfinite(row['wgs84_km'])
    groups = {g: describe([r for r in rows if r['group'] == g])
              for g in ('stratified', 'coincident', 'antipodal', 'near_antipodal')}
    groups['by_length'] = {str(km): describe([r for r in base if r['planned_wgs84_km'] == km])
                           for km in LENGTHS_KM}
    groups['by_latitude'] = {str(lat): describe([r for r in base if r['center_lat'] == lat])
                             for lat in LATITUDES}
    groups['source_validation'] = validation
    groups['sampling'] = {'center_latitudes':LATITUDES, 'center_longitudes':LONGITUDES,
                          'initial_azimuths':AZIMUTHS, 'test_lengths_km':LENGTHS_KM,
                          'radius_km':EARTH_RADIUS_KM,
                          'wgs84_a_m':Geodesic.WGS84.a,
                          'wgs84_f':Geodesic.WGS84.f}
    write_csv(args.out / 'spherical-vs-wgs84-routes.csv', rows)
    (args.out / 'spherical-vs-wgs84-summary.json').write_text(
        json.dumps(groups, indent=2) + '\n', encoding='utf-8')
    if not args.no_figures:
        write_figures(rows, args.out)
    report = {'source_validation': validation, 'stratified': groups['stratified'],
              'antipodal': groups['antipodal'], 'near_antipodal': groups['near_antipodal']}
    if args.verify_pyproj:
        report['pyproj_crosscheck'] = crosscheck_pyproj(rows)
    print(json.dumps(report, indent=2))


if __name__ == '__main__':
    main()