import numpy as np
from scipy.linalg import eigh

# -----------------------------
# MATERIAL / SECTION PROPERTIES
# -----------------------------

E = 7.75e9       # Young modulus of AS 30% reinforced [Pa]
I = 2.64e-7      # computed moment of area [m^4]

# -----------------------------
# ELEMENT LENGTHS [m]
# -----------------------------

L = np.array([
    56.1,
    54.8,
    54.8,
    54.8,
    54.8,
    54.8,
    54.8,
    54.8,
    54.2
]) * 1e-3   # mm to m

# -----------------------------
# LUMPED MASSES AT NODES [kg]
# -----------------------------

masses = np.array([
    0.078,
    0.038,
    0.038,
    0.038,
    0.038,
    0.038,
    0.038,
    0.038,
    0.038,
    0.054
])

# -----------------------------
# EFFECTIVE RADII FOR ROTATIONAL INERTIA [m]
# blade CG approx. diameter = (64 + 79)/2 = 71.5 mm
# radius = 35.75 mm
# -----------------------------

r_eff = np.ones(len(masses)) * 35.75e-3

# -----------------------------
# GLOBAL MATRICES
# -----------------------------

n_nodes = len(masses)
ndof = 2 * n_nodes

K = np.zeros((ndof, ndof))
M = np.zeros((ndof, ndof))

# -----------------------------
# MASS MATRIX
# -----------------------------

for i in range(n_nodes):
    # translational mass
    M[2*i, 2*i] = masses[i]

    # rotational mass moment of inertia
    J = masses[i] * r_eff[i]**2
    M[2*i+1, 2*i+1] = J

# -----------------------------
# STIFFNESS MATRIX ASSEMBLY
# -----------------------------

for e in range(len(L)):
    Le = L[e]

    Ke = E * I * np.array([
        [ 12/Le**3,   6/Le**2,  -12/Le**3,   6/Le**2],
        [  6/Le**2,    4/Le,    -6/Le**2,    2/Le   ],
        [-12/Le**3,  -6/Le**2,   12/Le**3,  -6/Le**2],
        [  6/Le**2,    2/Le,    -6/Le**2,    4/Le   ]
    ])

    dofs = [2*e, 2*e+1, 2*(e+1), 2*(e+1)+1]
    K[np.ix_(dofs, dofs)] += Ke

# -----------------------------
# BOUNDARY CONDITIONS
# simply supported ends:
# first and last vertical displacement fixed, rotations free
# -----------------------------

fixed = [0, 2*(n_nodes-1)]
free = np.setdiff1d(np.arange(ndof), fixed)

Kred = K[np.ix_(free, free)]
Mred = M[np.ix_(free, free)]

# -----------------------------
# SOLVE NATURAL FREQUENCIES
# -----------------------------

eigvals, eigvecs = eigh(Kred, Mred)

eigvals = eigvals[eigvals > 0]

omega = np.sqrt(eigvals)        # rad/s
freq = omega / (2*np.pi)        # Hz

print("Natural frequencies [Hz]:")
for i, f in enumerate(freq):
    print(f"Mode {i+1}: {f:.2f} Hz")

rpm_min = 2000
rpm_max = 2600

f_motor_min = rpm_min / 60
f_motor_max = rpm_max / 60

print("\nMotor operating frequency range:")
print(f"{f_motor_min:.2f} Hz to {f_motor_max:.2f} Hz")
