"""Original symmetric four-point loading and its elementary beam response."""
from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np


def main():
    span, offset, load = .160, .040, 40.0
    width, height, young = .010, .004, 70e9
    inertia = width * height**3 / 12
    moment = load * offset / 2
    curvature = moment / (young * inertia)
    assert np.isclose(moment, .8)
    assert np.isclose(moment * height / (2*inertia), 30e6)
    assert np.isclose(curvature * height / 2, 3e-4 / .7)

    fig = plt.figure(figsize=(10, 4.6), dpi=100, layout="constrained")
    fig.set_facecolor("white")
    grid = fig.add_gridspec(2, 2, width_ratios=[1.4, 1], height_ratios=[1, 1])
    fixture = fig.add_subplot(grid[0, 0])
    bending = fig.add_subplot(grid[1, 0])
    strain = fig.add_subplot(grid[:, 1])
    fixture.plot([0, 160], [0, 0], color="#555555", linewidth=5)
    for x in [0, 160]:
        fixture.scatter(x, -1.1, marker="^", color="#555555", s=90)
        fixture.annotate("", (x, -.8), (x, -3.6),
                         arrowprops=dict(arrowstyle="->", color="#1565c0", lw=2))
        fixture.text(x, -4.3, "20 N", ha="center", color="#1565c0")
    for x in [40, 120]:
        fixture.annotate("", (x, .4), (x, 4),
                         arrowprops=dict(arrowstyle="->", color="#b44e17", lw=2))
        fixture.text(x, 4.3, "20 N", ha="center", color="#b44e17")
    fixture.set(xlim=(-15, 175), ylim=(-6.3, 6), xticks=[0, 40, 120, 160],
                yticks=[], title="Four-point fixture: total load 40 N")
    fixture.set_xlabel("Position along beam (mm)")
    for side in ["top", "right", "left"]:
        fixture.spines[side].set_visible(False)

    x = np.linspace(0, span, 501)
    m = np.minimum(np.minimum(x, span-x), offset) * load/2
    assert np.allclose(m[(x > offset) & (x < span-offset)], moment)
    bending.plot(x*1e3, m, color="#1565c0", lw=2)
    bending.fill_between(x*1e3, 0, m, color="#1565c0", alpha=.12)
    bending.axvspan(40, 120, color="#d5e8ce", alpha=.4)
    bending.text(80, .43, "Pure bending\nV = 0", ha="center")
    bending.set(xlabel="Position along beam (mm)", ylabel="Moment (N m)",
                ylim=(0, 1), xlim=(0, 160), title="Constant central moment: 0.800 N m")
    bending.grid(alpha=.18)

    y = np.linspace(-height/2, height/2, 101)
    epsilon = -curvature*y
    strain.plot(epsilon*1e6, y*1e3, color="#b44e17", lw=2)
    strain.axhline(0, color="#777777", ls="--", lw=1)
    strain.axvline(0, color="#bbbbbb", lw=1)
    strain.text(0, .15, "Neutral axis", ha="center", fontsize=10)
    strain.scatter(epsilon[[0, -1]]*1e6, y[[0, -1]]*1e3, color="#b44e17")
    strain.set(xlabel="Longitudinal strain (microstrain)",
               ylabel="Height y from neutral axis (mm)",
               title="Central-section strain\nE = 70 GPa; b = 10 mm; h = 4 mm",
               xlim=(-550, 550), ylim=(-2.5, 2.5))
    strain.grid(alpha=.18)
    fig.savefig(Path(__file__).with_suffix(".png").name, facecolor="white")
    plt.close(fig)


if __name__ == "__main__":
    main()