# SPDX-License-Identifier: MIT

# earthscience — classify climates and storms, and calculate the atmosphere and
# the solid Earth.
#
# Two halves that share nothing but the subject.
#
# STRUCTURE. You supply three kinds of record as facts.
#
# Stations (`earthscience.Station`) by their coldest and warmest month, their
# annual mean temperature in degrees Celsius, their annual rainfall in
# millimetres and the share of it that falls in the high-sun half of the
# year. The rules conclude the Köppen group of each, A (tropical), B (arid),
# C (temperate), D (continental) or E (polar), the kinds of B (desert or
# steppe) and of E (tundra or ice cap), and why a station is not in the group
# you expected:
#
#     import earthscience
#     earthscience.station("manaus", 26.0, 28.0, 27.0, 2300.0, 0.5, 0.0)
#     earthscience.station("cairo", 14.0, 29.0, 22.0, 25.0, 0.5, 0.0)
#     print(earthscience.climate("manaus"), earthscience.climate("cairo"))
#     n = why_not(earthscience.Station.cairo.tropical)
#
# Storms (`Storm`) by their sustained wind in miles per hour. The rules
# conclude whether it is a hurricane, its Saffir-Simpson category from 1 to
# 5, and whether it is a major one:
#
#     earthscience.storm("a", 74)
#     print(earthscience.category("a"), earthscience.category("b"))
#
# Minerals (`Mineral`) by their hardness. `minerals()` writes a pair for each
# two of them, and the rules conclude which scratches which:
#
#     earthscience.minerals(["quartz", "calcite", "topaz"], [7.0, 3.0, 8.0])
#     print(earthscience.scratched_by("quartz"), earthscience.scratches("quartz", "topaz"))
#
# Each conclusion answers `why`, and `why_not` names the premise that failed.
# A corrected fact withdraws what rested on it: a station's rainfall rewritten
# moves it out of the arid group. Where two sources disagree about a fact,
# both values are reasoned about and `unsettled()` names what they disagree
# about. No read refuses but `dryness_threshold()` and `category()`, which read
# one value and refuse a cell holding two.
#
# CALCULATION. Value scripts that return a value and write nothing: the dew
# point by the Magnus formula, the heat index, the wind chill, the standard
# atmosphere's temperature and pressure by height, potential temperature,
# radiometric age, the energy of an earthquake and a magnitude from an energy,
# the ratio of two magnitudes' energies, the S-P interval as a distance, the
# Coriolis parameter and the height of the lifting condensation level.
#
# Each formula holds over a range and is not checked outside it; the comment
# of each script says where it stops.
#
# Sources, fetched 2026-11-08:
#   - Köppen's groups are those of
#     https://en.wikipedia.org/wiki/K%C3%B6ppen_climate_classification.
#     A: the coldest month at least 18 C. E: the warmest month under 10 C; ET
#     where it is 0 C or more, EF below. B: regions that are not E and have
#     annual rain below 20 times the annual mean in C plus 280 mm where 70% or
#     more of the rain falls in the high-sun half-year, 140 mm for 30% to 70%,
#     and nothing for under 30%; under half the threshold is desert (BW) and
#     the rest steppe (BS). C: the coldest month above 0 C (or -3 C in some
#     versions) and under 18 C; D: the coldest at or below it, and the
#     warmest at least 10 C. The page leaves the exact values 0 C and 10 C
#     between its words "above" and "below"; this module puts a coldest month
#     of exactly the limit in D and a warmest month of exactly 10 C in C or D.
#     It tests E first, then B, as the page's "not E" says; the page gives no
#     order of B against A, C and D, and this module puts B before them.
#   - Saffir-Simpson, sustained wind in mph: 1 is 74-95, 2 96-110, 3 111-129,
#     4 130-156, 5 157 or higher; 3 to 5 are the major hurricanes
#     (https://www.nhc.noaa.gov/aboutsshws.php).
#   - Mohs introduced his scale in 1812 (a harder material scratches a softer
#     one); its ten reference minerals are talc 1, gypsum 2, calcite 3,
#     fluorite 4, apatite 5, orthoclase 6, quartz 7, topaz 8, corundum 9 and
#     diamond 10 (https://en.wikipedia.org/wiki/Mohs_scale_of_mineral_hardness).
#   - Magnus: Td = c·g / (b − g) with g = ln(RH/100) + b·T / (c + T), b = 17.625
#     and c = 243.04 C, chosen to fit -40 C to +50 C
#     (https://en.wikipedia.org/wiki/Dew_point).
#   - Heat index: the Rothfusz regression and its nine coefficients, fitted for
#     80 F and 40% humidity or more, within 1.3 F
#     (https://en.wikipedia.org/wiki/Heat_index).
#   - Wind chill, 35.74 + 0.6215T − 35.75·V^0.16 + 0.4275·T·V^0.16 in F and
#     mph, for air at 50 F or colder and wind above 3 mph
#     (https://www.weather.gov/safety/cold-wind-chill-chart).
#   - Pressure by height, p = 1013.25 (1 − h / 44307.694)^5.25530 hPa, which that
#     page marks as needing a citation; sea level in the standard atmosphere
#     is 15 C and 101325 Pa and the troposphere cools 6.5 C per km
#     (https://en.wikipedia.org/wiki/Pressure_altitude and
#     https://en.wikipedia.org/wiki/International_Standard_Atmosphere).
#   - Potential temperature θ = T (1000 hPa / P)^0.286 for air
#     (https://en.wikipedia.org/wiki/Potential_temperature).
#   - Gutenberg and Richter: log10 E = 4.8 + 1.5 M in joules for the
#     surface-wave magnitude, so a whole step is 10^1.5, about 32 times the
#     energy (https://en.wikipedia.org/wiki/Moment_magnitude_scale).
#   - The Coriolis parameter f = 2Ω sin φ with Ω = 7.2921e-5 rad/s
#     (https://en.wikipedia.org/wiki/Coriolis_frequency).
#   - Espy's lifting condensation level, 125 (T − Td) metres, accurate to
#     about 1% (https://en.wikipedia.org/wiki/Lifted_condensation_level).
#   - Uranium-238 decays to lead-206 with a half-life of about 4.463 billion
#     years (https://en.wikipedia.org/wiki/Uranium-238); the script takes the
#     half-life as an argument.
#
# Floats are IEEE doubles. `math.log10`, `math.log`, `math.sin` and `**` return
# the correctly rounded value, so a result agrees with CPython wherever
# CPython's own answer is right; a magnitude from an energy no longer goes
# through log(x) / log(10).

import math

# ── Climates as facts ────────────────────────────────────────────────────────

class Station:
    # A weather station, by a name you choose: `coldest` and `warmest`, the
    # mean temperatures in C of its coldest and warmest month, `annual_temp`
    # the annual mean in C, `rain` the annual rainfall in mm, `summer_share`
    # the fraction from 0 to 1 of the rain that falls in the half of the year
    # when the sun is high, and `boundary`, the coldest-month temperature that
    # divides C from D, 0.0 or -3.0 as the version of the system you follow
    # has it. `station()` writes them. The rules conclude `addition`, 280, 140
    # or 0 mm by the share; the computed `threshold`, 20 times the annual mean
    # plus that, is the rainfall under which the station is dry, and `half`
    # and `dry` follow from it. Then the groups: `polar` (E), a
    # warmest month under 10 C, with the rules `tundra` and
    # `icecap`; `arid` (B), with `desert` and `steppe`; `tropical` (A); `temperate` (C); and `continental` (D). Every station
    # belongs to exactly one group.
    coldest: float
    warmest: float
    annual_temp: float
    rain: float
    summer_share: float
    boundary: float
    addition: float
    threshold: float = 20.0 * annual_temp + addition
    half: float = threshold / 2.0
    dry: bool = rain < threshold
    polar: bool = warmest < 10.0
    tundra: bool
    icecap: bool
    arid: bool
    desert: bool
    steppe: bool
    tropical: bool
    temperate: bool
    continental: bool

@rule
def threshold_summer_rain():
    for c in Station:
        if c.summer_share >= 0.7:
            c.addition = 280.0

@rule
def threshold_even_rain():
    for c in Station:
        if c.summer_share >= 0.3 and c.summer_share < 0.7:
            c.addition = 140.0

@rule
def threshold_winter_rain():
    for c in Station:
        if c.summer_share < 0.3:
            c.addition = 0.0

@rule
def tundra():
    for c in Station:
        if c.polar == True and c.warmest >= 0.0:
            c.tundra = True

@rule
def icecap():
    for c in Station:
        if c.polar == True and c.warmest < 0.0:
            c.icecap = True

@rule
def arid():
    for c in Station:
        if c.polar != True and c.dry == True:
            c.arid = True

@rule
def desert():
    for c in Station:
        if c.arid == True and c.rain < c.half:
            c.desert = True

@rule
def steppe():
    for c in Station:
        if c.arid == True and c.rain >= c.half:
            c.steppe = True

@rule
def tropical():
    for c in Station:
        if c.polar != True and c.dry != True and c.coldest >= 18.0:
            c.tropical = True

@rule
def temperate():
    for c in Station:
        if c.polar != True and c.dry != True and c.coldest < 18.0 and c.coldest > c.boundary:
            c.temperate = True

@rule
def continental():
    for c in Station:
        if c.polar != True and c.dry != True and c.coldest <= c.boundary:
            c.continental = True

# ── Storms as facts ──────────────────────────────────────────────────────────

class Storm:
    # A tropical cyclone, by a name you choose, and `wind_mph`, its sustained
    # wind in miles per hour. The rules conclude `hurricane` from 74, its
    # Saffir-Simpson `category`, and `major` from category 3.
    wind_mph: int
    hurricane: bool
    category: int
    major: bool

@rule
def hurricane():
    for s in Storm:
        if s.wind_mph >= 74:
            s.hurricane = True

@rule
def category_one():
    for s in Storm:
        if s.wind_mph >= 74 and s.wind_mph <= 95:
            s.category = 1

@rule
def category_two():
    for s in Storm:
        if s.wind_mph >= 96 and s.wind_mph <= 110:
            s.category = 2

@rule
def category_three():
    for s in Storm:
        if s.wind_mph >= 111 and s.wind_mph <= 129:
            s.category = 3

@rule
def category_four():
    for s in Storm:
        if s.wind_mph >= 130 and s.wind_mph <= 156:
            s.category = 4

@rule
def category_five():
    for s in Storm:
        if s.wind_mph >= 157:
            s.category = 5

@rule
def major():
    for s in Storm:
        if s.wind_mph >= 111:
            s.major = True

# ── Minerals as facts ────────────────────────────────────────────────────────

class Mineral:
    # A mineral, by a name you choose, and `hardness` on a scale where a
    # harder one scratches a softer one, the Mohs scale or a finer one.
    hardness: float

class Pair:
    # Two minerals in order, named `<first>__<second>`, so a mineral's name
    # holds no two underscores in a row: `first` and `second` name them and
    # `h_first` and `h_second` are their hardness. `minerals()` writes one for
    # each ordered two. `scratches` is True where the first is strictly
    # harder; two of equal hardness scratch neither.
    first: string
    second: string
    h_first: float
    h_second: float
    scratches: bool = h_first > h_second

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

def station(name: string, coldest: float, warmest: float, annual_temp: float, rain: float, summer_share: float, boundary: float):
    # A station. summer_share is the fraction of the rain in the high-sun
    # half-year, from 0 to 1; boundary is 0.0 or -3.0.
    update(Station(name).coldest, coldest, "stations")
    update(Station(name).warmest, warmest, "stations")
    update(Station(name).annual_temp, annual_temp, "stations")
    update(Station(name).rain, rain, "stations")
    update(Station(name).summer_share, summer_share, "stations")
    update(Station(name).boundary, boundary, "stations")

def storm(name: string, wind_mph: int):
    update(Storm(name).wind_mph, wind_mph, "storms")

def minerals(names: list<string>, hardness: list<float>):
    # Each mineral, and a pair for each ordered two of them, so n minerals
    # write n² pairs.
    for i in range(len(names)):
        m = names[i]
        update(Mineral(m).hardness, hardness[i], "minerals")
    for i in range(len(names)):
        for j in range(len(names)):
            if i != j:
                s = names[i] + "__" + names[j]
                update(Pair(s).first, names[i], "minerals")
                update(Pair(s).second, names[j], "minerals")
                update(Pair(s).h_first, hardness[i], "minerals")
                update(Pair(s).h_second, hardness[j], "minerals")

# ── Reading what the rules concluded ─────────────────────────────────────────
#
# The readers that return names or a truth read every value in force, so a
# fact two sources disagree about is reported with both values.
# `dryness_threshold` and `category` return one number and read it with
# `current()`, which refuses a cell holding two different values rather than
# choose one: a station whose annual mean or summer share two sources
# disagree about has two thresholds, and a storm whose wind two sources put in
# two categories has two. `unsettled()` names the fact.

def groups(name: string) -> list<string>:
    # The Köppen groups a station is in, from "A", "B", "C", "D" and "E".
    # Exactly one for a station whose facts are all written.
    out = []
    for f in ["A", "B", "C", "D", "E"]:
        if f == "A":
            vs = list_values(Station.tropical)
        elif f == "B":
            vs = list_values(Station.arid)
        elif f == "C":
            vs = list_values(Station.temperate)
        elif f == "D":
            vs = list_values(Station.continental)
        else:
            vs = list_values(Station.polar)
        for v in vs:
            parts = v.path.split(".")
            if parts[len(parts) - 2] == name and v.value == True and f not in out:
                out = out + [f]
    return out

def climate(name: string) -> string:
    # The group with its kind where it has one: "A", "BW" (desert), "BS"
    # (steppe), "C", "D", "ET" (tundra) or "EF" (ice cap); "" before the
    # facts are written.
    g = groups(name)
    if len(g) == 0:
        return ""
    if g[0] == "B":
        vs = list_values(Station.desert)
        for v in vs:
            parts = v.path.split(".")
            if parts[len(parts) - 2] == name:
                return "BW"
        return "BS"
    if g[0] == "E":
        vs = list_values(Station.icecap)
        for v in vs:
            parts = v.path.split(".")
            if parts[len(parts) - 2] == name:
                return "EF"
        return "ET"
    return g[0]

def dryness_threshold(name: string) -> float:
    # The rainfall in mm under which a station is dry, 0.0 before its facts
    # are written.
    t = current(Station.name.threshold)
    match t:
        case active_claim:
            match t.value:
                case real as x:
                    return x
                case default:
                    pass
        case empty:
            pass
    return 0.0

def category(name: string) -> int:
    # A storm's Saffir-Simpson category, 0 below hurricane force.
    c = current(Storm.name.category)
    match c:
        case active_claim:
            match c.value:
                case int as k:
                    return k
                case default:
                    pass
        case empty:
            pass
    return 0

def scratches(first: string, second: string) -> bool:
    # Whether the first mineral is strictly harder than the second.
    vs = list_values(Pair.scratches)
    for v in vs:
        parts = v.path.split(".")
        if parts[len(parts) - 2] == first + "__" + second and v.value == True:
            return True
    return False

def scratched_by(name: string) -> list<string>:
    # The minerals the named one scratches, sorted.
    out = []
    vs = list_values(Pair.scratches)
    for v in vs:
        parts = v.path.split(".")
        p = parts[len(parts) - 2]
        if p.startswith(name + "__") and v.value == True:
            x = p[len(name) + 2:]
            if x not in out:
                out = out + [x]
    return sorted(out)

def unsettled() -> list<string>:
    # Each station, storm and mineral figure that holds two different values
    # for a fact you write, as "station x rain" and so on, sorted.
    out = []
    for f in ["coldest", "warmest", "annual_temp", "rain", "summer_share", "boundary", "wind_mph", "hardness"]:
        if f == "coldest":
            vs = list_values(Station.coldest)
            label = "station "
        elif f == "warmest":
            vs = list_values(Station.warmest)
            label = "station "
        elif f == "annual_temp":
            vs = list_values(Station.annual_temp)
            label = "station "
        elif f == "rain":
            vs = list_values(Station.rain)
            label = "station "
        elif f == "summer_share":
            vs = list_values(Station.summer_share)
            label = "station "
        elif f == "boundary":
            vs = list_values(Station.boundary)
            label = "station "
        elif f == "wind_mph":
            vs = list_values(Storm.wind_mph)
            label = "storm "
        else:
            vs = list_values(Mineral.hardness)
            label = "mineral "
        seen = {}
        for v in vs:
            parts = v.path.split(".")
            s = label + parts[len(parts) - 2] + " " + parts[len(parts) - 1]
            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)

# ── Calculation: the atmosphere ──────────────────────────────────────────────

def dew_point(t: float, rh: float) -> float:
    # The dew point in C of air at t C and rh percent humidity, 0 < rh <= 100,
    # by the Magnus formula; the constants were fitted for -40 C to 50 C. At
    # 100% it returns the temperature, to within a rounding. A humidity outside
    # 0 < rh <= 100 refuses. Where b t passes the largest real, b t / (c + t) is
    # taken as b (t / (c + t)).
    if not (rh > 0.0 and rh <= 100.0):
        raise ValueError("a relative humidity is above 0 and at most 100")
    b = 17.625
    c = 243.04
    term = b * t / (c + t)
    if not math.isfinite(term) and math.isfinite(t):
        term = b * (t / (c + t))
    g = math.log(rh / 100.0) + term
    return c * g / (b - g)

def heat_index(t: float, rh: float) -> float:
    # The heat index in F for air at t F and rh percent, the Rothfusz
    # regression. It was fitted for 80 F and 40% or more and is within 1.3 F
    # there; below that it returns a value that means little.
    return -42.379 + 2.04901523 * t + 10.14333127 * rh - 0.22475541 * t * rh - 0.00683783 * t * t - 0.05481717 * rh * rh + 0.00122874 * t * t * rh + 0.00085282 * t * rh * rh - 0.00000199 * t * t * rh * rh

def wind_chill(t: float, v: float) -> float:
    # The wind chill in F for air at t F and a wind of v mph. It is defined for
    # air at 50 F or colder and a wind above 3 mph.
    w = v ** 0.16
    return 35.74 + 0.6215 * t - 35.75 * w + 0.4275 * t * w

def temperature_at(h: float) -> float:
    # The temperature in K at h metres in the standard atmosphere's
    # troposphere, 288.15 K at sea level and 6.5 K cooler per km, to 11,000 m.
    return 288.15 - 0.0065 * h

def pressure_at(h: float) -> float:
    # The pressure in hPa at h metres in the standard atmosphere,
    # 1013.25 (1 − h / 44307.694)^5.2553. It holds in the troposphere, to
    # about 11,000 m, and h is under 44,307.
    return 1013.25 * (1.0 - h / 44307.694) ** 5.2553

def potential_temperature(t_kelvin: float, p_hpa: float) -> float:
    # The temperature in K the air would have if brought dry-adiabatically to
    # 1000 hPa, T (1000 / p)^0.286. Where that is not a normal real while T and
    # p are, it is T e^(0.286 (ln 1000 - ln p)), so a pressure whose 1000 / p
    # passes the largest real still has its temperature.
    if t_kelvin == 0.0 and p_hpa > 0.0:
        return t_kelvin
    v = t_kelvin * (1000.0 / p_hpa) ** 0.286
    if (math.isfinite(v) and abs(v) >= 2.2250738585072014e-308) or not (p_hpa > 0.0 and math.isfinite(p_hpa)) or not math.isfinite(t_kelvin):
        return v
    return math.copysign(math.exp(math.log(abs(t_kelvin)) + 0.286 * (math.log(1000.0) - math.log(p_hpa))), t_kelvin)

def coriolis(lat_deg: float) -> float:
    # The Coriolis parameter in 1/s at a latitude in degrees, 2Ω sin φ with
    # Ω = 7.2921e-5 rad/s.
    return 2.0 * 7.2921e-5 * math.sin(math.radians(lat_deg))

def lcl_height(t: float, td: float) -> float:
    # The height in metres at which rising air reaches saturation, Espy's
    # 125 (T − Td), for t and the dew point td in the same unit; about 1%
    # accurate in ordinary air.
    return 125.0 * (t - td)

# ── Calculation: the solid Earth ─────────────────────────────────────────────

def radiometric_age(daughter_per_parent: float, half_life: float) -> float:
    # The age of a sample, in the unit of the half-life, from the ratio of
    # daughter atoms to parent atoms now: half_life · log2(1 + D/P). It
    # assumes no daughter at the start and none gained or lost since.
    return half_life * math.log2(1.0 + daughter_per_parent)

def quake_energy_joules(m: float) -> float:
    # The energy in joules radiated by an earthquake of surface-wave
    # magnitude m, 10^(1.5 m + 4.8).
    return 10.0 ** (1.5 * m + 4.8)

def magnitude_from_energy(joules: float) -> float:
    # The magnitude whose energy is the given joules: (log10 E − 4.8) / 1.5.
    # It inverts quake_energy_joules to within a rounding.
    return (math.log10(joules) - 4.8) / 1.5

def quake_ratio(m1: float, m2: float) -> float:
    # How many times more energy a magnitude m2 releases than m1,
    # 10^(1.5 (m2 − m1)): one step is 31.6, two are 1000.
    return 10.0 ** (1.5 * (m2 - m1))

def sp_distance_km(sp_seconds: float, vp: float, vs: float) -> float:
    # The distance in km to a quake from the interval in seconds between the
    # P and the S wave, for wave speeds vp and vs in km/s: the two travel the
    # same distance d, d/vs − d/vp = interval. vp − vs is taken of halves where
    # it passes the largest real, and where the distance is then not a normal
    # real, it is worked on the mantissas and the powers of two apart, as
    # math.frexp gives them.
    d = vp - vs
    if math.isinf(d) and math.isfinite(vp) and math.isfinite(vs):
        return sp_seconds * vp * vs / (0.5 * vp - 0.5 * vs) * 0.5
    v = sp_seconds * vp * vs / d
    if math.isfinite(v) and abs(v) >= 2.2250738585072014e-308:
        return v
    a = math.frexp(sp_seconds)
    b = math.frexp(vp)
    c = math.frexp(vs)
    e = math.frexp(d)
    return math.ldexp(a[0] * b[0] * c[0] / e[0], a[1] + b[1] + c[1] - e[1])