"""
SDOF receptance — closed form vs modal superposition
====================================================

A single-DOF mass-spring-damper has the analytical receptance

.. math::

    H(\\omega) \\;=\\; \\frac{1}{k - m\\,\\omega^{2} + i\\,\\omega\\,c}

This example drives the same SDOF system through
:func:`~femorph_solver.solvers.harmonic.solve_harmonic` and
overlays the modal-superposition output on the closed form.  With a
single retained mode, the superposition reproduces the closed form to
machine precision — the curves overlap exactly.

Tuning the example
------------------

* The eigenvalue ``ω_n = √(k / m)`` sets the resonance peak.
* The damping ratio ``ζ = c / (2 √(k m))`` sets the height.
* Reducing damping makes the peak taller; the closed-form and modal
  curves still match.

References
----------

* Ewins, D. J.  *Modal Testing: Theory, Practice and Application*,
  2nd ed., Research Studies Press (2000), §2.3 — receptance
  formulation used throughout this example.
* Craig, R. R. and Kurdila, A. J.  *Fundamentals of Structural
  Dynamics*, 2nd ed., Wiley (2006), §11.4 — modal superposition for
  forced harmonic response.
"""

from __future__ import annotations

import matplotlib.pyplot as plt
import numpy as np
import scipy.sparse as sp

from femorph_solver.solvers.harmonic import solve_harmonic

# %%
# 1-DOF mass-spring-damper
# ------------------------
m = 1.0  # kg
k = 4.0e4  # N/m → ω_n = 200 rad/s, f_n ≈ 31.83 Hz
zeta = 0.02
c = 2.0 * zeta * np.sqrt(k * m)

omega_n = np.sqrt(k / m)
f_n_hz = omega_n / (2.0 * np.pi)
print(f"Resonance: ω_n = {omega_n:.2f} rad/s   f_n = {f_n_hz:.3f} Hz   ζ = {zeta:.3f}")

K = sp.csr_array(np.array([[k]], dtype=np.float64))
M = sp.csr_array(np.array([[m]], dtype=np.float64))

# %%
# Excitation sweep
# ----------------
# Span 0.1 ω_n through 3 ω_n on a log grid so the resonance peak is
# well-resolved without piling samples at the high-frequency tail.
freq = np.geomspace(0.1 * f_n_hz, 3.0 * f_n_hz, 400)

# %%
# Modal-superposition response
# ----------------------------
F = np.array([1.0], dtype=np.complex128)  # unit forcing on the single DOF
res = solve_harmonic(
    K,
    M,
    F,
    frequencies=freq,
    n_modes=1,
    modal_damping_ratio=zeta,
)

# %%
# Closed-form receptance
# ----------------------
omega = 2.0 * np.pi * freq
H_closed = 1.0 / (k - m * omega**2 + 1j * omega * c)

# Sanity-check: with one retained mode the two should match to floating
# point precision.
np.testing.assert_allclose(res.displacement[:, 0], H_closed, atol=1e-12, rtol=1e-12)
print(
    "Modal-superposition matches closed form to "
    f"max |Δ| = {np.max(np.abs(res.displacement[:, 0] - H_closed)):.2e}"
)

# %%
# Overlay the magnitude / phase
# -----------------------------
fig, (ax_mag, ax_phase) = plt.subplots(2, 1, sharex=True, figsize=(8, 6))

ax_mag.loglog(freq, np.abs(H_closed), "k-", lw=2.5, label="closed form")
ax_mag.loglog(freq, np.abs(res.displacement[:, 0]), "C1--", lw=1.5, label="modal superposition")
ax_mag.axvline(f_n_hz, color="grey", linestyle=":", alpha=0.6, label=f"f_n = {f_n_hz:.2f} Hz")
ax_mag.set_ylabel("|H(ω)|  [m / N]")
ax_mag.set_title("SDOF receptance")
ax_mag.grid(True, which="both", alpha=0.3)
ax_mag.legend()

phase = np.angle(res.displacement[:, 0], deg=True)
phase_closed = np.angle(H_closed, deg=True)
ax_phase.semilogx(freq, phase_closed, "k-", lw=2.5, label="closed form")
ax_phase.semilogx(freq, phase, "C1--", lw=1.5, label="modal superposition")
ax_phase.axvline(f_n_hz, color="grey", linestyle=":", alpha=0.6)
ax_phase.set_xlabel("frequency [Hz]")
ax_phase.set_ylabel("phase [deg]")
ax_phase.grid(True, which="both", alpha=0.3)
ax_phase.legend()

fig.tight_layout()
plt.show()