Demo: Seismometer Response - EPS 164 Seismology Course Materials

Colab note: This notebook is designed to run on Google Colab. The first code cell installs dependencies.

This notebook demonstrates the response of a damped harmonic oscillator as a simple model of a seismometer.

It includes interactive widgets so you can vary:

natural frequency $f_0$
damping constant $h$
what the sensor measures (displacement, velocity, acceleration)
what ground motion you want to recover (displacement, velocity, acceleration)

The plots show the amplitude and phase response as functions of frequency.

Background¶

For a simple inertial seismometer,

\ddot z + 2\epsilon \dot z + \omega_0^2 z = -\ddot u,

(1)

where:

$u(t)$ is ground displacement,
$z(t)$ is displacement of the suspended mass relative to the frame,
$\omega_0 = 2\pi f_0$ is the natural angular frequency,
$h = \epsilon/\omega_0$ is the damping constant.

For harmonic motion, the basic response from ground displacement to sensor displacement is

H(\omega) = \frac{Z(\omega)}{U(\omega)} = \frac{\omega^2}{\omega_0^2 - 2 i \epsilon \omega - \omega^2}.

(2)

Each time derivative applied to the sensor motion multiplies the response by $(-i\omega)$ . Each time derivative applied to the ground motion divides the response by $(-i\omega)$ .

# Install dependencies (for Google Colab or missing packages)
import sys

# Check if running in Colab
try:
    import google.colab

    IN_COLAB = True
    print("Running in Google Colab")
except:
    IN_COLAB = False
    print("Running in local environment")

# Required packages for this notebook
required_packages = {
    'numpy': 'numpy',
    'matplotlib': 'matplotlib',
    'ipywidgets': 'ipywidgets'
}

missing_packages = []
for package, pip_name in required_packages.items():
    try:
        __import__(package)
        print(f"✓ {package} is already installed")
    except ImportError:
        missing_packages.append(pip_name)
        print(f"✗ {package} not found")

# Install missing packages
if missing_packages:
    print(f"\nInstalling missing packages: {', '.join(missing_packages)}")
    import subprocess

    subprocess.check_call([sys.executable, "-m", "pip", "install", "-q"] + missing_packages)
    print("✓ Installation complete!")
else:
    print("\n✓ All required packages are installed!")

# Enable widgets in Colab
if IN_COLAB:
    print("\nEnabling interactive widgets...")
    from google.colab import output

    output.enable_custom_widget_manager()
    print("✓ Widgets enabled!")

Running in local environment
✓ numpy is already installed
✓ matplotlib is already installed
✓ ipywidgets is already installed

✓ All required packages are installed!

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatLogSlider, Dropdown, Checkbox

plt.style.use('default')
plt.rcParams.update({
    'figure.facecolor': 'white',
    'axes.facecolor': 'white',
    'savefig.facecolor': 'white',
    'font.size': 12,
    'lines.linewidth': 2,
})

def response_function(freq_hz, f0_hz=1.0, h=1.0, sensor_quantity='displacement', ground_quantity='displacement'):
    """
    Compute the complex frequency response of a damped harmonic oscillator
    (simple inertial seismometer).

    This function returns the transfer function relating a chosen ground-motion
    quantity (displacement, velocity, or acceleration) to a chosen sensor-measured
    quantity.

    Parameters
    ----------
    freq_hz : array_like
        Frequencies (Hz) at which to evaluate the response.
    f0_hz : float, optional
        Natural (resonant) frequency of the instrument in Hz.
    h : float, optional
        Damping constant (h = epsilon / omega0). Critical damping occurs at h = 1.
    sensor_quantity : {'displacement', 'velocity', 'acceleration'}, optional
        Quantity measured by the sensor (mass motion).
    ground_quantity : {'displacement', 'velocity', 'acceleration'}, optional
        Ground-motion quantity used as input.

    Returns
    -------
    H : complex ndarray
        Complex frequency response (transfer function) such that:
            H(omega) = sensor_quantity / ground_quantity

    Notes
    -----
    The underlying equation of motion is:
        z̈ + 2ε ż + ω0^2 z = -ü

    The base response (sensor displacement / ground displacement) is:
        H0(ω) = ω^2 / (ω0^2 - 2iεω - ω^2)

    Time derivatives are applied in the frequency domain using factors of (-iω):
        velocity  → (-iω)
        acceleration → (-iω)^2

    Thus, this function generalizes the response to arbitrary combinations of
    ground and sensor quantities.
    """
    omega = 2 * np.pi * freq_hz
    omega0 = 2 * np.pi * f0_hz
    eps = h * omega0

    # Basic response: sensor displacement / ground displacement
    H0 = omega**2 / (omega0**2 - 2j * eps * omega - omega**2)

    deriv_order = {'displacement': 0, 'velocity': 1, 'acceleration': 2}
    ns = deriv_order[sensor_quantity]
    ng = deriv_order[ground_quantity]

    # Apply derivatives in frequency domain
    H = ((-1j * omega) ** ns) * H0 / (((-1j * omega) ** ng) + 0j)
    return H


def amplitude_phase(H):
    """
    Compute amplitude and phase of a complex frequency response.

    Parameters
    ----------
    H : complex ndarray
        Complex response function (e.g., output of response_function).

    Returns
    -------
    amp : ndarray
        Amplitude (magnitude) of the response, |H|.
    phase : ndarray
        Phase of the response in degrees.

    Notes
    -----
    The amplitude describes how the instrument scales signals at each frequency,
    while the phase describes the phase shift (lag or lead) introduced by the
    instrument.
    """
    amp = np.abs(H)
    phase = np.angle(H, deg=True)
    return amp, phase

def plot_response(
    f0_hz=1.0,
    h=1.0,
    sensor_quantity='displacement',
    ground_quantity='displacement',
    show_reference=True,
):
    freq = np.logspace(-2, 1.3, 1200)  # about 0.01 to 20 Hz
    H = response_function(
        freq,
        f0_hz=f0_hz,
        h=h,
        sensor_quantity=sensor_quantity,
        ground_quantity=ground_quantity,
    )
    amp, phase = amplitude_phase(H)

    fig = plt.figure(figsize=(8, 5))
    ax = fig.add_subplot(111)
    ax.loglog(freq, amp)
    if show_reference:
        ax.axvline(f0_hz, linestyle='--')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Amplitude response')
    ax.set_title(
        f'Amplitude: sensor {sensor_quantity} / ground {ground_quantity}\n'
        f'f0 = {f0_hz:.3g} Hz, h = {h:.3g}'
    )
    ax.grid(True, which='both', alpha=0.3)
    plt.show()

    fig = plt.figure(figsize=(8, 5))
    ax = fig.add_subplot(111)
    ax.semilogx(freq, phase)
    if show_reference:
        ax.axvline(f0_hz, linestyle='--')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Phase (degrees)')
    ax.set_title(
        f'Phase: sensor {sensor_quantity} / ground {ground_quantity}\n'
        f'f0 = {f0_hz:.3g} Hz, h = {h:.3g}'
    )
    ax.grid(True, which='both', alpha=0.3)
    plt.show()

    print('Interpretation:')
    print(f'- Natural frequency f0 = {f0_hz:.3g} Hz')
    print(f'- Damping h = {h:.3g} (critical damping is h = 1)')
    print(f'- Displayed response = sensor {sensor_quantity} / ground {ground_quantity}')
    if h < 1:
        print('- Because h < 1, expect a more pronounced resonant peak near f0.')
    elif h > 1:
        print('- Because h > 1, the response is more heavily damped and broader.')
    else:
        print('- h = 1 gives critical damping, often close to optimal for seismometers.')

Interactive response viewer¶

Use the widgets below to explore how the response changes.

interact(
    plot_response,
    f0_hz=FloatLogSlider(value=1.0, base=10, min=-1, max=1, step=0.01, description='f0 (Hz)'),
    h=FloatLogSlider(value=1.0, base=10, min=-1, max=0.7, step=0.01, description='h'),
    sensor_quantity=Dropdown(options=['displacement', 'velocity', 'acceleration'], value='displacement', description='Sensor'),
    ground_quantity=Dropdown(options=['displacement', 'velocity', 'acceleration'], value='displacement', description='Ground'),
    show_reference=Checkbox(value=True, description='Show f0 line'),
);

Compare damping values¶

This second widget fixes the quantity pair and lets you compare several damping choices directly.

def compare_damping(
    f0_hz=1.0,
    sensor_quantity='displacement',
    ground_quantity='displacement',
):
    freq = np.logspace(-2, 1.3, 1200)
    h_values = [0.25, 0.5, 1.0, 2.0, 4.0]

    fig = plt.figure(figsize=(8, 5))
    ax = fig.add_subplot(111)
    for h in h_values:
        H = response_function(freq, f0_hz, h, sensor_quantity, ground_quantity)
        amp, _ = amplitude_phase(H)
        ax.loglog(freq, amp, label=f'h = {h}')
    ax.axvline(f0_hz, linestyle='--')
    ax.set_xlabel('Frequency (Hz)')
    ax.set_ylabel('Amplitude response')
    ax.set_title(
        f'Effect of damping: sensor {sensor_quantity} / ground {ground_quantity}\n'
        f'f0 = {f0_hz:.3g} Hz'
    )
    ax.grid(True, which='both', alpha=0.3)
    ax.legend()
    plt.show()

interact(
    compare_damping,
    f0_hz=FloatLogSlider(value=1.0, base=10, min=-1, max=1, step=0.01, description='f0 (Hz)'),
    sensor_quantity=Dropdown(options=['displacement', 'velocity', 'acceleration'], value='displacement', description='Sensor'),
    ground_quantity=Dropdown(options=['displacement', 'velocity', 'acceleration'], value='displacement', description='Ground'),
);

Suggested things to try¶

Set sensor = displacement and ground = displacement. Then vary $h$ from 0.25 to 4.
Fix $h=1$ and compare ground displacement to ground acceleration.
Increase the natural frequency $f_0$ and see how the instrument shifts sensitivity toward higher frequencies.
Identify cases where the response is approximately flat over a broad frequency band.