"""Three-element Windkessel — single source of truth.

Defines the exact discrete scheme that BOTH the Excel sheet and the PDF
worked example replicate, so every number agrees to the digit.

Model (Westerhof 3-element Windkessel):
    Rc  = characteristic (proximal) resistance   [mmHg*s/mL]
    Rd  = distal (peripheral) resistance          [mmHg*s/mL]
    C   = total arterial compliance               [mL/mmHg]

Exact decomposition used here (identical to the full 3WK ODE):
    node pressure behind the compliance, 2-element WK driven by inflow Q:
        dP_c/dt = Q/C - P_c/(Rd*C)
    aortic pressure adds the characteristic-resistance jump:
        P(t) = P_c(t) + Rc*Q(t)

Discrete (explicit Euler, step dt), i = 0,1,2,...:
    P_c[i+1] = P_c[i] + dt*( Q[i]/C - P_c[i]/(Rd*C) )
    P[i]     = P_c[i] + Rc*Q[i]

Inflow waveform Q(t) over a cycle of period T: half-sine ejection of
duration Ts, zero during diastole.
"""
from __future__ import annotations

import math

# --- fixed inputs (shared by Excel + PDF) ---
HR = 75.0          # heart rate            [beats/min]
T = 60.0 / HR      # cardiac period        [s]           -> 0.8 s
TS = 0.30          # systolic ejection     [s]
SV = 70.0          # stroke volume         [mL]
MAP_TARGET = 92.0  # target mean pressure  [mmHg]
PP_TARGET = 40.0   # target pulse pressure [mmHg]
RC_FRACTION = 0.05 # Rc as a fraction of total resistance R_T

DT = 0.01          # time step             [s]
N_CYCLES = 6       # cycles simulated to reach periodic steady state

Q_PEAK = math.pi * SV / (2.0 * TS)   # half-sine peak so integral = SV
Q_MEAN = SV / T                      # mean inflow = cardiac output [mL/s]


def q_of_t(t: float) -> float:
    """Inflow Q at time t (s). Half-sine ejection then zero. Periodic in T."""
    phase = t % T
    if phase <= TS:
        return Q_PEAK * math.sin(math.pi * phase / TS)
    return 0.0


def resistances(map_target: float = MAP_TARGET) -> tuple[float, float, float]:
    """Return (R_T, Rc, Rd) that place mean pressure at map_target."""
    r_t = map_target / Q_MEAN
    rc = RC_FRACTION * r_t
    rd = r_t - rc
    return r_t, rc, rd


def simulate(C: float, map_target: float = MAP_TARGET,
             dt: float = DT, n_cycles: int = N_CYCLES):
    """Run the discrete scheme. Return dict of the LAST cycle's results."""
    _, rc, rd = resistances(map_target)
    n = int(round(n_cycles * T / dt))
    pc = map_target                      # start near mean -> fast convergence
    last_start = int(round((n_cycles - 1) * T / dt))
    p_last: list[float] = []
    for i in range(n + 1):
        t = i * dt
        q = q_of_t(t)
        p = pc + rc * q
        if i >= last_start:
            p_last.append(p)
        # advance node pressure
        pc = pc + dt * (q / C - pc / (rd * C))
    systolic = max(p_last)
    diastolic = min(p_last)
    return {
        "Rc": rc, "Rd": rd, "R_T": rc + rd, "C": C,
        "systolic": systolic, "diastolic": diastolic,
        "PP": systolic - diastolic, "MAP": sum(p_last) / len(p_last),
    }


def tune_compliance(pp_target: float = PP_TARGET,
                    map_target: float = MAP_TARGET) -> float:
    """Bisection on C so pulse pressure hits pp_target. (PP decreases with C.)"""
    lo, hi = 0.3, 5.0
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        pp = simulate(mid, map_target)["PP"]
        if pp > pp_target:   # too stiff -> need more compliance
            lo = mid
        else:
            hi = mid
    return 0.5 * (lo + hi)


if __name__ == "__main__":
    r_t, rc, rd = resistances()
    print(f"T={T:.3f}s  Q_peak={Q_PEAK:.2f} mL/s  Q_mean={Q_MEAN:.4f} mL/s")
    print(f"R_T={r_t:.4f}  Rc={rc:.4f}  Rd={rd:.4f} mmHg*s/mL")
    guess = simulate(1.5)
    print(f"\nGuess C=1.5: MAP={guess['MAP']:.2f}  sys={guess['systolic']:.2f}"
          f"  dia={guess['diastolic']:.2f}  PP={guess['PP']:.2f}")
    c_tuned = tune_compliance()
    tuned = simulate(c_tuned)
    print(f"\nTuned C={c_tuned:.4f}: MAP={tuned['MAP']:.2f}  sys={tuned['systolic']:.2f}"
          f"  dia={tuned['diastolic']:.2f}  PP={tuned['PP']:.2f}")