# SPDX-License-Identifier: MIT

# sound — state the rules of a noise dose and of the range of hearing, and
# calculate sound.
#
# THIS MODULE IS NOT OCCUPATIONAL-HEALTH ADVICE. It states published limits as
# data and arithmetic with published formulas. A limit is a figure on a page;
# whether a worker is protected is decided by a measurement and a hearing
# conservation programme, not by this module, and a dose over one is not a
# finding of hearing loss. It gives no diagnosis and no treatment. Ask an
# occupational-health professional.
#
# Two halves that share nothing but the subject.
#
# STRUCTURE. You supply three kinds of record as facts.
#
# Exposures (`sound.Exposure`), each a level in decibels (A-weighted, taken as
# given, never computed from a spectrum) held for some hours, with the limit it
# is judged against: a reference level for eight hours and an exchange rate,
# the decibels by which a level may rise while the permitted time halves.
# `exposure()` writes the permitted hours from that pair. The dose is the hours
# over the permitted hours, and the rule concludes an exposure whose dose is
# more than one is over its limit; a dose of exactly one is at the limit and is
# not over. An exposure over its limit is marked and kept, never dropped:
#
#     import sound
#     sound.exposure("e1", 95.0, 4.0, 85.0, 3.0)
#     sound.exposure("e2", 95.0, 4.0, 90.0, 5.0)
#     print(sound.over_limit_names(), sound.dose("e1"), sound.dose("e2"))
#     n = why_not(sound.Exposure.e2.over)
#
# Workers' days (`Worker`, with `Segment`s). A day of mixed exposures is judged
# by the sum of the doses of its segments, and the rule concludes a day whose
# sum is more than one:
#
#     sound.segment("w1", "s1", 95.0, 2.0, 85.0, 3.0)
#     sound.segment("w1", "s2", 88.0, 3.0, 85.0, 3.0)
#     print(sound.days_over(), sound.total_of("w1"))
#
# Tones (`Tone`), by frequency in hertz. The rules conclude whether a frequency
# is in the range of hearing, 20 to 20,000 Hz with both ends inside, below it
# (infrasound) or above it (ultrasound):
#
#     sound.tone("t1", 440.0)
#     print(sound.range_of("t1"))
#
# Each conclusion answers `why`, and `why_not` names the premise that failed
# (the first one, when several fail). A corrected fact withdraws what rested on
# it: an exposure rewritten with fewer hours is no longer over.
#
# Where two sources disagree about a fact, such as one exposure's hours, the
# rules reason about both values and the names they return hold both verdicts.
# `dose` returns one number and refuses a cell that holds two values.
# `unsettled()` names every fact you write that two sources give different
# values.
#
# A day's total is the exception, since a sum adds every value it reads: a
# segment whose hours two sources give as 4 and 3 would be counted under both,
# and a day under the limit by either source would be over it. So while a
# segment is bound to its worker's day twice, as two values of its hours,
# level or permitted hours bind it, the day is `disputed`: no total and no
# verdict is concluded for it, `disputed_days()` names it, `total_of` refuses
# it, and `why_not` on its total names the dispute. Forgetting the claim that
# is not to stand settles it. A segment whose worker two sources disagree on
# counts in each worker's day a value names, so each of those totals is one
# source's.
#
# CALCULATION. Value scripts that return a value and write nothing: the hours
# a limit permits, the level in decibels of a sound pressure, the level of two
# incoherent sources added, the level at a distance in free field, the speed
# of sound at a temperature, a wavelength, the Doppler shift and the 8-hour
# average level of a dose.
#
# Results are reals and are not rounded. `math.log10` is the correctly rounded
# logarithm, and CPython's can differ from it by one unit in the last place at
# some arguments (the level of a pressure of 1.4142135623730951 against a
# reference of 1 is 3.0102999566398125 here and 3.010299956639812 in CPython
# 3.11), so a decibel agrees with CPython's to that unit, not always to the
# last bit. Nothing here claims otherwise. A pressure of zero has no level and
# refuses, as does a distance of zero. Two levels are added in a form whose
# powers stay in range, so any two levels have a sum.
#
# What the rules and scripts do not cover: frequency weighting, octave bands,
# room acoustics and reverberation, sound in water, the action level and the
# hearing conservation programme of either limit, and any limit of another
# country.
#
# Sources, fetched 2026-11-09:
#   - The Bell System renamed its Transmission Unit, introduced in 1924, the
#     decibel in 1928, and the bel is named for Alexander Graham Bell. A power
#     quantity is 10 log10(P/P0) decibels and a root-power (field) quantity,
#     which a sound pressure is, 20 log10(F/F0); the reference pressure for
#     sound in air is 20 micropascals; a doubling of power is 3.0103 dB,
#     usually rounded to 3 dB (https://en.wikipedia.org/wiki/Decibel).
#   - OSHA 29 CFR 1910.95, Table G-16: 90 dBA for 8 hours, 92 for 6, 95 for 4,
#     97 for 3, 100 for 2, 102 for 1.5, 105 for 1, 110 for half an hour and 115
#     for a quarter of an hour or less; where a day has periods at different
#     levels, if the sum C1/T1 + C2/T2 + ... + Cn/Tn "exceeds unity" the mixed
#     exposure "should be considered to exceed the limit value", Cn being the
#     time at a level and Tn the time permitted at it
#     (https://www.osha.gov/laws-regs/regulations/standardnumber/1910/1910.95).
#     Its Appendix A gives the dose D = 100 C/T, for several levels D = 100
#     (C1/T1 + C2/T2 + ... + Cn/Tn), and the 8-hour time-weighted average TWA =
#     16.61 log10(D/100) + 90
#     (https://www.osha.gov/laws-regs/regulations/standardnumber/1910/1910.95AppA).
#   - NIOSH recommends that exposures be controlled below a level equivalent to
#     85 dBA for eight hours, with a 3 dBA exchange rate at which each rise of
#     3 dBA halves the allowed time, where OSHA's rule uses 5 dB
#     (https://en.wikipedia.org/wiki/Noise-induced_hearing_loss); NIOSH's own
#     page states the 85 dBA average over an eight-hour workday and no
#     exchange rate (https://www.cdc.gov/niosh/noise/about/index.html).
#     THE DURATION FORMULA THIS MODULE USES, 8 / 2^((L - R) / E) for a level L,
#     a reference level R and an exchange rate E, IS THIS MODULE'S OWN
#     DERIVATION from a reference level, eight hours and an exchange rate. It
#     is not a formula read from NIOSH or from OSHA, whose Appendix A shows its
#     own as an image; it reproduces Table G-16's 95 dBA for 4 hours and 100
#     dBA for 2 hours at 90 dBA and 5 dB, and the module's tests check that. The
#     table's other rows are rounded where the derivation is not: it gives 6.06
#     hours at 92 dBA where the table has 6, 3.03 at 97 where it has 3, and
#     1.52 at 102 where it has 1.5.
#   - The commonly stated range of human hearing is 20 to 20,000 Hz, and it
#     shrinks during life, usually beginning at around the age of eight with the
#     upper limit being reduced (https://en.wikipedia.org/wiki/Hearing_range).
#     That page gives no formal cut-offs for the names; "infrasound" and
#     "ultrasound" are used here for what lies below 20 Hz and above 20 kHz.
#   - For a point source the sound pressure falls by half when the distance
#     doubles, 6.02 dB per doubling, while the intensity falls to a quarter
#     (https://en.wikipedia.org/wiki/Inverse-square_law).
#   - The speed of sound in dry air is about 343 m/s at 20 C (the page's table
#     gives 343.2) and about 331 m/s at 0 C (331.3 theoretical), and c is
#     331.32 m/s x sqrt(1 + t/273.15) for a Celsius temperature t, which is the
#     constant used below; Newton computed the speed with isothermal and not
#     adiabatic assumptions, and Laplace later corrected it
#     (https://en.wikipedia.org/wiki/Speed_of_sound).
#   - Doppler first proposed the effect in 1842, in a treatise on the coloured
#     light of double stars, and Buys Ballot tested it for sound in 1845,
#     finding the pitch higher as the source approached and lower as it
#     receded; the observed frequency is f = (v +/- v_r) / (v -/+ v_s) f0, the
#     receiver's speed added when it moves toward the source and the source's
#     speed subtracted when it moves toward the receiver
#     (https://en.wikipedia.org/wiki/Doppler_effect).

import math

# ── Exposures as facts ───────────────────────────────────────────────────────

class Exposure:
    # A level held for some hours, by a name you choose. `ref` and `rate` are
    # the limit it is judged against, a reference level for eight hours and an
    # exchange rate, both in decibels; `allowed` is the permitted hours, which
    # `exposure()` writes from them because a computed field cannot take a
    # power. `dose` is `hours / allowed`, and `over` is concluded where it is
    # more than one.
    level: real
    hours: real
    ref: real
    rate: real
    allowed: real
    dose: real = hours / allowed
    over: bool

@rule
def over_limit():
    for e in Exposure:
        if e.dose > 1.0:
            e.over = True

# ── A worker's day as facts ──────────────────────────────────────────────────

class Segment:
    # One period of a worker's day at one level, by a name you choose.
    # `worker` names the day it belongs to, and `allowed` is the hours the
    # limit permits at that level, written by `segment()`.
    worker: string
    level: real
    hours: real
    allowed: real
    dose: real = hours / allowed

class Worker:
    # A worker's day, by a name you choose. `total` is the sum of the doses of
    # its segments, concluded by the rule; `over` is concluded where the total
    # is more than one. `label` is written so that the day exists.
    # The rules also conclude `claims` and `segments`, which `disputed`
    # compares: a segment bound to the day twice makes it `disputed`, and then
    # no total is concluded.
    label: string
    claims: int
    segments: int
    disputed: bool
    total: real
    over: bool

@rule
def day_claims():
    for w in Worker, s in Segment:
        if s.worker == w and s.dose > -1.0e300:
            w.claims = count(s)

@rule
def day_segments():
    for w in Worker, s in Segment:
        if s.worker == w and s.dose > -1.0e300:
            w.segments = distinct(s)

@rule
def day_disputed():
    # A sum adds every value it reads, so a segment bound to the day twice
    # would be counted under both values, which no source gives.
    for w in Worker:
        if w.claims > w.segments:
            w.disputed = True

@rule
def day_total():
    for w in Worker, s in Segment:
        if s.worker == w and absent(w.disputed):
            w.total = sum(s.dose)

@rule
def day_over():
    for w in Worker:
        if w.total > 1.0:
            w.over = True

# ── Tones as facts ───────────────────────────────────────────────────────────

class Tone:
    # A frequency in hertz, by a name you choose. `audible` is concluded for 20
    # to 20,000 Hz with both ends inside, `infrasound` below 20 and
    # `ultrasound` above 20,000.
    hz: real
    audible: bool
    infrasound: bool
    ultrasound: bool

@rule
def audible_tone():
    for t in Tone:
        if t.hz >= 20.0 and t.hz <= 20000.0:
            t.audible = True

@rule
def infrasound_tone():
    for t in Tone:
        if t.hz < 20.0:
            t.infrasound = True

@rule
def ultrasound_tone():
    for t in Tone:
        if t.hz > 20000.0:
            t.ultrasound = True

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

def permissible_hours(level: float, ref_level: float, exchange_rate: float) -> float:
    # The hours a day may hold at a level: 8 / 2^((level - ref_level) /
    # exchange_rate). It is this module's derivation, not a published formula;
    # at 85 dBA and 3 dB, or at 90 dBA and 5 dB, it gives 8 hours at the
    # reference level and half the time for each rise of the exchange rate.
    # The rise is taken of halves where level - ref_level passes the largest
    # real.
    d = level - ref_level
    if math.isinf(d) and math.isfinite(level) and math.isfinite(ref_level):
        return 8.0 / (2.0 ** ((0.5 * level - 0.5 * ref_level) / exchange_rate * 2.0))
    return 8.0 / (2.0 ** (d / exchange_rate))

def exposure(name: string, level: float, hours: float, ref_level: float, exchange_rate: float):
    # A level held for some hours, judged against a reference level and an
    # exchange rate. Write it again with other hours to correct it.
    update(Exposure(name).level, level, "exposure")
    update(Exposure(name).hours, hours, "exposure")
    update(Exposure(name).ref, ref_level, "exposure")
    update(Exposure(name).rate, exchange_rate, "exposure")
    update(Exposure(name).allowed, permissible_hours(level, ref_level, exchange_rate), "exposure")

def segment(worker: string, name: string, level: float, hours: float, ref_level: float, exchange_rate: float):
    # One period of a worker's day. A name is used once; the worker's day
    # exists from the first segment written for it.
    update(Worker(worker).label, worker, "segment")
    update(Segment(name).worker, worker, "segment")
    update(Segment(name).level, level, "segment")
    update(Segment(name).hours, hours, "segment")
    update(Segment(name).allowed, permissible_hours(level, ref_level, exchange_rate), "segment")

def tone(name: string, hertz: float):
    update(Tone(name).hz, hertz, "tone")

# ── Reading what the rules concluded ─────────────────────────────────────────
#
# The readers that return names, a range or a truth read every value in force,
# so a fact two sources disagree about gives each of its values a verdict and
# is reported with both. `dose` and `total_of` return one number and read it
# with `current()`, which refuses a cell holding two different values rather
# than choose one. A day's total is such a cell only when it holds two: it sums
# every dose a segment holds, so a segment with two doses adds both into one
# value, and `unsettled()` is what names the segment.

def holds(values: claim_list, name: string) -> bool:
    for v in values:
        parts = v.path.split(".")
        if parts[len(parts) - 2] == name and v.value == True:
            return True
    return False

def names_of(values: claim_list) -> list<string>:
    out = []
    for v in values:
        parts = v.path.split(".")
        x = parts[len(parts) - 2]
        if v.value == True and x not in out:
            out = out + [x]
    return sorted(out)

def over_limit_names() -> list<string>:
    # The exposures whose dose is more than one, sorted.
    vs = list_values(Exposure.over)
    return names_of(vs)

def days_over() -> list<string>:
    # The workers whose segments' doses sum to more than one, sorted.
    vs = list_values(Worker.over)
    return names_of(vs)

def disputed_days() -> list<string>:
    # The workers a segment is bound to twice, whose total waits until one of
    # its values is forgotten, sorted.
    vs = list_values(Worker.disputed)
    return names_of(vs)

def dose(name: string) -> float:
    # An exposure's dose, hours over permitted hours; 0.0 before the facts are
    # written.
    c = current(Exposure.name.dose)
    match c:
        case active_claim:
            match c.value:
                case real as x:
                    return x
                case default:
                    pass
        case empty:
            pass
    return 0.0

def total_of(name: string) -> float:
    # The sum of the doses of a worker's segments; 0.0 before the rule
    # concludes it. It refuses a disputed day, which has no total until the
    # dispute is settled.
    if name in disputed_days():
        raise ValueError("day " + name + " is disputed: a segment of it is bound to it twice; forget the claim that is not to stand")
    c = current(Worker.name.total)
    match c:
        case active_claim:
            match c.value:
                case real as x:
                    return x
                case default:
                    pass
        case empty:
            pass
    return 0.0

def range_of(name: string) -> string:
    # "audible", "infrasound" or "ultrasound" for a tone; "" before the facts
    # are written.
    a = list_values(Tone.audible)
    b = list_values(Tone.infrasound)
    c = list_values(Tone.ultrasound)
    if holds(a, name):
        return "audible"
    if holds(b, name):
        return "infrasound"
    if holds(c, name):
        return "ultrasound"
    return ""

def unsettled() -> list<string>:
    # Each exposure, segment, worker and tone that holds two different values
    # for a fact you write, as "exposure x hours", "segment x worker" and so
    # on, sorted. A fact is named only when its values differ.
    out = []
    for f in ["e.level", "e.hours", "e.ref", "e.rate", "e.allowed", "s.worker", "s.level", "s.hours", "s.allowed", "w.label", "t.hz"]:
        if f == "e.level":
            vs = list_values(Exposure.level)
            noun = "exposure "
        elif f == "e.hours":
            vs = list_values(Exposure.hours)
            noun = "exposure "
        elif f == "e.ref":
            vs = list_values(Exposure.ref)
            noun = "exposure "
        elif f == "e.rate":
            vs = list_values(Exposure.rate)
            noun = "exposure "
        elif f == "e.allowed":
            vs = list_values(Exposure.allowed)
            noun = "exposure "
        elif f == "s.worker":
            vs = list_values(Segment.worker)
            noun = "segment "
        elif f == "s.level":
            vs = list_values(Segment.level)
            noun = "segment "
        elif f == "s.hours":
            vs = list_values(Segment.hours)
            noun = "segment "
        elif f == "s.allowed":
            vs = list_values(Segment.allowed)
            noun = "segment "
        elif f == "w.label":
            vs = list_values(Worker.label)
            noun = "worker "
        else:
            vs = list_values(Tone.hz)
            noun = "tone "
        seen = {}
        for v in vs:
            parts = v.path.split(".")
            s = noun + 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: sound ───────────────────────────────────────────────────────

def level_of_pressure(pascals: float, reference_pascals: float) -> float:
    # The sound pressure level in decibels, 20 log10(p / p0); the reference
    # for sound in air is 0.00002 pascal. A pressure, and its reference, of
    # zero or below has no level and refuses with a math domain error. Where
    # p / p0 is not a normal real, it is 20 (log10 p - log10 p0), whose two
    # logarithms are then far apart.
    if not (pascals > 0.0 and reference_pascals > 0.0):
        raise ValueError("math domain error: a pressure and its reference are above 0")
    q = pascals / reference_pascals
    if not (math.isfinite(q) and q >= 2.2250738585072014e-308):
        return 20.0 * (math.log10(pascals) - math.log10(reference_pascals))
    return 20.0 * math.log10(q)

def add_levels(a: float, b: float) -> float:
    # The level of two incoherent sources of levels a and b decibels,
    # 10 log10(10^(a/10) + 10^(b/10)). Two equal sources add 3.0103 dB. It is
    # computed as the louder level plus 10 log10(1 + 10^(-d/10)), d the
    # difference, which is the same sum: no power passes the largest real, and
    # the answer is the exact one rounded more often than the sum of powers.
    louder = max(a, b)
    quieter = min(a, b)
    if louder == quieter:
        return louder + 10.0 * math.log10(2.0)
    return louder + 10.0 * math.log10(1.0 + 10.0 ** ((quieter - louder) / 10.0))

def at_distance(level_at_r1: float, r1: float, r2: float) -> float:
    # The level at distance r2 from a point source whose level at r1 is given,
    # in free field: level - 20 log10(r2 / r1), a fall of 6.02 dB for each
    # doubling of the distance. A distance of zero or below refuses with a math
    # domain error. Where r2 / r1 is not a normal real, its logarithm is
    # log10 r2 - log10 r1.
    if not (r1 > 0.0 and r2 > 0.0):
        raise ValueError("math domain error: a distance is above 0")
    q = r2 / r1
    if not (math.isfinite(q) and q >= 2.2250738585072014e-308):
        return level_at_r1 - 20.0 * (math.log10(r2) - math.log10(r1))
    return level_at_r1 - 20.0 * math.log10(q)

def speed_of_sound(celsius: float) -> float:
    # The speed of sound in dry air in metres a second, 331.32 x sqrt(1 +
    # t / 273.15). Below -273.15 C it refuses with a math domain error.
    return 331.32 * math.sqrt(1.0 + celsius / 273.15)

def wavelength(speed: float, hertz: float) -> float:
    # The wavelength in metres, speed over frequency.
    return speed / hertz

def doppler(source_hz: float, speed: float, receiver_toward: float, source_toward: float) -> float:
    # The observed frequency, f0 (c + vr) / (c - vs): `receiver_toward` is the
    # receiver's speed toward the source and `source_toward` the source's
    # toward the receiver, each negative for motion away. A source moving at
    # the speed of sound refuses with a division by zero. The sums are taken of
    # halves where one passes the largest real, and where f0 (c + vr) / (c - vs)
    # is then not a normal real, it is worked on the mantissas and the powers
    # of two apart, as math.frexp gives them.
    num = speed + receiver_toward
    den = speed - source_toward
    if (math.isinf(num) or math.isinf(den)) and math.isfinite(speed) and math.isfinite(receiver_toward) and math.isfinite(source_toward):
        num = 0.5 * speed + 0.5 * receiver_toward
        den = 0.5 * speed - 0.5 * source_toward
    v = source_hz * num / den
    if math.isfinite(v) and abs(v) >= 2.2250738585072014e-308:
        return v
    a = math.frexp(source_hz)
    b = math.frexp(num)
    c = math.frexp(den)
    return math.ldexp(a[0] * b[0] / c[0], a[1] + b[1] - c[1])

def twa_from_dose(dose_percent: float) -> float:
    # The 8-hour time-weighted average level of a dose written in percent, 16.61
    # log10(D / 100) + 90, as OSHA's Appendix A gives it.
    return 16.61 * math.log10(dose_percent / 100.0) + 90.0