#!/usr/bin/env python3
"""Spherical Mercator local scale, plus a reproducible schematic SVG.

Run: python3 mercator-area-lab.py [--svg mercator-area-diagram.svg]
The ratios apply to infinitesimal patches on an ideal sphere, not measured
country areas or exact WGS 84 ellipsoid ground areas.
"""

import argparse
import math
import platform
from pathlib import Path
from xml.etree import ElementTree as ET

RADIUS_M = 6_378_137.0
WEB_TILE_EDGE_DEG = math.degrees(math.atan(math.sinh(math.pi)))


def latitude_radians(degrees):
    if not math.isfinite(degrees) or not -90 < degrees < 90:
        raise ValueError("latitude must be finite and strictly between -90 and 90")
    return math.radians(degrees)


def mercator_y(degrees):
    phi = latitude_radians(degrees)
    return RADIUS_M * math.asinh(math.tan(phi))


def linear_scale(degrees):
    phi = latitude_radians(degrees)
    return 1.0 / math.cos(phi)


def local_area_scale(degrees):
    scale = linear_scale(degrees)
    return scale * scale


def check_boundaries():
    assert math.isclose(linear_scale(0), 1.0)
    assert math.isclose(local_area_scale(30), 4.0 / 3.0)
    assert math.isclose(linear_scale(60), 2.0)
    assert math.isclose(local_area_scale(60), 4.0)
    assert math.isclose(local_area_scale(-60), local_area_scale(60))
    assert math.isclose(mercator_y(-60), -mercator_y(60), abs_tol=1e-8)
    assert math.isclose(mercator_y(WEB_TILE_EDGE_DEG), math.pi * RADIUS_M,
                        abs_tol=1e-7)
    for degrees in (-90, 90, float("nan"), float("inf")):
        try:
            local_area_scale(degrees)
        except ValueError:
            pass
        else:
            raise AssertionError(f"invalid latitude accepted: {degrees}")


def make_svg():
    # Equal infinitesimal ground-area patches are drawn at one common map scale.
    # Their schematic side lengths grow by the spherical local linear factor.
    cells = [(0, 105, "#a8d08d"), (60, 305, "#f0c4a8"),
             (80, 545, "#94b8d8")]
    parts = [
        '<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 840 355" '
        'role="img" aria-labelledby="title desc">',
        '<title id="title">Spherical Mercator local area comparison</title>',
        '<desc id="desc">Equal tiny ground patches at 0, 60 and 80 degrees '
        'latitude occupy schematic map areas in the ratio 1 to 4 to 33.2.</desc>',
        '<rect width="840" height="355" rx="18" fill="#fdfcfa"/>',
        '<text x="32" y="38" font-family="system-ui,sans-serif" font-size="20" '
        'fill="#3d3a36">Equal tiny ground patches · one map scale</text>',
        '<line x1="45" y1="300" x2="800" y2="300" stroke="#7a756e" '
        'stroke-width="1.5"/>',
    ]
    for latitude, left, fill in cells:
        side = 40.0 * linear_scale(latitude)
        top = 300.0 - side
        ratio = local_area_scale(latitude)
        parts += [
            f'<rect x="{left:.1f}" y="{top:.1f}" width="{side:.1f}" '
            f'height="{side:.1f}" rx="4" fill="{fill}" '
            'stroke="#3d3a36" stroke-width="2"/>',
            f'<text x="{left + side / 2:.1f}" y="{top - 13:.1f}" '
            'text-anchor="middle" font-family="system-ui,sans-serif" '
            f'font-size="20" fill="#3d3a36">×{ratio:.1f}</text>',
            f'<text x="{left + side / 2:.1f}" y="331" text-anchor="middle" '
            'font-family="system-ui,sans-serif" font-size="18" '
            f'fill="#3d3a36">{latitude}°</text>',
        ]
    parts.append('</svg>')
    return '\n'.join(parts) + '\n'


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--svg", type=Path,
                        help="Write the schematic as an SVG file")
    options = parser.parse_args()
    check_boundaries()
    print(f"Python {platform.python_version()} on {platform.system()}")
    print("latitude  linear_scale  local_area_scale  projected_y_m")
    for latitude in (0, 30, 60, 80, 85):
        print(f"{latitude:>3}°      {linear_scale(latitude):>8.4f}"
              f"       {local_area_scale(latitude):>8.4f}"
              f"      {mercator_y(latitude):>12.2f}")
    print(f"square-tile edge: {WEB_TILE_EDGE_DEG:.8f}°"
          f" -> y={mercator_y(WEB_TILE_EDGE_DEG):.2f} m = πR")
    print("Boundary checks passed: equator, ±60°, 30°, tile edge, poles and non-finite input.")
    if options.svg:
        data = make_svg()
        ET.fromstring(data)  # The diagram must be well-formed XML.
        options.svg.write_text(data, encoding="utf-8")
        print(f"SVG written: {options.svg}")


if __name__ == "__main__":
    main()