import handcalcs.render          # registers the %%render cell magic
from handcalcs.render import render
import forallpeople as si
si.environment('structural', top_level=True)   # injects mm, m, kN, MPa, etc. into global scope
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib.gridspec import GridSpec

%%render params
# Section geometry
b_f = 600 * mm     # flange width
b_w = 250 * mm     # web width
h_f = 100 * mm     # flange thickness
h   = 450 * mm     # total section depth
h_w = h - h_f      # web height
c_nom = 35 * mm    # nominal cover to rebar
phi   = 20 * mm    # assumed rebar diameter
d     = h - c_nom - phi / 2    # effective depth to rebar centroid

%%render params
# Material properties
f_ck = 25 * MPa    # characteristic concrete compressive strength (cylinder)
f_yk = 500 * MPa   # characteristic steel yield strength
E_c  = 31 * GPa    # concrete modulus of elasticity (EC2 approx for C25)
E_s  = 200 * GPa   # steel modulus of elasticity

%%render
L   = 6.0 * m           # span
q   = 15.0 * kN / m     # uniformly distributed load
M_Ed = (q * L**2) / 8   # design bending moment at mid-span

%%render long
# Areas of sub-regions
A_f = b_f * h_f                # flange area
A_w = b_w * h_w                # web area
A_total = A_f + A_w            # total cross-section area

%%render
# Centroid distances from top fibre
y_f = h_f / 2                          # centroid of flange from top
y_w = h_f + h_w / 2                    # centroid of web from top
y_bar = (A_f * y_f + A_w * y_w) / A_total    # section centroid from top

%%render
# Second moment of area — parallel axis theorem for each part
I_f = (b_f * h_f**3) / 12 + A_f * (y_bar - y_f)**2   # flange
I_w = (b_w * h_w**3) / 12 + A_w * (y_bar - y_w)**2   # web
I   = I_f + I_w                                        # total

%%render
# Distance from neutral axis to extreme fibres
y_top = y_bar                  # top fibre to centroid
y_bot = h - y_bar              # bottom fibre to centroid

%%render symbolic
# Show the formula structure only — no numeric values yet
sigma_top_sym = M_Ed * y_top / I
sigma_bot_sym = M_Ed * y_bot / I
epsilon_top_sym = sigma_top_sym / E_c
epsilon_bot_sym = sigma_bot_sym / E_c

# @handcalc decorator — the programmatic API for using handcalcs outside of a magic cell
# Useful when you want to call a calculation as a reusable function
from handcalcs import handcalc
from IPython.display import Latex, display

@handcalc(jupyter_display=True)
def bending_stress(M, y, I):
    sigma = M * y / I
    return locals()

# Call it like a regular function — it prints the rendered LaTeX AND returns the result dict
results = bending_stress(M_Ed, y_bar, I)
print("sigma =", results["sigma"])

sigma = 4.511 MPa

%%render
# Bending stresses at extreme fibres
sigma_top = M_Ed * y_top / I    # compressive stress at top (positive = compression here)
sigma_bot = M_Ed * y_bot / I    # tensile stress at bottom

%%render
# Elastic strains in the uncracked section (illustration — not ULS design strains)
epsilon_top = sigma_top / E_c   # strain at top fibre (concrete in compression)
epsilon_bot = sigma_bot / E_c   # strain at bottom fibre (concrete in tension, uncracked)

%%render
# EC2 design strengths — partial factors gamma_c = 1.5 (concrete), gamma_s = 1.15 (steel)
gamma_c = 1.5
gamma_s = 1.15
f_cd = f_ck / gamma_c      # design compressive strength of concrete
f_yd = f_yk / gamma_s      # design yield strength of steel

%%render
# Moment capacity parameter K
# Using the web width b_w (conservative for a T-beam — see note above).
# A rigorous check would first locate the neutral axis: if x <= h_f the flange
# is in compression alone and b_f could be used, giving a smaller K.
K     = M_Ed / (b_w * d**2 * f_cd)
K_bal = 0.167    # balanced limit, EC2, no moment redistribution, fck <= 50 MPa

# Check: K must be <= K_bal for singly reinforced section
K_val = float(K)   # strip units for comparison (K is dimensionless)
print(f"K = {K_val:.4f}  |  K_bal = {K_bal}  |  {'OK — singly reinforced' if K_val <= K_bal else 'FAIL — need compression steel'}")

K = 0.0988  |  K_bal = 0.167  |  OK — singly reinforced

%%render
# Lever arm z (EC2 Eq. 6.18, capped at 0.95d)
z_calc = d * (0.5 + (0.25 - K_val / 1.134)**0.5)
z      = min(z_calc, 0.95 * d)

%%render
# Required tension reinforcement area
A_s_req = M_Ed / (f_yd * z)

%%render
# Minimum reinforcement check (EC2 Cl. 9.2.1.1)
f_ctm    = 0.3 * (f_ck / MPa)**(2/3) * MPa    # mean tensile strength
A_s_min  = max(0.26 * (f_ctm / f_yk) * b_w * d,  0.0013 * b_w * d)

# Final required area: governing value
A_s_design = max(A_s_req, A_s_min)
print(f"A_s,req  = {float(A_s_req * 1e6):.0f} mm²")
print(f"A_s,min  = {float(A_s_min * 1e6):.0f} mm²")
print(f"A_s,use  = {float(A_s_design * 1e6):.0f} mm²")

# Typical bar choice: 3 x H25 = 1473 mm², 4 x H20 = 1257 mm², 5 x H16 = 1005 mm²
A_s_prov_options = {
    "3 x H25": 3 * np.pi * (25e-3)**2 / 4,
    "4 x H20": 4 * np.pi * (20e-3)**2 / 4,
    "5 x H16": 5 * np.pi * (16e-3)**2 / 4,
}
print("\nBar options:")
for label, area in A_s_prov_options.items():
    status = "OK" if area * 1e6 >= float(A_s_req * 1e6) else "insufficient"
    print(f"  {label}: {area*1e6:.0f} mm²  — {status}")

A_s,req  = 424 mm²
A_s,min  = 135 mm²
A_s,use  = 424 mm²

Bar options:
  3 x H25: 1473 mm²  — OK
  4 x H20: 1257 mm²  — OK
  5 x H16: 1005 mm²  — OK

# Wong palette + shared numeric extracts for all plots
BLUE   = '#0072B2'
ORANGE = '#E69F00'
GREEN  = '#009E73'
RED    = '#D55E00'
GRAY   = '#666666'
PINK   = '#CC79A7'

import matplotlib
matplotlib.rcParams.update({'font.family': 'sans-serif', 'font.size': 9,
                             'axes.linewidth': 0.6, 'lines.linewidth': 1.4})

bf    = float(b_f / mm)
bw    = float(b_w / mm)
hf    = float(h_f / mm)
htot  = float(h / mm)
hw    = htot - hf
yb    = float(y_bar / mm)
d_val = float(d / mm)
s_top = float(sigma_top / MPa)
s_bot = float(sigma_bot / MPa)
A_req = float(A_s_req * 1e6)
A_min_val = float(A_s_min * 1e6)
A_prov_vals = {k: v * 1e6 for k, v in A_s_prov_options.items()}

fig1, ax1 = plt.subplots(figsize=(4.5, 5.5))
ax1.set_aspect('equal')
ax1.set_xlim(-60, bf + 90)
ax1.set_ylim(-35, htot + 50)
ax1.axis('off')
ax1.set_title('T-Section Geometry', fontsize=10, pad=10)

flange_x = (bf - bw) / 2
ax1.add_patch(patches.Rectangle((0, htot - hf), bf, hf, lw=0.8, edgecolor=GRAY, facecolor='#d1d5db'))
ax1.add_patch(patches.Rectangle((flange_x, 0), bw, hw, lw=0.8, edgecolor=GRAY, facecolor='#d1d5db'))

ax1.axhline(htot - yb, color=BLUE, lw=1.0, ls='--')
rebar_y = htot - d_val
for xi in [flange_x + 30, flange_x + bw / 2, flange_x + bw - 30]:
    ax1.plot(xi, rebar_y, 'o', color=RED, ms=6, zorder=5)

def dim_h(ax, y1, y2, x, top_label, bot_label):
    ax.annotate('', xy=(x, y2), xytext=(x, y1),
                arrowprops=dict(arrowstyle='<->', color=GRAY, lw=0.8))
    ax.text(x + 8, (y1 + y2) / 2, top_label, va='center', ha='left', fontsize=7.5, color=GRAY)

def dim_w(ax, x1, x2, y, label):
    ax.annotate('', xy=(x2, y), xytext=(x1, y),
                arrowprops=dict(arrowstyle='<->', color=GRAY, lw=0.8))
    ax.text((x1 + x2) / 2, y + 8, label, va='bottom', ha='center', fontsize=7.5, color=GRAY)

dim_h(ax1, 0, htot - hf, bf + 15, 'h_w = ' + str(int(hw)) + ' mm', '')
dim_h(ax1, htot - hf, htot, bf + 15, 'h_f = ' + str(int(hf)) + ' mm', '')
dim_w(ax1, 0, bf, htot + 22, 'b_f = ' + str(int(bf)) + ' mm')
dim_w(ax1, flange_x, flange_x + bw, -25, 'b_w = ' + str(int(bw)) + ' mm')
ax1.text(bf + 16, htot - yb + 4, 'y = ' + str(int(yb)) + ' mm (NA)', fontsize=7.5, color=BLUE)
ax1.text(flange_x + bw / 2 + 10, rebar_y, 'As', fontsize=9, color=RED, va='center', fontstyle='italic')

plt.savefig('fig1_section.png', dpi=150, bbox_inches='tight')
plt.show()

fig2, ax2 = plt.subplots(figsize=(4, 5))
ax2.set_title('Elastic Stress Distribution', fontsize=10, pad=10)

depths   = np.linspace(0, htot, 300)
stresses = np.interp(depths, [0, htot], [s_top, -s_bot])
comp_mask = stresses >= 0

ax2.fill_betweenx(depths, 0, stresses, where=comp_mask,  color=BLUE, alpha=0.3, label='Compression')
ax2.fill_betweenx(depths, 0, stresses, where=~comp_mask, color=RED,  alpha=0.3, label='Tension')
ax2.plot(stresses, depths, color=GRAY, lw=1.4)
ax2.axvline(0, color=GRAY, lw=0.6)
ax2.axhline(yb, color=BLUE, lw=0.8, ls='--', label='NA (y = ' + str(int(yb)) + ' mm)')

ax2.set_xlabel('Stress (MPa)', fontsize=9)
ax2.set_ylabel('Depth from top (mm)', fontsize=9)
ax2.set_ylim(htot, 0)
ax2.legend(fontsize=8, loc='lower left')
ax2.tick_params(labelsize=8)

ax2.annotate(str(round(s_top, 2)) + ' MPa', xy=(s_top, 5), xytext=(s_top - 2.0, 50),
             fontsize=7.5, color=BLUE, arrowprops=dict(arrowstyle='->', color=BLUE, lw=0.7))
ax2.annotate(str(round(s_bot, 2)) + ' MPa', xy=(-s_bot, htot - 5), xytext=(-s_bot - 2.5, htot - 55),
             fontsize=7.5, color=RED, arrowprops=dict(arrowstyle='->', color=RED, lw=0.7))

fig2.tight_layout()
plt.savefig('fig2_stress.png', dpi=150, bbox_inches='tight')
plt.show()

fig3, ax3 = plt.subplots(figsize=(5, 4))
ax3.set_title('Reinforcement Area', fontsize=10, pad=10)

bar_labels = ['As,req', 'As,min'] + list(A_prov_vals.keys())
bar_vals   = [A_req, A_min_val] + list(A_prov_vals.values())
bar_colors = [ORANGE, PINK] + [GREEN if v >= A_req else RED for v in A_prov_vals.values()]

bars = ax3.bar(bar_labels, bar_vals, color=bar_colors, edgecolor=GRAY, lw=0.6, width=0.55)
ax3.axhline(A_req, color=ORANGE, lw=1.0, ls='--', label='Required: ' + str(int(A_req)) + ' mm2')

for bar, val in zip(bars, bar_vals):
    ax3.text(bar.get_x() + bar.get_width() / 2, val + 12,
             str(int(val)), ha='center', va='bottom', fontsize=8)

ax3.set_ylabel('Area (mm2)', fontsize=9)
ax3.set_ylim(0, max(bar_vals) * 1.18)
ax3.tick_params(labelsize=8)
ax3.legend(fontsize=8)

fig3.tight_layout()
plt.savefig('fig3_rebar.png', dpi=150, bbox_inches='tight')
plt.show()

Library	What it does
🔢 handcalcs	Renders Python math as LaTeX. Write code, get calc sheet.
📐 forallpeople	SI units through every operation. Wrong units = instant error.

Magic	What it does
`%%render`	Default. Shows: symbol = expanded formula = numeric result
`%%render params`	Parameters only. No symbolic expansion — just name and value.
`%%render long`	Forces every expression onto its own line, even short ones.
`%%render symbolic`	Shows only the symbolic form — no numbers at all.
`render(expr)`	Inline function — render a single expression anywhere in a cell.

Comment	Effect
`# Label` at end of line	handcalcs uses your comment as the description label
Works on any line	Keeps rendered output readable without variable-name guesswork

Parameter	Value
Span $L$	6.0 m
Uniform load $q$	15 kN/m
Design moment $M_{Ed} = qL^2/8$	67.5 kN·m
Section centroid $\bar{y}$ (from top)	~184 mm
Second moment of area $I$	~$2.74 \times 10^9$ mm⁴
Effective depth $d$	405 mm
$f_{cd} = f_{ck}/\gamma_c$	16.67 MPa
$f_{yd} = f_{yk}/\gamma_s$	434.8 MPa
$K = M_{Ed}/(b_w d^2 f_{cd})$	0.099 < 0.167 ✅ singly reinforced
Lever arm $z$	~366 mm
$A_{s,req}$	424 mm²
$A_{s,min}$ (EC2 9.2.1.1)	135 mm²
Governing bar choice	3×H25 / 4×H20 / 5×H16 — all OK ✅

Feature	Where used
`%%render params`	Parameter tables (geometry, materials)
`%%render`	All intermediate calculations
`%%render long`	Section areas — forces each step to its own block
`%%render symbolic`	Stress formulas — formula only, no numbers
`render()`	Inline $K$ preview before formal calculation

🧮 Engineering Calculations in Python: No Excel Required¶

What we will calculate¶

⚙️ Setup¶

🧰 The Tools — Full Feature Tour¶

`handcalcs` — render modes¶

`forallpeople` — how units work¶

1. 🏗️ The Problem¶

2. 📐 Geometry and Material Properties¶

3. ⬇️ Loading and Bending Moment¶

4. 📏 Section Properties¶

5. 💪 Elastic Bending Stress and Strain (Uncracked Section)¶

6. 🔩 Required Tension Reinforcement (EC2 Simplified, ULS)¶

7. 📊 Plots¶

📐 Plot 1 — Section Geometry¶

💪 Plot 2 — Elastic Stress Distribution¶

🔩 Plot 3 — Required Reinforcement¶

8. 🏁 Summary¶

🧰 `handcalcs` features used in this notebook¶

📐 `forallpeople` features used¶

🤔 Was this worth it vs. Excel?¶

Feature	Example
Unit objects	`600 * mm`, `25 * MPa`, `15 * kN/m`
Dimensional arithmetic	`(kN/m) × m²` → `kN·m` auto-simplified
Dimensionless result	`Pa / Pa` → plain float, no stripping needed
`float(x / unit)`	Extract number for comparisons and numpy

🧮 Engineering Calculations in Python: No Excel Required¶

What we will calculate¶

⚙️ Setup¶

🧰 The Tools — Full Feature Tour¶

handcalcs — render modes¶

forallpeople — how units work¶

1. 🏗️ The Problem¶

2. 📐 Geometry and Material Properties¶

3. ⬇️ Loading and Bending Moment¶

4. 📏 Section Properties¶

5. 💪 Elastic Bending Stress and Strain (Uncracked Section)¶

6. 🔩 Required Tension Reinforcement (EC2 Simplified, ULS)¶

7. 📊 Plots¶

📐 Plot 1 — Section Geometry¶

💪 Plot 2 — Elastic Stress Distribution¶

🔩 Plot 3 — Required Reinforcement¶

8. 🏁 Summary¶

🧰 handcalcs features used in this notebook¶

📐 forallpeople features used¶

🤔 Was this worth it vs. Excel?¶

`handcalcs` — render modes¶

`forallpeople` — how units work¶

🧰 `handcalcs` features used in this notebook¶

📐 `forallpeople` features used¶