// diffeq.sy -- Ordinary and partial differential equation solvers
// as native Scyla. Euler, Heun, RK4, heat/wave/Laplace, and
// biological models (Michaelis-Menten, Hodgkin-Huxley, SIR, PK).
//
// Sources:
//   Burden RL, Faires JD. Numerical Analysis, 10th ed, Cengage, 2016.
//   Strogatz SH. Nonlinear Dynamics and Chaos, 2nd ed, CRC Press, 2015.
//   Keener J, Sneyd J. Mathematical Physiology, 2nd ed, Springer, 2009.
//   Murray JD. Mathematical Biology, 3rd ed, Springer, 2002.
//
// Marjorie McCubbins + Aether, April 2026
// Patent pending: 64/034,536

// =========================================================
// ODE SOLVERS -- Single Step Methods
// Since Scyla has no function pointers, user evaluates f(t,y)
// and passes the result. This is "call this, not go find this."
// =========================================================

// Forward Euler: y_{n+1} = y_n + dt * f(t_n, y_n)
// First-order, O(dt). The simplest ODE integrator.
// Euler L. Institutionum calculi integralis, 1768.
fn euler_step(y: Float, f_ty: Float, dt: Float) -> Float {
    return y + dt * f_ty
}

// Heun's method (improved Euler / explicit trapezoidal):
// Predictor: y_euler = y + dt * f(t, y)
// Corrector: y_{n+1} = y + (dt/2) * (f(t,y) + f(t+dt, y_euler))
// Second-order, O(dt^2). Burden 10th ed, Sec 5.2.
// f_ty = f at current point, f_ty_euler = f evaluated at Euler prediction
fn heun_step(y: Float, f_ty: Float, f_ty_euler: Float, dt: Float) -> Float {
    return y + (dt / 2.0) * (f_ty + f_ty_euler)
}

// Classical Runge-Kutta (RK4): fourth-order, O(dt^4).
// y_{n+1} = y_n + (dt/6) * (k1 + 2*k2 + 2*k3 + k4)
// User computes k1..k4 using the helper functions below.
// Runge C. Math Ann 46:167-178, 1895. Kutta W. Z Math Phys 46:435-453, 1901.
fn rk4_step(y: Float, k1: Float, k2: Float, k3: Float, k4: Float, dt: Float) -> Float {
    return y + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
}

// RK4 helpers: where to evaluate f for each stage
// k1 = f(t, y)                  -- evaluate at current point
// k2 = f(t + dt/2, y + dt*k1/2) -- evaluate at midpoint using k1
// k3 = f(t + dt/2, y + dt*k2/2) -- evaluate at midpoint using k2
// k4 = f(t + dt, y + dt*k3)     -- evaluate at endpoint using k3

fn rk4_k2_y(y: Float, k1: Float, dt: Float) -> Float {
    return y + dt * k1 / 2.0
}

fn rk4_k3_y(y: Float, k2: Float, dt: Float) -> Float {
    return y + dt * k2 / 2.0
}

fn rk4_k4_y(y: Float, k3: Float, dt: Float) -> Float {
    return y + dt * k3
}

fn rk4_t_mid(t: Float, dt: Float) -> Float {
    return t + dt / 2.0
}

fn rk4_t_end(t: Float, dt: Float) -> Float {
    return t + dt
}

// =========================================================
// SYSTEMS OF ODEs (2D)
// For dx/dt = f(t,x,y), dy/dt = g(t,x,y)
// =========================================================

// Euler step for x-component: x_{n+1} = x + dt * f
fn euler_step_2d_x(x: Float, f_val: Float, dt: Float) -> Float {
    return x + dt * f_val
}

// Euler step for y-component: y_{n+1} = y + dt * g
fn euler_step_2d_y(y: Float, g_val: Float, dt: Float) -> Float {
    return y + dt * g_val
}

// RK4 for either component of a 2D system (same formula)
fn rk4_step_2d(val: Float, k1: Float, k2: Float, k3: Float, k4: Float, dt: Float) -> Float {
    return val + (dt / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
}

// =========================================================
// STABILITY ANALYSIS
// =========================================================

// Forward Euler stability limit for dy/dt = lambda*y
// Stability region: |1 + dt*lambda| <= 1
// For real negative lambda: dt_max = 2/|lambda|
// Burden 10th ed, Sec 5.10.
fn euler_stability_limit(eigenvalue_real: Float) -> Float {
    if eigenvalue_real >= 0.0 { return 0.0 }
    return 0.0 - 2.0 / eigenvalue_real
}

// Stiffness ratio: |lambda_fast / lambda_slow|
// If ratio > 1000, the system is stiff and needs implicit methods.
// Strogatz 2nd ed, Ch 6.
fn stiffness_ratio(fast_eigenvalue: Float, slow_eigenvalue: Float) -> Float {
    return abs(fast_eigenvalue) / abs(slow_eigenvalue)
}

// =========================================================
// PDE METHODS -- Finite Differences
// =========================================================

// Heat equation (FTCS explicit): du/dt = alpha * d2u/dx2
// u_next = u_center + r * (u_left - 2*u_center + u_right)
// where r = alpha * dt / dx^2 (must be <= 0.5 for stability)
// Burden 10th ed, Sec 12.2.
fn heat_equation_explicit(
    u_left: Float,
    u_center: Float,
    u_right: Float,
    alpha: Float,
    dt: Float,
    dx: Float
) -> Float {
    let r = alpha * dt / (dx * dx)
    return u_center + r * (u_left - 2.0 * u_center + u_right)
}

// CFL number for heat equation: r = alpha * dt / dx^2
// Stable when r <= 0.5. Burden 10th ed, Sec 12.2.
fn heat_cfl(alpha: Float, dt: Float, dx: Float) -> Float {
    return alpha * dt / (dx * dx)
}

// Wave equation (explicit): d2u/dt2 = c^2 * d2u/dx2
// u_next = 2*u_center - u_prev + C^2 * (u_left - 2*u_center + u_right)
// where C = c*dt/dx (Courant number, must be <= 1)
// Burden 10th ed, Sec 12.4.
fn wave_equation_step(
    u_left: Float,
    u_center: Float,
    u_right: Float,
    u_prev_center: Float,
    c: Float,
    dt: Float,
    dx: Float
) -> Float {
    let courant_sq = (c * dt / dx) * (c * dt / dx)
    return 2.0 * u_center - u_prev_center + courant_sq * (u_left - 2.0 * u_center + u_right)
}

// Courant number for wave equation: C = c * dt / dx
// Stable when C <= 1. Burden 10th ed, Sec 12.4.
fn wave_cfl(c: Float, dt: Float, dx: Float) -> Float {
    return c * dt / dx
}

// Laplace equation (Jacobi iteration): d2u/dx2 + d2u/dy2 = 0
// u_new = (u_left + u_right + u_up + u_down) / 4
// Iterate until convergence. Burden 10th ed, Sec 12.1.
fn laplace_jacobi(u_left: Float, u_right: Float, u_up: Float, u_down: Float) -> Float {
    return (u_left + u_right + u_up + u_down) / 4.0
}

// Poisson equation (Jacobi): nabla^2(u) = f
// u_new = (u_left + u_right + u_up + u_down - dx^2 * f_center) / 4
fn poisson_jacobi(
    u_left: Float,
    u_right: Float,
    u_up: Float,
    u_down: Float,
    f_center: Float,
    dx: Float
) -> Float {
    return (u_left + u_right + u_up + u_down - dx * dx * f_center) / 4.0
}

// Gauss-Seidel for Laplace (uses already-updated neighbors)
// Same formula as Jacobi but converges ~2x faster
fn laplace_gauss_seidel(u_left_new: Float, u_right: Float, u_up_new: Float, u_down: Float) -> Float {
    return (u_left_new + u_right + u_up_new + u_down) / 4.0
}

// =========================================================
// BIOLOGICAL DIFFERENTIAL EQUATION MODELS
// =========================================================

// Logistic growth: dy/dt = r*y*(1 - y/K)
// r = intrinsic growth rate, K = carrying capacity
// Verhulst PF. Correspondance math et phys 10:113-121, 1838.
fn logistic_growth(y: Float, r: Float, k_capacity: Float) -> Float {
    return r * y * (1.0 - y / k_capacity)
}

// Logistic growth exact solution: y(t) = K / (1 + ((K-y0)/y0) * exp(-r*t))
fn logistic_exact(y0: Float, r: Float, k_capacity: Float, t: Float) -> Float {
    return k_capacity / (1.0 + ((k_capacity - y0) / y0) * exp(0.0 - r * t))
}

// Exponential decay rate: dy/dt = -lambda * y
fn exponential_decay(y: Float, lambda: Float) -> Float {
    return 0.0 - lambda * y
}

// Exponential decay exact: y(t) = y0 * exp(-lambda * t)
fn exponential_decay_exact(y0: Float, lambda: Float, t: Float) -> Float {
    return y0 * exp(0.0 - lambda * t)
}

// Half-life from decay constant: t_half = ln(2) / lambda
fn half_life(lambda: Float) -> Float {
    return 0.6931471805599453 / lambda
}

// Decay constant from half-life: lambda = ln(2) / t_half
fn decay_constant(t_half: Float) -> Float {
    return 0.6931471805599453 / t_half
}

// Michaelis-Menten enzyme kinetics: v = Vmax * S / (Km + S)
// Rate of product formation as a function of substrate concentration.
// Michaelis L, Menten ML. Biochem Z 49:333-369, 1913.
fn michaelis_menten(substrate: Float, vmax: Float, km: Float) -> Float {
    return vmax * substrate / (km + substrate)
}

// Hill equation: v = Vmax * S^n / (K^n + S^n)
// Cooperative binding. n > 1 = positive cooperativity.
// Hill AV. J Physiol 40:iv-vii, 1910.
fn hill_equation(substrate: Float, vmax: Float, k: Float, n: Float) -> Float {
    let s_n = pow(substrate, n)
    let k_n = pow(k, n)
    return vmax * s_n / (k_n + s_n)
}

// Hodgkin-Huxley gating variable: dx/dt = alpha*(1-x) - beta*x
// x = fraction of open gates (m, h, or n)
// alpha, beta = voltage-dependent rate constants
// Hodgkin AL, Huxley AF. J Physiol 117:500-544, 1952 (Nobel 1963).
fn hodgkin_huxley_gate(x: Float, alpha: Float, beta: Float) -> Float {
    return alpha * (1.0 - x) - beta * x
}

// HH steady-state: x_inf = alpha / (alpha + beta)
fn hh_steady_state(alpha: Float, beta: Float) -> Float {
    return alpha / (alpha + beta)
}

// HH time constant: tau = 1 / (alpha + beta)
fn hh_time_constant(alpha: Float, beta: Float) -> Float {
    return 1.0 / (alpha + beta)
}

// Lotka-Volterra prey: dx/dt = alpha*x - beta*x*y
// x = prey, y = predator, alpha = prey growth, beta = predation rate
// Lotka AJ. Elements of Physical Biology, 1925.
// Volterra V. Mem Accad Lincei 2:31-113, 1926.
fn lotka_volterra_prey(x: Float, y: Float, alpha: Float, beta: Float) -> Float {
    return alpha * x - beta * x * y
}

// Lotka-Volterra predator: dy/dt = delta*x*y - gamma*y
// delta = predator growth from predation, gamma = predator death
fn lotka_volterra_predator(x: Float, y: Float, delta: Float, gamma_r: Float) -> Float {
    return delta * x * y - gamma_r * y
}

// SIR epidemic model: dS/dt = -beta*S*I/N
// S = susceptible, I = infected, beta = transmission rate, N = total population
// Kermack WO, McKendrick AG. Proc R Soc Lond A 115:700-721, 1927.
fn sir_dS(s_pop: Float, i_pop: Float, beta: Float, n_pop: Float) -> Float {
    return 0.0 - beta * s_pop * i_pop / n_pop
}

// SIR: dI/dt = beta*S*I/N - gamma*I
fn sir_dI(s_pop: Float, i_pop: Float, beta: Float, gamma_r: Float, n_pop: Float) -> Float {
    return beta * s_pop * i_pop / n_pop - gamma_r * i_pop
}

// SIR: dR/dt = gamma*I
fn sir_dR(i_pop: Float, gamma_r: Float) -> Float {
    return gamma_r * i_pop
}

// Basic reproduction number R0 = beta / gamma
fn sir_r0(beta: Float, gamma_r: Float) -> Float {
    return beta / gamma_r
}

// One-compartment PK: Bateman equation
// C(t) = (Dose * ka) / (Vd * (ka - ke)) * (exp(-ke*t) - exp(-ka*t))
// ka = absorption rate, ke = elimination rate, Vd = volume of distribution
// Bateman H. Proc Cambridge Philos Soc 15:423-427, 1910.
fn pk_one_compartment(dose: Float, ka: Float, ke: Float, vd: Float, t: Float) -> Float {
    return (dose * ka) / (vd * (ka - ke)) * (exp(0.0 - ke * t) - exp(0.0 - ka * t))
}

// Time to peak concentration (Tmax) for one-compartment:
// Tmax = ln(ka/ke) / (ka - ke)
fn pk_tmax(ka: Float, ke: Float) -> Float {
    return ln(ka / ke) / (ka - ke)
}

// Steady-state concentration for constant infusion:
// Css = infusion_rate / clearance
fn pk_steady_state(infusion_rate: Float, clearance: Float) -> Float {
    return infusion_rate / clearance
}

// =========================================================
// ERROR ESTIMATION
// =========================================================

// Richardson extrapolation: combine two approximations at h and h/2
// Improved = (2^p * A(h/2) - A(h)) / (2^p - 1)
// p = order of the method. Burden 10th ed, Sec 4.2.
fn richardson_extrapolation(approx_h: Float, approx_h2: Float, p: Float) -> Float {
    let factor = pow(2.0, p)
    return (factor * approx_h2 - approx_h) / (factor - 1.0)
}

// Leading truncation error for Euler method: (dt/2) * |f''|
fn truncation_error_euler(d2y_dt2: Float, dt: Float) -> Float {
    return dt / 2.0 * abs(d2y_dt2)
}

// Leading truncation error for RK4: (dt^4/120) * |f^(5)|
fn truncation_error_rk4(d5y_dt5: Float, dt: Float) -> Float {
    return pow(dt, 4.0) / 120.0 * abs(d5y_dt5)
}

// =========================================================
// COUNTS
// =========================================================

fn diffeq_function_count() -> Int {
    return 42
}
