<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 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-z^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...
Step1:  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 pylab aka pl 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.
Step2:  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.
Step3:  Now, the real plotting
Step4:  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
<ASSISTANT_TASK:>
Python Code:
def x_2z_over_dst(z):
    w = 2*pi
    # beta = 1, wn =w
    wd = w*sqrt(1-z*z)
    # Clough Penzien p. 43
    A = z/sqrt(1-z*z)
    def f(t):
        return (cos(wd*t)+A*sin(wd*t))*exp(-z*w*t)-cos(w*t)
    return pl.vectorize(f)
Explanation: 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-z^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...
End of explanation
t = pl.linspace(0,20,1001)
print(t)
Explanation: 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 pylab aka pl 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.
End of explanation
zetas = (.02, .05, .10, .20)
print(zetas)
Explanation: 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.
End of explanation
for z in zetas:
    # call the function of zeta that returns
    # a function of time, assign the name bar_x to this function
    bar_x = x_2z_over_dst(z)
    # do the plotting...
    pl.plot(t,bar_x(t))
    pl.ylim((-1.0, 1.0))
    pl.title(r'$\zeta=%4.2f$'%(z,))
    pl.show()
Explanation: 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 (pl.show()) the plot to be produced.
End of explanation
t = pl.linspace(0,5,501)
for z in zetas:
    # call the function of zeta that returns
    # a function of time, assign the name bar_x to this function
    bar_x = x_2z_over_dst(z)
    # do the plotting...
    pl.plot(t,bar_x(t)/2/z, label=r'$\zeta=%4.2f$'%(z,))
pl.legend(ncol=5,loc='lower center', fancybox=1, shadow=1, framealpha=.95)
pl.grid()
Explanation: 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 pl.plot().
End of explanation
<END_TASK>