M = np.array(((2.0, 0.0), ( 0.0, 1.0)))
K = np.array(((3.0,-2.0), (-2.0, 2.0)))
p = np.array(( 0.0, 1.0)); w = 2.0

print_mat(M, pre='\\boldsymbol{M}=m\\,', fmt='%d')
print_mat(K, pre='\\boldsymbol{K}=k\\,', fmt='%d')
print_mat(p[:,None], pre=r'\boldsymbol{p}(t) = p_0\,', fmt='%d',
          post='\\sin(%d\\omega_0t)'%w, mt='B')

"""
Explanation: 2 DOF System
The data that describes our system and its loading:
End of explanation
"""


evals, Psi = eigh(K, M)
Mstar = Psi.T@M@Psi
Kstar = Psi.T@K@Psi
pstar = Psi.T@p

"""
Explanation: Computing the eigenvalues and the eigenvectors
Using the eigh library function, specialized for Hermitian (symmetric) matrices.
We compute also the starred matrices.  Note that the eigenvectors returned by eigh are normalized with respect to M, i.e., $M_i=m$ for every $i$.
End of explanation
"""


print_mat(evals[None,:], mt='p', pre=r'\omega^2_i=\omega^2_o\,')
print_mat(Psi, pre=r'\boldsymbol{\Psi}=')
print_mat(Mstar, pre=r'\boldsymbol{M}^\star=m\,')
print_mat(Kstar, pre=r'\boldsymbol{K}^\star=k\,')
print_mat(pstar[:,None], pre=r'\boldsymbol{p}^\star=p_o\,', mt='B')

"""
Explanation: The @ operator stands, in this context, for matrix multiplication.
We print the computed values.
End of explanation
"""


L  = np.sqrt(evals)
DAF = 1.0/(L**2-w**2)
beta = w/L
t = np.linspace(0,60,601)[:,None]
q = pstar*DAF*(np.sin(w*t)-beta*np.sin(L*t))

"""
Explanation: Modal Response
The modal equation of motion, divided by $m$
$$ \ddot q_i + \Lambda^2_i\omega_0^2 q_i = \frac{p_0}{k}\frac{k}{m}\,p^\star_i \sin2\omega_0t.$$
The particular integral $\xi_i(t)$ is given by
$$\xi_i = \Delta_{st} \frac{\omega_0^2}{(\Lambda^2_i-4)\omega_0^2}p^\star_i\sin2\omega_0t.$$
Eventually the modal response, computed in the hypotesis of initial rest conditions is
$$q_i = \frac{\Delta_{st}}{\Lambda^2_i-4}p^\star_i(\sin2\omega_0t-\beta\sin\Lambda_i\omega_0t).$$
End of explanation
"""


curves = plt.plot(t,q)
plt.legend(curves,['q1', 'q2'])
plt.title('Modal Response')
plt.xlabel('$\omega_0t$')
plt.ylabel('$q_i/\Delta_{st}$');

x = (Psi@q.T).T

curves = plt.plot(t, x)
plt.legend(curves,['x1', 'x2'])
plt.title('Structural Response')
plt.xlabel('$\omega_0t$')
plt.ylabel('$X_i/\Delta_{st}$');

"""
Explanation: The definition of time vector is a bit complicated...
End of explanation
"""


%matplotlib inline

import matplotlib.pyplot as plt
plt.style.use(['fivethirtyeight', '00_mplrc'])

import numpy as np ; from scipy.linalg import eigh
np.set_printoptions(suppress=False, linewidth=120)

from IPython.display import Latex

def print_mat(mat, pre='', post='', fmt='%.6f', mt='b'):
    display(Latex(
        r'$$' + pre + r'\begin{%smatrix}'%mt +
        r'\\'.join('&'.join(fmt%x for x in row) for row in mat) +
        r'\end{%smatrix}'%mt + post + r'$$'))

import sympy as sy
sy.init_printing(use_latex=1)
o = sy.symbols('Omega')
display(o)
sM = sy.Matrix(((2,0,),(0,1)))
sK = sy.Matrix(((3, -2),(-2,2)))
KooM = sK - o*sM
iKooM = KooM.inv()
sp = sy.Matrix(((0,),(1,)))
a,b=(iKooM*sp)
a.expand().simplify(), b.expand().simplify()

with plt.style.context('classic'):
    plot = sy.plot(a, b, (o, 0, 5), ylim=(-10,10), 
                   show=False, adaptive=False, nb_of_points=601)
    plot[0].line_color = 'black'; plot[0].label = '$x_1$'
    plot[1].line_color = 'red'  ; plot[1].label = '$x_2$'
    plot.xlabel = '$\\beta^2$'; plot.ylabel = r'$x_{i,\mathrm{ss}}/\Delta_\mathrm{st}$'
    plot.legend = True
    plot.show()

from IPython.display import HTML
HTML(open('00_custom.css').read())

"""
Explanation: The following code cell (that is executed before any other code cell by the notebook) loads libraries (or functions from libraries) and determines the style of plots and of the notebook itself.  Besides the cell defines a function to format conveniently matrices and vectors.
End of explanation
"""


