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

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import numpy as np
from numpy import pi, sin, cos, exp, sqrt
import style ; style.clean()


# # Resonant excitation

# We want to study the behaviour of an undercritically damped SDOF system when it is
# subjected to a harmonic force $p(t) = p_o \sin\omega_nt$, i.e., when the excitation frequency equals the free vibration frequency of the system.
# 
# Of course, $\beta=1$, $D(\beta,\zeta)|_{\beta=1}=\displaystyle\frac{1}{2\zeta}$
# and $\theta=\pi/2$, hence $$\xi(t)=\Delta_{st}\,\frac{1}{2\zeta}\cos\omega_nt.$$
# 
# Starting from rest conditions, we have
# 
# $$\frac{x(t)}{\Delta_{st}} = \exp(-\zeta\omega_n t)\left(
# -\frac{\omega_n}{2\omega_D}\sin(\omega_n t)
# -\frac{1}{2\zeta}\cos(\omega_n t)\right) + \frac{1}{2\zeta}\cos(\omega_n t)$$
# 
# and, multiplying both sides by $2\zeta$
# 
# \begin{align*}
# x(t)\frac{2\zeta}{\Delta_{st}} = \bar{x}(t)& = 
#     \exp(-\zeta\omega_n t)\left(
#     -\zeta\frac{\omega_n}{\omega_D}\sin(\omega_n t)
#     -\cos(\omega_n t)\right) + \cos(\omega_n t)\\
# & = \exp(-\zeta\omega_n t)\left(
#     -\frac{\zeta}{\sqrt{1-\zeta^2}}\sin(\omega_n t)
#     -\cos(\omega_n t)\right) + \cos(\omega_n t).
# \end{align*}
# 
# We have now a normalized function of time that grows, oscillating, from 0 to 1,
# where the free parameters are just $\omega_n$ and $\zeta$.
# 
# To go further, we set arbitrarily $\omega_n=2\pi$ (our plots will be nicer...)
# and have just a dependency on $t$ and $\zeta$.
# 
# Eventually, we define a function of $\zeta$ that returns a function of $t$ only,
# here it is...

# In[2]:


def x_normalized(t, z):
    wn = w = 2*pi
    wd = wn*sqrt(1-z*z)
    # Clough Penzien p. 43
    A = z/sqrt(1-z*z)
    return (-cos(wd*t)-A*sin(wd*t))*exp(-z*wn*t) + cos(w*t)


# Above we compute some constants that depend on $\zeta$,
# i.e., the damped frequency and the coefficient in
# front of the sine term, then  we define a function of time
# in terms of these constants and of $\zeta$ itself.
# 
# Because we are going to use this function with a vector argument,
# the last touch is to _vectorize_ the function just before returning it
# to the caller.

# ## Plotting our results

# We start by using a function defined in the `numpy` aka `np` module to
# generate a vector whose entries are 1001 equispaced real numbers, starting from
# zero and up to 20, inclusive of both ends, and assigning the name `t` to this vector.

# In[3]:


t = np.linspace(0,20,1001)
print(t)


# We want to see what happens for different values of $\zeta$, so we create
# a list of values and assign the name `zetas` to this list.

# In[4]:


zetas = (.02, .05, .10, .20)
print(zetas)


# Now, the real plotting: 
# 
# - `z` takes in turn each of the values in `zetas`,
#   - then we generate a function of time for the current `z`
#   - we generate a plot with a line that goes through the point
#     (a(0),b(0)), (a(1),b(1)), (a(2),b(2)), ...
#     where, in our case, a is the vector `t` and b is the vector
#     returned from the vectorized function `bar_x`
#   - we make a slight adjustement to the extreme values of the y-axis
#     of the plot
#   - we give a title to the plot
#   - we FORCE (`plt.show()`) the plot to be produced.

# In[5]:


for z in zetas:
    plt.plot(t, x_normalized(t, z))
    plt.ylim((-1.0, 1.0))
    plt.title(r'$\zeta=%4.2f$'%(z,))
    plt.show()


# ## Wait a minute!

# So, after all this work, we have that the greater the damping, the smaller the
# number of cycles that's needed to reach the maximum value of the response...
# 
# Yes, it's exactly like that, and there is a reason. Think of it.
# 
# .
# 
# .
# 
# .
# 
# .
# 
# .
# 
# .
# 
# .
# 

# We have normalized the response functions to have always a maximum absolute
# value of one, but in effect the max values __are different__, and a heavily damped 
# system needs less cycles to reach steady-state because the maximum value is much,
# much smaller.
# 
# Let's plot the unnormalized (well, there's still the $\Delta_{st}$ normalization)
# responses.
# 
# Note the differences with above:
# 
# - we focus on a shorter interval of time and, in each step
#  - we don't add a title
#  - we don't force the creation of a distinct plot in each cycle,
#  - we add a `label` to each curve
# 
# at the end of the cycle, 
# 
# - we ask for the generation of a `legend` that uses the labels 
#   we specified to generate a, well, a legend for the curves
# - we ask to plot all the properly labeled curves using `plt.plot()`.

# In[7]:


t = np.linspace(0,5,501)
for z in zetas:
    plt.plot(t, x_normalized(t, z)/(2*z), label=r'$\zeta=%4.2f$'%(z,))
plt.legend(ncol=5,loc='lower center', fancybox=1, shadow=1, framealpha=.95)
plt.grid()