```python
from sympy.physics.units import *
from sympy import *

# Rounding:
import decimal
from decimal import Decimal as DX
def iso_round(obj, pv, rounding=decimal.ROUND_HALF_EVEN):
    import sympy
    """
    Rounding acc. to DIN EN ISO 80000-1:2013-08
    place value = Rundestellenwert
    """
    assert pv in set([
        # place value   #  round to:
        100,            #  3rd last digit before decimal
        10,             #  2nd last
        1,              #  last
        0.1,            #  1st digit after decimal
        0.01,           #  2nd
        0.001,          #  3rd
        0.0001,         #  4th
        0.00001,        #  5th
        0.000001,       #  6th
        0.0000001,      #  7th
        0.00000001,     #  8th
        0.000000001,    #  9th
        0.0000000001,   # 10th
        ])
    try:
        tmp = DX(str(float(obj)))
        obj = tmp.quantize(DX(str(pv)), rounding=rounding)
    except:
        for i in range(len(obj)):
            tmp = DX(str(float(obj[i])))
            obj[i] = tmp.quantize(DX(str(pv)), rounding=rounding)
    return obj

# LateX:
kwargs = {}
kwargs["mat_str"] = "bmatrix"
kwargs["mat_delim"] = ""
# kwargs["symbol_names"] = {FB: "F^{\mathsf B}", }

# Units:
(k, M, G ) = ( 10**3, 10**6, 10**9 )
(mm, cm, deg) = ( m/1000, m/100, pi/180)
Newton = kg*m/s**2
Pa     = Newton/m**2
MPa    = M*Pa
GPa    = G*Pa
kN     = k*Newton

# ---

c1, c2, m1, m2 = var("c1, c2, m1, m2", positive=True)

c, mass = var("c, mass", positive=True)

K = Matrix([
[c1 + c2,  -c2],
[-c2,  c2],
])

M = Matrix([
[m1, 0],
[0, m2],
])

sub_list = [
    (m1,mass),
    (m2,2*mass),
    (c1,c),
    (c2,c),
]

# ξ = λ²
xi = var("xi")
w = var("omega")

pprint("\nK + λ² M:")
A = K + xi*M
pprint(A)

pprint("\nK + λ²M for (m1=m, m2=2m, c1=c, c2=c):")
A = A.subs(sub_list)
pprint(A)
# exit()
pprint("\nCharacteristic equation:")
eq = Eq(det(A))
pprint(eq)

sol_xi = solve(eq,xi)

# Init:
eigenfreqs = [0,0]
eigenvects = [Matrix([1,0]), Matrix([1,0])]

x1, x2 = var("x1, x2")
w = Matrix([x1, x2])

pprint("\n")
for i in range(len(sol_xi)):
    pprint("\n--- (Value, Vector)-Pair Number "+str(i+1)+" ---")
    pprint("\nValue λ² / (c/m):")
    xii = sol_xi[i]
    eigenfreqs[i] = sqrt(-xii)
    tmp = xii / (c/mass)
    pprint(tmp)
    Ai = A.subs(xi,xii)
    pprint("\nVector:")
    eq = Eq(Ai*w,Matrix([0,0]))
    sol = solve(eq, x2)
    pprint(sol)
    tmp = sol[x2]
    tmp = tmp.subs(x1, 1)
    eigenvects[i][1] = tmp

pprint("\nEigenfrequencies / sqrt(c/m) and resp. Eigenvector for x1=1:")
for i in range(2):
    pprint("\n--- (Value, Vector)-Pair Number "+str(i+1)+" ---")
    tmp = eigenfreqs[i] / sqrt(c/mass)
    pprint(tmp)
    tmp = iso_round(tmp,0.01)
    print("≈ ", end="")
    pprint(tmp)
    pprint("\n")
    pprint(eigenvects[i])

sub_list_ex = [
    (c,1*Newton/m),
    (mass, 1*kg),
    ]

# Set Me!
first = True  # first eigenmode
first = False # second eigenmode

if first:
    pprint("\nFirst Eigenmode:")
    tmp = eigenfreqs[0]
    fac = eigenvects[0][1]
else:
    pprint("\nSecond Eigenmode:")
    tmp = eigenfreqs[1]
    fac = eigenvects[1][1]

pprint("\nAngular Frequency / (1/s):")
tmp = tmp.subs(sub_list)
tmp = tmp.subs(sub_list_ex)
tmp /= (1/s)
w_ex = N(tmp,15)
tmp = iso_round(tmp,0.01)
pprint(tmp)
pprint("\nTime Period / s:")
tmp = 2*pi/w_ex
tmp = iso_round(tmp,0.01)
pprint(tmp)


from pylab import *
t_in_s = linspace(0,15,100)
cos_arg = w_ex  * t_in_s
x1_in_m = array([cos(float(x)) for x in cos_arg])
x2_in_m = fac * x1_in_m

tmp = plt.axis()
tmp = plt.grid()
tmp = plt.plot(t_in_s, x1_in_m, label=r"$x_1\,\, /  \,\, \mathrm{m}$")
tmp = plt.plot(t_in_s, x2_in_m, label=r"$x_2\,\, /  \,\, \mathrm{m}$")
tmp = plt.xlabel(r"$t  \,\, /  \,\, \mathrm{s}$")
tmp = plt.legend()
plt.show()


#
# K + λ² M:
# ⎡c₁ + c₂ + m₁⋅ξ     -c₂   ⎤
# ⎢                         ⎥
# ⎣     -c₂        c₂ + m₂⋅ξ⎦
#
# K + λ²M for (m1=m, m2=2m, c1=c, c2=c):
# ⎡2⋅c + mass⋅ξ       -c     ⎤
# ⎢                          ⎥
# ⎣     -c       c + 2⋅mass⋅ξ⎦
#
# Characteristic equation:
#    2
# - c  + (c + 2⋅mass⋅ξ)⋅(2⋅c + mass⋅ξ) = 0
#
#
# --- (Value, Vector)-Pair Number 1 ---
#
# Value λ² / (c/m):
#   5   √17
# - ─ + ───
#   4    4
#
# Vector:
# ⎧    x₁⋅(3 + √17)⎫
# ⎨x₂: ────────────⎬
# ⎩         4      ⎭
#
# --- (Value, Vector)-Pair Number 2 ---
#
# Value λ² / (c/m):
#   5   √17
# - ─ - ───
#   4    4
#
# Vector:
# ⎧    x₁⋅(-√17 + 3)⎫
# ⎨x₂: ─────────────⎬
# ⎩          4      ⎭
#
# Eigenfrequencies / sqrt(c/m) and resp. Eigenvector for x1=1:
#
# --- (Value, Vector)-Pair Number 1 ---
#   __________
# ╲╱ -√17 + 5
# ────────────
#      2
# ≈ 0.47
#
# ⎡   1   ⎤
# ⎢       ⎥
# ⎢3   √17⎥
# ⎢─ + ───⎥
# ⎣4    4 ⎦
# ≈
# ⎡1.0 ⎤
# ⎢    ⎥
# ⎣1.78⎦
#
# --- (Value, Vector)-Pair Number 2 ---
#   _________
# ╲╱ √17 + 5
# ───────────
#      2
# ≈ 1.51
#
# ⎡    1    ⎤
# ⎢         ⎥
# ⎢  √17   3⎥
# ⎢- ─── + ─⎥
# ⎣   4    4⎦
# ≈
# ⎡ 1.0 ⎤
# ⎢     ⎥
# ⎣-0.28⎦
#
# Example:
#
# Angular Frequency / (1/s):
# 0.47
#
# Time Period / s:
# 13.42

```
