"""TopEdu project 9: When the textbook beam formula breaks

A beam fixed at one end sags under its own weight. The textbook (Euler-Bernoulli) answer for the tip sag is
    w = q L^4 / (8 E I)      with q = weight per length, E the stiffness of the material, I = h^4 / 12 for a square.
It assumes a long, slender beam. A short beam also bends by shear, which the formula leaves out. Timoshenko's
theory adds it:  w = q L^4 / (8 E I) + q L^2 / (2 kappa G A),  with G = E / (2 (1 + nu)) and kappa = 5/6.

The full 3D answer needs the online lab (elasticity in 3D is too heavy for a starter). This file does the parts
you can do on a laptop:
Stage 1: solve the beam equation  E I w'''' = q  numerically and check it against the formula. CHECK should say PASS.
Then it compares both theories with the lab's 3D results that we measured when we tested this project.
Stage 2: measure real rulers, fit the L^4 law, and find where the formula stops working.

Run:  python 09_cantilever_beam.py  [--show]      or     python 09_cantilever_beam.py my_rulers.csv
Needs: numpy and matplotlib
"""
import csv
import os
import sys

import numpy as np
import matplotlib

if "--show" not in sys.argv:
    matplotlib.use("Agg")
import matplotlib.pyplot as plt

OUT = os.path.join(os.path.dirname(os.path.abspath(__file__)), "output")
os.makedirs(OUT, exist_ok=True)

NU = 0.3   # Poisson ratio

# Tip sag divided by the textbook formula, from the online lab's 3D elasticity solver (9 points across the section,
# nu = 0.3). These are the numbers in TopEdu's test report; rerun them in the lab with a finer grid before quoting.
LAB_3D = {14: 0.996, 10: 0.998, 5: 1.017, 3: 1.073}


def beam_sag(length=1.0, n=201, q=1.0, ei=1.0):
    """Finite differences for E I w'''' = q on 0 < x < L. Clamped at x = 0 (w = 0, w' = 0), free at x = L
    (no bending moment w'' = 0, no shear force w''' = 0). Returns x and w, with w measured downward."""
    x = np.linspace(0, length, n)
    h = x[1] - x[0]
    size = n + 4                       # two ghost points at each end carry the end conditions
    M = np.zeros((size, size)); rhs = np.zeros(size)
    for i in range(n):                 # the beam equation at every real point
        row = i + 2
        M[row, row - 2:row + 3] = np.array([1, -4, 6, -4, 1]) * ei / h ** 4
        rhs[row] = q
    M[0, 2] = 1                                    # w(0) = 0
    M[1, [1, 3]] = [-1, 1]                         # w'(0) = 0
    M[n + 2, [n, n + 1, n + 2]] = [1, -2, 1]       # w''(L) = 0
    M[n + 3, [n - 1, n, n + 2, n + 3]] = [-1, 2, -2, 1]   # w'''(L) = 0
    w = np.linalg.solve(M, rhs)
    return x, w[2:n + 2]


def timoshenko_ratio(slenderness, nu=NU, kappa=5 / 6):
    """Timoshenko tip sag divided by the textbook formula, for a square section: 1 + E/(3 kappa G) (h/L)^2."""
    return 1 + 2 * (1 + nu) / (3 * kappa) / slenderness ** 2


def read_rulers(path):
    """A CSV with columns length_cm and sag_mm (one row per overhang you measured)."""
    with open(path, newline="") as handle:
        rows = list(csv.DictReader(handle))
    return np.array([float(r["length_cm"]) for r in rows]), np.array([float(r["sag_mm"]) for r in rows])


def main():
    x, w = beam_sag()
    formula = 1 / 8                  # q L^4 / (8 E I) with q = L = E I = 1
    error = abs(w[-1] - formula) / formula
    print(f"beam equation: tip sag {w[-1]:.5f}, formula {formula:.5f}, error {100 * error:.3f}%")
    print(f"CHECK numerical beam against q L^4 / (8 E I) -> {'PASS' if error < 0.01 else 'FAIL'} (limit 1%)")

    print("tip sag / textbook formula, by length / depth:")
    print("  L/h   Timoshenko   lab 3D (our test)")
    for s in sorted(LAB_3D, reverse=True):
        print(f"  {s:3d}   {timoshenko_ratio(s):.3f}        {LAB_3D[s]:.3f}")
    s5 = np.sqrt(2 * (1 + NU) / (3 * 5 / 6) / 0.05)
    print(f"Timoshenko says the textbook formula is 5% low once L/h < {s5:.1f}; does the 3D solver agree?")

    fig, ax = plt.subplots(figsize=(5.5, 3.6))
    s = np.linspace(2, 15, 100)
    ax.plot(s, timoshenko_ratio(s), label="Timoshenko (adds shear)")
    ax.plot(list(LAB_3D), list(LAB_3D.values()), "o", label="lab, 3D elasticity (coarse grid)")
    ax.axhline(1, color="gray", lw=0.6); ax.axhline(1.05, color="gray", lw=0.6, ls="--")
    ax.set_xlabel("length / depth"); ax.set_ylabel("sag / textbook formula"); ax.legend()
    fig.tight_layout(); fig.savefig(os.path.join(OUT, "09_ratio.png"), dpi=150)

    rulers = [arg for arg in sys.argv[1:] if arg.endswith(".csv")]
    if rulers:
        lengths, sags = read_rulers(rulers[0])
        power, _ = np.polyfit(np.log(lengths), np.log(sags), 1)
        print(f"your rulers: sag grows like length^{power:.2f} (the formula says 4 for self-weight, 3 for a weight at the tip)")
    print(f"plots saved in {OUT}")

    # ---- YOUR TURN ----
    # 1. Rulers: clamp a ruler to a desk with 10, 15, 20, 25 and 30 cm hanging over, measure the tip sag, and save a
    #    CSV with columns length_cm,sag_mm. Run:  python 09_cantilever_beam.py my_rulers.csv
    #    Predict the power first. Coins at the tip change it: why?
    # 2. From your fit, work out E for the ruler (weigh it to get q). Compare with a published value for the material.
    # 3. In the online lab, rerun the 3D beam at L/h = 3 and 5 with a finer grid. How much do the numbers above move?
    #    The clamped end has a stress singularity, so report how your answer depends on the grid.
    # 4. Change nu to 0.05 and 0.4. Does the 5% threshold move? Check one value in the lab.
    if "--show" in sys.argv:
        plt.show()


if __name__ == "__main__":
    main()