"""Original equal-material model, not measured fibre data; save to cwd.

Python 3.14.4, NumPy 2.3.5, Matplotlib 3.10.7.
Axial loading assumes equal extension of every filament. Bending assumes
independent neutral axes, free inter-filament sliding and equal curvature.
"""
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import numpy as np

area = 1.0  # square millimetres of material, excluding voids
young = 70_000.0  # N/mm^2, stipulated illustrative material
length = 100.0  # millimetres
radius = np.sqrt(area / np.pi)
counts = np.logspace(0, 3, 301)
small_radii = radius / np.sqrt(counts)
axial = counts * young * np.pi * small_radii**2 / length
bending = counts * young * np.pi * small_radii**4 / 4
assert np.allclose(axial, 700)
assert np.allclose(bending / bending[0], 1 / counts)
assert np.isclose(bending[0], young * area**2 / (4 * np.pi))

plt.rcParams.update({"font.size": 10, "figure.facecolor": "white"})
fig = plt.figure(figsize=(9.6, 4.8), dpi=100)
grid = fig.add_gridspec(1, 3, width_ratios=[1, 1, 1.8], left=0.055,
                       right=0.98, top=0.79, bottom=0.27, wspace=0.26)
solid, bundle, scaling = [fig.add_subplot(grid[0, i]) for i in range(3)]
fig.suptitle("Same material area: axial and bending responses scale differently", fontsize=14, y=0.96)
solid.add_patch(Circle((0, 0), radius, facecolor="#d9e9f3", edgecolor="#226c99", lw=1.4))
r = radius / 10
for row in range(10):
    for col in range(10):
        bundle.add_patch(Circle(((col - 4.5) * 2.15*r, (row - 4.5) * 2.15*r), r,
                                facecolor="#f3dfd0", edgecolor="#ae581e", lw=0.7))
for ax in [solid, bundle]:
    ax.set(xlim=(-0.7, 0.7), ylim=(-0.7, 0.7), xticks=[], yticks=[])
    ax.set_aspect("equal")
    ax.spines[:].set_visible(False)
solid.set_title("One solid circular rod", fontsize=10)
bundle.set_title("100 sliding filaments", fontsize=10)
solid.set_xlabel("Material area: 1 mm²")
bundle.set_xlabel("Total material area: 1 mm²")
scaling.loglog(counts, axial/axial[0], color="#226c99", lw=2,
               label="Axial stiffness: 1")
scaling.loglog(counts, bending/bending[0], color="#ae581e", lw=2,
               label="Independent bending: 1/N")
scaling.scatter([100], [0.01], color="#ae581e", s=28, zorder=4)
scaling.set(xlabel="Number of filaments N", ylabel="Ratio to the solid rod",
            xlim=(1, 1000), ylim=(0.0008, 1.7))
scaling.legend(loc="lower left", frameon=False, fontsize=9)
scaling.grid(which="major", alpha=0.2)
scaling.spines[["top", "right"]].set_visible(False)
fig.text(0.5, 0.13, "E, material area and length fixed. Bending model permits sliding and neglects friction.",
         ha="center", fontsize=10)
fig.text(0.5, 0.065, "Drawings exclude voids from material area; the two load cases use their stated constraints.",
         ha="center", fontsize=10)
fig.savefig(Path("unbonded-fibre-bundle-scaling.png"), facecolor="white")
plt.close(fig)