"""Reproduce fr.aime labs' short-term, singly reinforced rectangular example.

EN 1992-1-1:2004, section 7.3.4. Units: mm, N, MPa.
No axial force, compression reinforcement, creep, shrinkage or imposed strain.
This reproduces the article arithmetic; it is not a general design tool.
References and scope: https://www.fraimelabs.com/blog/crack-width-calculation-eurocode-2
Run with Python 3; no external dependencies required.
"""
import math

B, H, D = 300.0, 500.0, 450.0
COVER, DIAMETER, BAR_COUNT = 40.0, 20.0, 4
ES, ECM, FCT_EFF = 200000.0, 33000.0, 2.9
KT, K1, K2, K3, K4 = 0.6, 0.8, 0.5, 3.4, 0.425


def calculate(moment_knm):
    area_s = BAR_COUNT * math.pi * DIAMETER**2 / 4
    alpha_e = ES / ECM
    x = (-alpha_e * area_s + math.sqrt((alpha_e * area_s)**2
         + 2 * B * alpha_e * area_s * D)) / B
    i_cr = B * x**3 / 3 + alpha_e * area_s * (D - x)**2
    moment = moment_knm * 1e6
    sigma_s = alpha_e * moment * (D - x) / i_cr
    # Independent equilibrium/lever-arm expression for the same steel stress.
    sigma_equilibrium = moment / (area_s * (D - x / 3))
    assert math.isclose(sigma_s, sigma_equilibrium, rel_tol=1e-12)
    h_eff = min(2.5 * (H - D), (H - x) / 3, H / 2)
    area_eff = B * h_eff
    rho_eff = area_s / area_eff
    centre_spacing = (B - 2 * (COVER + DIAMETER / 2)) / (BAR_COUNT - 1)
    assert centre_spacing <= 5 * (COVER + DIAMETER / 2)
    sr_max = K3 * COVER + K1 * K2 * K4 * DIAMETER / rho_eff
    raw_strain = (sigma_s - KT * FCT_EFF / rho_eff
                  * (1 + alpha_e * rho_eff)) / ES
    strain_floor = 0.6 * sigma_s / ES
    delta_strain = max(raw_strain, strain_floor)
    return dict(moment_knm=moment_knm, area_s_mm2=area_s,
                alpha_e=alpha_e, x_mm=x, i_cr_mm4=i_cr,
                steel_stress_mpa=sigma_s, h_eff_mm=h_eff,
                area_eff_mm2=area_eff, rho_eff=rho_eff,
                bar_centres_mm=centre_spacing, sr_max_mm=sr_max,
                raw_strain=raw_strain, strain_floor=strain_floor,
                delta_strain=delta_strain, wk_mm=sr_max * delta_strain)


if __name__ == "__main__":
    for moment_knm in (120.0, 170.0):
        print(f"\nShort-term illustrative case: M = {moment_knm:g} kN.m")
        for name, value in calculate(moment_knm).items():
            print(f"{name}: {value:.10g}")