calcs/wood-stud-column/calc.py

from __future__ import annotations

import json
import math
import sys
import argparse
from pathlib import Path

try:
    import yaml
except ImportError:
    raise SystemExit("Install dependencies: python -m pip install -r requirements.txt")

try:
    from pint import DimensionalityError, UndefinedUnitError, UnitRegistry
except ImportError:
    raise SystemExit("Install dependencies: python -m pip install -r requirements.txt")


HERE = Path(__file__).resolve().parent

ureg = UnitRegistry()
ureg.define("kip = 1000 * force_pound")
ureg.define("ksi = kip / inch ** 2")
ureg.define("psi = force_pound / inch ** 2")


def quantity(value, unit: str, name: str) -> float:
    try:
        q = ureg.Quantity(value).to(unit)
    except (DimensionalityError, UndefinedUnitError, TypeError, ValueError) as exc:
        raise ValueError(f"{name}: expected {unit}, got {value!r}") from exc
    magnitude = float(q.magnitude)
    if magnitude <= 0:
        raise ValueError(f"{name} must be positive")
    return magnitude


def factor(value, name: str) -> float:
    magnitude = float(value)
    if magnitude <= 0:
        raise ValueError(f"{name} must be positive")
    return magnitude


# NDS 2018 constants for sawn-lumber columns (NDS 3.7.1.5, Eqs. 3.7-1 and 3.7-2).
KC_EULER = 0.822  # Euler buckling coefficient for a rectangular section about its
                 # weak axis (pi^2 / 12); appears directly in Eq. 3.7-2.
C_SAWN = 0.8      # Column-stability denominator coefficient c for sawn lumber
                 # (Eq. 3.7-1). c = 0.85 for round timber, 0.9 for glulam/SCL/CLT.


def _column_stability_factor(Fc_star: float, Emin_prime: float, Le: float, d: float):
    """NDS 3.7.1.5, Eq. 3.7-1 — column stability factor Cp for one buckling axis."""
    FcE = KC_EULER * Emin_prime / (Le / d) ** 2
    ratio = FcE / Fc_star
    Cp = (1.0 + ratio) / (2.0 * C_SAWN) - math.sqrt(((1.0 + ratio) / (2.0 * C_SAWN)) ** 2 - ratio / C_SAWN)
    if not math.isfinite(Cp) or Cp <= 0:
        raise ValueError(f"Cp = {Cp:.4f} is not finite/positive; column stability invalid (NDS 3.7.1.4)")
    return min(Cp, 1.0), FcE, ratio


def compute(inp: dict) -> dict:
    # Geometry
    b_in = quantity(inp["width"], "in", "width")
    d_in = quantity(inp["depth"], "in", "depth")
    L_in = quantity(inp["length"], "in", "length")
    Ke = factor(inp["Ke"], "Ke")

    # Material (NDS Supplement reference values)
    Fc_psi = quantity(inp["Fc"], "psi", "Fc")
    E_psi = quantity(inp["E"], "psi", "E")
    Emin_psi = quantity(inp["E_min"], "psi", "E_min")
    species_grade = str(inp["species_grade"])

    # Adjustment factors (NDS Table 4.3.1)
    CD = factor(inp["CD"], "CD")
    CM = factor(inp["CM"], "CM")
    Ct = factor(inp["Ct"], "Ct")
    CF = factor(inp["CF"], "CF")
    Ci = factor(inp["Ci"], "Ci")

    # Checked axial demand (ASD)
    P_lbf = quantity(inp["P"], "lbf", "P")

    A_in2 = b_in * d_in
    Le_in = Ke * L_in
    d_least = min(b_in, d_in)
    d_great = max(b_in, d_in)

    # Adjusted reference compression (without Cp), NDS 3.7.1.5
    Fc_star = Fc_psi * CD * CM * Ct * CF * Ci
    Emin_prime = Emin_psi * CM * Ct * Ci

    # Column stability factor about each axis; the weak axis governs.
    Cp_weak, FcE_weak, ratio_weak = _column_stability_factor(Fc_star, Emin_prime, Le_in, d_least)
    Cp_strong, FcE_strong, ratio_strong = _column_stability_factor(Fc_star, Emin_prime, Le_in, d_great)
    governing_axis = "weak" if Cp_weak <= Cp_strong else "strong"
    Cp = min(Cp_weak, Cp_strong)
    FcE = FcE_weak if governing_axis == "weak" else FcE_strong

    # Adjusted compression design value, NDS Table 4.3.1
    Fc_prime = Fc_star * Cp
    fc_psi = P_lbf / A_in2
    Pa_lbf = Fc_prime * A_in2

    def q(value: float) -> float:
        return round(value, 6)

    values = {
        "P_lbf": q(P_lbf),
        "L_in": q(L_in), "b_in": q(b_in), "d_in": q(d_in),
        "A_in2": q(A_in2), "Ke": q(Ke), "Le_in": q(Le_in),
        "d_least_in": q(d_least), "d_great_in": q(d_great),
        "Le_d_weak": q(Le_in / d_least), "Le_d_strong": q(Le_in / d_great),
        "species_grade": species_grade,
        "Fc_psi": q(Fc_psi), "E_psi": q(E_psi), "Emin_psi": q(Emin_psi),
        "CD": q(CD), "CM": q(CM), "Ct": q(Ct), "CF": q(CF), "Ci": q(Ci),
        "kc_euler": q(KC_EULER), "c_sawn": q(C_SAWN),
        "Fc_star_psi": q(Fc_star), "Emin_prime_psi": q(Emin_prime),
        "FcE_weak_psi": q(FcE_weak), "FcE_strong_psi": q(FcE_strong),
        "ratio_weak": q(ratio_weak), "ratio_strong": q(ratio_strong),
        "Cp_weak": q(Cp_weak), "Cp_strong": q(Cp_strong),
        "Cp": q(Cp), "governing_axis": governing_axis,
        "Fc_prime_psi": q(Fc_prime), "fc_psi": q(fc_psi), "Pa_lbf": q(Pa_lbf),
    }

    checks = {
        "axial": {"demand": q(P_lbf), "capacity": q(Pa_lbf), "ok": P_lbf <= Pa_lbf},
    }

    return {
        "tool": "wood_stud_column",
        "version": "0.1",
        "project": inp.get("project", ""),
        "prepared_by": inp.get("prepared_by", ""),
        "values": values,
        "checks": checks,
    }


def main(argv: list[str] | None = None) -> int:
    parser = argparse.ArgumentParser(description="Calculate the wood stud column axial capacity from a YAML input file.")
    parser.add_argument("--input", "-i", type=Path, default=HERE / "input.yaml", help="YAML input path")
    parser.add_argument("--output", "-o", type=Path, help="JSON output path; defaults beside the input")
    parser.add_argument("--stdout", action="store_true", help="Write the complete JSON result to stdout instead of a file")
    args = parser.parse_args(argv)

    input_path = args.input
    output_path = args.output or input_path.with_name("results.json")
    with input_path.open(encoding="utf-8") as handle:
        result = compute(yaml.safe_load(handle))
    serialized = json.dumps(result, indent=2) + "\n"
    if args.stdout:
        sys.stdout.write(serialized)
    else:
        output_path.write_text(serialized, encoding="utf-8")
        print(output_path)
    return 0


if __name__ == "__main__":
    raise SystemExit(main())