{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import pyvista\npyvista.OFF_SCREEN = True\npyvista.set_jupyter_backend('static')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# SDOF receptance \u2014 closed form vs modal superposition\n\nA single-DOF mass-spring-damper has the analytical receptance\n\n\\begin{align}H(\\omega) \\;=\\; \\frac{1}{k - m\\,\\omega^{2} + i\\,\\omega\\,c}\\end{align}\n\nThis example drives the same SDOF system through\n:func:`~femorph_solver.solvers.harmonic.solve_harmonic` and\noverlays the modal-superposition output on the closed form.  With a\nsingle retained mode, the superposition reproduces the closed form to\nmachine precision \u2014 the curves overlap exactly.\n\n## Tuning the example\n\n* The eigenvalue ``\u03c9_n = \u221a(k / m)`` sets the resonance peak.\n* The damping ratio ``\u03b6 = c / (2 \u221a(k m))`` sets the height.\n* Reducing damping makes the peak taller; the closed-form and modal\n  curves still match.\n\n## References\n\n* Ewins, D. J.  *Modal Testing: Theory, Practice and Application*,\n  2nd ed., Research Studies Press (2000), \u00a72.3 \u2014 receptance\n  formulation used throughout this example.\n* Craig, R. R. and Kurdila, A. J.  *Fundamentals of Structural\n  Dynamics*, 2nd ed., Wiley (2006), \u00a711.4 \u2014 modal superposition for\n  forced harmonic response.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from __future__ import annotations\n\nimport matplotlib.pyplot as plt\nimport numpy as np\nimport scipy.sparse as sp\n\nfrom femorph_solver.solvers.harmonic import solve_harmonic"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## 1-DOF mass-spring-damper\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "m = 1.0  # kg\nk = 4.0e4  # N/m \u2192 \u03c9_n = 200 rad/s, f_n \u2248 31.83 Hz\nzeta = 0.02\nc = 2.0 * zeta * np.sqrt(k * m)\n\nomega_n = np.sqrt(k / m)\nf_n_hz = omega_n / (2.0 * np.pi)\nprint(f\"Resonance: \u03c9_n = {omega_n:.2f} rad/s   f_n = {f_n_hz:.3f} Hz   \u03b6 = {zeta:.3f}\")\n\nK = sp.csr_array(np.array([[k]], dtype=np.float64))\nM = sp.csr_array(np.array([[m]], dtype=np.float64))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Excitation sweep\nSpan 0.1 \u03c9_n through 3 \u03c9_n on a log grid so the resonance peak is\nwell-resolved without piling samples at the high-frequency tail.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "freq = np.geomspace(0.1 * f_n_hz, 3.0 * f_n_hz, 400)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Modal-superposition response\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "F = np.array([1.0], dtype=np.complex128)  # unit forcing on the single DOF\nres = solve_harmonic(\n    K,\n    M,\n    F,\n    frequencies=freq,\n    n_modes=1,\n    modal_damping_ratio=zeta,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Closed-form receptance\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "omega = 2.0 * np.pi * freq\nH_closed = 1.0 / (k - m * omega**2 + 1j * omega * c)\n\n# Sanity-check: with one retained mode the two should match to floating\n# point precision.\nnp.testing.assert_allclose(res.displacement[:, 0], H_closed, atol=1e-12, rtol=1e-12)\nprint(\n    \"Modal-superposition matches closed form to \"\n    f\"max |\u0394| = {np.max(np.abs(res.displacement[:, 0] - H_closed)):.2e}\"\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Overlay the magnitude / phase\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, (ax_mag, ax_phase) = plt.subplots(2, 1, sharex=True, figsize=(8, 6))\n\nax_mag.loglog(freq, np.abs(H_closed), \"k-\", lw=2.5, label=\"closed form\")\nax_mag.loglog(freq, np.abs(res.displacement[:, 0]), \"C1--\", lw=1.5, label=\"modal superposition\")\nax_mag.axvline(f_n_hz, color=\"grey\", linestyle=\":\", alpha=0.6, label=f\"f_n = {f_n_hz:.2f} Hz\")\nax_mag.set_ylabel(\"|H(\u03c9)|  [m / N]\")\nax_mag.set_title(\"SDOF receptance\")\nax_mag.grid(True, which=\"both\", alpha=0.3)\nax_mag.legend()\n\nphase = np.angle(res.displacement[:, 0], deg=True)\nphase_closed = np.angle(H_closed, deg=True)\nax_phase.semilogx(freq, phase_closed, \"k-\", lw=2.5, label=\"closed form\")\nax_phase.semilogx(freq, phase, \"C1--\", lw=1.5, label=\"modal superposition\")\nax_phase.axvline(f_n_hz, color=\"grey\", linestyle=\":\", alpha=0.6)\nax_phase.set_xlabel(\"frequency [Hz]\")\nax_phase.set_ylabel(\"phase [deg]\")\nax_phase.grid(True, which=\"both\", alpha=0.3)\nax_phase.legend()\n\nfig.tight_layout()\nplt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}