#!/usr/bin/env python3
"""Companion to Engineering Foundations v1.0, DevinaHQ, 2026-09-20.

Run with Python 3: python3 oscillator-check.py
Standard library only. Prints results; does not create files.

Hypothetical parameters, NOT experimental data:
    m = 1 kg; k = 100 N/m; c = 0 or 2 N s/m
    x(0) = 0.01 m; v(0) = 0 m/s; external force = 0
Model: m*x'' + c*x' + k*x = 0, displacement from static equilibrium.
Compares RK4 against an analytic underdamped solution over 0 <= t <= 5 s.
Analytic agreement verifies computation; it does not validate a physical model.
"""

import math
import platform

MASS = 1.0
STIFFNESS = 100.0
INITIAL_X = 0.01
INITIAL_V = 0.0
DURATION = 5.0


def exact(t, damping):
    alpha = damping / (2 * MASS)
    omega_sq = STIFFNESS / MASS - alpha**2
    if omega_sq <= 0:
        raise ValueError("This analytic solution requires an underdamped case.")
    omega = math.sqrt(omega_sq)
    a = INITIAL_X
    b = (INITIAL_V + alpha * a) / omega
    cos_t, sin_t = math.cos(omega * t), math.sin(omega * t)
    decay = math.exp(-alpha * t)
    x = decay * (a * cos_t + b * sin_t)
    v = decay * ((-alpha * a + omega * b) * cos_t
                 + (-alpha * b - omega * a) * sin_t)
    return x, v


def derivative(state, damping):
    x, v = state
    return v, -(damping * v + STIFFNESS * x) / MASS


def shifted(state, slope, step):
    return tuple(a + step * b for a, b in zip(state, slope))


def rk4(state, step, damping):
    k1 = derivative(state, damping)
    k2 = derivative(shifted(state, k1, step / 2), damping)
    k3 = derivative(shifted(state, k2, step / 2), damping)
    k4 = derivative(shifted(state, k3, step), damping)
    return tuple(a + step * (b + 2*c + 2*d + e) / 6
                 for a, b, c, d, e in zip(state, k1, k2, k3, k4))


def energy(state):
    x, v = state
    return 0.5 * MASS * v**2 + 0.5 * STIFFNESS * x**2


def compare(step, damping):
    count = round(DURATION / step)
    if not math.isclose(count * step, DURATION):
        raise ValueError("Time step must divide the simulation duration.")
    state = (INITIAL_X, INITIAL_V)
    initial_energy = energy(state)
    max_x_error = 0.0
    max_energy_error = 0.0
    for i in range(count + 1):
        expected = exact(i * step, damping)
        max_x_error = max(max_x_error, abs(state[0] - expected[0]))
        max_energy_error = max(max_energy_error,
                               abs(energy(state) - energy(expected)))
        if i < count:
            state = rk4(state, step, damping)
    return max_x_error, max_energy_error / initial_energy, energy(state)


def main():
    print("Engineering Foundations v1.0 — hypothetical simulation")
    print("Python:", platform.python_version())
    print("No experimental data; analytical comparison only.")
    for damping in (0.0, 2.0):
        omega_n = math.sqrt(STIFFNESS / MASS)
        zeta = damping / (2 * math.sqrt(STIFFNESS * MASS))
        omega_d = omega_n * math.sqrt(1 - zeta**2)
        print("\nc = %.1f N s/m, zeta = %.2f, period = %.6f s"
              % (damping, zeta, 2 * math.pi / omega_d))
        print("dt_s    max_x_error_m   max_energy_error/E0   final_energy_J")
        errors = []
        for step in (0.02, 0.01, 0.005):
            x_error, relative_energy_error, final_energy = compare(step, damping)
            errors.append(x_error)
            print("%.3f   %.8e   %.8e        %.8e"
                  % (step, x_error, relative_energy_error, final_energy))
        if not errors[0] > errors[1] > errors[2]:
            raise RuntimeError("Expected decreasing displacement errors.")
        if errors[-1] >= 1e-7:
            raise RuntimeError("Fine-step displacement error exceeded 1e-7 m.")
        print("PASS: errors decrease; fine-step error below 1e-7 m.")
    print("\nNext step for validation: independently collected measurements.")


if __name__ == "__main__":
    main()