#!/usr/bin/env python
# coding: utf-8

# # Dynamic Testing

# We want to measure the dynamical characteristics of a SDOF building system,
# i.e., its mass, its damping coefficient and its elastic stiffness.
# 
# To this purpose, we demonstrate that is sufficient to measure the steady-state
# response of the SDOF when subjected to a number of harmonic excitations with
# different frequencies.
# 
# The steady-state response is characterized by its amplitude, $ρ$ and the phase
# delay, $θ$, as in $x_{SS}(t) = ρ \sin(ωt − θ)$.
# 
# A SDOF structural system is excited by a vibrodyne that exerts a harmonic force
# $p(t) = p_o\sin ωt$, with $p_o = 2.224\,{}$kN at different frequencies, and we can
# measure the steady-state response parameters for two different input frequencies,
# as detailed in the following table.
# 
# <style type="text/css">
# .tg  {border-collapse:collapse;border-spacing:0;text-align:center;}
# .tg td{font-size:14px;padding:10px 5px;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;text-align:center;}
# .tg th{font-size:14px;font-weight:normal;padding:10px 5px;border-style:solid;border-width:1px;overflow:hidden;word-break:normal;text-align:center;}
# </style>
# <center>
# <table class="tg">
#   <tr>
#     <th class="tg-031e">$i$</th>
#     <th class="tg-031e">$ω_i$ (rad/s)</th>
#     <th class="tg-031e">$ρ_i$ (μm)</th>
#     <th class="tg-031e">$θ_i$ (deg) </th>
#     <th class="tg-031e">$\cos θ_i$</th>
#     <th class="tg-031e">$\sin θ_i$</th>
#   </tr>
#   <tr>
#     <td class="tg-031e">1</td>
#     <td class="tg-031e">16.0</td>
#     <td class="tg-031e">183.0</td>
#     <td class="tg-031e">15.0</td>
#     <td class="tg-031e">0.9660</td>
#     <td class="tg-031e">0.2590</td>
#   </tr>
#   <tr>
#     <td class="tg-031e">2</td>
#     <td class="tg-031e">25.0</td>
#     <td class="tg-031e">368.0</td>
#     <td class="tg-031e">55.0</td>
#     <td class="tg-031e">0.5740</td>
#     <td class="tg-031e">0.8190</td>
#   </tr>
# </table>
# </center>

# ## Determination of $k$ and $m$

# We start from the equation for steady-state response amplitude,
# 
# $$\rho=\frac{p_o}{k}\frac{1}{\sqrt{(1-\beta^2)^2+(2\zeta\beta)^2}}$$
# 
# where we collect $(1-\beta^2)^2$ in the radicand in the right member,
# 
# $$\rho=\frac{p_o}{k}\frac{1}{1-\beta^2}\frac{1}{\sqrt{1+[2\zeta\beta/(1-\beta^2)]^2}}$$
# 
# but the equation for the phase angle,
# $\tan\vartheta=\frac{2\zeta\beta}{1-\beta^2}$, can be substituted in
# the radicand, so that, using simple trigonometric identities, we find that
# 
# $$\rho=\frac{p_o}{k}\frac{1}{1-\beta^2}\frac{1}{\sqrt{1+\tan^2\vartheta}}=
# \frac{p_o}{k}\frac{\cos\vartheta}{1-\beta^2}.$$
# 
# With $k(1-\beta^2)=k-k\frac{\omega^2}{k/m}=k-\omega^2m$   and using a
# simple rearrangement, we eventually have
# 
# 
# <center>
# $\displaystyle{k-\omega^2m=\frac{p_o}{\rho}\cos\vartheta.}$
# </center>

# In[ ]:


from scipy import matrix, sqrt, pi, cos, sin, set_printoptions

p0 = 2224.0 # converted from kN to Newton
rho1 = 183E-6 ; rho2 = 368E-6 # converted from μm to m
w1 = 16.0 ; w2 = 25.0 
th1 = 15.0 ; th2 = 55.0
d2r = pi/180.
cos1 = cos(d2r*th1) ; cos2 = cos(d2r*th2)
sin1 = sin(d2r*th1) ; sin2 = sin(d2r*th2)


# In[6]:


# the unknowns are k and m
# coefficient matrix, row i is 1, omega_i^2
coeff = matrix(((1, -w1**2),(1, -w2**2)))
# kt i.e., know term, cos(theta_i)/rho_i * p_0
kt = matrix((cos1/rho1,cos2/rho2)).T*p0
print(coeff)
print(kt)


# In[7]:


k_and_m = coeff.I*kt
k, m = k_and_m[0,0], k_and_m[1,0]
wn2, wn = k/m, sqrt(k/m)
print('       k          m         wn2           wn')
print(k, m, wn2, wn)


# ## Determination of $\zeta$

# Using the previously established relationship for $\cos\vartheta$, we
# can write $\cos\vartheta=k(1-\beta^2)\frac{\rho}{p_o}$, from the
# equation of the phase angle (see above), we can write $\cos\vartheta =
# \frac{1-\beta^2}{2\zeta\beta}\sin\vartheta$, and finally
# 
# $$\frac{\rho k}{p_o}=\frac{\sin\vartheta}{2\zeta\beta},
# \quad\text{hence}\quad
# \zeta=\frac{p_o}{\rho k}\frac{\sin\vartheta}{2\beta}$$
# 
# Lets write some code that gives us our two wstimates

# In[8]:


z1 = p0*sin1/rho1/k/2/(w1/wn)
z2 = p0*sin2/rho2/k/2/(w2/wn)
print(z1*100, z2*100)


# ## Experimental approximation

# We have seen that our two estimates for $\zeta$ are sligtly different, this is due to
# the fact that every measurement is slightly approximated...
# 
# One can partly obviate using a number of measurements larger, or even much larger, than
# the number of parameters s/he's trying to determine.
# 
# In our case, for each set of $M$ observations (an observation for us is a set of three values, $\omega_i, \rho_i, \theta_i$) we can write a linear equation in the $N=2$ unknowns $k$ and $m$ so that, with $b_i = p_o \cos\theta_i/\rho_i$, we can write
# 
# $$Ax-b=0$$

# $Ax-b=0$ is a set of $M$ equations in $N<M$ variables that cannot have an exact solution but we can introduce a vector of errors (or _residuals_ ) with respect to an approximate solution  $y \approx x$,
# 
# $$e = A\,y - b$$
# 
# and assume (reasonably) that the _best solution_ is the one that minimizes the scalar norm of the residual
# 
# $$\left|e\right|=(y^TA^T- b^T)(A\,y - b) = y^TA^TA\,y - 2 y^TA^Tb +b^Tb$$
# 
# Minimizing the norm is equivalent to take the derivatives of the above expression with respect to each
# of the $y_i$ unknowns and equate every derivative to zero,
# 
# $$\frac{\partial |e|}{\partial y_i} = 0,\quad\text{for }i=1,\ldots,N$$
# 
# and the set of equations above can be written using matrix notation: 
# 
# $$A^TA\,y - A^Tb = 0.$$
# 
# Solving the above equatio gives $y  =(A^TA)^{-1}A^Tb$, and $y$ is the vector of parameters that minimizes the norm
# of the residual.
# 
# <span style="color:blue">The usual name of the procedure that I have sketched is, of course, _"Least Squares Parameter Estimation"_.</span>
# 
# We'll deal again with the derivative of a quadratic form in the future...