Function bodies 335 total

    def summary(self) -> str:
        return (
            f"BearingModel ({self.bearing_type}): "
            f"{len(self.nodes)} nodes, {len(self.elements)} elements, "
            f"{len(self.materials)} materials | "
            f"top={self.top_nodes}, bot={self.bottom_nodes}"
        )

class TagAllocator:
    """Sequential tag allocator for OpenSees objects."""

    def __init__(self, start: int = 1):
        self._next = start

    def next(self, count: int = 1) -> int | list[int]:
        if count == 1:
            tag = self._next
            self._next += 1
            return tag
        tags = list(range(self._next, self._next + count))
        self._next += count
        return tags

    @property
    def current(self) -> int:
        return self._next

    def next(self, count: int = 1) -> int | list[int]:
        if count == 1:
            tag = self._next
            self._next += 1
            return tag
        tags = list(range(self._next, self._next + count))
        self._next += count
        return tags

def _elastomeric_stiffness(cfg: ElastomericConfig) -> dict:
    """Compute elastomeric bearing stiffness properties.

    Per AASHTO LRFD 14.7.5 and 14.7.6.

    Returns dict with Kh, Kv, shape_factor, rotation_capacity.
    """
    A = cfg.length_in * cfg.width_in  # plan area (in²)
    h_rt = cfg.total_rubber_thickness_in  # total rubber thickness
    G = cfg.shear_modulus_ksi

    # Shape factor S = loaded plan area / (perimeter × layer thickness)
    # Per AASHTO 14.7.5.1-1
    perimeter = 2 * (cfg.length_in + cfg.width_in)
    t_layer = cfg.layer_thickness_in
    S = A / (perimeter * t_layer) if (perimeter * t_layer) > 0 else 1.0

    # Horizontal stiffness: Kh = G × A / h_rt
    Kh = G * A / h_rt if h_rt > 0 else 1.0e6

    # Effective compressive modulus (shape factor dependent)
    # Ec ≈ 3G(1 + 2κS²) for steel-reinforced bearings
    # κ = material constant ≈ 0.93 for plain pads, ≈ 1.0 for reinforced
    kappa = 1.0
    Ec = 3.0 * G * (1.0 + 2.0 * kappa * S ** 2)

    # Vertic

def build_elastomeric(
    cfg: ElastomericConfig,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build elastomeric bearing model.

    Creates zeroLength element with Elastic material for shear
    and ENT (compression-only) for vertical.
    """
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="elastomeric", compression_only=True)

    props = _elastomeric_stiffness(cfg)
    model.properties = props

    # --- Nodes ---
    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # --- Materials ---
    # Horizontal shear (elastic)
    mat_h = tags.next()
    model.materials.append({
        "tag": mat_h, "type": "Elastic",
        "name": "elastomeric_shear",
        "k_kip_per_in": props["Kh_

def build_pot_fixed(
    cfg: PotBearingConfig,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build fixed pot bearing — rigid in both horizontal directions."""
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="pot_fixed", compression_only=True)

    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # Fixed horizontal: very high stiffness
    mat_h = tags.next()
    model.materials.append({
        "tag": mat_h, "type": "Elastic",
        "name": "pot_fixed_horizontal",
        "k_kip_per_in": cfg.stiffness_kip_per_in,
    })

    # Vertical: compression only
    mat_v = tags.next()
    model.materials.append({
        "tag": mat_v, "type": "ENT",
        "name": "pot_vertical",
        "k

def build_pot_guided(
    cfg: PotBearingConfig,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build guided pot bearing — free in guided direction, fixed in other."""
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="pot_guided", compression_only=True)

    guide_dir = cfg.guide_direction or 1  # default: free longitudinal

    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # Fixed direction: very stiff
    mat_fixed = tags.next()
    model.materials.append({
        "tag": mat_fixed, "type": "Elastic",
        "name": "pot_guided_fixed",
        "k_kip_per_in": cfg.stiffness_kip_per_in,
    })

    # Guided direction: very low stiffness (essentially free)
    mat_free = tags.next()
  

def _ptfe_friction(ptfe_type: str, temperature: str = "reference") -> dict:
    """Look up PTFE friction coefficients.

    Args:
        ptfe_type: Key into PTFE_FRICTION table.
        temperature: "reference" (~68°F), "cold" (~20°F), "hot" (~100°F).

    Returns:
        Dict with mu_slow, mu_fast.
    """
    if ptfe_type not in PTFE_FRICTION:
        ptfe_type = "glass_filled"  # default

    table = PTFE_FRICTION[ptfe_type]

    if temperature == "cold":
        return {
            "mu_slow": table["mu_slow_cold"],
            "mu_fast": table["mu_fast_cold"],
        }
    elif temperature == "hot":
        # Hot = lower friction (interpolate 20% reduction from reference)
        return {
            "mu_slow": table["mu_slow"] * 0.80,
            "mu_fast": table["mu_fast"] * 0.80,
        }
    else:
        return {
            "mu_slow": table["mu_slow"],
            "mu_fast": table["mu_fast"],
        }

def build_ptfe_sliding(
    cfg: PTFEConfig,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build PTFE sliding bearing with velocity-dependent friction.

    Generates upper and lower bound cases for temperature effects.
    Uses flatSliderBearing element concept modeled with zeroLength
    and friction material.
    """
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="ptfe_sliding", compression_only=True)

    # Reference friction
    mu_ref = _ptfe_friction(cfg.ptfe_type, "reference")
    mu_cold = _ptfe_friction(cfg.ptfe_type, "cold")
    mu_hot = _ptfe_friction(cfg.ptfe_type, "hot")

    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # Friction material (reference temperature)
    

def build_fp_single(
    cfg: FPSingleConfig,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build single friction pendulum bearing.

    Self-centering period T = 2π√(R/g).
    Uses SingleFPBearing element concept.
    """
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="fp_single", compression_only=True)

    # Self-centering period
    T = 2.0 * math.pi * math.sqrt(cfg.radius_in / G_ACCEL)

    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # Friction material
    mat_fric = tags.next()
    model.materials.append({
        "tag": mat_fric, "type": "frictionModel",
        "name": "fp_friction",
        "mu": cfg.mu,
    })

    # Vertical: compression only
    mat_v = tags.next()
  

def build_fp_triple(
    cfg: FPTripleConfig,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build triple friction pendulum bearing.

    Multi-stage behavior with 4 friction surfaces and 4 radii.
    Uses TripleFrictionPendulum element concept.
    """
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="fp_triple", compression_only=True)

    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # Friction materials per surface
    fric_tags = []
    for i, (mu, R) in enumerate(zip(
        [cfg.mu1, cfg.mu2, cfg.mu3, cfg.mu4],
        [cfg.R1, cfg.R2, cfg.R3, cfg.R4],
    )):
        mat_tag = tags.next()
        model.materials.append({
            "tag": mat_tag, "type": "frictionModel",
   

def build_integral(
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build integral (monolithic) connection.

    Uses equalDOF constraint — no relative displacement.
    No bearing element, just a constraint between super and sub nodes.
    """
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="integral", compression_only=False)

    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # equalDOF constraint: all DOFs tied
    model.constraints.append({
        "type": "equalDOF",
        "master": bot_node,
        "slave": top_node,
        "dofs": [1, 2, 3, 4, 5, 6],
    })

    model.properties = {
        "connection": "monolithic",
        "relative_displacement": 0.0,
    }

    # No temp

def build_rocker_roller(
    cfg: RockerRollerConfig,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
) -> BearingModel:
    """Build rocker/roller bearing (vintage steel).

    Uses ENT (compression-only) vertical + near-free horizontal.
    Flags uplift vulnerability.
    """
    tags = TagAllocator(tag_start)
    model = BearingModel(bearing_type="rocker_roller", compression_only=True)

    bot_node = tags.next()
    model.nodes.append({"tag": bot_node, "x": x, "y": y, "z": z})
    model.bottom_nodes.append(bot_node)

    top_node = tags.next()
    model.nodes.append({"tag": top_node, "x": x, "y": y, "z": z})
    model.top_nodes.append(top_node)

    # Horizontal: near-free (roller) or with some resistance (rocker)
    if cfg.roller_diameter_in > 0:
        # Roller: very free
        mat_h = tags.next()
        model.materials.append({
            "tag": mat_h, "type": "Elastic",
            "name": "roller_horizontal",
            "k_kip_per_in"

def layout_bearings(
    num_girders: int,
    girder_spacing_ft: float,
    bearing_type: str,
    support_type: str = "fixed",
    x: float = 0.0,
    y: float = 0.0,
    z_center: float = 0.0,
    tag_start: int = 1,
    **kwargs: Any,
) -> list[BearingModel]:
    """Generate bearings for a multi-girder bridge at one support line.

    Places one bearing per girder line at regular transverse spacing.

    Args:
        num_girders: Number of girders.
        girder_spacing_ft: Girder spacing (ft).
        bearing_type: BearingType value string.
        support_type: "fixed" or "expansion" (affects pot bearing choice).
        x: Longitudinal coordinate (inches).
        y: Vertical coordinate (inches).
        z_center: Transverse center coordinate (inches).
        tag_start: Starting tag.
        **kwargs: Forwarded to bearing config constructors.

    Returns:
        List of BearingModel, one per girder line.
    """
    spacing_in = girder_spacing_ft * FT_TO_IN
    total_width 

def create_bearing(
    bearing_type: str,
    x: float = 0.0,
    y: float = 0.0,
    z: float = 0.0,
    tag_start: int = 1,
    **kwargs: Any,
) -> BearingModel:
    """Create a bearing model from type string and parameters.

    Main entry point for the NLB pipeline.

    Args:
        bearing_type: One of BearingType values.
        x, y, z: Bearing location (inches).
        tag_start: Starting tag.
        **kwargs: Parameters forwarded to specific builder.

    Returns:
        BearingModel.

    Raises:
        ValueError: If bearing_type is not recognized.
    """
    bt = bearing_type.lower().replace(" ", "_").replace("-", "_")

    if bt == "elastomeric":
        cfg = kwargs.get("config", ElastomericConfig(**{
            k: v for k, v in kwargs.items()
            if k in ElastomericConfig.__dataclass_fields__
        }))
        return build_elastomeric(cfg, x, y, z, tag_start)

    elif bt == "pot_fixed":
        cfg = kwargs.get("config", PotBearingConfig(**{
         

class DistributionFactors:
    """Live-load distribution factors for a single girder line.

    Governing values (``gM_int``, etc.) are the maximum of the 1-lane and
    2+-lane cases.  The AASHTO formulas already embed the multiple-presence
    factor (m=1.20 for 1 lane, m=1.00 for 2 lanes); the governing factor is
    therefore simply ``max(g_1lane, g_2plus)``.

    The lever-rule values (exterior, 1 lane) explicitly include m=1.20.
    """

    # --- Governing (max of 1-lane and 2+-lane) ----------------------------
    gM_int: float   # Interior girder moment DF
    gM_ext: float   # Exterior girder moment DF
    gV_int: float   # Interior girder shear DF
    gV_ext: float   # Exterior girder shear DF

    # --- Individual cases -------------------------------------------------
    gM_int_1: float   # Interior moment, 1 lane loaded
    gM_int_2: float   # Interior moment, 2+ lanes loaded
    gM_ext_1: float   # Exterior moment, 1 lane (lever rule × m=1.20)
    gM_ext_2: float   # E

def compute_kg(n: float, I_girder: float, A_girder: float, eg: float) -> float:
    """Longitudinal stiffness parameter per AASHTO 4.6.2.2.1-1.

    Args:
        n:        Modular ratio, Es / Ec (dimensionless).
        I_girder: Moment of inertia of the steel girder alone (in⁴).
        A_girder: Cross-sectional area of the steel girder alone (in²).
        eg:       Distance from girder centroid to centroid of deck (in).

    Returns:
        Kg in in⁴.
    """
    return n * (I_girder + A_girder * eg ** 2)

def _lever_rule_exterior(S: float, de: float, m: float = 1.20) -> float:
    """Compute the exterior-girder DF using the lever rule.

    AASHTO 4.6.2.2.2d-1 note: "Lever Rule" for one design lane loaded.

    Lane placement (AASHTO 3.6.1.3.1):
      - First wheel at 2 ft from the interior face of the traffic barrier.
      - Second wheel 6 ft further toward interior (standard 6-ft wheel spacing).

    Geometry (all distances measured from exterior girder, positive toward
    interior of bridge):
      - ``de`` = distance from interior face of barrier to exterior web of
        exterior girder (AASHTO sign convention: positive when girder is
        inboard of barrier).
      - Wheel 1 at  (de - 2) ft from exterior girder.
      - Wheel 2 at  (de + 4) ft from exterior girder.

    The lever rule treats the deck as a simply-supported beam spanning *S* ft
    between the exterior and first interior girder.

    Args:
        S:  Girder spacing (ft).
        de: Distance from exterior gir

def compute_distribution_factors(
    S: float,
    L: float,
    ts: float,
    Kg: float,
    de: float = 2.0,
    Nb: int = 5,
    num_lanes: int = 0,
    roadway_width: float = 0.0,
    strict_roa: bool = False,
) -> DistributionFactors:
    """Compute AASHTO LRFD live-load distribution factors.

    Valid for Type (a) cross-sections: steel I-girder with composite concrete
    deck (AASHTO Tables 4.6.2.2.2b-1, 4.6.2.2.2d-1, 4.6.2.2.3a-1, 4.6.2.2.3c-1).

    The multiple-presence factor is **embedded** in both the 1-lane and 2+-lane
    tabulated formulas (m = 1.20 and 1.00 respectively), so the governing factor
    is simply ``max(g_1lane, g_2plus)`` with no additional scaling.

    For exterior girders with one lane, the lever rule is used directly and
    m = 1.20 is applied explicitly (it is not embedded in the formula).

    Args:
        S:             Girder spacing (ft), range 3.5–16.0.
        L:             Span length (ft), range 20–240.
        ts:            Structural 

class SoilType(str, Enum):
    SOFT_CLAY = "soft_clay"
    STIFF_CLAY = "stiff_clay"
    SAND = "sand"
    ROCK = "rock"

class FoundationType(str, Enum):
    DRILLED_SHAFT = "drilled_shaft"
    DRIVEN_PILE_GROUP = "driven_pile_group"
    SPREAD_FOOTING = "spread_footing"
    PILE_BENT = "pile_bent"

class SoilLayer:
    """Single soil layer in a profile.

    Attributes:
        soil_type: Classification (soft_clay, stiff_clay, sand, rock).
        top_depth_ft: Depth to top of layer from ground surface (ft).
        thickness_ft: Layer thickness (ft).
        su_ksf: Undrained shear strength (ksf) — clays.
        phi_deg: Friction angle (degrees) — sands.
        gamma_pcf: Total unit weight (pcf).
        N_spt: SPT blow count (blows/ft).
        eps50: Strain at 50% strength — clays (default from su).
        k_py_pci: Initial modulus of subgrade reaction for p-y (pci) — sand.
        qu_ksf: Unconfined compressive strength (ksf) — rock.
    """
    soil_type: str
    top_depth_ft: float
    thickness_ft: float
    su_ksf: float = 0.0
    phi_deg: float = 0.0
    gamma_pcf: float = 120.0
    N_spt: int = 0
    eps50: float | None = None
    k_py_pci: float | None = None
    qu_ksf: float = 0.0

    @property
    def bot_depth_ft(self) -> float:
        return self.top_depth_ft

    def get_eps50(self) -> float:
        """Return eps50; estimate from su if not provided (FHWA Table)."""
        if self.eps50 is not None:
            return self.eps50
        # Matlock recommended values based on su
        su = self.su_ksf
        if su <= 0.5:
            return 0.02
        elif su <= 1.0:
            return 0.01
        elif su <= 2.0:
            return 0.007
        elif su <= 4.0:
            return 0.005
        else:
            return 0.004

    def get_k_py(self) -> float:
        """Initial modulus of subgrade reaction (pci) for sand p-y.

        Per API RP 2GEO Table for submerged sand.
        """
        if self.k_py_pci is not None:
            return self.k_py_pci
        # Approximate from friction angle (pci)
        phi = self.phi_deg
        if phi <= 25:
            return 20.0
        elif phi <= 30:
            return 35.0
        elif phi <= 35:
            return 55.0
        elif phi <= 40:
            return 90.0
        else:
            return 115.0

class SiteProfile:
    """Complete site geotechnical profile."""
    layers: list[SoilLayer]
    gwt_depth_ft: float = 10.0  # groundwater table depth from surface
    scour: dict | None = None   # {water_crossing: bool, depth_ft: float}

    @property
    def gwt_depth_in(self) -> float:
        return self.gwt_depth_ft * FT_TO_IN

    def layer_at_depth(self, depth_ft: float) -> SoilLayer:
        """Return the soil layer at a given depth."""
        for layer in self.layers:
            if layer.top_depth_ft <= depth_ft < layer.bot_depth_ft:
                return layer
        # Below last layer — return last
        return self.layers[-1]

    def effective_vertical_stress(self, depth_ft: float) -> float:
        """Compute effective vertical stress at depth (ksi).

        Integrates gamma through layers, subtracts pore pressure below GWT.
        """
        sigma_v = 0.0
        prev_depth = 0.0
        for layer in sorted(self.layers, key=lambda l: l.top_depth_ft):
           

    def layer_at_depth(self, depth_ft: float) -> SoilLayer:
        """Return the soil layer at a given depth."""
        for layer in self.layers:
            if layer.top_depth_ft <= depth_ft < layer.bot_depth_ft:
                return layer
        # Below last layer — return last
        return self.layers[-1]

    def effective_vertical_stress(self, depth_ft: float) -> float:
        """Compute effective vertical stress at depth (ksi).

        Integrates gamma through layers, subtracts pore pressure below GWT.
        """
        sigma_v = 0.0
        prev_depth = 0.0
        for layer in sorted(self.layers, key=lambda l: l.top_depth_ft):
            top = max(layer.top_depth_ft, prev_depth)
            bot = min(layer.bot_depth_ft, depth_ft)
            if top >= bot:
                continue
            dz = (bot - top) * FT_TO_IN
            sigma_v += layer.gamma_pci * dz
            prev_depth = bot
            if bot >= depth_ft:
                break

        # Pore pressure below GWT
        if depth_ft > self.gwt_depth_ft:
            u = GAMMA_WATER_PCI * (depth_ft - self.gwt_depth_ft) * FT_TO_IN
            sigma_v -= u

        return max(sigma_v, 0.0)

class FoundationModel:
    """Complete foundation model output.

    All tags are integers suitable for direct use in OpenSees commands.
    All coordinates and forces are in kip-inch-second units.
    """
    nodes: list[dict] = field(default_factory=list)
    elements: list[dict] = field(default_factory=list)
    springs: list[dict] = field(default_factory=list)
    materials: list[dict] = field(default_factory=list)
    boundary_conditions: list[dict] = field(default_factory=list)
    top_node: int = 0
    base_node: int = 0
    capacity: dict = field(default_factory=dict)
    cases: dict = field(default_factory=dict)
    opensees_commands: list[str] = field(default_factory=list)

    def summary(self) -> str:
        return (
            f"FoundationModel: {len(self.nodes)} nodes, "
            f"{len(self.elements)} elements, {len(self.springs)} springs, "
            f"{len(self.materials)} materials | "
            f"top_node={self.top_node}, base_node={self.base_node}"
        

    def summary(self) -> str:
        return (
            f"FoundationModel: {len(self.nodes)} nodes, "
            f"{len(self.elements)} elements, {len(self.springs)} springs, "
            f"{len(self.materials)} materials | "
            f"top_node={self.top_node}, base_node={self.base_node}"
        )

class PYCurve:
    """p-y curve computation for lateral soil resistance.

    References:
        - Matlock (1970): Soft clay
        - Reese et al. (1975): Stiff clay
        - API RP 2GEO: Sand
        - FHWA-NHI-10-016 Ch. 9: Drilled shafts
    """

    @staticmethod
    def soft_clay_matlock(
        depth_in: float,
        su_ksi: float,
        eps50: float,
        diameter_in: float,
        gamma_pci: float,
        sigma_v_ksi: float | None = None,
    ) -> dict:
        """Matlock (1970) soft clay p-y curve.

        pu increases from 3*su*D at surface to 9*su*D at depth.
        Transition depth: J*su*D / (gamma*D + J*su) where J ≈ 0.5 (soft) or 0.25.

        Args:
            depth_in: Depth below ground (in).
            su_ksi: Undrained shear strength (ksi).
            eps50: Strain at 50% strength.
            diameter_in: Shaft/pile diameter (in).
            gamma_pci: Effective unit weight (pci).
            sigma_v_ksi: Effective overburden (ksi), computed if No

    def soft_clay_matlock(
        depth_in: float,
        su_ksi: float,
        eps50: float,
        diameter_in: float,
        gamma_pci: float,
        sigma_v_ksi: float | None = None,
    ) -> dict:
        """Matlock (1970) soft clay p-y curve.

        pu increases from 3*su*D at surface to 9*su*D at depth.
        Transition depth: J*su*D / (gamma*D + J*su) where J ≈ 0.5 (soft) or 0.25.

        Args:
            depth_in: Depth below ground (in).
            su_ksi: Undrained shear strength (ksi).
            eps50: Strain at 50% strength.
            diameter_in: Shaft/pile diameter (in).
            gamma_pci: Effective unit weight (pci).
            sigma_v_ksi: Effective overburden (ksi), computed if None.

        Returns:
            dict with pu (kip/in), y50 (in), curve points [(y, p), ...].
        """
        D = diameter_in
        J = 0.5  # Matlock recommended for soft clay

        if sigma_v_ksi is None:
            sigma_v_ksi = gamma_pci * depth_in

      

    def stiff_clay_reese(
        depth_in: float,
        su_ksi: float,
        eps50: float,
        diameter_in: float,
        gamma_pci: float,
        sigma_v_ksi: float | None = None,
    ) -> dict:
        """Reese et al. (1975) stiff clay p-y curve (no free water).

        Args:
            depth_in: Depth below ground (in).
            su_ksi: Undrained shear strength (ksi).
            eps50: Strain at 50% strength.
            diameter_in: Shaft/pile diameter (in).
            gamma_pci: Effective unit weight (pci).

        Returns:
            dict with pu, y50, curve points.
        """
        D = diameter_in

        if sigma_v_ksi is None:
            sigma_v_ksi = gamma_pci * depth_in

        # Ultimate resistance — Reese formulation
        # Shallow wedge failure
        ca = su_ksi  # adhesion ≈ su for stiff clay
        pu_wedge = (2.0 * ca * D
                    + sigma_v_ksi * D
                    + 2.83 * su_ksi * depth_in)
        # Deep flow-around fail

    def sand_api(
        depth_in: float,
        phi_deg: float,
        diameter_in: float,
        gamma_pci: float,
        k_py_pci: float,
    ) -> dict:
        """API RP 2GEO sand p-y curve.

        pu from passive pressure coefficients; initial modulus k*z.

        Args:
            depth_in: Depth below ground (in).
            phi_deg: Friction angle (degrees).
            diameter_in: Shaft/pile diameter (in).
            gamma_pci: Effective unit weight (pci).
            k_py_pci: Initial modulus of subgrade reaction (pci).

        Returns:
            dict with pu, k_initial, curve points.
        """
        D = diameter_in
        phi = math.radians(phi_deg)

        # Coefficients per API
        beta_angle = math.pi / 4.0 + phi / 2.0
        Kp = math.tan(beta_angle) ** 2
        Ka = math.tan(math.pi / 4.0 - phi / 2.0) ** 2
        K0 = 1.0 - math.sin(phi)

        alpha_a = phi / 2.0

        # Shallow failure (wedge)
        C1 = (Kp - Ka) * math.tan(phi) + K0

    def rock(
        depth_in: float,
        qu_ksi: float,
        diameter_in: float,
    ) -> dict:
        """Simplified rock p-y curve.

        pu = qu * D (unconfined compressive strength × diameter).
        Linear-elastic up to yield, then constant.

        Args:
            depth_in: Depth below ground (in).
            qu_ksi: Unconfined compressive strength (ksi).
            diameter_in: Shaft/pile diameter (in).

        Returns:
            dict with pu, curve points.
        """
        D = diameter_in
        pu = qu_ksi * D

        # Initial stiffness: 100*qu*D / D = 100*qu per unit length
        k_rock = 100.0 * qu_ksi
        y_yield = pu / k_rock if k_rock > 0 else 0.01

        points = [
            (0.0, 0.0),
            (y_yield * 0.5, pu * 0.5),
            (y_yield, pu),
            (y_yield * 5.0, pu),
            (y_yield * 20.0, pu),
        ]

        return {"pu": pu, "y_yield": y_yield, "points": points}

class TZCurve:
    """t-z curve computation for skin friction resistance.

    References:
        - API RP 2GEO
        - AASHTO 10th Ed. Section 10.8
        - FHWA-NHI-10-016 Ch. 13
    """

    @staticmethod
    def clay_alpha(
        su_ksi: float,
        diameter_in: float,
        spacing_in: float,
    ) -> dict:
        """Alpha method for clay skin friction.

        tult = alpha * su (per unit length = alpha * su * pi * D * dz)
        alpha from API: alpha = 0.5*(su/sigma_v')^-0.5 for su/sigma_v' <= 1
                        alpha = 0.5*(su/sigma_v')^-0.25 for su/sigma_v' > 1
        Simplified: alpha ≈ min(1.0, 0.55 - 0.1*(su_ksf - 1.5)) per FHWA

        Args:
            su_ksi: Undrained shear strength (ksi).
            diameter_in: Shaft diameter (in).
            spacing_in: Spring spacing (in).

        Returns:
            dict with tult (kip/in), z_peak (in), curve points [(z, t)].
        """
        su_ksf = su_ksi / KSF_TO_KSI

        # Alpha per AASHTO/FHWA

    def clay_alpha(
        su_ksi: float,
        diameter_in: float,
        spacing_in: float,
    ) -> dict:
        """Alpha method for clay skin friction.

        tult = alpha * su (per unit length = alpha * su * pi * D * dz)
        alpha from API: alpha = 0.5*(su/sigma_v')^-0.5 for su/sigma_v' <= 1
                        alpha = 0.5*(su/sigma_v')^-0.25 for su/sigma_v' > 1
        Simplified: alpha ≈ min(1.0, 0.55 - 0.1*(su_ksf - 1.5)) per FHWA

        Args:
            su_ksi: Undrained shear strength (ksi).
            diameter_in: Shaft diameter (in).
            spacing_in: Spring spacing (in).

        Returns:
            dict with tult (kip/in), z_peak (in), curve points [(z, t)].
        """
        su_ksf = su_ksi / KSF_TO_KSI

        # Alpha per AASHTO/FHWA simplified
        if su_ksf <= 1.5:
            alpha = 0.55
        else:
            alpha = max(0.25, 0.55 - 0.1 * (su_ksf - 1.5))
        alpha = min(1.0, alpha)

        # Unit skin friction
        fs_ksi

    def sand_beta(
        sigma_v_ksi: float,
        phi_deg: float,
        diameter_in: float,
        spacing_in: float,
        depth_in: float,
    ) -> dict:
        """Beta method for sand skin friction.

        tult = beta * sigma_v' (per unit area)
        beta from AASHTO Table 10.8.3.5.2b-1 (simplified).

        Args:
            sigma_v_ksi: Effective vertical stress (ksi).
            phi_deg: Friction angle (degrees).
            diameter_in: Shaft diameter (in).
            spacing_in: Spring spacing (in).
            depth_in: Depth below ground (in).

        Returns:
            dict with tult (kip/in), z_peak (in), curve points.
        """
        # Beta from AASHTO — increases with friction angle, decreases with depth
        depth_ft = depth_in / FT_TO_IN
        N_corr = depth_ft  # simplified depth correction

        if phi_deg <= 25:
            beta_base = 0.25
        elif phi_deg <= 30:
            beta_base = 0.35
        elif phi_deg <= 35:
          

class QZCurve:
    """Q-z curve computation for end bearing resistance.

    References:
        - Reese & O'Neill (1988): Drilled shafts
        - AASHTO 10th Ed. Section 10.8.3.5
    """

    @staticmethod
    def clay(su_ksi: float, diameter_in: float) -> dict:
        """Clay end bearing Q-z curve.

        qult = Nc * su, Nc ≈ 9 per AASHTO/Reese & O'Neill.

        Args:
            su_ksi: Undrained shear strength at tip (ksi).
            diameter_in: Shaft diameter (in).

        Returns:
            dict with Qult (kip), z_peak (in), curve points [(z, Q)].
        """
        Nc = 9.0
        area = math.pi * diameter_in ** 2 / 4.0
        qult = Nc * su_ksi
        Qult = qult * area  # Total end bearing (kip)

        z_peak = 0.05 * diameter_in  # 5% of diameter per Reese & O'Neill

        points = [
            (0.0, 0.0),
            (0.002 * diameter_in, 0.15 * Qult),
            (0.013 * diameter_in, 0.50 * Qult),
            (0.042 * diameter_in, 0.85 * Qult),
       

    def clay(su_ksi: float, diameter_in: float) -> dict:
        """Clay end bearing Q-z curve.

        qult = Nc * su, Nc ≈ 9 per AASHTO/Reese & O'Neill.

        Args:
            su_ksi: Undrained shear strength at tip (ksi).
            diameter_in: Shaft diameter (in).

        Returns:
            dict with Qult (kip), z_peak (in), curve points [(z, Q)].
        """
        Nc = 9.0
        area = math.pi * diameter_in ** 2 / 4.0
        qult = Nc * su_ksi
        Qult = qult * area  # Total end bearing (kip)

        z_peak = 0.05 * diameter_in  # 5% of diameter per Reese & O'Neill

        points = [
            (0.0, 0.0),
            (0.002 * diameter_in, 0.15 * Qult),
            (0.013 * diameter_in, 0.50 * Qult),
            (0.042 * diameter_in, 0.85 * Qult),
            (z_peak, Qult),
            (0.10 * diameter_in, Qult),
        ]

        return {"Qult": Qult, "z_peak": z_peak, "Nc": Nc, "points": points}

    def sand(sigma_v_ksi: float, phi_deg: float, diameter_in: float) -> dict:
        """Sand end bearing Q-z curve.

        qult = Nq * sigma_v', Nq from AASHTO tables.

        Args:
            sigma_v_ksi: Effective vertical stress at tip (ksi).
            phi_deg: Friction angle at tip (degrees).
            diameter_in: Shaft diameter (in).

        Returns:
            dict with Qult (kip), z_peak (in), Nq, curve points.
        """
        # Nq from AASHTO (simplified Meyerhof/Berezantzev)
        if phi_deg <= 25:
            Nq = 8.0
        elif phi_deg <= 28:
            Nq = 12.0
        elif phi_deg <= 30:
            Nq = 20.0
        elif phi_deg <= 33:
            Nq = 30.0
        elif phi_deg <= 36:
            Nq = 50.0
        elif phi_deg <= 40:
            Nq = 80.0
        else:
            Nq = 120.0

        area = math.pi * diameter_in ** 2 / 4.0

        # Limit qult to 100 ksf ≈ 0.694 ksi per AASHTO
        qult = min(Nq * sigma_v_ksi, 0.694)
        Qult

    def rock(qu_ksi: float, diameter_in: float) -> dict:
        """Rock end bearing Q-z curve.

        qult = 2.5 * qu per FHWA-NHI-10-016.

        Args:
            qu_ksi: Unconfined compressive strength (ksi).
            diameter_in: Shaft diameter (in).

        Returns:
            dict with Qult (kip), z_peak (in), curve points.
        """
        area = math.pi * diameter_in ** 2 / 4.0
        qult = 2.5 * qu_ksi
        Qult = qult * area

        z_peak = 0.01 * diameter_in  # Rock — very stiff

        points = [
            (0.0, 0.0),
            (z_peak * 0.25, 0.40 * Qult),
            (z_peak * 0.50, 0.70 * Qult),
            (z_peak, Qult),
            (z_peak * 5.0, Qult),
        ]

        return {"Qult": Qult, "z_peak": z_peak, "points": points}

class TagAllocator:
    """Simple sequential tag allocator for OpenSees objects."""

    def __init__(self, start: int = 1):
        self._next = start

    def next(self, count: int = 1) -> int | list[int]:
        if count == 1:
            tag = self._next
            self._next += 1
            return tag
        tags = list(range(self._next, self._next + count))
        self._next += count
        return tags

    @property
    def current(self) -> int:
        return self._next

    def next(self, count: int = 1) -> int | list[int]:
        if count == 1:
            tag = self._next
            self._next += 1
            return tag
        tags = list(range(self._next, self._next + count))
        self._next += count
        return tags

def build_circular_section(
    tag_alloc: TagAllocator,
    diameter_in: float,
    fc_ksi: float,
    fy_ksi: float = 60.0,
    n_bars: int = 12,
    bar_size: str = "#10",
    cover_in: float = 3.0,
) -> tuple[list[dict], int]:
    """Build OpenSees fiber section for circular RC section.

    Uses Concrete01 (confined/unconfined) and Steel02 materials.

    Returns:
        (materials_list, section_tag)
    """
    materials = []

    # Concrete — Mander confined model (simplified)
    # Unconfined concrete
    fc = fc_ksi
    ec0 = 0.002
    fcu = 0.0  # spalling
    ecu = 0.005

    unconf_tag = tag_alloc.next()
    materials.append({
        "tag": unconf_tag,
        "type": "Concrete01",
        "params": [-fc, -ec0, -fcu, -ecu],
        "description": f"Unconfined concrete f'c={fc} ksi",
    })

    # Confined concrete (simplified Mander — ~1.3x f'c for typical confinement)
    fcc = 1.3 * fc
    ecc0 = 0.004
    fccu = 0.2 * fcc
    eccu = 0.015

    conf_tag = tag_alloc.next

def _spring_depths(
    length_in: float,
    spacing_in: float = 12.0,
    scour_depth_in: float = 0.0,
) -> list[float]:
    """Generate spring depths along shaft, skipping scour zone.

    Args:
        length_in: Total embedded length (in).
        spacing_in: Spring spacing (in), default 12 (1 ft).
        scour_depth_in: Depth of scour below ground (in). Springs above
            this depth are omitted.

    Returns:
        List of depths (in) from ground surface where springs are placed.
    """
    depths = []
    d = spacing_in / 2.0  # First spring at half-spacing
    while d < length_in:
        if d >= scour_depth_in:
            depths.append(d)
        d += spacing_in
    return depths

def _build_py_spring(
    layer: SoilLayer,
    depth_in: float,
    diameter_in: float,
    spacing_in: float,
    site: SiteProfile,
    multiplier: float = 1.0,
) -> dict:
    """Build a single p-y spring at given depth.

    Selects formulation based on soil type. Returns curve data with
    multiplier applied for upper/lower bound analysis.
    """
    sigma_v = site.effective_vertical_stress(depth_in / FT_TO_IN)

    if layer.soil_type in (SoilType.SOFT_CLAY, "soft_clay"):
        curve = PYCurve.soft_clay_matlock(
            depth_in=depth_in,
            su_ksi=layer.su_ksi,
            eps50=layer.get_eps50(),
            diameter_in=diameter_in,
            gamma_pci=layer.gamma_pci,
            sigma_v_ksi=sigma_v,
        )
    elif layer.soil_type in (SoilType.STIFF_CLAY, "stiff_clay"):
        curve = PYCurve.stiff_clay_reese(
            depth_in=depth_in,
            su_ksi=layer.su_ksi,
            eps50=layer.get_eps50(),
            diameter_in=diameter_in,
        

def _build_tz_spring(
    layer: SoilLayer,
    depth_in: float,
    diameter_in: float,
    spacing_in: float,
    site: SiteProfile,
    multiplier: float = 1.0,
) -> dict:
    """Build a single t-z spring at given depth."""
    if layer.soil_type in (SoilType.SOFT_CLAY, SoilType.STIFF_CLAY, "soft_clay", "stiff_clay"):
        curve = TZCurve.clay_alpha(
            su_ksi=layer.su_ksi,
            diameter_in=diameter_in,
            spacing_in=spacing_in,
        )
    elif layer.soil_type in (SoilType.SAND, "sand"):
        sigma_v = site.effective_vertical_stress(depth_in / FT_TO_IN)
        curve = TZCurve.sand_beta(
            sigma_v_ksi=sigma_v,
            phi_deg=layer.phi_deg,
            diameter_in=diameter_in,
            spacing_in=spacing_in,
            depth_in=depth_in,
        )
    elif layer.soil_type in (SoilType.ROCK, "rock"):
        # Rock skin friction — use alpha method with high su
        su_equiv = layer.qu_ksi / 2.0  # su ≈ qu/2
        curve = TZCurv

def build_drilled_shaft(
    params: dict,
    site: SiteProfile,
    tag_alloc: TagAllocator | None = None,
    spring_spacing_ft: float = 1.0,
) -> FoundationModel:
    """Build complete drilled shaft foundation model.

    Args:
        params: {
            diameter_ft: float,        # Shaft diameter (ft)
            length_ft: float,          # Embedded length (ft)
            above_grade_ft: float,     # Extension above ground (ft), default 0
            fc_ksi: float,             # Concrete strength (ksi), default 4.0
            fy_ksi: float,             # Rebar yield strength (ksi), default 60
            n_bars: int,               # Number of longitudinal bars
            bar_size: str,             # e.g. "#10"
            cover_in: float,           # Clear cover (in), default 3
        }
        site: SiteProfile with soil layers and GWT.
        tag_alloc: Tag allocator (created if None).
        spring_spacing_ft: Spring spacing in ft (default 1.0).

    Returns:
        

def build_spread_footing(
    params: dict,
    site: SiteProfile,
    tag_alloc: TagAllocator | None = None,
) -> FoundationModel:
    """Build spread footing with Winkler spring model.

    Uses compression-only (ENT material) zeroLength springs on a grid.

    Args:
        params: {
            length_ft: float,          # Footing length in plan (ft)
            width_ft: float,           # Footing width in plan (ft)
            depth_ft: float,           # Embedment depth (ft)
            thickness_ft: float,       # Footing thickness (ft), default 3
            subgrade_modulus_kcf: float, # kcf (default 100)
            n_springs_l: int,          # Springs along length (default 5)
            n_springs_w: int,          # Springs along width (default 5)
        }
        site: SiteProfile.
        tag_alloc: Tag allocator.

    Returns:
        FoundationModel with translational + rotational stiffness.
    """
    if tag_alloc is None:
        tag_alloc = TagAllocator(start=2000)

def _get_p_multipliers(spacing_over_D: float, n_rows: int) -> list[float]:
    """Get p-multipliers for pile group per AASHTO 10.7.2.4.

    Interpolates between tabulated values.

    Args:
        spacing_over_D: Center-to-center spacing / pile diameter.
        n_rows: Number of rows in loading direction.

    Returns:
        List of p-multiplier per row [lead, 2nd, 3rd, ... trailing].
    """
    # Breakpoints
    breakpoints = [3.0, 4.0, 5.0, 6.0, 8.0]
    values_table = [
        [0.80, 0.40, 0.30, 0.30],
        [0.85, 0.55, 0.45, 0.40],
        [0.90, 0.65, 0.55, 0.50],
        [0.95, 0.75, 0.65, 0.60],
        [1.00, 0.90, 0.80, 0.75],
    ]

    s = max(3.0, min(8.0, spacing_over_D))

    # Find bracketing indices
    for i in range(len(breakpoints) - 1):
        if breakpoints[i] <= s <= breakpoints[i + 1]:
            frac = (s - breakpoints[i]) / (breakpoints[i + 1] - breakpoints[i])
            interp = [
                values_table[i][j] + frac * (values_table[i + 1][j