calculator.py

"""Deterministic backend: underdamped second-order step response.

Model: y(s)/u(s) = K * wn^2 / (s^2 + 2*zeta*wn*s + wn^2), step of size A at t = 0,
zero initial conditions, 0 < zeta < 1. Normalised response:
    c(t) = 1 - exp(-zeta*wn*t) / sqrt(1 - zeta^2) * sin(wd*t + theta)
with wd = wn*sqrt(1 - zeta^2) and theta = atan(sqrt(1 - zeta^2) / zeta) = acos(zeta).
Physical output is y(t) = y_final * c(t), where y_final = K*A (or F/k for a mass-spring-damper).
"""
import math

import numpy as np

VERSION = "1.0"

# Valid input ranges, enforced before any calculation runs
RANGES = {
    "zeta": (0.05, 0.95),        # underdamped only, away from the singular ends
    "wn": (0.1, 500.0),          # rad/s
    "mass": (0.01, 10_000.0),    # kg
    "stiffness": (1.0, 1e7),     # N/m
    "damping": (0.0, 1e6),       # N*s/m, zeta range is checked after conversion
    "force": (-1e6, 1e6),        # N, nonzero
    "amplitude": (-1e6, 1e6),    # step size in input units, nonzero
    "gain": (-1e6, 1e6),         # DC gain K, nonzero
    "os_limit": (0.1, 99.0),     # percent
    "ts_limit": (1e-4, 1e5),     # s
}
BANDS = {"2%": 0.02, "5%": 0.05}  # settling band as a fraction of the final value


class InputError(ValueError):
    """Raised with a user-facing message when inputs are out of scope."""


def _check(name, value, label, unit=""):
    """Validate one numeric input against RANGES and return it as float."""
    if value is None:
        raise InputError(f"Enter a value for {label}.")
    value = float(value)
    if not math.isfinite(value):
        raise InputError(f"{label} must be a finite number.")
    lo, hi = RANGES[name]
    if not lo <= value <= hi:
        raise InputError(f"{label} must be between {lo:g} and {hi:g} {unit}".strip() + ".")
    return value


def normalised_response(t, zeta, wn):
    """Closed-form unit step response c(t) of the standard second-order system."""
    root = math.sqrt(1.0 - zeta**2)
    wd = wn * root                          # damped natural frequency
    theta = math.acos(zeta)                 # phase, equals atan(root / zeta)
    return 1.0 - np.exp(-zeta * wn * t) / root * np.sin(wd * t + theta)


def _bisect(f, a, b, iters=80):
    """Root of f on [a, b] assuming a sign change; plain bisection for determinism."""
    fa = f(a)
    for _ in range(iters):
        m = 0.5 * (a + b)
        fm = f(m)
        if (fm > 0) == (fa > 0):
            a, fa = m, fm                   # root lies in the right half
        else:
            b = m                           # root lies in the left half
    return 0.5 * (a + b)


def _first_crossing(level, zeta, wn, t_max, n=20_000):
    """First time c(t) reaches a level, refined by bisection."""
    t = np.linspace(0.0, t_max, n)
    c = normalised_response(t, zeta, wn)
    idx = int(np.argmax(c >= level))        # first grid point at or above the level
    f = lambda x: float(normalised_response(np.array(x), zeta, wn)) - level
    return _bisect(f, t[max(idx - 1, 0)], t[idx])


def _settling_time(band, zeta, wn):
    """Last time |c(t) - 1| exceeds the band, searched up to the envelope bound."""
    root = math.sqrt(1.0 - zeta**2)
    sigma = zeta * wn
    t_env = -math.log(band * root) / sigma  # envelope is inside the band after this
    t = np.linspace(0.0, t_env, 40_000)
    err = np.abs(normalised_response(t, zeta, wn) - 1.0)
    outside = np.nonzero(err > band)[0]
    last = int(outside[-1])                 # last grid sample still outside the band
    f = lambda x: abs(float(normalised_response(np.array(x), zeta, wn)) - 1.0) - band
    return _bisect(f, t[last], t[min(last + 1, len(t) - 1)]), t_env


def _q(label, value, unit, formula=""):
    """One labelled quantity in the structured record."""
    return {"label": label, "value": float(value), "unit": unit, "formula": formula}


def compute(mode, wn=None, zeta=None, amplitude=1.0, gain=1.0,
            mass=None, damping=None, stiffness=None, force=None,
            band_name="2%", os_limit=None, ts_limit=None):
    """Validate inputs, compute every metric, and return a structured record.

    mode is "normalised" (wn, zeta, amplitude, gain) or "mass-spring-damper"
    (mass, damping, stiffness, force). Design limits are optional; None skips a check.
    """
    if band_name not in BANDS:
        raise InputError("Choose a settling band of 2% or 5%.")
    band = BANDS[band_name]
    inputs = {}

    if mode == "normalised":
        wn = _check("wn", wn, "Natural frequency ωn", "rad/s")
        zeta = _check("zeta", zeta, "Damping ratio ζ")
        amplitude = _check("amplitude", amplitude, "Step size A")
        gain = _check("gain", gain, "DC gain K")
        if amplitude == 0 or gain == 0:
            raise InputError("Step size and DC gain must be nonzero, otherwise the output never moves.")
        y_final, out_unit, out_name = gain * amplitude, "", "output"  # dimensionless output
        inputs["wn"] = _q("Natural frequency ωn", wn, "rad/s")
        inputs["zeta"] = _q("Damping ratio ζ", zeta, "")
        inputs["amplitude"] = _q("Step size A", amplitude, "")
        inputs["gain"] = _q("DC gain K", gain, "")
    elif mode == "mass-spring-damper":
        mass = _check("mass", mass, "Mass m", "kg")
        stiffness = _check("stiffness", stiffness, "Spring stiffness k", "N/m")
        damping = _check("damping", damping, "Damping coefficient c", "N·s/m")
        force = _check("force", force, "Step force F", "N")
        if force == 0:
            raise InputError("Step force must be nonzero, otherwise the mass never moves.")
        wn = math.sqrt(stiffness / mass)                     # rad/s
        zeta = damping / (2.0 * math.sqrt(stiffness * mass))
        lo, hi = RANGES["zeta"]
        if not lo <= zeta <= hi:
            kind = "overdamped or critically damped" if zeta >= 1 else "outside the supported range"
            raise InputError(
                f"These values give ζ = {zeta:.3f}, which is {kind}. This calculator covers "
                f"underdamped systems with {lo} ≤ ζ ≤ {hi}; adjust c, k, or m.")
        if not RANGES["wn"][0] <= wn <= RANGES["wn"][1]:
            raise InputError(f"These values give ωn = {wn:.3g} rad/s, outside 0.1 to 500 rad/s.")
        y_final, out_unit, out_name = force / stiffness * 1000.0, "mm", "displacement"
        inputs["mass"] = _q("Mass m", mass, "kg")
        inputs["damping"] = _q("Damping coefficient c", damping, "N·s/m")
        inputs["stiffness"] = _q("Spring stiffness k", stiffness, "N/m")
        inputs["force"] = _q("Step force F", force, "N")
    else:
        raise InputError("Choose an input mode.")

    if os_limit is not None:
        os_limit = _check("os_limit", os_limit, "Overshoot limit", "%")
    if ts_limit is not None:
        ts_limit = _check("ts_limit", ts_limit, "Settling time limit", "s")

    # Core closed-form quantities
    root = math.sqrt(1.0 - zeta**2)
    sigma = zeta * wn                                    # decay rate of the envelope
    wd = wn * root                                       # damped natural frequency
    theta = math.acos(zeta)                              # phase angle in the formula
    t_peak = math.pi / wd                                # first peak
    os_frac = math.exp(-zeta * math.pi / root)           # fractional overshoot
    t_rise_0_100 = (math.pi - theta) / wd                # first time c(t) = 1
    ts_textbook = (4.0 if band == 0.02 else 3.0) / sigma  # classic 4/σ or 3/σ rule

    # Numerical quantities, still deterministic
    ts_exact, t_env = _settling_time(band, zeta, wn)
    t10 = _first_crossing(0.1, zeta, wn, t_peak)
    t90 = _first_crossing(0.9, zeta, wn, t_peak)

    # Self-check: sampled peak should match the closed-form peak time and height
    grid = np.linspace(0.0, 2.0 * t_peak, 200_001)
    c = normalised_response(grid, zeta, wn)
    i_max = int(np.argmax(c))
    peak_time_err = abs(grid[i_max] - t_peak)
    peak_height_err = abs(c[i_max] - (1.0 + os_frac))

    derived = {}
    if mode == "mass-spring-damper":
        # Only physical inputs need wn and zeta derived; normalised inputs already list them
        derived["wn"] = _q("Natural frequency ωn", wn, "rad/s", "√(k/m)")
        derived["zeta"] = _q("Damping ratio ζ", zeta, "", "c / (2√(km))")
    derived.update({
        "wd": _q("Damped frequency ωd", wd, "rad/s", "ωn√(1−ζ²)"),
        "sigma": _q("Envelope decay rate σ", sigma, "1/s", "ζωn"),
        "theta": _q("Phase angle θ", math.degrees(theta), "deg", "atan(√(1−ζ²)/ζ)"),
        "period": _q("Damped period Td", 2 * math.pi / wd, "s", "2π/ωd"),
    })
    metrics = {
        "final": _q(f"Final {out_name}", y_final, out_unit, "K·A" if mode == "normalised" else "F/k"),
        "overshoot": _q("Percent overshoot", 100 * os_frac, "%", "100·exp(−ζπ/√(1−ζ²))"),
        "peak_value": _q(f"Peak {out_name}", y_final * (1 + os_frac), out_unit, "final·(1 + OS)"),
        "peak_time": _q("Peak time Tp", t_peak, "s", "π/ωd"),
        "rise_10_90": _q("Rise time, 10% to 90%", t90 - t10, "s", "numerical, bisection"),
        "rise_0_100": _q("Rise time, 0% to 100%", t_rise_0_100, "s", "(π−θ)/ωd"),
        "ts_exact": _q(f"Settling time, {band_name} band", ts_exact, "s", "numerical, last band exit"),
        "ts_textbook": _q("Settling time, textbook estimate", ts_textbook, "s",
                          "4/(ζωn)" if band == 0.02 else "3/(ζωn)"),
    }

    checks = []
    if os_limit is not None:
        # Minimum damping ratio that meets the overshoot limit
        ln_os = math.log(os_limit / 100.0)
        zeta_req = -ln_os / math.sqrt(math.pi**2 + ln_os**2)
        ok = 100 * os_frac <= os_limit
        checks.append({
            "name": "Overshoot", "value": 100 * os_frac, "limit": os_limit, "unit": "%",
            "result": "PASS" if ok else "FAIL",
            "requirement": _q("ζ needed for this overshoot limit", zeta_req, "", "−ln(OS)/√(π²+ln²OS)"),
        })
    if ts_limit is not None:
        # Envelope-based ωn that settles within the limit at the current ζ
        wn_req = -math.log(band * root) / (zeta * ts_limit)
        ok = ts_exact <= ts_limit
        checks.append({
            "name": "Settling time", "value": ts_exact, "limit": ts_limit, "unit": "s",
            "result": "PASS" if ok else "FAIL",
            "requirement": _q("ωn needed for this settling limit (envelope bound)", wn_req, "rad/s",
                              "−ln(band·√(1−ζ²)) / (ζ·Ts,max)"),
        })

    warnings = []
    if zeta < 0.2:
        warnings.append("Very light damping: many oscillations before settling; small parameter errors change results noticeably.")
    if zeta > 0.8:
        warnings.append("Close to critical damping: overshoot is tiny and the textbook settling estimate is least accurate here.")
    if abs(ts_exact - ts_textbook) / ts_exact > 0.25:
        warnings.append("The textbook settling estimate differs from the exact value by more than 25%.")

    return {
        "calculator": "Underdamped second-order step response",
        "version": VERSION,
        "mode": mode,
        "output_name": out_name,
        "output_unit": out_unit,
        "settling_band": band_name,
        "inputs": inputs,
        "derived": derived,
        "metrics": metrics,
        "checks": checks,
        "poles": {"real": -sigma, "imag": wd, "unit": "rad/s"},
        "assumptions": [
            "Linear, time-invariant second-order system with no zeros.",
            "Ideal step applied at t = 0 with zero initial position and velocity.",
            "Underdamped only: 0.05 ≤ ζ ≤ 0.95.",
            "Settling time is the last exit from the band around the final value.",
        ],
        "warnings": warnings,
        "self_check": {"peak_time_error_s": peak_time_err, "peak_height_error": peak_height_err},
        "_curve": {"zeta": zeta, "wn": wn, "y_final": y_final, "band": band,
                   "t_end": max(1.3 * ts_exact, 3 * t_peak)},
    }


def fmt(value, sig=4):
    """Readable number: sig figs, no scientific notation for everyday sizes."""
    if value == 0:
        return "0"
    mag = math.floor(math.log10(abs(value)))
    if -4 <= mag < 6:
        decimals = max(sig - 1 - mag, 0)
        text = f"{value:.{decimals}f}"
        return text.rstrip("0").rstrip(".") if "." in text else text
    return f"{value:.{sig - 1}e}"


def to_text(record):
    """Human-readable version of the record, the exact text the LLM receives."""
    lines = [f"{record['calculator']} (v{record['version']}, {record['mode']} inputs)", "", "Inputs:"]
    for q in record["inputs"].values():
        lines.append(f"  {q['label']}: {fmt(q['value'])} {q['unit']}".rstrip())
    lines.append("Derived parameters:")
    for q in record["derived"].values():
        lines.append(f"  {q['label']}: {fmt(q['value'])} {q['unit']}".rstrip())
    p = record["poles"]
    lines.append(f"  Poles: {fmt(p['real'])} ± j{fmt(p['imag'])} rad/s")
    lines.append("Step response metrics:")
    for q in record["metrics"].values():
        lines.append(f"  {q['label']}: {fmt(q['value'])} {q['unit']}".rstrip())
    if record["checks"]:
        lines.append("Design checks:")
        for ch in record["checks"]:
            req = ch["requirement"]
            lines.append(f"  {ch['name']}: {fmt(ch['value'])} {ch['unit']} vs limit {fmt(ch['limit'])} "
                         f"{ch['unit']} -> {ch['result']}; {req['label']}: {fmt(req['value'])} {req['unit']}".rstrip())
    else:
        lines.append("Design checks: none requested")
    if record["warnings"]:
        lines.append("Warnings:")
        lines += [f"  {w}" for w in record["warnings"]]
    return "\n".join(lines)


def to_rows(record):
    """Table rows (quantity, value, unit, formula) for the numerical results panel."""
    rows = []
    for group in ("metrics", "derived"):
        for q in record[group].values():
            rows.append([q["label"], fmt(q["value"]), q["unit"], q["formula"]])
    p = record["poles"]
    rows.append(["Poles", f"{fmt(p['real'])} ± j{fmt(p['imag'])}", "rad/s", "−ζωn ± jωd"])
    return rows