# How Bridges Work 4, Exercise 2: a Pratt truss with the deck on top
# Joints A..E on the bottom chord at x = 0, 4, 8, 12, 16 m; F, G, H on the top
# chord at (4, 3), (8, 3), (12, 3). Members: bottom AB BC CD DE; top FG GH; end
# posts AF and HE; verticals BF CG DH; diagonals FC and HC (sloping down towards
# the centre). A pin at A, a roller at E, and 50 kN downward at each of F, G, H.
# Print every member's force and say which are zero-force members. Which
# members are in tension and which in compression?

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": (12, 0), "E": (16, 0),
          "F": (4, 3), "G": (8, 3), "H": (12, 3)}
members = [("A", "B"), ("B", "C"), ("C", "D"), ("D", "E"), ("F", "G"), ("G", "H"),
           ("A", "F"), ("H", "E"), ("B", "F"), ("C", "G"), ("D", "H"), ("F", "C"), ("H", "C")]
loads = {"F": (0, -50), "G": (0, -50), "H": (0, -50)}
supports = [("A", "x"), ("A", "y"), ("E", "y")]
res = solve_truss(joints, members, loads, supports)
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:
# AB 100.0 T
# BC 100.0 T
# CD 100.0 T
# DE 100.0 T
# FG -133.33 C
# GH -133.33 C
# AF -125.0 C
# HE -125.0 C
# BF 0.0 -
# CG -50.0 C
# DH 0.0 -
# FC 41.67 T
# HC 41.67 T
# reactions [0.0, 75.0, 75.0]
#
# Top chords FG, GH and end posts AF, HE are in compression; the bottom chord
# AB..DE and the diagonals FC and HC are in tension; the vertical CG is in
# compression. That is the Pratt pattern with the deck on top: long diagonals
# in tension, short verticals in compression. BF and DH are zero-force members:
# at joints B and D three members meet, two of them in a straight line (the
# bottom chord) and there is no load, so the third must carry nothing. Zero-force
# members still matter: they brace the chords against buckling and carry load
# when a vehicle crosses over them.