SRS - the shock response spectrum
=================================

This and other notebooks are available here:
https://github.com/twmacro/pyyeti/tree/master/docs/tutorials.

First, do some imports:

.. code:: ipython3

    import numpy as np
    import matplotlib.pyplot as plt
    from pyyeti import srs

Some settings specifically for the jupyter notebook.

.. code:: ipython3

    %matplotlib inline
    plt.rcParams['figure.figsize'] = [6.4, 4.8]
    plt.rcParams['figure.dpi'] = 150.

Compute an SRS of a half cycle square wave input
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The `SRS <../modules/generated/pyyeti.srs.srs.html#pyyeti.srs.srs>`__
routine has many options, most of which are set to sensible defaults.
For this demo, we’ll mostly use the defaults. Note that parallel
processing is on by default for non-Windows systems; on Windows, the
default is to not use parallel processing because it only seems to slow
things down – I’m not certain why. Consult the documentation to learn
more about the options available for computing SRS’s.

.. code:: ipython3

    sig = np.zeros(1000)
    sig[10:200] = 1.
    sr = 1000.
    freq = np.arange(.1, 100, .1)
    Q = 25
    sh = srs.srs(sig, sr, freq, Q)

Plot time signal and resulting SRS:

.. code:: ipython3

    t = np.arange(len(sig))/sr
    
    plt.subplot(211)
    plt.plot(t, sig)
    plt.ylim(0, 1.2)
    plt.title('Signal')
    plt.xlabel('Time (s)')
    plt.ylabel('Base Acceleration (g)')
    plt.grid(True)
    
    plt.subplot(212)
    plt.plot(freq, sh)
    plt.title('SRS Q={:g}'.format(Q))
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Absolute Acceleration (g)')
    plt.grid(True)
    plt.tight_layout()



.. image:: temp_srs_files/temp_srs_8_0.png


You can also get the response histories. This output is a dict with
members: ‘t’, ‘hist’, ‘sr’. The ‘hist’ value is a 3D array to accomodate
multiple input signals: ``(len(t) x nsignals x len(freq))``. The
following plots a couple different frequencies for illustration:

.. code:: ipython3

    sh, resp = srs.srs(sig, sr, freq, Q, getresp=1)
    resp['hist'].shape




.. code-block:: none

    (2000, 1, 999)



.. code:: ipython3

    for j in np.searchsorted(freq, [20., 60.]):
        plt.plot(resp['t'], resp['hist'][:, 0, j],
                 label='{:g} Hz Oscillator'.format(freq[j]))
    plt.legend(loc='best')
    plt.xlabel('Time (s)')
    plt.ylabel('Acceleration (g)')
    plt.grid(True)



.. image:: temp_srs_files/temp_srs_11_0.png


--------------

Compare the different roll-off methods.
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

The recommendation is to use either the “lanczos” or “fft” methods.

.. code:: ipython3

    sr = 200
    t = np.arange(0, 5, 1/sr)
    sig = np.sin(2*np.pi*15*t) + 3*np.sin(2*np.pi*85*t)
    Q = 50
    frq = np.linspace(5, 100, 476)
    for meth in ['none', 'linear', 'prefilter', 'lanczos', 'fft']:
        sh = srs.srs(sig, sr, frq, Q, rolloff=meth)
        plt.plot(frq, sh, label=meth)
    plt.legend(loc='best')
    ttl = '85 Hz peak should approach 150'
    plt.title(ttl)
    plt.xlabel('Frequency (Hz)')
    plt.grid(True)



.. image:: temp_srs_files/temp_srs_13_0.png


**Warning:** if your points-per-cycle (``ppc``) value is too low, the
rolloff method will be irrelevant. The exception to this rule is that
the “prefilter” method does not depend on ``ppc``. To demonstrate, the
following example sets the ``ppc`` value to 2:

.. code:: ipython3

    for meth in ['none', 'linear', 'prefilter', 'lanczos', 'fft']:
        sh = srs.srs(sig, sr, frq, Q, rolloff=meth, ppc=2)
        plt.plot(frq, sh, label=meth)
    plt.legend(loc='best')
    ttl = '85 Hz peak should approach 150'
    plt.title(ttl)
    plt.xlabel('Frequency (Hz)')
    plt.grid(True)



.. image:: temp_srs_files/temp_srs_15_0.png