# How Bridges Work 4, Exercise 1: the Warren truss under a different load
# Use the truss of the chapter: joints A, B, C on the bottom chord at x = 0, 4,
# 8 m, and D, E on the top chord at (2, 3) and (6, 3); members AB, BC, DE, AD,
# DB, BE, EC. Support A is a pin, support C a roller. Solve it for
# (a) 60 kN downward at the top joint E, and
# (b) 40 kN at B and 60 kN at E together.
# Print each member's force (tension positive), T or C, and the reactions.
# Check that the vertical reactions add up to the total load.

import math

def solve_truss(joints, members, loads, supports):
    """Method of joints as one linear system. Tension is positive."""
    names = list(joints)
    n = len(members) + len(supports)
    rows = {}
    for j in names:
        rows[(j, "x")] = [0.0] * n
        rows[(j, "y")] = [0.0] * n
    for k, (a, b) in enumerate(members):
        (xa, ya), (xb, yb) = joints[a], joints[b]
        length = math.hypot(xb - xa, yb - ya)
        cx, cy = (xb - xa) / length, (yb - ya) / length
        rows[(a, "x")][k] += cx
        rows[(a, "y")][k] += cy
        rows[(b, "x")][k] -= cx
        rows[(b, "y")][k] -= cy
    for k, (j, d) in enumerate(supports):
        rows[(j, d)][len(members) + k] = 1.0
    keys = list(rows)
    a = [rows[key] + [-(loads.get(key[0], (0, 0))[0 if key[1] == "x" else 1])] for key in keys]
    for c in range(n):                                   # Gaussian elimination
        p = max(range(c, n), key=lambda r: abs(a[r][c]))
        a[c], a[p] = a[p], a[c]
        for r in range(n):
            if r != c:
                f = a[r][c] / a[c][c]
                a[r] = [x - f * y for x, y in zip(a[r], a[c])]
    return [a[i][n] / a[i][i] for i in range(n)]

joints = {"A": (0, 0), "B": (4, 0), "C": (8, 0), "D": (2, 3), "E": (6, 3)}
members = [("A", "B"), ("B", "C"), ("D", "E"), ("A", "D"), ("D", "B"), ("B", "E"), ("E", "C")]
supports = [("A", "x"), ("A", "y"), ("C", "y")]

for label, loads in (("a", {"E": (0, -60)}), ("b", {"B": (0, -40), "E": (0, -60)})):
    res = solve_truss(joints, members, loads, supports)
    print("case", label)
    for (p, q), f in zip(members, res):
        print(p + q, round(f, 2), "T" if f > 1e-9 else ("C" if f < -1e-9 else "-"))
    print("reactions", [round(r, 2) for r in res[len(members):]])

# Output when run:
# case a
# AB 10.0 T
# BC 30.0 T
# DE -20.0 C
# AD -18.03 C
# DB 18.03 T
# BE -18.03 C
# EC -54.08 C
# reactions [0.0, 15.0, 45.0]
# case b
# AB 23.33 T
# BC 43.33 T
# DE -46.67 C
# AD -42.06 C
# DB 42.06 T
# BE 6.01 T
# EC -78.12 C
# reactions [0.0, 35.0, 65.0]
#
# Checks: case (a) reactions 15 + 45 = 60 kN, the load; case (b) 35 + 65 = 100
# kN, the total load. The reaction nearer the load is larger (C is 2 m from the
# 60 kN load at E). Member BE changes sign between the cases: it is compressed
# (-18.03 kN) in (a) and stretched (+6.01 kN) in (b), because the 40 kN at B now
# pulls it down. In a truss the members' roles depend on where the load is, and
# a moving train or lorry can flip a member between tension and compression;
# designers must check both.