"""Plot an original simply supported tensioned-beam model, not experimental data.

Writes vibroscope-bending-bias.png in the current working directory. The actual
clamped-grip vibroscope problem has different boundary conditions.
"""

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def main():
    for ratio in (0.0, 0.1, 1.0, 2.0):
        intercept = 4e6
        gradient = ratio * intercept
        first = np.sqrt(intercept + gradient)
        second = 2 * np.sqrt(intercept + 4 * gradient)
        np.testing.assert_allclose(
            (second**2 / 4 - first**2) / 3, gradient, atol=1e-8
        )
        np.testing.assert_allclose(
            (4 * first**2 - second**2 / 4) / 3, intercept
        )
        np.testing.assert_allclose(intercept / first**2, 1 / (1 + ratio))
        assert 2 <= second / first < 4

    plt.rcParams.update({"font.size": 11, "axes.titlesize": 11})
    fig, axes = plt.subplots(1, 2, figsize=(8.8, 4.4), dpi=100, facecolor="white")
    fig.subplots_adjust(left=0.08, right=0.975, top=0.78, bottom=0.17, wspace=0.33)
    fig.suptitle("Hinged fibre model: frequency and area bias", fontsize=15, y=0.97)
    fig.text(0.5, 0.88, r"Simply supported ends; $r=\pi^2 B/(T L^2)$; $B=EI$",
             ha="center", fontsize=12)
    ratio = np.linspace(0, 2, 601)
    area = 1 / (1 + ratio)
    modes = 2 * np.sqrt((1 + 4 * ratio) / (1 + ratio))
    axes[0].plot(ratio, area, color="#0072b2", lw=2.4)
    axes[0].set(xlabel="Bending-to-tension parameter r", ylabel=r"Area ratio $\hat A/A$",
                title="Uncorrected flexible-string estimate", xlim=(0, 2), ylim=(0.25, 1.05))
    axes[0].plot(0.1, 1 / 1.1, "o", color="#0072b2")
    axes[0].annotate("r = 0.10:\n9.09% below true area", xy=(0.1, 1 / 1.1),
                     xytext=(0.55, 0.88), fontsize=10,
                     arrowprops={"arrowstyle": "->", "color": "#555555"})
    axes[1].plot(ratio, modes, color="#d55e00", lw=2.4)
    axes[1].axhline(4, ls="--", color="#888888", lw=0.9)
    axes[1].set(xlabel="Bending-to-tension parameter r", ylabel=r"Mode ratio $f_2/f_1$",
                title="Diagnostic from two frequencies", xlim=(0, 2), ylim=(1.9, 4.15))
    axes[1].text(1.97, 3.94, "large-r limit", ha="right", va="top", fontsize=9)
    for ax in axes:
        ax.grid(alpha=0.23)
        ax.set_facecolor("white")
        ax.spines[["top", "right"]].set_visible(False)
    fig.savefig(Path.cwd() / "vibroscope-bending-bias.png", dpi=100,
                facecolor="white", transparent=False)
    plt.close(fig)


if __name__ == "__main__":
    main()