Function bodies 335 total

def build_driven_pile_group(
    params: dict,
    site: SiteProfile,
    tag_alloc: TagAllocator | None = None,
    spring_spacing_ft: float = 1.0,
) -> FoundationModel:
    """Build driven pile group foundation model.

    Individual piles modeled with p-y/t-z/Q-z springs (like shafts).
    Pile cap modeled as rigid body with rigidLink constraints.
    Group effects via AASHTO p-multipliers.

    Args:
        params: {
            pile_type: str,            # "HP14x73", "PIPE16", "CONC18", etc.
            n_rows: int,               # Rows in loading direction
            n_cols: int,               # Columns perpendicular
            spacing_ft: float,         # Center-to-center spacing (ft)
            length_ft: float,          # Embedded pile length (ft)
            cap_length_ft: float,      # Pile cap L (ft)
            cap_width_ft: float,       # Pile cap W (ft)
            cap_thickness_ft: float,   # Pile cap thickness (ft)
        }
        site: SiteProfile.
        tag_a

def build_pile_bent(
    params: dict,
    site: SiteProfile,
    tag_alloc: TagAllocator | None = None,
    spring_spacing_ft: float = 1.0,
) -> FoundationModel:
    """Build pile bent (trestle) foundation model.

    Piles extend through soil with p-y springs and above ground.
    Simple pile bent: no cap, piles connect directly to superstructure.
    Capped pile bent: cap beam ties pile heads.

    Args:
        params: {
            pile_type: str,            # "HP14x73", "PIPE16", etc.
            n_piles: int,              # Number of piles in bent
            spacing_ft: float,         # Center-to-center spacing (ft)
            embedded_length_ft: float, # Below ground (ft)
            exposed_length_ft: float,  # Above ground to cap/deck (ft)
            batter_deg: float,         # Batter angle (degrees), 0 = vertical
            has_cap: bool,             # Whether to include cap beam
        }
        site: SiteProfile.
        tag_alloc: Tag allocator.
        spring_spaci

def create_foundation(
    foundation_type: str,
    params: dict,
    site_profile: dict | SiteProfile,
) -> FoundationModel:
    """Create complete foundation model with OpenSees commands.

    This is the main entry point for the foundation tool. It makes
    engineering decisions internally based on foundation type and site data.

    Args:
        foundation_type: One of "drilled_shaft", "driven_pile_group",
            "spread_footing", "pile_bent".
        params: Foundation-specific parameters (see individual builders).
        site_profile: Either a SiteProfile object or a dict with:
            {
                "layers": [
                    {"soil_type": "soft_clay", "top_depth_ft": 0,
                     "thickness_ft": 20, "su_ksf": 1.0, "gamma_pcf": 110},
                    ...
                ],
                "gwt_depth_ft": 10,
                "scour": {"water_crossing": true, "depth_ft": 6}
            }

    Returns:
        FoundationModel with nodes, elements,

class CaseForces:
    """Unfactored element forces for a single load case.

    Attributes:
        case_name:      Unique identifier, e.g. ``"DC1_deck_slab"``.
        case_type:      AASHTO load type code: ``"DC"``, ``"DW"``, ``"LL"``,
                        ``"WS"``, ``"WL"``, ``"TU"``, ``"EQ"``, ``"BR"``.
        element_forces: Mapping of element tag → force dict.
                        Each force dict may contain any subset of:
                        ``"Mz_i"`` (moment at i-end, kip-ft or kip-in),
                        ``"Vy_i"`` (shear at i-end),
                        ``"N_i"``  (axial at i-end).
                        Missing keys are treated as zero.
    """

    case_name: str
    case_type: str  # "DC","DW","LL","WS","WL","TU","EQ","BR"
    element_forces: dict = field(default_factory=dict)

class ForceEnvelope:
    """Factored force envelope for a single element.

    Stores the maximum and minimum factored force effects across all
    applicable AASHTO load combinations, together with the controlling
    combination name.

    Attributes:
        element_tag:   OpenSees element tag (integer key).
        Mz_max:        Maximum (most positive) factored bending moment.
        Mz_min:        Minimum (most negative) factored bending moment.
        Vy_max:        Maximum (most positive) factored shear.
        Vy_min:        Minimum (most negative) factored shear.
        N_max:         Maximum (most tensile/positive) factored axial force.
        N_min:         Minimum (most compressive/negative) factored axial force.
        Mz_max_combo:  Name of the controlling load combination for Mz_max.
        Mz_min_combo:  Name of the controlling load combination for Mz_min.
        Vy_max_combo:  Name of the controlling load combination for Vy_max.
        Vy_min_combo:  Name of 

    def to_dict(self) -> dict:
        """Serialize to plain dict (for JSON / report output)."""
        return {
            "element_tag": self.element_tag,
            "Mz_max": self.Mz_max,
            "Mz_min": self.Mz_min,
            "Vy_max": self.Vy_max,
            "Vy_min": self.Vy_min,
            "N_max": self.N_max,
            "N_min": self.N_min,
            "Mz_max_combo": self.Mz_max_combo,
            "Mz_min_combo": self.Mz_min_combo,
            "Vy_max_combo": self.Vy_max_combo,
            "Vy_min_combo": self.Vy_min_combo,
            "N_max_combo": self.N_max_combo,
            "N_min_combo": self.N_min_combo,
        }

def _resolve_limit_states(
    limit_states: Optional[list[str]],
) -> dict[str, dict[str, object]]:
    """Return the factor table filtered to the requested limit states.

    Parameters
    ----------
    limit_states:
        List of limit state keys (e.g. ``["Strength_I", "Service_II"]``).
        If ``None``, all limit states in :data:`LIMIT_STATE_FACTORS` are used.

    Returns
    -------
    dict
        Subset of :data:`LIMIT_STATE_FACTORS` matching the request.
    """
    if limit_states is None:
        return LIMIT_STATE_FACTORS

    resolved: dict[str, dict[str, object]] = {}
    for key in limit_states:
        if key in LIMIT_STATE_FACTORS:
            resolved[key] = LIMIT_STATE_FACTORS[key]
        else:
            logger.warning("Unknown limit state '%s' — skipped.", key)
    return resolved

def _expand_combos(
    factor_table: dict[str, dict[str, object]],
) -> list[tuple[str, dict[str, float]]]:
    """Expand the factor table into concrete combo variants.

    For limit states where DC (or DW/TU) has ``(max_factor, min_factor)``
    tuples, two variants are generated:
      * ``"<LimitState>_max"`` — all tuples resolved to their [0] (max) value
      * ``"<LimitState>_min"`` — all tuples resolved to their [1] (min) value

    For limit states without any tuple entries a single variant with the
    original name is produced.

    Parameters
    ----------
    factor_table:
        Mapping of limit state name → {case_type: factor_or_tuple_or_None}.

    Returns
    -------
    list of (combo_name, {case_type: resolved_factor})
        Ready-to-use combo definitions.
    """
    combos: list[tuple[str, dict[str, float]]] = []

    for ls_name, type_factors in factor_table.items():
        has_dual = any(
            isinstance(v, tuple) for v in type_factors.values() if v 

def _group_cases_by_type(
    case_forces: list[CaseForces],
) -> dict[str, list[CaseForces]]:
    """Group CaseForces objects by their case_type."""
    groups: dict[str, list[CaseForces]] = {}
    for cf in case_forces:
        groups.setdefault(cf.case_type, []).append(cf)
    return groups

def _collect_element_tags(case_forces: list[CaseForces]) -> set:
    """Collect all unique element tags across all load cases."""
    tags: set = set()
    for cf in case_forces:
        tags.update(cf.element_forces.keys())
    return tags

def _compute_combo_force(
    elem_tag: int,
    type_factors: dict[str, float],
    cases_by_type: dict[str, list[CaseForces]],
) -> tuple[float, float, float]:
    """Compute factored (Mz, Vy, N) for one element under one combo variant.

    Parameters
    ----------
    elem_tag:      Element tag to evaluate.
    type_factors:  {case_type: factor} for this combo variant.
    cases_by_type: Grouped CaseForces by type.

    Returns
    -------
    (Mz, Vy, N) — summed factored forces.
    """
    Mz = 0.0
    Vy = 0.0
    N = 0.0

    for case_type, factor in type_factors.items():
        for cf in cases_by_type.get(case_type, []):
            forces = cf.element_forces.get(elem_tag, {})
            Mz += factor * forces.get("Mz_i", 0.0)
            Vy += factor * forces.get("Vy_i", 0.0)
            N += factor * forces.get("N_i", 0.0)

    return Mz, Vy, N

def compute_factored_envelopes(
    case_forces: list[CaseForces],
    limit_states: Optional[list[str]] = None,
) -> dict[int, ForceEnvelope]:
    """Compute AASHTO LRFD factored force envelopes for all elements.

    For each element found in ``case_forces`` and for each applicable load
    combination, sums ``γᵢ × unfactored_force_i`` over all contributing
    load cases.  Tracks the algebraic maximum and minimum of each force
    component together with the controlling combination name.

    DC (and DW, TU when applicable) are evaluated with *both* their maximum
    and minimum load factors to capture reversed-effect scenarios (e.g.
    hogging vs. sagging at continuous span supports).

    The LL "envelope" is naturally obtained by summing all CaseForces objects
    whose ``case_type == "LL"`` — callers should pass the LL cases that
    represent the governing truck/tandem positions.

    Parameters
    ----------
    case_forces:
        List of :class:`CaseForces` objects, one p

def compute_dcr_from_envelopes(
    envelopes: dict[int, ForceEnvelope],
    section_capacities: "dict | list",
) -> list[dict]:
    """Compute demand-to-capacity ratios (DCR) from factored envelopes.

    Parameters
    ----------
    envelopes:
        Output of :func:`compute_factored_envelopes`.
    section_capacities:
        Either:

        * A **dict** keyed by element tag →
          ``{"Mn": float, "Vn": float, "Pn": float}``,
          where Mn is nominal moment capacity, Vn shear, Pn axial.
          Any key may be omitted — missing capacities are simply not
          checked.

        * A **list** of dicts each containing an ``"element_tag"`` key
          plus ``"Mn"``, ``"Vn"``, ``"Pn"`` as above.

    Returns
    -------
    list[dict]
        One entry per element in ``envelopes``.  Each dict contains:

        * ``"element_tag"``
        * ``"DCR_Mz"``, ``"DCR_Mz_max"``, ``"DCR_Mz_min"``  (if Mn given)
        * ``"DCR_Vy"``, ``"DCR_Vy_max"``, ``"DCR_Vy_min"``  (if Vn

class LoadCase:
    """A single load case for OpenSees analysis.

    Attributes:
        name:        Unique identifier, e.g. ``"DC1_deck_slab"``.
        category:    ``'standard'`` or ``'adversarial'``.
        load_type:   AASHTO code (DC, DW, LL, EQ, WS, WL, BR, TU, TG, P, SC).
        description: Human-readable description for reports.
        loads:       List of load dicts, each with a ``type`` key:
                     ``'distributed'``, ``'point'``, ``'spectrum'``,
                     ``'thermal'``, ``'modification'``.
    """
    name: str
    category: str
    load_type: str
    description: str
    loads: list[dict] = field(default_factory=list)

    def to_dict(self) -> dict:
        return asdict(self)

class LoadCombination:
    """A factored load combination per AASHTO Table 3.4.1-1.

    Attributes:
        name:        e.g. ``"Strength_I_max"``.
        limit_state: e.g. ``"Strength I"``.
        factors:     Mapping of load case name → load factor.
    """
    name: str
    limit_state: str
    factors: dict = field(default_factory=dict)

    def to_dict(self) -> dict:
        return asdict(self)

class LoadModel:
    """Complete load model: all cases, combinations, and distribution factors.

    This is the canonical output of :func:`generate_loads`.
    """
    cases: list[LoadCase] = field(default_factory=list)
    combinations: list[LoadCombination] = field(default_factory=list)
    distribution_factors: dict = field(default_factory=dict)
    adversarial_cases: list[LoadCase] = field(default_factory=list)
    adversarial_combos: list[LoadCombination] = field(default_factory=list)
    moving_load_positions: dict = field(default_factory=dict)

    @property
    def total_combinations(self) -> int:
        return len(self.combinations) + len(self.adversarial_combos)

    def to_dict(self) -> dict:
        return {
            "cases": [c.to_dict() for c in self.cases],
            "combinations": [c.to_dict() for c in self.combinations],
            "distribution_factors": self.distribution_factors,
            "adversarial_cases": [c.to_dict() for c in self.adversarial_cases],

    def to_dict(self) -> dict:
        return {
            "cases": [c.to_dict() for c in self.cases],
            "combinations": [c.to_dict() for c in self.combinations],
            "distribution_factors": self.distribution_factors,
            "adversarial_cases": [c.to_dict() for c in self.adversarial_cases],
            "adversarial_combos": [c.to_dict() for c in self.adversarial_combos],
            "moving_load_positions": self.moving_load_positions,
            "total_combinations": self.total_combinations,
        }

class BridgeGeometry:
    """Bridge geometry required for load generation.

    All dimensions in feet (input) — converted to inches internally where needed.

    Attributes:
        span_ft:            Span length (ft). For multi-span, use list.
        girder_spacing_ft:  Center-to-center girder spacing (ft).
        deck_thickness_in:  Concrete deck thickness (inches).
        num_girders:        Number of girders.
        deck_width_ft:      Total deck width curb-to-curb (ft).
        num_barriers:       Number of traffic barriers (typically 2).
        barrier_weight_klf: Weight per barrier (klf). Default 0.40.
        haunch_thickness_in: Haunch (pad) thickness between deck and girder (in).
        haunch_width_in:    Haunch width (in).
        girder_weight_plf:  Steel girder self-weight (lb/ft). If 0, estimated.
        girder_depth_in:    Steel girder depth (inches).
        num_lanes:          Number of design lanes (auto-computed if 0).
        overhang_ft:        Deck overh

def _compute_dead_loads(geom: BridgeGeometry) -> list[LoadCase]:
    """Compute DC1, DC2, DC3, and DW dead load cases.

    All output loads are in kip/in (distributed) for OpenSees.
    """
    cases: list[LoadCase] = []

    # --- DC1: Structural components ---

    # Deck slab: thickness × trib width × γ_concrete
    # Trib width per girder ≈ girder_spacing for interior, half-spacing + overhang for exterior
    deck_weight_klf = (
        (geom.deck_thickness_in / 12.0)
        * geom.girder_spacing_ft
        * GAMMA_CONCRETE / 1000.0
    )
    deck_kli = deck_weight_klf / 12.0  # kip/in

    cases.append(LoadCase(
        name="DC1_deck_slab",
        category="standard",
        load_type="DC",
        description=(
            f"Concrete deck slab: {geom.deck_thickness_in}\" thick × "
            f"{geom.girder_spacing_ft}' trib width × {GAMMA_CONCRETE} pcf = "
            f"{deck_weight_klf:.3f} klf per interior girder"
        ),
        loads=[{
            "type": "distribut

def compute_distribution_factors(
    girder_spacing_ft: float,
    span_ft: float,
    deck_thickness_in: float,
    girder_depth_in: float,
    num_girders: int,
    structure_type: str = "steel",
) -> dict:
    """Compute AASHTO LRFD live-load distribution factors.

    Per AASHTO Tables 4.6.2.2.2b-1 and 4.6.2.2.3a-1 for steel I-girder
    bridges (Type "a" cross-section).

    Parameters are in mixed units matching AASHTO formulas:
    spacing in ft, span in ft, thickness in inches.

    Returns a dict with keys:
        moment_interior_1, moment_interior_2p,
        moment_exterior_1, moment_exterior_2p,
        shear_interior_1, shear_interior_2p,
        shear_exterior_1, shear_exterior_2p,
        moment_interior, moment_exterior,
        shear_interior, shear_exterior
    """
    S = girder_spacing_ft
    L = span_ft
    ts = deck_thickness_in

    # Stiffness parameter Kg (AASHTO 4.6.2.2.1-1)
    # Kg = n * (I + A * eg²)
    # For steel girder with concrete deck: n ≈ 8 (modul

def _simple_span_moment(axles: list[dict], span_ft: float) -> float:
    """Max simple-span moment (kip-ft) for a set of axles using influence lines.

    Places the resultant of the axle group at midspan and checks for maximum
    moment under each axle.
    """
    if span_ft <= 0:
        return 0.0

    total_w = sum(a["weight_kip"] for a in axles)
    if total_w <= 0:
        return 0.0

    # Resultant position from first axle
    x_r = sum(a["weight_kip"] * a["position_ft"] for a in axles) / total_w

    best_moment = 0.0
    # Try placing each axle at or near midspan
    for target_axle in axles:
        # Place this axle at midspan offset to maximize moment
        offset = span_ft / 2.0 - target_axle["position_ft"] + (target_axle["position_ft"] - x_r) / 2.0

        for axle in axles:
            x = axle["position_ft"] + offset
            if x < 0 or x > span_ft:
                continue
            # Moment at axle location from simple beam
            ra = total_w * (span

def _simple_span_shear(axles: list[dict], span_ft: float) -> float:
    """Max simple-span shear (kip) at the support for a set of axles."""
    if span_ft <= 0:
        return 0.0

    best_shear = 0.0
    # Place first axle at left support for max left reaction
    for i, lead_axle in enumerate(axles):
        shift = -lead_axle["position_ft"]
        reaction = 0.0
        for axle in axles:
            x = axle["position_ft"] + shift
            if x < 0 or x > span_ft:
                continue
            reaction += axle["weight_kip"] * (span_ft - x) / span_ft
        best_shear = max(best_shear, reaction)

    return best_shear

def _compute_live_loads(geom: BridgeGeometry) -> list[LoadCase]:
    """Generate HL-93 live load cases."""
    cases: list[LoadCase] = []
    L = geom.span_ft

    # --- Design Truck ---
    truck_moment = _simple_span_moment(DESIGN_TRUCK_AXLES, L)
    truck_shear = _simple_span_shear(DESIGN_TRUCK_AXLES, L)

    # With impact
    truck_moment_im = truck_moment * (1 + IM_TRUCK)
    truck_shear_im = truck_shear * (1 + IM_TRUCK)

    cases.append(LoadCase(
        name="LL_HL93_truck",
        category="standard",
        load_type="LL",
        description=(
            f"HL-93 Design Truck (8k-32k-32k, 14' spacing) on {L}' span. "
            f"M={truck_moment:.1f} k-ft, V={truck_shear:.1f} k (before IM). "
            f"IM={IM_TRUCK*100:.0f}%: M={truck_moment_im:.1f} k-ft, V={truck_shear_im:.1f} k"
        ),
        loads=[{
            "type": "point",
            "vehicle": "HL93_truck",
            "axles": DESIGN_TRUCK_AXLES,
            "max_moment_kft": round(truck_moment, 2),
 

def _compute_permit_loads(geom: BridgeGeometry) -> list[LoadCase]:
    """Generate permit vehicle load cases."""
    cases: list[LoadCase] = []

    for key, vehicle in PERMIT_VEHICLES.items():
        axles = vehicle["axles"]
        moment = _simple_span_moment(axles, geom.span_ft)
        shear = _simple_span_shear(axles, geom.span_ft)
        moment_im = moment * (1 + IM_TRUCK)
        shear_im = shear * (1 + IM_TRUCK)

        cases.append(LoadCase(
            name=f"P_{key}",
            category="standard",
            load_type="P",
            description=(
                f"Permit: {vehicle['name']} ({vehicle['total_kip']}k) on {geom.span_ft}' span. "
                f"M={moment_im:.1f} k-ft (w/ IM), V={shear_im:.1f} k"
            ),
            loads=[{
                "type": "point",
                "vehicle": vehicle["name"],
                "axles": axles,
                "total_weight_kip": vehicle["total_kip"],
                "max_moment_kft": round(moment, 2),
    

def _compute_thermal_loads(
    geom: BridgeGeometry,
    thermal_profile: Optional[dict] = None,
    state: str = "",
) -> list[LoadCase]:
    """Generate uniform temperature (TU) and gradient (TG) load cases."""
    cases: list[LoadCase] = []

    # --- TU: Uniform temperature change ---
    if thermal_profile:
        delta_t = thermal_profile.get("delta_t", 120.0)
        t_min = thermal_profile.get("t_min", -10)
        t_max = thermal_profile.get("t_max", 110)
    else:
        delta_t = 120.0
        t_min = -10
        t_max = 110

    # Setting temperature assumed at 60°F
    t_set = 60.0
    delta_t_rise = t_max - t_set
    delta_t_fall = t_set - t_min

    # Coefficient of thermal expansion
    alpha = 6.5e-6 if geom.structure_type == "steel" else 5.5e-6  # per °F

    cases.append(LoadCase(
        name="TU_rise",
        category="standard",
        load_type="TU",
        description=(
            f"Uniform temperature rise: +{delta_t_rise}°F (T_set={t_set}°F → T_max={t_m

def _get_kz(height_ft: float) -> float:
    """Interpolate velocity pressure exposure coefficient Kz (Exposure C)."""
    if height_ft <= 15:
        return 0.85
    if height_ft >= 90:
        return 1.24
    # Linear interpolation
    keys = sorted(_KZ_TABLE_EXP_C.keys())
    for i in range(len(keys) - 1):
        if keys[i] <= height_ft <= keys[i + 1]:
            z1, z2 = keys[i], keys[i + 1]
            k1, k2 = _KZ_TABLE_EXP_C[z1], _KZ_TABLE_EXP_C[z2]
            return k1 + (k2 - k1) * (height_ft - z1) / (z2 - z1)
    return 1.0

def _compute_wind_loads(
    geom: BridgeGeometry,
    wind_v: int = 115,
    bridge_height_ft: float = 30.0,
) -> list[LoadCase]:
    """Generate wind load cases (WS and WL)."""
    cases: list[LoadCase] = []

    # --- WS: Wind on structure ---
    kz = _get_kz(bridge_height_ft)
    # Velocity pressure (psf)
    qz = WIND_PRESSURE_COEFF * kz * wind_v ** 2
    # Wind pressure on exposed surface
    p_ws = qz * WIND_GUST_FACTOR * WIND_CD_GIRDER  # psf

    # Exposed depth = girder depth + barrier (typ 3.5 ft)
    exposed_depth_ft = (geom.girder_depth_in or 48.0) / 12.0 + 3.5  # default 48" if None
    # Distributed lateral load (klf)
    ws_klf = p_ws * exposed_depth_ft / 1000.0
    ws_kli = ws_klf / 12.0

    cases.append(LoadCase(
        name="WS_wind_on_structure",
        category="standard",
        load_type="WS",
        description=(
            f"Wind on structure: V={wind_v} mph, Kz={kz:.3f}, qz={qz:.1f} psf, "
            f"p={p_ws:.1f} psf on {exposed_depth_ft:.1f}' expose

def _compute_braking(geom: BridgeGeometry) -> list[LoadCase]:
    """Compute braking force (BR) per AASHTO 3.6.4."""
    # 25% of truck axle weights
    truck_total = sum(a["weight_kip"] for a in DESIGN_TRUCK_AXLES)
    br_truck = BR_TRUCK_FRACTION * truck_total

    # 5% of (truck + lane on full span)
    lane_total = DESIGN_LANE_KLF * geom.span_ft
    br_combo = BR_COMBO_FRACTION * (truck_total + lane_total)

    br_force = max(br_truck, br_combo)

    return [LoadCase(
        name="BR_braking",
        category="standard",
        load_type="BR",
        description=(
            f"Braking force: max(25%×{truck_total:.0f}k = {br_truck:.1f}k, "
            f"5%×({truck_total:.0f}+{lane_total:.1f}) = {br_combo:.1f}k) = {br_force:.1f}k. "
            f"Applied longitudinally at deck level."
        ),
        loads=[{
            "type": "point",
            "force_kip": round(br_force, 2),
            "direction": "longitudinal",
            "height": "deck_level",
            "br_tr

def _compute_seismic(
    geom: BridgeGeometry,
    seismic_profile: Optional[dict] = None,
) -> list[LoadCase]:
    """Generate seismic load cases from site response spectrum."""
    cases: list[LoadCase] = []

    if seismic_profile is None:
        seismic_profile = {
            "sds": 0.267, "sd1": 0.160, "pga": 0.10,
            "sdc": "B", "site_class": "D",
        }

    sds = seismic_profile.get("sds", 0.267)
    sd1 = seismic_profile.get("sd1", 0.160)
    pga = seismic_profile.get("pga", 0.10)
    sdc = seismic_profile.get("sdc", "B")

    # Build AASHTO design response spectrum points
    # T0 = 0.2 * SD1/SDS, Ts = SD1/SDS
    if sds > 0:
        ts = sd1 / sds
        t0 = 0.2 * ts
    else:
        ts = 1.0
        t0 = 0.2

    spectrum_points = []
    # Ramp from PGA to SDS
    spectrum_points.append({"T": 0.0, "Sa": pga})
    spectrum_points.append({"T": t0, "Sa": sds})
    # Constant plateau
    spectrum_points.append({"T": ts, "Sa": sds})
    # 1/T descent
    for t 

def _compute_scour(water_crossing: bool = False) -> list[LoadCase]:
    """Generate scour modification flag."""
    if not water_crossing:
        return []

    return [LoadCase(
        name="SC_scour",
        category="standard",
        load_type="SC",
        description=(
            "Scour: modifies foundation springs — remove soil resistance above "
            "computed scour depth (Q100 for Strength, Q500 for Extreme Event). "
            "Not a direct load — flag passed to foundation tool."
        ),
        loads=[{
            "type": "modification",
            "target": "foundation_springs",
            "action": "remove_above_scour_depth",
            "design_flood": "Q100",
            "check_flood": "Q500",
        }],
    )]

def _generate_combinations(cases: list[LoadCase]) -> list[LoadCombination]:
    """Generate all AASHTO load combinations from the case list.

    For load types with max/min factors (DC, DW, TU), generates separate
    max and min envelope combinations.
    """
    combos: list[LoadCombination] = []

    # Build lookup: load_type → [case names]
    cases_by_type: dict[str, list[str]] = {}
    for c in cases:
        lt = c.load_type
        # Map P (permit) to LL for combination purposes
        if lt == "P":
            lt = "LL"
        if lt == "TG":
            lt = "TU"  # gradient grouped with thermal
        cases_by_type.setdefault(lt, []).append(c.name)

    for limit_state, factors in _COMBO_TABLE.items():
        # Determine which envelope variants to generate
        has_dual = any(isinstance(v, tuple) for v in factors.values() if v is not None)

        if has_dual:
            # Generate max and min envelope variants
            for suffix, idx in [("max", 0), ("min", 1)]

def _generate_adversarial_cases(geom: BridgeGeometry) -> list[LoadCase]:
    """Generate adversarial (red-team) load cases.

    These represent conditions that standard practice typically doesn't model
    but can cause failures in real bridges.
    """
    cases: list[LoadCase] = []

    # ------------------------------------------------------------------
    # 1. Construction loads
    # ------------------------------------------------------------------
    cases.append(LoadCase(
        name="ADV_crane_on_overhang",
        category="adversarial",
        load_type="CONST",
        description=(
            "Construction: Crane on overhang during deck pour. "
            "30-50k point load at deck edge, partial structure (no composite action)."
        ),
        loads=[{
            "type": "point",
            "force_kip": 40.0,
            "position": "deck_edge",
            "position_from_ext_girder_ft": geom.overhang_ft,
            "composite_action": False,
            "sce

def _generate_adversarial_combos(
    standard_cases: list[LoadCase],
    adversarial_cases: list[LoadCase],
) -> list[LoadCombination]:
    """Generate adversarial load combinations.

    Each adversarial case is combined with relevant standard loads at
    appropriate load levels.
    """
    combos: list[LoadCombination] = []

    # Build DC/DW case name lists
    dc_cases = [c.name for c in standard_cases if c.load_type == "DC"]
    dw_cases = [c.name for c in standard_cases if c.load_type == "DW"]
    ll_cases = [c.name for c in standard_cases
                if c.load_type == "LL" and "governing" in c.name]
    eq_cases = [c.name for c in standard_cases if c.load_type == "EQ"]
    tu_cases = [c.name for c in standard_cases if c.load_type == "TU"]

    for adv in adversarial_cases:
        factors: dict[str, float] = {}

        # Always include dead loads
        for dc in dc_cases:
            factors[dc] = 1.25
        for dw in dw_cases:
            factors[dw] = 1.50

       

def generate_moving_load_positions(
    span_lengths: list[float],
    num_positions: int = 20,
) -> dict:
    """Generate moving load positions for HL-93 influence-line analysis.

    For each position along the bridge, computes axle locations for
    design truck, design tandem, and lane load.

    Args:
        span_lengths: List of span lengths in feet.
        num_positions: Number of truck front-axle positions to evaluate.

    Returns:
        Dict with keys:
            ``truck_positions``: list of truck position dicts
            ``tandem_positions``: list of tandem position dicts
            ``lane_load``: uniform lane load dict
    """
    total_length_ft = sum(span_lengths)

    # Generate evenly-spaced front-axle positions along bridge
    step = total_length_ft / (num_positions + 1)
    positions_ft = [step * (i + 1) for i in range(num_positions)]

    truck_positions = []
    for pos in positions_ft:
        # Design truck: 8k at front, 32k at 14ft, 32k at 28ft behind
  

def generate_load_combination_script(
    load_cases: list[LoadCase],
    combinations: list[LoadCombination],
    element_tags: list[int] | None = None,
    node_map: dict | None = None,
    girder_spacing_ft: float = 8.0,
) -> str:
    """Generate OpenSees Python code for load combination analysis.

    Produces a script block that:
      1. Defines a timeSeries + pattern for each load case
      2. For each combination, applies factored loads
      3. Runs static analysis
      4. Extracts element forces into a results dict
      5. Computes envelopes across all combinations

    Args:
        load_cases:       List of :class:`LoadCase` objects.
        combinations:     List of :class:`LoadCombination` objects.
        element_tags:     Element tags to extract forces from.
                          If None, uses ``ops.getEleTags()``.
        node_map:         Dict mapping load application descriptions to node tags.
                          If None, loads are applied to all beam el

def extract_force_envelopes(
    element_tags: list[int],
    combination_results: dict,
) -> dict:
    """Extract controlling force envelopes from load combination results.

    For each element, finds the maximum and minimum forces across all
    load combinations.

    Args:
        element_tags:        List of element tags to process.
        combination_results: Dict from combination analysis, keyed by
                             combo name, each with ``element_forces`` sub-dict.

    Returns:
        Dict keyed by element tag::

            {element_tag: {
                'max_moment': float,
                'min_moment': float,
                'max_shear': float,
                'min_shear': float,
                'max_axial': float,
                'min_axial': float,
                'controlling_combo_moment': str,
                'controlling_combo_shear': str,
                'controlling_combo_axial': str,
            }}
    """
    envelopes = {}

    for ele_tag in eleme

def generate_moving_load_script(
    span_lengths: list[float],
    num_positions: int = 20,
    girder_spacing_ft: float = 8.0,
) -> str:
    """Generate OpenSees Python code for HL-93 moving load analysis.

    Traverses the bridge with design truck and tandem at discrete positions,
    concurrent with lane load, to find the worst-case live load effects.

    Args:
        span_lengths:     List of span lengths (ft).
        num_positions:    Number of positions to evaluate.
        girder_spacing_ft: For distribution factor computation.

    Returns:
        Multi-line Python string for insertion into OpenSees script.
    """
    total_length_ft = sum(span_lengths)
    total_length_in = total_length_ft * 12.0
    step_in = total_length_in / (num_positions + 1)

    # Build axle data for inline use
    truck_axles_in = [
        (a["weight_kip"] * (1 + IM_TRUCK), a["position_ft"] * 12.0)
        for a in DESIGN_TRUCK_AXLES
    ]
    tandem_axles_in = [
        (a["weight_kip"] * (1 +

def generate_loads(
    geom: BridgeGeometry,
    site_profile: Optional[dict] = None,
    bridge_height_ft: float = 30.0,
) -> LoadModel:
    """Generate complete load model for a bridge.

    This is the canonical entry point for the loads tool.  It:
      1. Computes dead loads (DC, DW) from geometry
      2. Computes live loads (HL-93) with impact and distribution factors
      3. Generates permit vehicle loads
      4. Computes thermal (TU, TG), wind (WS, WL), braking (BR)
      5. Generates seismic (EQ) response spectrum
      6. Checks for scour (SC) modification
      7. Assembles AASHTO load combinations per Table 3.4.1-1
      8. Generates adversarial (red-team) load cases and combinations

    Args:
        geom:             :class:`BridgeGeometry` with bridge dimensions.
        site_profile:     Dict from :class:`~nlb.tools.site_recon.SiteProfile.to_dict`
                          or None (uses conservative defaults).
        bridge_height_ft: Height of bridge deck above g

class Finding:
    """A single red-team finding.

    Attributes:
        severity:          ``'CRITICAL'``, ``'WARNING'``, or ``'NOTE'``.
        vector:            Which attack vector found this.
        element:           Element tag (if applicable), else None.
        location:          Human-readable location ("Span 2, 0.4L").
        description:       What was found.
        dcr:               Demand/capacity ratio (if applicable).
        controlling_combo: Load combination that governs.
        recommendation:    What to do about it.
        precedent:         Historical failure match (if any).
    """
    severity: str
    vector: str
    element: int | None
    location: str
    description: str
    dcr: float | None
    controlling_combo: str
    recommendation: str
    precedent: str | None = None

    def to_dict(self) -> dict:
        return asdict(self)

class CascadeChain:
    """A progressive collapse chain.

    Attributes:
        trigger_element: The element whose failure starts the chain.
        chain:           Ordered list of (element, dcr) tuples in cascade.
        causes_collapse: True if the chain leads to global instability.
        description:     Human-readable cascade narrative.
    """
    trigger_element: int
    chain: list[tuple[int, float]] = field(default_factory=list)
    causes_collapse: bool = False
    description: str = ""

class SensitivityResult:
    """Result of varying one parameter in the sensitivity sweep.

    Attributes:
        parameter:      Name of the varied parameter.
        base_dcr:       DCR at baseline parameter value.
        low_dcr:        DCR at parameter -20%.
        high_dcr:       DCR at parameter +20%.
        delta_dcr:      Maximum absolute change in DCR.
        classification: ``'DOMINANT'``, ``'MODERATE'``, or ``'INSENSITIVE'``.
    """
    parameter: str
    base_dcr: float
    low_dcr: float
    high_dcr: float
    delta_dcr: float
    classification: str

class HistoryMatch:
    """A match against the bridge failure database.

    Attributes:
        failure_name: Name of the historical failure.
        year:         Year of the failure.
        score:        Similarity score (higher = more similar).
        lesson:       Key engineering lesson from the failure.
        matching_factors: What matched (list of factor descriptions).
    """
    failure_name: str
    year: int
    score: int
    lesson: str
    matching_factors: list[str] = field(default_factory=list)

class RedTeamReport:
    """Complete red-team analysis report.

    Attributes:
        findings:            All findings across all attack vectors.
        risk_rating:         ``'GREEN'``, ``'YELLOW'``, or ``'RED'``.
        summary:             1-paragraph executive summary.
        attack_vectors_run:  List of attack vector names executed.
        total_load_cases:    Total number of load cases analyzed.
        total_combinations:  Total number of load combinations checked.
        analysis_time_sec:   Wall-clock time for the full red-team run.
        cascade_chains:      Progressive collapse chains (from vector 2).
        sensitivity_results: Tornado diagram data (from vector 4).
        history_matches:     Failure database matches (from vector 7).
        robustness_results:  Component removal results (from vector 6).
    """
    findings: list[Finding] = field(default_factory=list)
    risk_rating: str = "GREEN"
    summary: str = ""
    attack_vectors_run: list[str] = field

    def to_dict(self) -> dict:
        d = asdict(self)
        return d

def load_failure_database() -> list[dict]:
    """Load the bridge failure database from JSON.

    Returns:
        List of failure records with keys: name, year, type, spans,
        material, failure_mode, cause, lesson, details.
    """
    try:
        with open(_FAILURES_JSON, "r") as f:
            data = json.load(f)
        if isinstance(data, list):
            return data
        return []
    except (FileNotFoundError, json.JSONDecodeError) as exc:
        logger.warning("Failed to load failure database: %s", exc)
        return []

def dcr_scanner(
    element_results: list[dict],
    critical_threshold: float = 1.0,
    warning_threshold: float = 0.85,
    note_threshold: float = 0.70,
) -> list[Finding]:
    """Scan all elements for demand/capacity ratio exceedances.

    Classifies each element as CRITICAL (DCR > 1.0), WARNING (> 0.85),
    or NOTE (> 0.70). Elements below 0.70 are not flagged.

    Args:
        element_results: List of dicts, each with:
            - ``element``: int — element tag
            - ``location``: str — human-readable location
            - ``dcr``: float — demand/capacity ratio
            - ``force_type``: str — "moment", "shear", or "axial"
            - ``controlling_combo``: str — load combination name
            - ``demand``: float (optional) — demand value
            - ``capacity``: float (optional) — capacity value
        critical_threshold: DCR threshold for CRITICAL. Default 1.0.
        warning_threshold:  DCR threshold for WARNING. Default 0.85.
        note_thresho

def failure_cascade(
    element_results: list[dict],
    analyze_fn=None,
    dcr_threshold: float = 1.0,
    progressive_collapse_combo: str = "DC + 0.5×LL",
) -> tuple[list[Finding], list[CascadeChain]]:
    """Analyze progressive collapse by removing failed elements.

    For each element with DCR above the threshold, conceptually removes it
    (sets stiffness to near-zero) and re-analyzes to find chain reactions.

    Args:
        element_results:  Element DCR results (same format as dcr_scanner).
        analyze_fn:       Callable(removed_elements: list[int]) -> list[dict].
                          Re-runs analysis with specified elements removed and
                          returns updated element_results. If None, uses a
                          simplified estimation model.
        dcr_threshold:    DCR threshold to trigger cascade analysis.
        progressive_collapse_combo: Name of the load combo used.

    Returns:
        Tuple of (findings, cascade_chains).
    """
 

def _estimate_cascade(
    trigger: dict,
    all_results: list[dict],
    chain: CascadeChain,
) -> None:
    """Simplified cascade estimation without re-analysis.

    Estimates demand redistribution based on element adjacency and
    current utilization. Elements already near capacity are most
    vulnerable to cascade.
    """
    trigger_elem = trigger.get("element", 0)
    trigger_dcr = trigger.get("dcr", 0.0)

    # Estimate redistribution: nearby elements pick up ~30% more load
    redistribution_factor = 0.30

    for er in all_results:
        elem = er.get("element")
        if elem == trigger_elem:
            continue
        current_dcr = er.get("dcr", 0.0)
        if current_dcr is None:
            continue

        # Simple adjacency: elements within ±2 tags are "adjacent"
        if elem is not None and trigger_elem is not None:
            if abs(elem - trigger_elem) <= 2:
                new_dcr = current_dcr * (1.0 + redistribution_factor)
                if new_dc

def construction_vulnerability(
    stage_results: list[dict],
    site_constraints: dict | None = None,
) -> list[Finding]:
    """Analyze vulnerability during construction stages.

    Args:
        stage_results: List of dicts with:
            - ``stage_name``: str — name of the construction stage
            - ``stage_number``: int — sequential stage number
            - ``max_dcr``: float — maximum DCR across all elements in this stage
            - ``max_dcr_element``: int — element with highest DCR
            - ``max_dcr_location``: str — location of highest DCR element
            - ``requires_temp_support``: bool — whether temp supports needed
            - ``description``: str — description of what happens in this stage
        site_constraints: Dict of site constraints, e.g.:
            - ``no_equipment_in_water``: bool
            - ``restricted_access``: list[str]
            - ``night_work_only``: bool

    Returns:
        List of Finding objects for construction vuln

def sensitivity_sweep(
    base_dcr: float,
    analyze_fn=None,
    parameter_results: list[dict] | None = None,
) -> tuple[list[Finding], list[SensitivityResult]]:
    """Vary key parameters ±20% and classify sensitivity.

    Can operate in two modes:
    1. With ``analyze_fn``: calls function for each parameter variation
    2. With ``parameter_results``: uses pre-computed results

    Args:
        base_dcr:          DCR at baseline parameter values.
        analyze_fn:        Callable(param_key: str, factor: float) -> float.
                           Returns max DCR with the parameter scaled by factor.
                           factor < 1.0 means reduced, > 1.0 means increased.
        parameter_results: Pre-computed list of dicts with:
            - ``parameter``: str — parameter name
            - ``low_dcr``: float — DCR at -20%
            - ``high_dcr``: float — DCR at +20%

    Returns:
        Tuple of (findings, sensitivity_results).
    """
    findings: list[Finding]