# vakwerk.py - staafkrachten van een statisch bepaald vakwerk (evenwicht per knoop, als stelsel: hoofdstuk 5)
import math
import numpy as np


def vakwerk(knopen, staven, steunen, lasten):
    """knopen {naam: (x, y)} in m; staven [(p, q)]; steunen [(knoop, richting)]; lasten {knoop: (Fx, Fy)} in kN.
    Geeft de staafkrachten (trek positief) en de reacties."""
    rij = {k: 2 * i for i, k in enumerate(knopen)}           # rij 2i: evenwicht in x, rij 2i+1: in y
    m, r, k = len(staven), len(steunen), len(knopen)
    if m + r != 2 * k:
        raise ValueError(f"m + r = {m + r} en 2k = {2 * k}: niet statisch bepaald")
    A = np.zeros((2 * k, m + r))
    for j, (p, q) in enumerate(staven):
        e = np.subtract(knopen[q], knopen[p])
        e = e / np.linalg.norm(e)                             # eenheidsvector van p naar q
        A[rij[p]:rij[p] + 2, j] += e                          # een trekkracht trekt p naar q ...
        A[rij[q]:rij[q] + 2, j] -= e                          # ... en q naar p
    for j, (knoop, richting) in enumerate(steunen):
        A[rij[knoop]:rij[knoop] + 2, m + j] += richting
    if np.linalg.matrix_rank(A) < 2 * k:
        raise ValueError("singuliere matrix: een mechanisme of een kritische vorm")
    F = np.zeros(2 * k)
    for knoop, f in lasten.items():
        F[rij[knoop]:rij[knoop] + 2] += f
    x = np.linalg.solve(A, -F)                                # A x + F = 0 in elke knoop
    return x[:m], x[m:]


# het dakspant van de woonmodule: Fink, 6,0 m, 30 graden, spanten om de 0,60 m
h = math.tan(math.radians(30))
knopen = {"A": (0, 0), "D": (1.5, 1.5 * h), "C": (3, 3 * h), "E": (4.5, 1.5 * h), "B": (6, 0), "F": (2, 0), "G": (4, 0)}
staven = [("A", "D"), ("D", "C"), ("C", "E"), ("E", "B"), ("A", "F"), ("F", "G"), ("G", "B"),
          ("D", "F"), ("F", "C"), ("C", "G"), ("G", "E")]
steunen = [("A", (1, 0)), ("A", (0, 1)), ("B", (0, 1))]      # scharnier in A, rol in B

# rekenwaarden 1,35 G + 1,5 S per spant, in kN per meter horizontaal
q_boven = 1.35 * 0.60 / math.cos(math.radians(30)) * 0.60 + 1.5 * 0.40 * 0.60
q_onder = 1.35 * 0.30 * 0.60
lasten = {"D": (0, -1.5 * q_boven), "C": (0, -1.5 * q_boven), "E": (0, -1.5 * q_boven),
          "F": (0, -2.0 * q_onder), "G": (0, -2.0 * q_onder),
          "A": (0, -0.75 * q_boven - q_onder), "B": (0, -0.75 * q_boven - q_onder)}

krachten, reacties = vakwerk(knopen, staven, steunen, lasten)
for (p, q), s in zip(staven, krachten):
    print(f"staaf {p}{q}: {s:6.2f} kN  {'trek' if s > 0 else 'druk'}")
print(f"reacties: Ay = {reacties[1]:.2f} kN, By = {reacties[2]:.2f} kN (Ax = 0: geen horizontale last)")