###################### GeometryBuilder #########################
from shapely.geometry import Polygon, MultiPolygon
from pyproj import Transformer
import numpy as np
import math

NM_TO_M = 1852

def get_aeqd(lat, lon):
    proj = f"+proj=aeqd +lat_0={lat} +lon_0={lon}"
    fwd = Transformer.from_crs("EPSG:4326", proj, always_xy=True).transform
    inv = Transformer.from_crs(proj, "EPSG:4326", always_xy=True).transform
    return fwd, inv


def arc_points(center, start, end, radius_nm, clockwise):
    fwd, inv = get_aeqd(center.lat, center.lon)

    cx, cy = fwd(center.lon, center.lat)
    sx, sy = fwd(start.lon, start.lat)
    ex, ey = fwd(end.lon, end.lat)

    a1 = math.atan2(sy-cy, sx-cx)
    a2 = math.atan2(ey-cy, ex-cx)

    if clockwise:
        if a2 > a1: a2 -= 2*math.pi
    else:
        if a2 < a1: a2 += 2*math.pi

    r = radius_nm * NM_TO_M
    angles = np.linspace(a1, a2, 64)

    pts = []
    for a in angles:
        x = cx + r*math.cos(a)
        y = cy + r*math.sin(a)
        pts.append(inv(x, y))
    return pts


def circle_points(center, radius_nm):
    fwd, inv = get_aeqd(center.lat, center.lon)
    cx, cy = fwd(center.lon, center.lat)

    r = radius_nm * NM_TO_M
    angles = np.linspace(0, 2*math.pi, 128)

    return [inv(cx + r*math.cos(a), cy + r*math.sin(a)) for a in angles]


def build_geometry(ast: NotamAST):
    polys = []

    for area in ast.areas:
        pts = []

        for seg in area.segments:
            if isinstance(seg, LineSeg):
                pts.append((seg.start.lon, seg.start.lat))
                pts.append((seg.end.lon, seg.end.lat))

            elif isinstance(seg, ArcSeg):
                pts.extend(arc_points(
                    seg.center, seg.start, seg.end,
                    seg.radius_nm, seg.clockwise
                ))

            elif isinstance(seg, CircleSeg):
                pts.extend(circle_points(seg.center, seg.radius_nm))

        if pts and pts[0] != pts[-1]:
            pts.append(pts[0])

        if pts:
            polys.append(Polygon(pts))

    return polys[0] if len(polys) == 1 else MultiPolygon(polys)