const state = {
    m: 1.0,
    k: 100,
    c: 5.0,
    F0: 10,
    p: 5.0,
    u0: 0.0,
    v0: 0.0,
    tMax: 10.0,
    dt: 0.001
};

function derivatives(t, u, v) {
    const u_prime = v;
    const v_prime = (state.F0 * Math.sin(state.p * t) - state.c * v - state.k * u) / state.m;
    return { u_prime, v_prime };
}

function solve() {
    let t = 0;
    let u = state.u0;
    let v = state.v0;
    let maxU = 0;

    const steps = Math.ceil(state.tMax / state.dt);
    for (let i = 0; i <= steps; i++) {
        const k1 = derivatives(t, u, v);
        const t2 = t + state.dt / 2;
        const u2 = u + k1.u_prime * state.dt / 2;
        const v2 = v + k1.v_prime * state.dt / 2;
        const k2 = derivatives(t2, u2, v2);
        const t3 = t + state.dt / 2;
        const u3 = u + k2.u_prime * state.dt / 2;
        const v3 = v + k2.v_prime * state.dt / 2;
        const k3 = derivatives(t3, u3, v3);
        const t4 = t + state.dt;
        const u4 = u + k3.u_prime * state.dt;
        const v4 = v + k3.v_prime * state.dt;
        const k4 = derivatives(t4, u4, v4);

        u += (state.dt / 6) * (k1.u_prime + 2 * k2.u_prime + 2 * k3.u_prime + k4.u_prime);
        v += (state.dt / 6) * (k1.v_prime + 2 * k2.v_prime + 2 * k3.v_prime + k4.v_prime);
        t += state.dt;

        if (Math.abs(u) > maxU) maxU = Math.abs(u);
    }
    return { finalU: u, maxU: maxU };
}

function runTest(name, params) {
    Object.assign(state, params);
    const res = solve();
    console.log(`Test: ${name}`);
    console.log(`  Params: ${JSON.stringify(params)}`);
    console.log(`  Result: Max |u| = ${res.maxU.toFixed(5)}, Final u = ${res.finalU.toFixed(5)}`);
    return res;
}

console.log("--- Verification Suite ---");

// 1. Zero Values
runTest("Zero Values", { m: 1, k: 100, c: 5, F0: 0, p: 0, u0: 0, v0: 0 });

// 2. Resonance (omega_n = 10)
// Low damping to see large amplitude
runTest("Resonance", { m: 1, k: 100, c: 1, F0: 10, p: 10, u0: 0, v0: 0 });

// 3. Underdamped Free Vibration (zeta = 0.25)
runTest("Underdamped Free", { m: 1, k: 100, c: 5, F0: 0, p: 0, u0: 1, v0: 0 });

// 4. Critically Damped Free Vibration (zeta = 1.0)
// c = 2 * sqrt(1 * 100) = 20
runTest("Critically Damped Free", { m: 1, k: 100, c: 20, F0: 0, p: 0, u0: 1, v0: 0 });

// 5. Overdamped Free Vibration (zeta = 2.0)
// c = 40
runTest("Overdamped Free", { m: 1, k: 100, c: 40, F0: 0, p: 0, u0: 1, v0: 0 });