# SPDX-License-Identifier: MIT

# geography — reason over regions, maps and places, and calculate distances,
# directions, projections and areas on the Earth.
#
# Two halves that share nothing but the subject.
#
# STRUCTURE. You supply three kinds of record as facts.
#
# Regions (`geography.Region`), each with the region it lies in directly and
# a level, such as "country", "county" or "municipality". The rules conclude
# every region a region lies within, to any depth, and how many regions each
# one holds:
#
#     import geography
#     geography.region("norway", "", "country")
#     geography.region("vestland", "norway", "county")
#     geography.region("bergen", "vestland", "municipality")
#     print(geography.within("bergen"), geography.regions_in("norway", "municipality"))
#
# A map: borders between regions (`geography.Border`) and the colour the map
# gives each region. The rules conclude each pair of neighbours, and each pair
# of neighbours given one colour, a clash. `dsatur()` chooses colours so that
# no neighbours share one, and `paint()` writes them:
#
#     geography.border("vestland", "rogaland")
#     geography.colour("vestland", "green")
#     geography.colour("rogaland", "green")
#     print(geography.clashes())
#     geography.paint(["green", "yellow", "red", "blue"])
#     print(geography.clashes())
#
# Places by latitude and longitude (`geography.Place`) and boxes of latitude
# and longitude (`geography.Box`). The rules conclude which places lie in
# which boxes, a box across the antimeridian included:
#
#     geography.place("suva", -18.14, 178.44)
#     geography.box("fiji", -21.0, -12.0, 170.0, -173.0)
#     print(geography.places_in("fiji"))
#
# Each conclusion answers `why`, and withdraws itself when a fact it rests on
# is corrected: a region moved to another parent leaves every region above
# its old one. Where two sources disagree about a fact, both values are
# reasoned about, no query refuses, and `unsettled()` names what they
# disagree about. A region holds each region within it once, however many
# paths lead there; a region whose parent two sources disagree on lies within
# each region a value names, and each of those holds it.
#
# CALCULATION. Value scripts that return a value and write nothing: the
# great-circle distance, bearing and destination on a sphere of the Earth's
# mean radius; the distance and azimuths between two points, and the point a
# distance and azimuth away, on the WGS84 ellipsoid by Vincenty's method;
# the Mercator projection and its scale; and the area of a triangle on the
# sphere by L'Huilier's theorem. Angles are in degrees and distances in
# metres.

import math

# ── Regions as facts ─────────────────────────────────────────────────────────

class Region:
    # A region, by a name you choose: `parent`, the name of the region it lies
    # in directly, unwritten for one that lies in none you record; `level`,
    # what kind of region it is, in any word you use; and `colour`, the
    # colour a map gives it. The rules conclude `holds`, how many regions lie
    # within it at any depth, each counted once however many paths lead to it.
    parent: string
    level: string
    colour: string
    holds: int

class Within:
    # The relation `Within.region.outer`, True for a region's parent, its
    # parent's parent, and so on to the outermost region recorded.
    pass

class Border:
    # A border between the regions named `a` and `b`, named `<a>_<b>`.
    a: string
    b: string

class Neighbour:
    # The relation `Neighbour.x.y`, True both ways for each border.
    pass

class Clash:
    # The relation `Clash.x.y`, True both ways for two neighbours the map
    # gives one colour. It is declared so that a map with no clash reads as
    # empty rather than as an unknown kind.
    pass

# ── Places as facts ──────────────────────────────────────────────────────────

class Place:
    # A place, by a name you choose, at latitude `lat` and longitude `lon`,
    # in degrees, north and east positive.
    lat: real
    lon: real

class Box:
    # A box of latitude and longitude, by a name you choose: from `south` to
    # `north` and from `west` to `east`, in degrees. A box whose western edge
    # lies east of its eastern one crosses the antimeridian, 180°, and holds
    # the longitudes from `west` round to `east`.
    south: real
    north: real
    west: real
    east: real

class Inside:
    # The relation `Inside.place.box`, True for each place within a box, its
    # edges included.
    pass

# ── Standing rules: regions ──────────────────────────────────────────────────

@rule
def within_parent():
    for r in Region, p in Region:
        if r.parent == p:
            Within.r.p = True

@rule
def within_further():
    for r in Region, p in Region, q in Region:
        if Within.r.p == True and p.parent == q:
            Within.r.q = True

@rule
def holds():
    for o in Region, r in Region:
        if Within.r.o == True:
            o.holds = distinct(r)

# ── Standing rules: a map ────────────────────────────────────────────────────

@rule
def neighbour_ab():
    for b in Border, x in Region, y in Region:
        if b.a == x and b.b == y:
            Neighbour.x.y = True

@rule
def neighbour_ba():
    for b in Border, x in Region, y in Region:
        if b.a == x and b.b == y:
            Neighbour.y.x = True

@rule
def clash():
    for x in Region, y in Region:
        if Neighbour.x.y == True and x.colour == y.colour:
            Clash.x.y = True

# ── Standing rules: places in boxes ──────────────────────────────────────────

@rule
def inside():
    for p in Place, b in Box:
        if b.west <= b.east and p.lat >= b.south and p.lat <= b.north and p.lon >= b.west and p.lon <= b.east:
            Inside.p.b = True

@rule
def inside_before_antimeridian():
    # A box across 180°: the longitudes from its western edge to 180°.
    for p in Place, b in Box:
        if b.west > b.east and p.lat >= b.south and p.lat <= b.north and p.lon >= b.west:
            Inside.p.b = True

@rule
def inside_after_antimeridian():
    # A box across 180°: the longitudes from −180° to its eastern edge.
    for p in Place, b in Box:
        if b.west > b.east and p.lat >= b.south and p.lat <= b.north and p.lon <= b.east:
            Inside.p.b = True

# ── Writing facts ────────────────────────────────────────────────────────────

def region(name: string, parent: string, level: string):
    # A region; pass "" for a parent not recorded.
    update(Region(name).level, level, "regions")
    if parent != "":
        update(Region(name).parent, parent, "regions")

def border(a: string, b: string):
    # A border between two regions.
    s = a + "_" + b
    update(Border(s).a, a, "map")
    update(Border(s).b, b, "map")

def colour(name: string, c: string):
    # The colour a map gives a region.
    update(Region(name).colour, c, "map")

def place(name: string, lat: real, lon: real):
    update(Place(name).lat, lat, "places")
    update(Place(name).lon, lon, "places")

def box(name: string, south: real, north: real, west: real, east: real):
    update(Box(name).south, south, "places")
    update(Box(name).north, north, "places")
    update(Box(name).west, west, "places")
    update(Box(name).east, east, "places")

# ── Reading what the rules concluded ─────────────────────────────────────────
#
# Each reader reads every value in force, so a fact two sources disagree
# about is reported with both values rather than refused.

def regions() -> list<string>:
    # Every region recorded, sorted.
    out = []
    vs = list_values(Region.level)
    for v in vs:
        parts = v.path.split(".")
        r = parts[len(parts) - 2]
        if r not in out:
            out = out + [r]
    return sorted(out)

def within(name: string) -> list<string>:
    # Every region a region lies within, sorted.
    out = []
    for r in regions():
        w = current(Within.name.r)
        match w:
            case empty:
                pass
            case default:
                out = out + [r]
    return sorted(out)

def regions_in(name: string, at: string) -> list<string>:
    # Every region within a region at the level `at`, or at any level for "",
    # sorted.
    levels: dict<string, string> = {}
    vs = list_values(Region.level)
    for v in vs:
        parts = v.path.split(".")
        levels[parts[len(parts) - 2]] = str(v.value)
    out = []
    vs = list_values(Within.name)
    for v in vs:
        parts = v.path.split(".")
        r = parts[len(parts) - 2]
        if (at == "" or levels[r] == at) and r not in out:
            out = out + [r]
    return sorted(out)

def neighbours(name: string) -> list<string>:
    # Every region that shares a border with a region, sorted.
    out = []
    for r in regions():
        w = current(Neighbour.name.r)
        match w:
            case empty:
                pass
            case default:
                out = out + [r]
    return sorted(out)

def clashes() -> list<string>:
    # Every pair of neighbours the map gives one colour, as "x y" with x
    # before y, sorted.
    out = []
    for y in regions():
        vs = list_values(Clash.y)
        for v in vs:
            parts = v.path.split(".")
            x = parts[len(parts) - 2]
            s = x + " " + y
            if y < x:
                s = y + " " + x
            if s not in out:
                out = out + [s]
    return sorted(out)

def places_in(name: string) -> list<string>:
    # Every place within a box, sorted.
    out = []
    vs = list_values(Inside.name)
    for v in vs:
        parts = v.path.split(".")
        p = parts[len(parts) - 2]
        if p not in out:
            out = out + [p]
    return sorted(out)

def unsettled() -> list<string>:
    # Each region, border, place and box that holds two different values for
    # a fact you write, as "region x parent", "place x lat" and so on,
    # sorted.
    out = []
    for f in ["parent", "level", "colour", "a", "b", "lat", "lon", "south", "north", "west", "east"]:
        if f == "parent":
            vs = list_values(Region.parent)
            label = "region "
        elif f == "level":
            vs = list_values(Region.level)
            label = "region "
        elif f == "colour":
            vs = list_values(Region.colour)
            label = "region "
        elif f == "a":
            vs = list_values(Border.a)
            label = "border "
        elif f == "b":
            vs = list_values(Border.b)
            label = "border "
        elif f == "lat":
            vs = list_values(Place.lat)
            label = "place "
        elif f == "lon":
            vs = list_values(Place.lon)
            label = "place "
        elif f == "south":
            vs = list_values(Box.south)
            label = "box "
        elif f == "north":
            vs = list_values(Box.north)
            label = "box "
        elif f == "west":
            vs = list_values(Box.west)
            label = "box "
        else:
            vs = list_values(Box.east)
            label = "box "
        seen = {}
        for v in vs:
            parts = v.path.split(".")
            s = label + parts[len(parts) - 2] + " " + f
            x = str(v.value)
            if s in seen and seen[s] != x and s not in out:
                out = out + [s]
            seen[s] = x
    return sorted(out)

# ── Colouring a map ──────────────────────────────────────────────────────────

def dsatur() -> dict<string, int>:
    # A colour for every region named in a border, numbered from 0, by
    # Brélaz's DSatur: each step colours the region whose neighbours already
    # show the most colours, the most neighbours breaking a tie and then the
    # name, with the lowest colour none of them shows. It reads the borders
    # and writes nothing.
    adj: dict<string, list<string>> = {}
    vs = list_values(Border.a)
    ws = list_values(Border.b)
    ends: dict<string, list<string>> = {}
    for v in vs:
        parts = v.path.split(".")
        ends[parts[len(parts) - 2]] = [str(v.value)]
    for w in ws:
        parts = w.path.split(".")
        b = parts[len(parts) - 2]
        if b in ends:
            ends[b] = ends[b] + [str(w.value)]
    for b in ends:
        e = ends[b]
        if len(e) == 2 and e[0] != e[1]:
            for k in range(2):
                x = e[k]
                y = e[1 - k]
                if x not in adj:
                    adj[x] = []
                if y not in adj[x]:
                    adj[x] = adj[x] + [y]
    out: dict<string, int> = {}
    names = sorted([x for x in adj])
    for step in range(len(names)):
        best = ""
        bs = -1
        bd = -1
        for x in names:
            if x not in out:
                seen = []
                for y in adj[x]:
                    if y in out and out[y] not in seen:
                        seen = seen + [out[y]]
                if len(seen) > bs or (len(seen) == bs and len(adj[x]) > bd):
                    best = x
                    bs = len(seen)
                    bd = len(adj[x])
        used = [out[y] for y in adj[best] if y in out]
        c = 0
        for k in range(len(names) + 1):
            if c in used:
                c = c + 1
        out[best] = c
    return out

def paint(palette: list<string>):
    # Write the colours `dsatur()` chooses, colour k as palette[k]. The
    # palette needs as many colours as the map does.
    chosen = dsatur()
    for x in chosen:
        update(Region(x).colour, palette[chosen[x]], "map")

# ── Calculation: the sphere ──────────────────────────────────────────────────

def mean_radius() -> real:
    # The Earth's mean radius in metres, (2a + b) / 3 of the WGS84
    # ellipsoid, rounded as the IUGG gives it.
    return 6371008.8

def central_angle(lat1: real, lon1: real, lat2: real, lon2: real) -> real:
    # The angle between two points seen from the sphere's centre, in radians,
    # by the haversine formula, which keeps its precision for points close
    # together. Each coordinate is taken modulo 360 degrees first, which
    # math.fmod does exactly, so one past a turn loses no digits to radians;
    # where the haversine falls below the smallest normal real, its root is
    # the hypot of the two terms' roots, which squares nothing.
    p1 = math.radians(math.fmod(lat1, 360.0))
    p2 = math.radians(math.fmod(lat2, 360.0))
    dp = p2 - p1
    dl = math.radians(math.fmod(lon2, 360.0) - math.fmod(lon1, 360.0))
    cc = math.cos(p1) * math.cos(p2)
    h = math.sin(dp / 2) * math.sin(dp / 2) + cc * math.sin(dl / 2) * math.sin(dl / 2)
    if h < 2.2250738585072014e-308 and cc >= 0.0:
        return 2 * math.asin(math.hypot(math.sin(dp / 2), math.sqrt(cc) * math.sin(dl / 2)))
    return 2 * math.asin(math.sqrt(h))

def haversine(lat1: real, lon1: real, lat2: real, lon2: real) -> real:
    # The great-circle distance between two points on the mean sphere.
    return mean_radius() * central_angle(lat1, lon1, lat2, lon2)

def bearing(lat1: real, lon1: real, lat2: real, lon2: real) -> real:
    # The initial bearing of the great circle from the first point to the
    # second, clockwise from north, from 0 to 360. Each coordinate is taken
    # modulo 360 degrees first, as in central_angle.
    p1 = math.radians(math.fmod(lat1, 360.0))
    p2 = math.radians(math.fmod(lat2, 360.0))
    dl = math.radians(math.fmod(lon2, 360.0) - math.fmod(lon1, 360.0))
    y = math.sin(dl) * math.cos(p2)
    x = math.cos(p1) * math.sin(p2) - math.sin(p1) * math.cos(p2) * math.cos(dl)
    t = math.degrees(math.atan2(y, x))
    if t < 0:
        t = t + 360
    return t

def destination(lat: real, lon: real, heading: real, distance: real) -> list<real>:
    # The point reached by setting out from a point on a bearing and keeping
    # to the great circle for a distance: [lat, lon], lon from −180 to 180.
    # Each angle is taken modulo 360 degrees first, as in central_angle.
    d = distance / mean_radius()
    p1 = math.radians(math.fmod(lat, 360.0))
    t = math.radians(math.fmod(heading, 360.0))
    p2 = math.asin(math.sin(p1) * math.cos(d) + math.cos(p1) * math.sin(d) * math.cos(t))
    l2 = math.radians(math.fmod(lon, 360.0)) + math.atan2(math.sin(t) * math.sin(d) * math.cos(p1), math.cos(d) - math.sin(p1) * math.sin(p2))
    out = math.degrees(l2)
    if out > 180:
        out = out - 360
    if out < -180:
        out = out + 360
    return [math.degrees(p2), out]

def spherical_excess(lat1: real, lon1: real, lat2: real, lon2: real, lat3: real, lon3: real) -> real:
    # The spherical excess of the triangle with three points as its corners,
    # in radians, by L'Huilier's theorem from its sides: the amount by which
    # its angles sum past π. Where the product of the four tangents is not a
    # normal real, its root is taken on the mantissas and the powers of two
    # math.frexp gives, so a small triangle keeps its excess.
    a = central_angle(lat2, lon2, lat3, lon3)
    b = central_angle(lat1, lon1, lat3, lon3)
    c = central_angle(lat1, lon1, lat2, lon2)
    s = (a + b + c) / 2
    t1 = math.tan(s / 2)
    t2 = math.tan((s - a) / 2)
    t3 = math.tan((s - b) / 2)
    t4 = math.tan((s - c) / 2)
    t = t1 * t2 * t3 * t4
    if (math.isfinite(t) and abs(t) >= 2.2250738585072014e-308) or t1 == 0.0 or t2 == 0.0 or t3 == 0.0 or t4 == 0.0:
        if t <= 0:
            return 0.0
        return 4 * math.atan(math.sqrt(t))
    f1 = math.frexp(t1)
    f2 = math.frexp(t2)
    f3 = math.frexp(t3)
    f4 = math.frexp(t4)
    m = f1[0] * f2[0] * f3[0] * f4[0]
    e = f1[1] + f2[1] + f3[1] + f4[1]
    if m <= 0:
        return 0.0
    if e % 2 != 0:
        m = m * 2.0
        e = e - 1
    return 4 * math.atan(math.ldexp(math.sqrt(m), e // 2))

def triangle_area(lat1: real, lon1: real, lat2: real, lon2: real, lat3: real, lon3: real) -> real:
    # The area of that triangle on the mean sphere, in square metres.
    r = mean_radius()
    return spherical_excess(lat1, lon1, lat2, lon2, lat3, lon3) * r * r

def mercator(lat: real) -> real:
    # The Mercator projection's northing for a latitude, in units of the
    # sphere's radius: ln(tan(45° + φ/2)).
    return math.log(math.tan(math.pi / 4 + math.radians(lat) / 2))

def mercator_inverse(y: real) -> real:
    # The latitude a Mercator northing comes from.
    return math.degrees(2 * math.atan(math.exp(y)) - math.pi / 2)

def mercator_scale(lat: real) -> real:
    # How much the projection enlarges a short distance at a latitude, 1/cos φ.
    return 1 / math.cos(math.radians(lat))

def dms(x: real) -> list<real>:
    # An angle as [degrees, minutes, seconds], each carrying the angle's sign.
    s = 1.0
    if x < 0:
        s = -1.0
    v = abs(x)
    d = float(math.floor(v))
    m = float(math.floor((v - d) * 60))
    return [s * d, s * m, s * ((v - d) * 60 - m) * 60]

# ── Calculation: the ellipsoid ───────────────────────────────────────────────

def vincenty_inverse(lat1: real, lon1: real, lat2: real, lon2: real) -> list<real>:
    # The distance between two points on the WGS84 ellipsoid, by Vincenty's
    # iteration of 1975: [metres, the azimuth from the first point to the
    # second, the azimuth from the second back to the first, iterations]. Near antipodal points the
    # iteration need not converge; after 200 rounds it stops and returns
    # [0.0, 0.0, 0.0, −1.0]. Each coordinate is taken modulo 360 degrees
    # first, as in central_angle, and where the square of the sine of the
    # arc falls below the smallest normal real, the sine is a hypot.
    a = 6378137.0
    f = 1 / 298.257223563
    b = a * (1 - f)
    l = math.radians(math.fmod(lon2, 360.0) - math.fmod(lon1, 360.0))
    u1 = math.atan((1 - f) * math.tan(math.radians(math.fmod(lat1, 360.0))))
    u2 = math.atan((1 - f) * math.tan(math.radians(math.fmod(lat2, 360.0))))
    su1 = math.sin(u1)
    cu1 = math.cos(u1)
    su2 = math.sin(u2)
    cu2 = math.cos(u2)
    lam = l
    done = -1.0
    ss = 0.0
    cs = 1.0
    sig = 0.0
    c2a = 1.0
    c2m = 0.0
    sl = 0.0
    cl = 1.0
    for k in range(200):
        sl = math.sin(lam)
        cl = math.cos(lam)
        py = cu2 * sl
        px = cu1 * su2 - su1 * cu2 * cl
        q = py * py + px * px
        ss = math.sqrt(q)
        if q < 2.2250738585072014e-308:
            ss = math.hypot(py, px)
        if ss == 0:
            return [0.0, 0.0, 0.0, k + 1.0]
        cs = su1 * su2 + cu1 * cu2 * cl
        sig = math.atan2(ss, cs)
        sa = cu1 * cu2 * sl / ss
        c2a = 1 - sa * sa
        c2m = 0.0
        if c2a != 0:
            c2m = cs - 2 * su1 * su2 / c2a
        c = f / 16 * c2a * (4 + f * (4 - 3 * c2a))
        prev = lam
        lam = l + (1 - c) * f * sa * (sig + c * ss * (c2m + c * cs * (-1 + 2 * c2m * c2m)))
        if abs(lam - prev) < 1e-12:
            done = k + 1.0
            break
    if done < 0:
        return [0.0, 0.0, 0.0, -1.0]
    u = c2a * (a * a - b * b) / (b * b)
    aa = 1 + u / 16384 * (4096 + u * (-768 + u * (320 - 175 * u)))
    bb = u / 1024 * (256 + u * (-128 + u * (74 - 47 * u)))
    ds = bb * ss * (c2m + bb / 4 * (cs * (-1 + 2 * c2m * c2m) - bb / 6 * c2m * (-3 + 4 * ss * ss) * (-3 + 4 * c2m * c2m)))
    s = b * aa * (sig - ds)
    az1 = math.degrees(math.atan2(cu2 * sl, cu1 * su2 - su1 * cu2 * cl))
    az2 = math.degrees(math.atan2(cu1 * sl, -su1 * cu2 + cu1 * su2 * cl)) + 180
    if az1 < 0:
        az1 = az1 + 360
    if az2 >= 360:
        az2 = az2 - 360
    return [s, az1, az2, done]

def vincenty_direct(lat1: real, lon1: real, az1: real, s: real) -> list<real>:
    # The point a distance s along the geodesic from a point at an azimuth,
    # on the WGS84 ellipsoid, by Vincenty's direct method: [lat, lon, the
    # azimuth from there back to the start]. The latitude and the azimuth are
    # taken modulo 360 degrees first, as in central_angle.
    a = 6378137.0
    f = 1 / 298.257223563
    b = a * (1 - f)
    t = math.radians(math.fmod(az1, 360.0))
    sa1 = math.sin(t)
    ca1 = math.cos(t)
    tu1 = (1 - f) * math.tan(math.radians(math.fmod(lat1, 360.0)))
    cu1 = 1 / math.sqrt(1 + tu1 * tu1)
    su1 = tu1 * cu1
    s1 = math.atan2(tu1, ca1)
    sa = cu1 * sa1
    c2a = 1 - sa * sa
    u = c2a * (a * a - b * b) / (b * b)
    aa = 1 + u / 16384 * (4096 + u * (-768 + u * (320 - 175 * u)))
    bb = u / 1024 * (256 + u * (-128 + u * (74 - 47 * u)))
    sig = s / (b * aa)
    c2m = 0.0
    ss = 0.0
    cs = 1.0
    for k in range(200):
        c2m = math.cos(2 * s1 + sig)
        ss = math.sin(sig)
        cs = math.cos(sig)
        ds = bb * ss * (c2m + bb / 4 * (cs * (-1 + 2 * c2m * c2m) - bb / 6 * c2m * (-3 + 4 * ss * ss) * (-3 + 4 * c2m * c2m)))
        prev = sig
        sig = s / (b * aa) + ds
        if abs(sig - prev) < 1e-12:
            break
    c2m = math.cos(2 * s1 + sig)
    ss = math.sin(sig)
    cs = math.cos(sig)
    x = su1 * ss - cu1 * cs * ca1
    p2 = math.atan2(su1 * cs + cu1 * ss * ca1, (1 - f) * math.sqrt(sa * sa + x * x))
    lam = math.atan2(ss * sa1, cu1 * cs - su1 * ss * ca1)
    c = f / 16 * c2a * (4 + f * (4 - 3 * c2a))
    l = lam - (1 - c) * f * sa * (sig + c * ss * (c2m + c * cs * (-1 + 2 * c2m * c2m)))
    az2 = math.degrees(math.atan2(sa, -x)) + 180
    if az2 >= 360:
        az2 = az2 - 360
    lon = lon1 + math.degrees(l)
    if lon > 180:
        lon = lon - 360
    if lon < -180:
        lon = lon + 360
    return [math.degrees(p2), lon, az2]