Script: Engineering Reference Tables Generator

#!/usr/bin/env python3
"""
Generate engineering reference tables for the Recovery Library.

Produces tables-engineering.md (Doc #017) with:
- Beam deflection tables for common NZ timber sections (Radiata Pine)
- Pipe flow tables via Hazen-Williams equation
- Timber span tables for floor joists and rafters (NZS 3603 principles)
- Wire/cable ampacity for copper conductors (AS/NZS 3008 principles)
- ISO metric thread dimensions (M3–M24, coarse and fine pitch)
- Metric bolt proof loads (Grade 4.8, 8.8, 10.9, 12.9)

All structural/hydraulic values computed from first principles using the
formulae stated. Electrical and thread values use hardcoded reference data
derived from the cited standards.

Usage:
    python generate_engineering_tables.py [output_file]
    Default output: ../tables-engineering.md
"""

import math
import sys
import os
from datetime import datetime


# ---------------------------------------------------------------------------
# Material / section constants
# ---------------------------------------------------------------------------

# Radiata Pine modulus of elasticity (NZS 3603 Table 2.1, SG8 / MSG8)
E_PINE_GPA = 8.0          # GPa
E_PINE_PA  = E_PINE_GPA * 1e9  # Pa = N/m²

# Standard NZ timber sections: (breadth_mm, depth_mm)
# Sections listed as depth × breadth in common NZ usage; listed here as
# (breadth, depth) so depth is the strong-axis dimension.
NZ_SECTIONS = [
    (45,  90),
    (45, 140),
    (45, 190),
    (45, 240),
    (45, 290),
]


def second_moment(b_mm, d_mm):
    """Second moment of area I for a rectangular section (mm⁴ → m⁴)."""
    I_mm4 = (b_mm * d_mm**3) / 12.0
    return I_mm4 * 1e-12  # convert mm⁴ → m⁴


# ---------------------------------------------------------------------------
# Document header
# ---------------------------------------------------------------------------

def header():
    return f"""---
title: "Engineering Reference Tables"
subtitle: "Recovery Library — Doc #017"
date: "{datetime.now().strftime('%Y-%m-%d')}"
---

# Engineering Reference Tables

*Recovery Library — Doc #017*

**Scope:** Structural, hydraulic, electrical, and mechanical reference tables for
field engineering without specialist software. Values are computed from the formulae
stated or taken from the cited standards. Any calculation can be re-derived from
the governing equations given in each section preamble.

**Material basis:** Timber data uses Radiata Pine (SG8/MSG8 grade) per NZS 3603.
Hydraulic data uses Hazen-Williams (H-W). Electrical data follows AS/NZS 3008.1.1
principles. Thread data follows ISO 68-1 and ISO 261/262.

**Conservative use:** These tables give working-load deflections and capacities.
Apply appropriate safety factors per the relevant design standard. In the absence
of a licensed engineer, use a minimum factor of 2.0 on calculated loads for
life-safety applications.

---

"""


# ---------------------------------------------------------------------------
# Section 1 — Beam Deflection Tables
# ---------------------------------------------------------------------------

SPAN_M = list(range(2, 9))   # 2 m to 8 m inclusive

BEAM_CONFIGS = [
    ("Simply Supported — Central Point Load",
     "δ = PL³ / (48EI)",
     "Deflection per kN of point load (mm/kN)",
     lambda P, w, L, E, I: (P * L**3) / (48 * E * I) * 1000,  # mm per N→kN scaling below
     "point"),

    ("Simply Supported — Uniformly Distributed Load (UDL)",
     "δ = 5wL⁴ / (384EI)",
     "Deflection per kN/m of UDL (mm per kN/m)",
     lambda P, w, L, E, I: (5 * w * L**4) / (384 * E * I) * 1000,
     "udl"),

    ("Cantilever — End Point Load",
     "δ = PL³ / (3EI)",
     "Deflection per kN of point load (mm/kN)",
     lambda P, w, L, E, I: (P * L**3) / (3 * E * I) * 1000,
     "point"),

    ("Cantilever — Uniformly Distributed Load (UDL)",
     "δ = wL⁴ / (8EI)",
     "Deflection per kN/m of UDL (mm per kN/m)",
     lambda P, w, L, E, I: (w * L**4) / (8 * E * I) * 1000,
     "udl"),
]


def beam_deflection_tables():
    """Beam deflection tables for common NZ Radiata Pine sections."""
    out  = "## 1. Beam Deflection Tables\n\n"
    out += ("Radiata Pine, E = 8.0 GPa (NZS 3603 SG8/MSG8). "
            "Sections listed as depth × breadth (standard NZ notation).\n\n"
            "**Reading the tables:** For a point-load table, multiply the tabulated "
            "value by your actual load in kN to obtain mid-span deflection in mm. "
            "For a UDL table, multiply by load in kN/m.\n\n"
            "**Span/deflection ratio check:** L/360 is a common serviceability limit "
            "for floors (e.g. 5 m span → 13.9 mm limit). L/240 is typical for roofs.\n\n")

    E = E_PINE_PA

    for title, formula, col_label, calc_fn, load_type in BEAM_CONFIGS:
        out += f"### 1.{BEAM_CONFIGS.index((title, formula, col_label, calc_fn, load_type)) + 1}  {title}\n\n"
        out += f"Formula: {formula}  \n"
        out += f"Tabulated: {col_label}\n\n"

        # Build column headers — one column per section
        section_labels = [f"{d}×{b}" for b, d in NZ_SECTIONS]
        # Also show I values
        I_values = [second_moment(b, d) for b, d in NZ_SECTIONS]

        header_row  = "| Span (m) | " + " | ".join(section_labels) + " |\n"
        sep_row     = "|" + "---|" * (1 + len(NZ_SECTIONS)) + "\n"
        out += header_row + sep_row

        # Sub-header showing I values
        I_strs = [f"I={I*1e6:.1f}×10⁻⁶ m⁴" for I in I_values]
        out += "| *I (m⁴)* | " + " | ".join(I_strs) + " |\n"

        for L in SPAN_M:
            row = f"| {L} m"
            for (b, d), I in zip(NZ_SECTIONS, I_values):
                # For point load: P = 1000 N = 1 kN; deflection formula gives metres,
                # calc_fn already returns mm per kN (or mm per kN/m for UDL).
                if load_type == "point":
                    P = 1000.0   # 1 kN in Newtons
                    w = 0.0
                else:
                    P = 0.0
                    w = 1000.0   # 1 kN/m in N/m

                delta_mm_per_unit = calc_fn(P, w, L, E, I)
                row += f" | {delta_mm_per_unit:.2f}"
            row += " |\n"
            out += row

        out += "\n"

    out += (
        "> **Source:** Beam deflection formulae from Roark's Formulas for Stress "
        "and Strain (Young & Budynas, 8th ed., Table 3). "
        "Section properties from NZS 3603:1993 and standard metric timber sizes.\n\n"
    )
    out += "---\n\n"
    return out


# ---------------------------------------------------------------------------
# Section 2 — Pipe Flow Tables (Hazen-Williams)
# ---------------------------------------------------------------------------

# Hazen-Williams: Q = 0.849 × C × A × R^0.63 × S^0.54
# For a full circular pipe: R = D/4 (hydraulic radius), A = π D²/4
# Q in m³/s, D in m, S = slope (dimensionless, rise/run)
# Re-expressed per common engineering form below.

PIPE_DIAMETERS_MM = [25, 32, 40, 50, 65, 80, 100, 150, 200]
PIPE_GRADIENTS    = [("1:100", 1/100), ("1:200", 1/200), ("1:500", 1/500), ("1:1000", 1/1000)]
HW_C_VALUES       = [
    (150, "New smooth plastic (PVC/PE)"),
    (130, "New steel / galvanised"),
    (100, "Old steel / concrete"),
    ( 80, "Old cast iron / corroded steel"),
]


def hw_flow(D_m, C, S):
    """
    Hazen-Williams flow for a full circular pipe.
    Returns Q in L/s.
    D_m : internal diameter (metres)
    C   : Hazen-Williams roughness coefficient
    S   : hydraulic gradient (dimensionless)
    """
    A = math.pi * D_m**2 / 4.0       # cross-sectional area m²
    R = D_m / 4.0                     # hydraulic radius for full pipe, m
    Q_m3s = 0.849 * C * A * (R**0.63) * (S**0.54)
    return Q_m3s * 1000.0             # convert m³/s → L/s


def pipe_flow_tables():
    """Hazen-Williams pipe flow tables."""
    out  = "## 2. Pipe Flow Tables (Hazen-Williams)\n\n"
    out += ("Formula: Q = 0.849 × C × A × R^0.63 × S^0.54  \n"
            "where Q = flow (m³/s), C = roughness coefficient, "
            "A = pipe cross-section area (m²), R = hydraulic radius = D/4 (m), "
            "S = hydraulic gradient (m/m, i.e. head loss per unit length).  \n"
            "**Flow values below are in litres per second (L/s) for full-bore flow.**\n\n"
            "Velocity can be estimated as V = Q / A. "
            "Recommended maximum design velocity: 1.5–3.0 m/s for water supply; "
            "0.6–1.0 m/s minimum for self-cleansing in sewers.\n\n")

    for C, C_desc in HW_C_VALUES:
        out += f"### 2.{HW_C_VALUES.index((C, C_desc)) + 1}  C = {C} — {C_desc}\n\n"

        # Table: rows = diameters, cols = gradients
        grad_labels = [g[0] for g in PIPE_GRADIENTS]
        out += "| Diameter (mm) | " + " | ".join(grad_labels) + " |\n"
        out += "|" + "---|" * (1 + len(PIPE_GRADIENTS)) + "\n"
        out += "| *Internal dia* | " + " | ".join([f"*S={g[1]:.5f}*" for g in PIPE_GRADIENTS]) + " |\n"

        for D_mm in PIPE_DIAMETERS_MM:
            D_m = D_mm / 1000.0
            row = f"| {D_mm}"
            for _, S in PIPE_GRADIENTS:
                q = hw_flow(D_m, C, S)
                if q < 0.1:
                    row += f" | {q:.3f}"
                elif q < 10.0:
                    row += f" | {q:.2f}"
                else:
                    row += f" | {q:.1f}"
            row += " |\n"
            out += row

        out += "\n"

    # Velocity reference sub-table
    out += "### Pipe Velocity Reference\n\n"
    out += ("Full-bore velocity V = Q/A. At 1 L/s through each pipe size:\n\n")
    out += "| Diameter (mm) | Area (cm²) | V at 1 L/s (m/s) |\n"
    out += "|---|---|---|\n"
    for D_mm in PIPE_DIAMETERS_MM:
        D_m = D_mm / 1000.0
        A_m2 = math.pi * D_m**2 / 4.0
        A_cm2 = A_m2 * 1e4
        V = 0.001 / A_m2   # 1 L/s = 0.001 m³/s
        out += f"| {D_mm} | {A_cm2:.2f} | {V:.3f} |\n"

    out += "\n"
    out += ("> **Source:** Hazen-Williams equation per Mays, L.W. (ed.), "
            "*Water Distribution Systems Handbook* (McGraw-Hill, 2000), §2.4. "
            "C-values per Chadwick, A. et al., *Hydraulics in Civil and "
            "Environmental Engineering* (5th ed., 2013), Table 8.1.\n\n")
    out += "---\n\n"
    return out


# ---------------------------------------------------------------------------
# Section 3 — Timber Span Tables
# ---------------------------------------------------------------------------

# Simplified span calculation based on NZS 3603 principles.
# Maximum span governed by serviceability (deflection L/360) or by
# bending stress, whichever is more restrictive.
#
# For a simply supported beam under UDL w (kN/m):
#   Bending:      M = wL²/8  ≤  f_b × Z       → L = sqrt(8 f_b Z / w)
#   Deflection:   δ = 5wL⁴/(384EI) ≤ L/limit  → L = (384EI / (5w × limit))^(1/3)
#
# NZS 3603 Table 2.4 — SG8 Radiata Pine bending strength:
#   f_b = 18.0 MPa (characteristic), design f_b' ≈ 14.4 MPa with k factors
#   We use f_b_design = 14.0 MPa as a round conservative figure.
#
# Load cases:
#   Floor joists:  dead + live = 1.5 kPa
#   Rafters:       dead + live = 0.5 kPa
#
# The tributary load on a joist = pressure × joist spacing (m) → kN/m.

F_B_DESIGN_MPA = 14.0    # Design bending strength, MPa
DEFLECTION_LIMIT_FLOOR  = 360   # L/360
DEFLECTION_LIMIT_ROOF   = 240   # L/240

FLOOR_SPACINGS_M  = [0.40, 0.60]
RAFTER_SPACINGS_M = [0.60, 0.90]

FLOOR_SECTIONS = [(45, d) for d in [140, 190, 240, 290]]
RAFTER_SECTIONS = [(45, d) for d in [90, 140, 190, 240, 290]]


def max_span_bending(b_mm, d_mm, w_kN_m):
    """Maximum span limited by bending (m)."""
    Z_mm3 = (b_mm * d_mm**2) / 6.0
    Z_m3  = Z_mm3 * 1e-9
    M_max = F_B_DESIGN_MPA * 1e6 * Z_m3   # N·m
    # wL²/8 = M_max  → L = sqrt(8M/w)
    w_N_m = w_kN_m * 1000.0
    if w_N_m <= 0:
        return float('inf')
    L = math.sqrt(8.0 * M_max / w_N_m)
    return L


def max_span_deflection(b_mm, d_mm, w_kN_m, limit_ratio):
    """Maximum span limited by deflection (m)."""
    I = second_moment(b_mm, d_mm)   # m⁴
    E = E_PINE_PA
    w_N_m = w_kN_m * 1000.0
    if w_N_m <= 0:
        return float('inf')
    # δ = 5wL⁴/(384EI) ≤ L/limit
    # → L³ ≤ 384EI/(5w × limit)  → L = (384EI/(5w·limit))^(1/3)
    L = (384.0 * E * I / (5.0 * w_N_m * limit_ratio)) ** (1.0 / 3.0)
    return L


def governed_span(b_mm, d_mm, w_kN_m, limit_ratio):
    """Return the governing (minimum) span in m, rounded down to nearest 50 mm."""
    L_bend = max_span_bending(b_mm, d_mm, w_kN_m)
    L_defl = max_span_deflection(b_mm, d_mm, w_kN_m, limit_ratio)
    L = min(L_bend, L_defl)
    # Round DOWN to nearest 50 mm (0.05 m) — conservative
    L_rounded = math.floor(L / 0.05) * 0.05
    return L_rounded, ("bending" if L_bend < L_defl else "deflection")


def timber_span_tables():
    out  = "## 3. Timber Span Tables\n\n"
    out += ("Material: Radiata Pine, MSG8/SG8 grade, NZS 3603:1993.  \n"
            "Design bending strength f'b = 14.0 MPa (characteristic 18.0 MPa × "
            "modification factors k₁=0.8 for long-term load duration, "
            "k₄=1.0 for dry service, k₈=1.0 for single member — conservative).  \n"
            "E = 8.0 GPa.  \n"
            "Spans governed by lesser of bending (M ≤ f'b·Z, wL²/8) or deflection "
            "(L/360 floor, L/240 roof). Values rounded DOWN to nearest 50 mm.\n\n"
            "**Load basis:** Floor = 1.5 kPa (dead + live). "
            "Roof/rafter = 0.5 kPa (dead + live, light roofing).  \n"
            "Tributary load per joist or rafter (kN/m) = pressure × spacing.\n\n")

    # 3.1 Floor joists
    out += "### 3.1  Floor Joists — Design Load 1.5 kPa\n\n"
    for spacing in FLOOR_SPACINGS_M:
        w = 1.5 * spacing   # kN/m on the joist
        out += f"#### Joist spacing: {int(spacing*1000)} mm centres  (w = {w:.3f} kN/m)\n\n"
        out += "| Section (d×b mm) | Max Span (m) | Governed by |\n"
        out += "|---|---|---|\n"
        for b, d in FLOOR_SECTIONS:
            L, gov = governed_span(b, d, w, DEFLECTION_LIMIT_FLOOR)
            out += f"| {d}×{b} | {L:.2f} | {gov} |\n"
        out += "\n"

    # 3.2 Rafters
    out += "### 3.2  Rafters — Design Load 0.5 kPa\n\n"
    for spacing in RAFTER_SPACINGS_M:
        w = 0.5 * spacing   # kN/m on the rafter
        out += f"#### Rafter spacing: {int(spacing*1000)} mm centres  (w = {w:.3f} kN/m)\n\n"
        out += "| Section (d×b mm) | Max Span (m) | Governed by |\n"
        out += "|---|---|---|\n"
        for b, d in RAFTER_SECTIONS:
            L, gov = governed_span(b, d, w, DEFLECTION_LIMIT_ROOF)
            out += f"| {d}×{b} | {L:.2f} | {gov} |\n"
        out += "\n"

    out += ("> **Note:** These spans are indicative for initial sizing. "
            "Final design must account for lateral restraint, notching, "
            "end fixity, load sharing (k₈ > 1.0 for multiple members), "
            "and connection capacity per NZS 3603:1993. "
            "Spans may be extended by 5–15% when load sharing applies.\n\n"
            "> **Source:** NZS 3603:1993 *Timber Structures Standard*, "
            "Tables 2.1, 2.4, B1–B6. Deflection limits per NZS 3604:2011 §6.\n\n")
    out += "---\n\n"
    return out


# ---------------------------------------------------------------------------
# Section 4 — Wire/Cable Ampacity
# ---------------------------------------------------------------------------

# Current-carrying capacity for copper conductors.
# Values derived from AS/NZS 3008.1.1:2017 Table 7 (thermoplastic insulation,
# 75°C rating, ambient 40°C for conduit, 30°C ambient for clipped/free air,
# single circuit basis).
#
# These are reference values. Derating for grouping, high ambient temperature,
# and thermal insulation must be applied in installation design.

# (conductor_mm2, in_conduit_A, clipped_direct_A, free_air_A)
AMPACITY_DATA = [
    ( 1.0,  10,  13,  15),
    ( 1.5,  13,  17,  20),
    ( 2.5,  18,  24,  27),
    ( 4.0,  24,  32,  36),
    ( 6.0,  31,  41,  46),
    (10.0,  42,  57,  63),
    (16.0,  57,  76,  85),
    (25.0,  75, 101, 112),
]


def ampacity_tables():
    out  = "## 4. Wire and Cable Ampacity — Copper Conductors\n\n"
    out += ("Current-carrying capacity (amperes) for copper conductors with "
            "thermoplastic (PVC/XLPE) insulation, 75°C conductor temperature rating.  \n"
            "Reference standard: AS/NZS 3008.1.1:2017, Table 7.  \n\n"
            "**Installation methods:**  \n"
            "- *In conduit:* conductors in enclosed conduit, 40°C ambient  \n"
            "- *Clipped direct:* surface-mounted, 30°C ambient  \n"
            "- *Free air:* spaced away from surface, 30°C ambient  \n\n"
            "**Voltage drop** (approx, 230V single-phase): "
            "V_drop = I × R × L / 1000 where R (mΩ/m) from table below.\n\n")

    out += "| Size (mm²) | In Conduit (A) | Clipped Direct (A) | Free Air (A) | R at 75°C (mΩ/m) |\n"
    out += "|---|---|---|---|---|\n"

    # Resistance values: R_75 = ρ_Cu × (1 + 0.00393×(75-20)) / A
    # ρ_Cu at 20°C = 17.2 nΩ·m
    rho_20 = 17.2e-9   # Ω·m
    alpha  = 0.00393   # /°C for copper

    for csa, i_conduit, i_clip, i_air in AMPACITY_DATA:
        A_m2 = csa * 1e-6
        R_20 = rho_20 / A_m2        # Ω/m
        R_75 = R_20 * (1 + alpha * (75 - 20))
        R_mOhm_m = R_75 * 1000      # mΩ/m
        out += f"| {csa:.1f} | {i_conduit} | {i_clip} | {i_air} | {R_mOhm_m:.2f} |\n"

    out += "\n"
    out += "### 4.1  Derating Factors\n\n"
    out += ("Apply derating factors by multiplying tabulated current by all applicable factors.\n\n")
    out += "| Condition | Derating Factor |\n"
    out += "|---|---|\n"
    derating = [
        ("2 circuits grouped (touching)", "0.80"),
        ("3 circuits grouped", "0.70"),
        ("4–6 circuits grouped", "0.60"),
        ("7–9 circuits grouped", "0.50"),
        ("Ambient 45°C (vs 30°C base)", "0.91"),
        ("Ambient 50°C", "0.82"),
        ("Ambient 55°C", "0.71"),
        ("Ambient 60°C", "0.58"),
        ("Fully enclosed in thermal insulation", "0.50"),
        ("Partially covered by thermal insulation", "0.75"),
    ]
    for cond, factor in derating:
        out += f"| {cond} | {factor} |\n"

    out += "\n"
    out += "### 4.2  Common Voltage Drop Budget (230 V, ≤5% limit = 11.5 V)\n\n"
    out += ("Maximum single-phase circuit length (m) for 5% voltage drop at full load:  \n"
            "L_max = 11.5 × 1000 / (I × R_mΩm × 2)   (factor 2 for return conductor)\n\n")
    out += "| Size (mm²) | 10 A load | 16 A load | 20 A load | 32 A load |\n"
    out += "|---|---|---|---|---|\n"

    load_currents = [10, 16, 20, 32]
    for csa, i_conduit, i_clip, i_air in AMPACITY_DATA:
        A_m2 = csa * 1e-6
        R_20 = rho_20 / A_m2
        R_75 = R_20 * (1 + alpha * (75 - 20))
        R_mOhm_m = R_75 * 1000
        row = f"| {csa:.1f}"
        for I in load_currents:
            L_max = 11500 / (I * R_mOhm_m * 2)
            row += f" | {L_max:.0f} m"
        row += " |\n"
        out += row

    out += "\n"
    out += ("> **Source:** AS/NZS 3008.1.1:2017 *Electrical Installations — "
            "Selection of Cables*, Table 7 (current capacity) and Table 30 "
            "(resistance). Derating factors from Tables 22–25 of the same standard.\n\n")
    out += "---\n\n"
    return out


# ---------------------------------------------------------------------------
# Section 5 — ISO Metric Thread Dimensions
# ---------------------------------------------------------------------------

# ISO 68-1 defines the thread profile. Relationships:
#   Minor diameter (external) = D - 1.2269 × P
#   Minor diameter (internal) = D - 1.0825 × P
#   Pitch diameter            = D - 0.6495 × P
#   Tap drill (75% thread)    ≈ D - P   (gives ~75% thread engagement)
#
# Standard coarse pitches per ISO 261; fine pitches per ISO 262.

THREADS = [
    # (M_size, pitch_coarse, pitch_fine_or_None)
    ( 3,  0.50, None),
    ( 4,  0.70, None),
    ( 5,  0.80, None),
    ( 6,  1.00, 0.75),
    ( 8,  1.25, 1.00),
    (10,  1.50, 1.25),
    (12,  1.75, 1.25),
    (14,  2.00, 1.50),
    (16,  2.00, 1.50),
    (18,  2.50, 1.50),
    (20,  2.50, 1.50),
    (22,  2.50, 1.50),
    (24,  3.00, 2.00),
]


def thread_dims(D, P):
    """Return (pitch_dia, minor_dia_ext, minor_dia_int, tap_drill) in mm."""
    pitch_dia    = D - 0.6495 * P
    minor_ext    = D - 1.2269 * P   # external (bolt) minor diameter
    minor_int    = D - 1.0825 * P   # internal (nut) minor diameter
    tap_drill    = D - P            # nominal 75% thread engagement
    return pitch_dia, minor_ext, minor_int, tap_drill


def thread_tables():
    out  = "## 5. ISO Metric Thread Dimensions\n\n"
    out += ("Based on ISO 68-1 (thread profile), ISO 261 (coarse series), "
            "ISO 262 (fine series).  \n\n"
            "Formulae (all dimensions in mm):  \n"
            "- Pitch diameter  d₂ = D − 0.6495 × P  \n"
            "- Minor dia (bolt) d₁ = D − 1.2269 × P  \n"
            "- Minor dia (nut)  D₁ = D − 1.0825 × P  \n"
            "- Tap drill (75% thread) ≈ D − P  \n\n"
            "For class 6g/6H (standard commercial tolerance), add tolerance "
            "from ISO 965-1.\n\n")

    out += "### 5.1  Coarse Pitch Series\n\n"
    out += "| Size | Pitch (mm) | Pitch Dia d₂ (mm) | Minor (bolt) d₁ (mm) | Minor (nut) D₁ (mm) | Tap Drill (mm) |\n"
    out += "|---|---|---|---|---|---|\n"
    for D, P_c, _ in THREADS:
        pd, m_ext, m_int, td = thread_dims(D, P_c)
        out += f"| M{D} | {P_c:.2f} | {pd:.3f} | {m_ext:.3f} | {m_int:.3f} | {td:.2f} |\n"

    out += "\n### 5.2  Fine Pitch Series\n\n"
    out += "| Size | Pitch (mm) | Pitch Dia d₂ (mm) | Minor (bolt) d₁ (mm) | Minor (nut) D₁ (mm) | Tap Drill (mm) |\n"
    out += "|---|---|---|---|---|---|\n"
    for D, _, P_f in THREADS:
        if P_f is None:
            continue
        pd, m_ext, m_int, td = thread_dims(D, P_f)
        out += f"| M{D} | {P_f:.2f} | {pd:.3f} | {m_ext:.3f} | {m_int:.3f} | {td:.2f} |\n"

    out += "\n### 5.3  Thread Engagement and Strength Notes\n\n"
    out += ("- Minimum thread engagement in steel: 1.0 × D  \n"
            "- Minimum thread engagement in aluminium: 1.5 × D  \n"
            "- Minimum thread engagement in cast iron / grey iron: 1.25 × D  \n"
            "- Minimum thread engagement in plastics: 2.0 × D  \n"
            "- Standard nut height (coarse) ≈ 0.8 × D (grade 8 nuts ≈ 1.0 × D)  \n"
            "- Self-locking torque (nyloc inserts): add ~20–30% to standard torque\n\n")

    out += ("> **Source:** ISO 68-1:2013 *General purpose screw threads — Basic profile*. "
            "ISO 261:1998 *ISO general-purpose metric screw threads — General plan*. "
            "Tap drill guidance per Machinery's Handbook (31st ed., p. 1924).\n\n")
    out += "---\n\n"
    return out


# ---------------------------------------------------------------------------
# Section 6 — Bolt Proof Loads
# ---------------------------------------------------------------------------

# ISO 898-1 metric bolt mechanical properties.
# Proof load = proof stress × tensile stress area
#
# Tensile stress area (ISO 898-1 Annex A):
#   A_s = π/4 × ((d₂ + d₁)/2)² = π/4 × (D - 0.9382×P)²
#   where d₂ = pitch diameter, d₁ = minor diameter
#
# Proof stresses per ISO 898-1:2013 Table 5:
#   Grade 4.8: σ_p = 310 MPa
#   Grade 8.8: σ_p = 600 MPa (≤M16), 580 MPa (>M16) — use 580 conservative
#   Grade 10.9: σ_p = 830 MPa
#   Grade 12.9: σ_p = 970 MPa
#
# Tightening torque (indicative, friction coeff μ = 0.15, standard nut):
#   T ≈ 0.2 × F_p × D  (Nut factor K ≈ 0.2)

BOLT_SIZES = [
    # (M, P_coarse)
    ( 6, 1.00),
    ( 8, 1.25),
    (10, 1.50),
    (12, 1.75),
    (14, 2.00),
    (16, 2.00),
    (18, 2.50),
    (20, 2.50),
    (22, 2.50),
    (24, 3.00),
]

# (grade_label, proof_stress_MPa)
BOLT_GRADES = [
    ("4.8",  310),
    ("8.8",  580),
    ("10.9", 830),
    ("12.9", 970),
]


def tensile_stress_area(D, P):
    """ISO 898-1 tensile stress area in mm²."""
    d_s = D - 0.9382 * P   # effective diameter for stress area
    A_s = math.pi / 4.0 * d_s**2
    return A_s


def bolt_proof_load(D, P, sigma_p_MPa):
    """Proof load in kN."""
    A_s = tensile_stress_area(D, P)   # mm²
    F_p = sigma_p_MPa * A_s           # N (MPa × mm² = N)
    return F_p / 1000.0               # kN


def bolt_torque(F_p_kN, D_mm, K=0.20):
    """Indicative tightening torque in N·m."""
    F_N = F_p_kN * 1000.0
    T   = K * F_N * (D_mm / 1000.0)   # N·m
    return T


def bolt_strength_tables():
    out  = "## 6. Metric Bolt Proof Loads and Tightening Torques\n\n"
    out += ("Per ISO 898-1:2013 *Mechanical properties of fasteners — "
            "Bolts, screws and studs*.  \n\n"
            "**Proof load** = proof stress × tensile stress area.  \n"
            "Tensile stress area: A_s = π/4 × (D − 0.9382P)² (ISO 898-1 Annex A).  \n\n"
            "**Tightening torque** is indicative for lightly oiled or as-received "
            "uncoated steel (nut factor K ≈ 0.20). "
            "T ≈ K × F_proof × D. Actual torque varies with friction; "
            "always use measured or published torque for critical joints.\n\n"
            "**Proof stresses used (MPa):**  \n"
            "- Grade 4.8: 310 MPa  \n"
            "- Grade 8.8: 580 MPa (conservative; spec is 600 for M≤16, 580 for M>16)  \n"
            "- Grade 10.9: 830 MPa  \n"
            "- Grade 12.9: 970 MPa  \n\n")

    out += "### 6.1  Tensile Stress Areas\n\n"
    out += "| Size | Pitch (mm) | Stress Area A_s (mm²) |\n"
    out += "|---|---|---|\n"
    for D, P in BOLT_SIZES:
        A_s = tensile_stress_area(D, P)
        out += f"| M{D} | {P:.2f} | {A_s:.1f} |\n"
    out += "\n"

    out += "### 6.2  Proof Loads (kN)\n\n"
    grade_labels = [g[0] for g in BOLT_GRADES]
    out += "| Size | " + " | ".join([f"Grade {g}" for g in grade_labels]) + " |\n"
    out += "|" + "---|" * (1 + len(BOLT_GRADES)) + "\n"
    for D, P in BOLT_SIZES:
        row = f"| M{D}"
        for _, sigma_p in BOLT_GRADES:
            Fp = bolt_proof_load(D, P, sigma_p)
            row += f" | {Fp:.1f}"
        row += " |\n"
        out += row
    out += "\n"

    out += "### 6.3  Indicative Tightening Torques (N·m), K = 0.20\n\n"
    out += "| Size | " + " | ".join([f"Grade {g}" for g in grade_labels]) + " |\n"
    out += "|" + "---|" * (1 + len(BOLT_GRADES)) + "\n"
    for D, P in BOLT_SIZES:
        row = f"| M{D}"
        for _, sigma_p in BOLT_GRADES:
            Fp = bolt_proof_load(D, P, sigma_p)
            T  = bolt_torque(Fp, D)
            row += f" | {T:.0f}"
        row += " |\n"
        out += row
    out += "\n"

    out += "### 6.4  Bolt Identification — Head Markings\n\n"
    out += "| Grade | Yield/Tensile (MPa) | Typical use |\n"
    out += "|---|---|---|\n"
    markings = [
        ("4.8", "320 / 420", "General light structural, non-critical joints"),
        ("8.8", "640 / 800", "Standard structural bolts; most structural connections"),
        ("10.9","940 / 1040","High-strength; flanged connections, critical machinery"),
        ("12.9","1100 / 1220","Very high strength; socket-head cap screws, aerospace"),
    ]
    for grade, ys_ts, use in markings:
        out += f"| {grade} | {ys_ts} | {use} |\n"

    out += "\n"
    out += "### 6.5  Nut and Washer Selection\n\n"
    out += ("Match nut grade to bolt grade:  \n"
            "- Grade 4 nuts → Grade 4.8 bolts  \n"
            "- Grade 8 nuts → Grade 8.8 bolts  \n"
            "- Grade 10 nuts → Grade 10.9 bolts  \n"
            "- Grade 12 nuts → Grade 12.9 bolts  \n\n"
            "Hardened washers (Grade 300HV) required under head and nut for "
            "Grade 8.8 and above in steel structures. Plain washers acceptable "
            "for Grade 4.8 in timber and light construction.\n\n")

    out += ("> **Source:** ISO 898-1:2013 *Mechanical properties of fasteners — "
            "Bolts, screws and studs*. Torque estimation from Shigley's "
            "Mechanical Engineering Design (10th ed., §8-9, Nut Factor Table 8-15).\n\n")
    out += "---\n\n"
    return out


# ---------------------------------------------------------------------------
# Footer
# ---------------------------------------------------------------------------

def footer():
    return (
        "## Notes on Using These Tables\n\n"
        "1. **Check units carefully.** Most errors in structural calculations come "
        "from unit inconsistency. Verify that load units (kN vs N, kPa vs Pa) are "
        "carried correctly through every calculation.\n\n"
        "2. **Deflection ≠ strength.** A beam may pass a deflection limit at a "
        "shorter span than its bending capacity permits, or vice versa. Both "
        "criteria must be checked independently.\n\n"
        "3. **Hydraulic gradients.** The pipe flow tables assume the pipe is "
        "flowing full. Partially full flow (gravity drainage) requires Manning's "
        "equation or a hydraulic elements chart. For pressurised systems, add head "
        "losses at fittings (typically 10–30% additional equivalent length).\n\n"
        "4. **Electrical safety.** Ampacity tables assume a single circuit in "
        "good condition. Overcurrent protection must be rated ≤ the cable ampacity "
        "after all derating factors are applied. When in doubt, use the next larger "
        "cable size.\n\n"
        "5. **Thread lubrication.** Tightening torques change significantly with "
        "lubrication: dry threads may require 20–30% more torque than oiled threads "
        "for the same clamp force; zinc-plated threads are typically between dry "
        "and oiled.\n\n"
        "6. **Professional review.** These tables are for initial sizing, "
        "field estimation, and educational use. Life-safety structures "
        "(bridges, occupied buildings, pressure vessels) require review by a "
        "licensed Chartered Professional Engineer (CPEng) or equivalent.\n\n"
        "---\n\n"
        f"*Generated by `scripts/generate_engineering_tables.py` — "
        f"Recovery Library Doc #017. "
        f"Date: {datetime.now().strftime('%Y-%m-%d')}.*\n"
    )


# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------

def main():
    output_file = sys.argv[1] if len(sys.argv) > 1 else os.path.join(
        os.path.dirname(os.path.abspath(__file__)), "..", "tables-engineering.md"
    )

    content  = header()
    content += beam_deflection_tables()
    content += pipe_flow_tables()
    content += timber_span_tables()
    content += ampacity_tables()
    content += thread_tables()
    content += bolt_strength_tables()
    content += footer()

    with open(output_file, 'w', encoding='utf-8') as f:
        f.write(content)

    lines = content.count('\n')
    print(f"Written to {output_file}")
    print(f"Total lines: {lines}")


if __name__ == "__main__":
    main()