#!/usr/bin/env python3
"""
Limit-state design of a singly reinforced rectangular RC beam to IS 456:2000.

Simply supported or cantilever beam under uniformly distributed load.
Steps follow the usual hand-calculation sequence:

    1. Trial section            cl 23.2.1 (span/depth ratios)
    2. Effective span           cl 22.2
    3. Loads                    self weight + imposed + finishes, factored x 1.5 (Table 18)
    4. Moment and shear         wL^2/8, wL/2   or   wL^2/2, wL
    5. Depth check              Mu,lim = k fck b d^2     (Annex G, G-1.1)
    6. Tension steel            Ast from Annex G-1.1 (b), min/max cl 26.5.1.1
    7. Shear                    Table 19, Table 20, cl 40.4, cl 26.5.1.5 / 26.5.1.6
    8. Deflection               cl 23.2.1 with Fig 4 modification factor

Units: mm, N, N/mm^2 internally; kN/m and kN.m at the interface.

    python3 beam_design.py            interactive
    python3 beam_design.py --demo     built-in example
    python3 beam_design.py --json     print the result dict as JSON (with --demo or piped input)
    python3 beam_design.py --png out.png   save the detailing drawing instead of showing it

Rithika A, 2023 (original competition entry); revised 2026.
"""
from __future__ import annotations

import argparse
import json
import math
import sys

# ---------------------------------------------------------------- code data

# Annex G-1.1 (c): Mu,lim = k fck b d^2 ; xu,max/d
K_LIM = {250: 0.148, 415: 0.138, 500: 0.133}
XU_MAX = {250: 0.53, 415: 0.48, 500: 0.46}

# Table 19: design shear strength of concrete tau_c (N/mm^2) vs 100 Ast/(b d)
PT_TABLE = [0.15, 0.25, 0.50, 0.75, 1.00, 1.25, 1.50, 1.75, 2.00, 2.25, 2.50, 2.75, 3.00]
TAU_C = {
    20: [0.28, 0.36, 0.48, 0.56, 0.62, 0.67, 0.72, 0.75, 0.79, 0.81, 0.82, 0.82, 0.82],
    25: [0.29, 0.36, 0.49, 0.57, 0.64, 0.70, 0.74, 0.78, 0.82, 0.85, 0.88, 0.90, 0.92],
    30: [0.29, 0.37, 0.50, 0.59, 0.66, 0.71, 0.76, 0.80, 0.84, 0.88, 0.91, 0.94, 0.96],
    35: [0.29, 0.37, 0.50, 0.59, 0.67, 0.73, 0.78, 0.82, 0.86, 0.90, 0.93, 0.96, 0.99],
    40: [0.30, 0.38, 0.51, 0.60, 0.68, 0.74, 0.79, 0.84, 0.88, 0.92, 0.95, 0.98, 1.01],
}
# Table 20: maximum shear stress tau_c,max (N/mm^2)
TAU_C_MAX = {20: 2.8, 25: 3.1, 30: 3.5, 35: 3.7, 40: 4.0}

BASIC_L_D = {"ss": 20.0, "cantilever": 7.0}   # cl 23.2.1 (a)
GAMMA_F = 1.5                                  # Table 18, DL + LL
UNIT_WEIGHT = 25.0                             # kN/m^3, reinforced concrete


def interp(x: float, xs: list[float], ys: list[float]) -> float:
    """Linear interpolation, clamped to the table ends."""
    if x <= xs[0]:
        return ys[0]
    if x >= xs[-1]:
        return ys[-1]
    for i in range(1, len(xs)):
        if x <= xs[i]:
            t = (x - xs[i - 1]) / (xs[i] - xs[i - 1])
            return ys[i - 1] + t * (ys[i] - ys[i - 1])
    return ys[-1]


def bar_area(dia: float) -> float:
    return math.pi * dia * dia / 4.0


# ---------------------------------------------------------------- design steps

def suggest_depth(beam_type: str, clear_span: float) -> float:
    """Trial effective depth from the basic span/depth ratio, rounded up to 25 mm."""
    d = clear_span / BASIC_L_D[beam_type]
    return math.ceil(d / 25.0) * 25.0


def effective_span(beam_type: str, clear_span: float, d: float, support_width: float) -> float:
    """cl 22.2 (a) and (c)."""
    if beam_type == "ss":
        return min(clear_span + d, clear_span + support_width)
    return clear_span + 0.5 * d


def loads(b: float, D: float, imposed: float, finishes: float) -> dict:
    """kN/m. Self weight from the gross section, no hidden allowances."""
    self_weight = b * D * 1e-6 * UNIT_WEIGHT
    w = self_weight + imposed + finishes
    return {"self_weight": self_weight, "w": w, "wu": GAMMA_F * w}


def moment_shear(beam_type: str, wu: float, leff: float) -> tuple[float, float]:
    """Returns (Mu in N.mm, Vu in N). wu in kN/m == N/mm, leff in mm."""
    if beam_type == "ss":
        return wu * leff ** 2 / 8.0, wu * leff / 2.0
    return wu * leff ** 2 / 2.0, wu * leff


def flexure(Mu: float, b: float, d: float, D: float, fck: float, fy: float, bar_dia: float) -> dict:
    k = K_LIM[int(fy)]
    Mu_lim = k * fck * b * d * d
    d_req = math.sqrt(Mu / (k * fck * b))
    out = {"k": k, "Mu_lim": Mu_lim, "d_req": d_req, "singly": Mu <= Mu_lim}
    if not out["singly"]:
        out["Ast_req"] = None
        return out
    # Annex G-1.1 (b)
    Ast = 0.5 * fck / fy * (1.0 - math.sqrt(1.0 - 4.6 * Mu / (fck * b * d * d))) * b * d
    Ast_min = 0.85 * b * d / fy            # cl 26.5.1.1 (a)
    Ast_max = 0.04 * b * D                 # cl 26.5.1.1 (b)
    Ast_design = max(Ast, Ast_min)
    a1 = bar_area(bar_dia)
    n = math.ceil(Ast_design / a1)
    Ast_prov = n * a1
    # rough fit check in one layer: side cover 25 + 8 mm stirrup each side, clear gap >= max(dia, 25)  cl 26.3.2
    gap = max(bar_dia, 25.0)
    width_needed = 2 * (25.0 + 8.0) + n * bar_dia + (n - 1) * gap
    out.update({
        "Ast_req": Ast, "Ast_min": Ast_min, "Ast_max": Ast_max, "Ast_design": Ast_design,
        "bar_dia": bar_dia, "n_bars": n, "Ast_prov": Ast_prov,
        "pt": 100.0 * Ast_prov / (b * d),
        "fits_one_layer": width_needed <= b, "width_needed": width_needed,
        "steel_ok": Ast_prov <= Ast_max,
    })
    return out


def shear(Vu: float, b: float, d: float, fck: float, fy: float, pt: float, stirrup_dia: float, legs: int = 2) -> dict:
    tau_v = Vu / (b * d)
    tau_c = interp(min(max(pt, 0.15), 3.0), PT_TABLE, TAU_C[int(fck)])
    tau_c_max = TAU_C_MAX[int(fck)]
    fy_s = min(fy, 415.0)                  # cl 40.4: stirrup fy not taken > 415
    Asv = legs * bar_area(stirrup_dia)
    out = {"tau_v": tau_v, "tau_c": tau_c, "tau_c_max": tau_c_max, "fy_s": fy_s, "Asv": Asv,
           "section_ok": tau_v <= tau_c_max}
    if not out["section_ok"]:
        out.update({"regime": "redesign", "Vus": None, "Sv_calc": None, "Sv_max": None, "Sv_prov": None})
        return out
    if tau_v <= tau_c:
        Vus = 0.0
        Sv_calc = 0.87 * fy_s * Asv / (0.4 * b)              # cl 26.5.1.6 minimum stirrups
        regime = "minimum"
    else:
        Vus = Vu - tau_c * b * d
        Sv_calc = 0.87 * fy_s * Asv * d / Vus                 # cl 40.4 (a)
        regime = "designed"
    Sv_max = min(0.75 * d, 300.0)                             # cl 26.5.1.5
    Sv_prov = max(50.0, math.floor(min(Sv_calc, Sv_max) / 10.0) * 10.0)
    out.update({"regime": regime, "Vus": Vus, "Sv_calc": Sv_calc, "Sv_max": Sv_max, "Sv_prov": Sv_prov})
    return out


def deflection(beam_type: str, leff: float, d: float, fy: float, Ast_req: float, Ast_prov: float, pt: float) -> dict:
    basic = BASIC_L_D[beam_type]
    if beam_type == "ss" and leff > 10000.0:
        basic *= 10000.0 / leff                                # cl 23.2.1 (b)
    fs = 0.58 * fy * Ast_req / Ast_prov                        # Fig 4 note
    # closed form of Fig 4 (SP 24 / Pillai & Menon), capped at 2.0
    kt = 1.0 / (0.225 + 0.00322 * fs + 0.625 * math.log10(pt))
    kt = min(max(kt, 0.0), 2.0)
    allowed = basic * kt
    actual = leff / d
    return {"basic": basic, "fs": fs, "kt": kt, "ld_allowed": allowed, "ld_actual": actual, "ok": actual <= allowed}


def design(inp: dict) -> dict:
    """Run the full design. Returns a flat dict of every number and check used on the calc sheet."""
    beam_type = inp["beam_type"]
    b = float(inp["b"]); l = float(inp["clear_span"]); d = float(inp["d"])
    cover = float(inp.get("cover", 50)); D = d + cover
    support = float(inp.get("support_width", b))
    fck = float(inp["fck"]); fy = float(inp["fy"])
    bar_dia = float(inp.get("bar_dia", 20)); st_dia = float(inp.get("stirrup_dia", 8))
    if int(fy) not in K_LIM:
        raise ValueError(f"fy must be one of {sorted(K_LIM)}")
    if int(fck) not in TAU_C:
        raise ValueError(f"fck must be one of {sorted(TAU_C)}")

    r = {"beam_type": beam_type, "b": b, "clear_span": l, "d": d, "cover": cover, "D": D,
         "support_width": support, "fck": fck, "fy": fy, "stirrup_dia": st_dia,
         "d_suggested": suggest_depth(beam_type, l)}
    r["leff"] = effective_span(beam_type, l, d, support)
    r.update(loads(b, D, float(inp["imposed"]), float(inp.get("finishes", 0))))
    r["imposed"] = float(inp["imposed"]); r["finishes"] = float(inp.get("finishes", 0))
    Mu, Vu = moment_shear(beam_type, r["wu"], r["leff"])
    r["Mu"] = Mu; r["Vu"] = Vu
    fl = flexure(Mu, b, d, D, fck, fy, bar_dia)
    r.update(fl)
    r["depth_ok"] = d >= fl["d_req"]
    if not fl["singly"]:
        r["status"] = "doubly"
        return r
    sh = shear(Vu, b, d, fck, fy, fl["pt"], st_dia)
    r.update(sh)
    if not sh["section_ok"]:
        r["status"] = "shear_fail"
        return r
    de = deflection(beam_type, r["leff"], d, fy, fl["Ast_req"], fl["Ast_prov"], fl["pt"])
    r.update({"basic_ld": de["basic"], "fs": de["fs"], "kt": de["kt"],
              "ld_allowed": de["ld_allowed"], "ld_actual": de["ld_actual"], "deflection_ok": de["ok"]})
    r["status"] = "ok" if (r["depth_ok"] and fl["steel_ok"] and de["ok"]) else "check"
    return r


# ---------------------------------------------------------------- reporting

def fmt(x, nd=1):
    return "-" if x is None else f"{x:,.{nd}f}"


def report(r: dict) -> str:
    kind = "Simply supported" if r["beam_type"] == "ss" else "Cantilever"
    L = []
    P = L.append
    P(f"{kind} beam  |  M{int(r['fck'])}  Fe{int(r['fy'])}")
    P("=" * 64)
    P("1. SECTION  (cl 23.2.1)")
    P(f"   b = {fmt(r['b'],0)} mm   d = {fmt(r['d'],0)} mm   cover = {fmt(r['cover'],0)} mm   D = {fmt(r['D'],0)} mm")
    P(f"   trial d from span/{int(BASIC_L_D[r['beam_type']])} = {fmt(r['d_suggested'],0)} mm")
    P("2. EFFECTIVE SPAN  (cl 22.2)")
    P(f"   Leff = {fmt(r['leff'],0)} mm")
    P("3. LOADS")
    P(f"   self weight = {fmt(r['self_weight'],2)} kN/m   imposed = {fmt(r['imposed'],2)}   finishes = {fmt(r['finishes'],2)}")
    P(f"   w = {fmt(r['w'],2)} kN/m   wu = 1.5 w = {fmt(r['wu'],2)} kN/m")
    P("4. MOMENT AND SHEAR")
    P(f"   Mu = {fmt(r['Mu']/1e6,2)} kN.m   Vu = {fmt(r['Vu']/1e3,2)} kN")
    P("5. DEPTH CHECK  (Annex G-1.1)")
    P(f"   Mu,lim = {r['k']} fck b d^2 = {fmt(r['Mu_lim']/1e6,2)} kN.m   d,req = {fmt(r['d_req'],1)} mm   "
      f"{'OK' if r['depth_ok'] else 'INCREASE d'}")
    if r["status"] == "doubly":
        P("   Mu > Mu,lim: section would be over-reinforced. Increase depth or design as doubly reinforced.")
        return "\n".join(L)
    P("6. TENSION STEEL  (Annex G-1.1, cl 26.5.1.1)")
    P(f"   Ast,req = {fmt(r['Ast_req'],0)} mm^2   Ast,min = {fmt(r['Ast_min'],0)}   Ast,max = {fmt(r['Ast_max'],0)}")
    P(f"   provide {r['n_bars']} no. {int(r['bar_dia'])} mm bars = {fmt(r['Ast_prov'],0)} mm^2   pt = {fmt(r['pt'],2)} %"
      + ("" if r["fits_one_layer"] else f"   (needs {fmt(r['width_needed'],0)} mm width: use two layers or wider beam)"))
    P("7. SHEAR  (Tables 19, 20; cl 40.4; cl 26.5.1.5/6)")
    P(f"   tau_v = {fmt(r['tau_v'],3)}   tau_c = {fmt(r['tau_c'],3)}   tau_c,max = {fmt(r['tau_c_max'],2)} N/mm^2")
    if r["status"] == "shear_fail":
        P("   tau_v > tau_c,max: section inadequate in shear. Increase b or d.")
        return "\n".join(L)
    if r["regime"] == "minimum":
        P(f"   tau_v <= tau_c: minimum stirrups, Sv = 0.87 fy Asv / (0.4 b) = {fmt(r['Sv_calc'],0)} mm")
    else:
        P(f"   Vus = Vu - tau_c b d = {fmt(r['Vus']/1e3,2)} kN   Sv = 0.87 fy Asv d / Vus = {fmt(r['Sv_calc'],0)} mm")
    P(f"   Sv,max = min(0.75 d, 300) = {fmt(r['Sv_max'],0)} mm   ->  2-legged {int(r['stirrup_dia'])} mm stirrups @ {fmt(r['Sv_prov'],0)} mm c/c")
    P("8. DEFLECTION  (cl 23.2.1, Fig 4)")
    P(f"   fs = {fmt(r['fs'],1)} N/mm^2   kt = {fmt(r['kt'],2)}   (L/d)allowed = {fmt(r['basic_ld'],0)} x kt = {fmt(r['ld_allowed'],2)}   "
      f"(L/d)actual = {fmt(r['ld_actual'],2)}   {'OK' if r['deflection_ok'] else 'NOT OK'}")
    P("=" * 64)
    P("RESULT: " + {"ok": "section is adequate", "check": "review the items marked above"}[r["status"]])
    return "\n".join(L)


# ---------------------------------------------------------------- drawing

def draw(r: dict, png: str | None = None) -> None:
    """Cross-section and elevation drawn from the design result."""
    import matplotlib
    if png:
        matplotlib.use("Agg")
    import matplotlib.pyplot as plt
    from matplotlib.patches import Rectangle, Circle

    b, D, d = r["b"], r["D"], r["d"]
    n, dia = r.get("n_bars", 0), r.get("bar_dia", 20)
    sv, sd = r.get("Sv_prov") or 0, r["stirrup_dia"]
    L = r["leff"]
    c = 25.0  # side cover to stirrup

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), gridspec_kw={"width_ratios": [1, 3]})
    # ---- section
    ax1.add_patch(Rectangle((0, 0), b, D, fill=False, lw=2))
    ax1.add_patch(Rectangle((c, c), b - 2 * c, D - 2 * c, fill=False, lw=1.2, ls="--"))
    if n:
        xs = [c + sd + dia / 2] if n == 1 else [c + sd + dia / 2 + i * (b - 2 * c - 2 * sd - dia) / (n - 1) for i in range(n)]
        for x in xs:
            ax1.add_patch(Circle((x, c + sd + dia / 2), dia / 2, color="C1"))
    for x in (c + sd + 6, b - c - sd - 6):
        ax1.add_patch(Circle((x, D - c - sd - 6), 6, color="C0"))
    ax1.set_xlim(-60, b + 60); ax1.set_ylim(-60, D + 60); ax1.set_aspect("equal"); ax1.axis("off")
    ax1.set_title(f"Section {int(b)} x {int(D)}\n{n}-{int(dia)}Ø bottom, 2-12Ø top, {int(sd)}Ø @ {int(sv)}")
    # ---- elevation
    ax2.add_patch(Rectangle((0, 0), L, D, fill=False, lw=2))
    ax2.plot([c, L - c], [c + sd + dia / 2] * 2, color="C1", lw=3)
    ax2.plot([c, L - c], [D - c - sd - 6] * 2, color="C0", lw=2)
    if sv:
        x = sv / 2
        while x < L:
            ax2.plot([x, x], [c, D - c], color="0.4", lw=1)
            x += sv
    if r["beam_type"] == "ss":
        ax2.plot([0], [-15], marker="^", ms=12, color="k"); ax2.plot([L], [-15], marker="^", ms=12, color="k")
    else:
        ax2.add_patch(Rectangle((-80, -60), 80, D + 120, hatch="///", fill=False))
    ax2.set_xlim(-120, L + 120); ax2.set_ylim(-120, D + 150); ax2.set_aspect("equal"); ax2.axis("off")
    ax2.set_title(f"Elevation, Leff = {int(L)} mm, stirrups {int(sd)} mm @ {int(sv)} c/c")
    fig.suptitle("Beam detailing (from design)")
    fig.tight_layout()
    if png:
        fig.savefig(png, dpi=150)
    else:
        plt.show()


# ---------------------------------------------------------------- cli

DEMO = {"beam_type": "ss", "b": 230, "clear_span": 5000, "support_width": 230, "imposed": 10,
        "finishes": 0, "fck": 20, "fy": 500, "d": 450, "cover": 50, "bar_dia": 20, "stirrup_dia": 8}


def ask(prompt, cast=float, default=None):
    s = input(f"{prompt}{'' if default is None else f' [{default}]'}: ").strip()
    if s == "" and default is not None:
        return default
    return cast(s)


def interactive() -> dict:
    t = ask("Beam type (ss / cantilever)", str, "ss").lower()
    beam_type = "cantilever" if t.startswith("c") else "ss"
    l = ask("Clear span (mm)")
    inp = {
        "beam_type": beam_type, "clear_span": l,
        "b": ask("Beam width b (mm)", float, 230),
        "support_width": ask("Support / wall width (mm)", float, 230),
        "imposed": ask("Imposed load (kN/m)"),
        "finishes": ask("Finishes / partitions (kN/m)", float, 0),
        "fck": ask("Concrete grade fck (20/25/30/35/40)", int, 20),
        "fy": ask("Steel grade fy (250/415/500)", int, 500),
    }
    inp["d"] = ask("Effective depth d (mm)", float, suggest_depth(beam_type, l))
    inp["cover"] = ask("Effective cover (mm)", float, 50)
    inp["bar_dia"] = ask("Main bar diameter (mm)", float, 20)
    inp["stirrup_dia"] = ask("Stirrup diameter (mm)", float, 8)
    return inp


def main(argv=None):
    ap = argparse.ArgumentParser(description=__doc__.split("\n")[1])
    ap.add_argument("--demo", action="store_true", help="run the built-in example")
    ap.add_argument("--json", action="store_true", help="print result as JSON")
    ap.add_argument("--input", help="JSON file (or '-' for stdin) with the input dict")
    ap.add_argument("--png", help="save the drawing to this file")
    ap.add_argument("--no-draw", action="store_true")
    a = ap.parse_args(argv)

    if a.demo:
        inp = dict(DEMO)
    elif a.input:
        inp = json.load(sys.stdin if a.input == "-" else open(a.input))
    else:
        inp = interactive()
    r = design(inp)
    if a.json:
        print(json.dumps(r, indent=1))
    else:
        print(report(r))
    if not a.no_draw and not a.json:
        draw(r, a.png)


if __name__ == "__main__":
    main()