Body Waves: P and S Waves

Body Waves: P and S Waves#

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

Learning Objectives:

  • Derive the wave equation for P and S waves

  • Understand particle motion polarization

  • Visualize wave propagation in elastic media

Prerequisites: Vector calculus (divergence, curl), Partial differential equations

Reference: Shearer, Chapter 3 (Seismic Waves)

Notebook Outline:

This notebook explores the theoretical foundations of body waves (P and S waves) through computational exercises:

  1. Toy problem: Plane-wave polarization to verify P vs S motion

Hide code cell source

 
# 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")

# Install required packages if needed
required_packages = {
    'numpy': 'numpy',
    'matplotlib': 'matplotlib',
    'scipy': 'scipy',
    'obspy': 'obspy'
}

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")

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!")

Hide code cell source

 
import numpy as np
import matplotlib.pyplot as plt

from numpy.linalg import norm, eig

Some functions

def central_diff(x, dt):
    """
    Central difference derivative for equally spaced samples.
    Returns derivative with same length, using one-sided at ends.
    """
    dx = np.zeros_like(x)
    dx[1:-1] = (x[2:] - x[:-2])/(2*dt)
    dx[0] = (x[1]-x[0])/dt
    dx[-1] = (x[-1]-x[-2])/dt
    return dx

def moving_window(x, i0, i1):
    """Safe slicing."""
    i0 = max(i0, 0)
    i1 = min(i1, len(x))
    return x[i0:i1]

def polarization_metrics(Z, N, E):
    """
    Polarization via covariance eigendecomposition.
    Returns eigenvalues (sorted desc), eigenvectors (columns), rectilinearity.
    """
    X = np.vstack([Z, N, E]).T
    X = X - X.mean(axis=0, keepdims=True)
    C = (X.T @ X) / X.shape[0]
    w, V = eig(C)
    idx = np.argsort(w)[::-1]
    w = np.real(w[idx])
    V = np.real(V[:, idx])
    # Common rectilinearity proxy: 1 - (λ2+λ3)/(2*λ1)
    rect = 1.0 - (w[1] + w[2])/(2*w[0] + 1e-30)
    return w, V, rect

def angle_from_vertical(v):
    """
    Incidence angle from vertical (Z axis).
    v is a 3-vector (Z,N,E) direction.
    """
    v = v / (norm(v) + 1e-30)
    # Z axis = [1,0,0]
    return np.degrees(np.arccos(np.clip(v[0], -1, 1)))

def baz_from_NE(v):
    """
    Back-azimuth estimate from horizontal projection of v in (N,E).
    Returns degrees in [0,360).
    """
    v = v / (norm(v) + 1e-30)
    n, e = v[1], v[2]
    baz = (np.degrees(np.arctan2(e, n)) + 360) % 360
    return baz

1. Toy problem: plane-wave polarization (P vs S)#

1.1 Build a monochromatic plane wave (Chapter 3.4–3.5)#

We generate a plane wave in 1-D (propagating in +x), but we treat displacement as a 3-component vector.

# Sampling
dt = 0.01
t = np.arange(0, 6, dt)

# Wave parameters
f = 1.0
w = 2*np.pi*f
cP = 6.0  # km/s (toy)
cS = 3.5  # km/s (toy)

# Choose an "x" location for a synthetic seismogram at a receiver
x = 30.0  # km
phaseP = w*(t - x/cP)
phaseS = w*(t - x/cS)

A = 1.0

# P wave polarized along x (map x -> N for convenience later)
# We'll use components (Z, N, E). For toy: let propagation be "N".
Zp = np.zeros_like(t)
Np = A*np.cos(phaseP)
Ep = np.zeros_like(t)

# S wave polarized transverse (along E)
Zs = np.zeros_like(t)
Ns = np.zeros_like(t)
Es = A*np.cos(phaseS)

# Plot
plt.figure()
plt.plot(t, Np, label="P: N component")
plt.plot(t, Es, label="S: E component")
plt.xlim(0, 6)
plt.xlabel("Time (s)")
plt.ylabel("Displacement (arb.)")
plt.legend()
plt.title("Toy plane waves: longitudinal P vs transverse S")
plt.show()

1.2 Test 1 (concept): “P is longitudinal, S is transverse”#

Automated test: For a “pure” P, almost all energy should be in N; for a “pure” S, almost all energy should be in E.

def energy(x): return np.sum(x**2)

P_energy = np.array([energy(Zp), energy(Np), energy(Ep)])
S_energy = np.array([energy(Zs), energy(Ns), energy(Es)])

print("P energy fractions (Z,N,E):", P_energy/P_energy.sum())
print("S energy fractions (Z,N,E):", S_energy/S_energy.sum())

assert (P_energy[1] / P_energy.sum()) > 0.99, "P not sufficiently longitudinal"
assert (S_energy[2] / S_energy.sum()) > 0.99, "S not sufficiently transverse"
print("PASS: Toy polarization test.")

2. Toy problem: div/curl intuition (P ↔ div, S ↔ curl) — and its limits#

Chapter 3 shows that in a homogeneous medium the equation of motion can be written (Eq. 3.19–3.27) and decoupled by taking divergence and curl:

P: wave equation for ∇·u

S: wave equation for ∇×u

At a single receiver, you do not measure spatial derivatives directly. So we do two things:

  • Build a tiny synthetic spatial grid to compute numerical div/curl.

  • Explicitly note what cannot be done with a single seismogram.

2.1 Synthetic 3-D displacement field on a small grid#

We build a plane P wave with u = ∇φ and an S wave with u = ∇×ψ (Helmholtz decomposition, Eq. 3.28). Then compute numerical divergence and curl and verify:

P: curl ≈ 0, div ≠ 0

S: div ≈ 0, curl ≠ 0

# Small spatial grid (km)
dx = 0.5
x = np.arange(-10, 10+dx, dx)
y = np.arange(-10, 10+dx, dx)
z = np.arange(-10, 10+dx, dx)

X, Y, Z = np.meshgrid(x, y, z, indexing='ij')

# Snapshot time
t0 = 2.0

# P-wave scalar potential φ = cos(ω(t - x/cP))
phi = np.cos(w*(t0 - X/cP))

# uP = ∇φ (approx via finite differences)
uPx = np.gradient(phi, dx, axis=0)
uPy = np.gradient(phi, dx, axis=1)
uPz = np.gradient(phi, dx, axis=2)

# S-wave vector potential ψ with only z-component -> produces shear in x-y plane
# Let ψz = cos(ω(t - x/cS)); ψx=ψy=0
psiz = np.cos(w*(t0 - X/cS))
psix = np.zeros_like(psiz)
psiy = np.zeros_like(psiz)

# uS = ∇×ψ
# curl(ψ) = [∂y ψz - ∂z ψy, ∂z ψx - ∂x ψz, ∂x ψy - ∂y ψx]
dpsiz_dy = np.gradient(psiz, dx, axis=1)
dpsiz_dx = np.gradient(psiz, dx, axis=0)
uSx = dpsiz_dy - np.gradient(psiy, dx, axis=2)
uSy = np.gradient(psix, dx, axis=2) - dpsiz_dx
uSz = np.gradient(psiy, dx, axis=0) - np.gradient(psix, dx, axis=1)

# divergence
div_uP = (np.gradient(uPx, dx, axis=0) +
          np.gradient(uPy, dx, axis=1) +
          np.gradient(uPz, dx, axis=2))

div_uS = (np.gradient(uSx, dx, axis=0) +
          np.gradient(uSy, dx, axis=1) +
          np.gradient(uSz, dx, axis=2))

# curl magnitude
curl_uP_x = np.gradient(uPz, dx, axis=1) - np.gradient(uPy, dx, axis=2)
curl_uP_y = np.gradient(uPx, dx, axis=2) - np.gradient(uPz, dx, axis=0)
curl_uP_z = np.gradient(uPy, dx, axis=0) - np.gradient(uPx, dx, axis=1)
curlmag_uP = np.sqrt(curl_uP_x**2 + curl_uP_y**2 + curl_uP_z**2)

curl_uS_x = np.gradient(uSz, dx, axis=1) - np.gradient(uSy, dx, axis=2)
curl_uS_y = np.gradient(uSx, dx, axis=2) - np.gradient(uSz, dx, axis=0)
curl_uS_z = np.gradient(uSy, dx, axis=0) - np.gradient(uSx, dx, axis=1)
curlmag_uS = np.sqrt(curl_uS_x**2 + curl_uS_y**2 + curl_uS_z**2)

print("P field: mean |div| =", np.mean(np.abs(div_uP)), "mean |curl| =", np.mean(curlmag_uP))
print("S field: mean |div| =", np.mean(np.abs(div_uS)), "mean |curl| =", np.mean(curlmag_uS))

2.2 Test 2 (concept): P is (approximately) irrotational, S is (approximately) divergence-free#

# We compare ratios to avoid scale dependence
ratio_P = np.mean(curlmag_uP) / (np.mean(np.abs(div_uP)) + 1e-30)
ratio_S = np.mean(np.abs(div_uS)) / (np.mean(curlmag_uS) + 1e-30)

print("ratio_P = <|curl|>/<|div|> (should be small):", ratio_P)
print("ratio_S = <|div|>/<|curl|> (should be small):", ratio_S)

assert ratio_P < 0.05, "P not sufficiently irrotational in this discretization"
assert ratio_S < 0.05, "S not sufficiently divergence-free in this discretization"
print("PASS: div/curl separation on a synthetic grid.")

2.3 Reflection: what you can/cannot do with a single station#

You can test polarization (particle motion direction) very well from 3C recordings (Ch. 3.5).

You cannot compute true spatial divergence/curl of displacement from one station; that requires an array or a wavefield model.

3. Real data: download 3-component seismograms with ObsPy (FDSN)#

We now pull real waveform data and repeat the polarization analysis on P and S windows.

3.1 Fetch a catalog event + station#

You can choose any event/station pair. The default below is designed to be editable.

from obspy.clients.fdsn import Client
from obspy import UTCDateTime
from obspy.geodetics import gps2dist_azimuth
from obspy.taup import TauPyModel

client = Client("IRIS")  # or "USGS", "RESIF", etc.

# --- Pick an event time and station (edit these freely) ---
origin_time = UTCDateTime("2021-11-28T10:52:14")  # example; change as needed
oo = client.get_events(starttime=origin_time - 10, endtime=origin_time + 10,minmagnitude=7)[0]
print(oo)

net, sta, loc, cha = "C1", "VA04", "*", "BH?"  # 3C broadband
sta = client.get_stations(network=net, station=sta, location=loc,
                          channel=cha, starttime=origin_time,
                          endtime=origin_time + 3600,
                          level="response")

ev_lat = oo.origins[0].latitude
ev_lon = oo.origins[0].longitude
ev_depth_km = oo.origins[0].depth / 1000.0  # m -> km
st_lat = sta[0].stations[0].latitude
st_lon = sta[0].stations[0].longitude
# Distance/azimuth
dist_m, az, baz = gps2dist_azimuth(ev_lat, ev_lon, st_lat, st_lon)
dist_deg = dist_m / (1000.0 * 111.19)  # rough deg; TauP wants degrees. We'll compute precisely below too.
dist_deg = dist_m / 1000.0 / 111.19

print("Approx distance (deg):", dist_deg, "az:", az, "baz:", baz)

3.2 Predict P and S arrivals (TauP)#

model = TauPyModel(model="iasp91")
arrivals = model.get_travel_times(source_depth_in_km=ev_depth_km,
                                  distance_in_degree=dist_deg)

for a in arrivals[:18]:
    print(a.name, "t(s)=", a.time)

tP = [a.time for a in arrivals if a.name == "P"][0]
tS = [a.time for a in arrivals if a.name == "S"][0]
print("Predicted P:", tP, "s ; S:", tS, "s")

3.3 Download waveforms around the arrivals#

t0 = origin_time + (tP - 300)  # start 5 min before P

st = client.get_waveforms(network=net, station=sta, location="*", channel=cha,
                          starttime=t0, endtime=t0 + 3600, attach_response=True)

st.merge(fill_value="interpolate")
st.detrend("linear")
st.detrend("demean")
st.filter("bandpass", freqmin=0.02, freqmax=1.0, corners=4, zerophase=True)
st.taper(0.02)

print(st)
st.plot()

3.4 Rotate to ZNE (if needed), then to RT (radial/transverse)#

For polarization and P/S separation, RT rotation is often helpful:

P energy tends to be strong on Z and R

SH tends to be strong on T

SV tends to mix Z and R

# Ensure we have Z, N, E traces
# Many networks already provide BHZ,BHN,BHE.
# If you have BH1/BH2, you'd need instrument azimuth metadata and rotate.
stZ = st.select(channel="*Z")[0]
stN = st.select(channel="*N")[0]
stE = st.select(channel="*E")[0]

# Rotate NE -> RT using back-azimuth
from obspy.signal.rotate import rotate_ne_rt

N = stN.data.astype(float)
E = stE.data.astype(float)
R, T = rotate_ne_rt(N, E, baz)

# Time axis
npts = stZ.stats.npts
dt = stZ.stats.delta
times = np.arange(npts)*dt
abs_times = stZ.stats.starttime + times

# Plot Z, R, T
plt.figure(figsize=(10,6))
plt.plot(times, stZ.data/np.max(np.abs(stZ.data)), label="Z (norm)")
plt.plot(times, R/np.max(np.abs(R)), label="R (norm)")
plt.plot(times, T/np.max(np.abs(T)), label="T (norm)")
plt.axvline(tP - 300, linestyle="--")
plt.axvline(tS - 300, linestyle="--")
plt.xlabel("Time since window start (s)")
plt.legend()
plt.title("3C components (normalized), with predicted P and S markers")
plt.show()

4. Real data polarization analysis (P window vs S window)#

We quantify polarization in two windows:

P window: around predicted P arrival

S window: around predicted S arrival

4.1 Define windows#

# Window definitions relative to our downloaded start (t0)
P_center = tP - 300
S_center = tS - 300

# Editable window lengths (seconds)
P_pre, P_post = 5, 25
S_pre, S_post = 5, 35

iP0 = int((P_center - P_pre)/dt)
iP1 = int((P_center + P_post)/dt)
iS0 = int((S_center - S_pre)/dt)
iS1 = int((S_center + S_post)/dt)

ZP = stZ.data[iP0:iP1]
NP = N[iP0:iP1]
EP = E[iP0:iP1]

ZS = stZ.data[iS0:iS1]
NS = N[iS0:iS1]
ES = E[iS0:iS1]

print("P window length (s):", (iP1-iP0)*dt, "S window length (s):", (iS1-iS0)*dt)

4.2 Compute polarization metrics (covariance eigenanalysis)#

wP, VP, rectP = polarization_metrics(ZP, NP, EP)
wS, VS, rectS = polarization_metrics(ZS, NS, ES)

vP = VP[:,0]  # principal polarization direction (Z,N,E)
vS = VS[:,0]

print("P eigenvalues:", wP, "rectilinearity:", rectP)
print("S eigenvalues:", wS, "rectilinearity:", rectS)

print("P principal direction (Z,N,E):", vP, "incidence angle from vertical:", angle_from_vertical(vP))
print("S principal direction (Z,N,E):", vS, "incidence angle from vertical:", angle_from_vertical(vS))

Interpretation:

P should be relatively rectilinear (high rectilinearity), with principal axis often near the ray direction (mix of Z and radial).

S can be more complex (SV+SH mixing) but often less rectilinear in noisy/converted contexts.

4.3 Test 3 (concept): P is “more rectilinear than S” (often, but not guaranteed)#

This is a probabilistic test (data quality matters), so we set a mild expectation and make the student diagnose failures.

print("rectP:", rectP, "rectS:", rectS)

# Mild expectation; adjust if needed
if rectP <= rectS:
    print("NOTE: P not more rectilinear than S for this record. Diagnose:")
    print("- Is the P window too long/short?")
    print("- Is bandpass appropriate?")
    print("- Are you near a nodal direction or strongly scattered site?")
else:
    print("PASS (typical): P rectilinearity exceeds S rectilinearity.")

5. P/S separation tests using component behavior (receiver-side)#

Chapter 3’s divergence/curl separation is a field statement in homogeneous media (Eq. 3.20–3.26), but with single-station data we use observable proxies:

5.1 Proxy A: Z/R vs T energy ratios (a practical separation heuristic)#

Compute energy fractions in the P and S windows:

P tends to have higher Z+R energy.

SH tends to show strongly on T around S.

# Build R,T windows too
RP = R[iP0:iP1]; TP = T[iP0:iP1]
RS = R[iS0:iS1]; TS = T[iS0:iS1]

def frac_T(Z, R, T):
    Ez, Er, Et = np.sum(Z**2), np.sum(R**2), np.sum(T**2)
    return Et/(Ez+Er+Et+1e-30), (Ez+Er)/(Ez+Er+Et+1e-30)

fTP, fZR_P = frac_T(ZP, RP, TP)
fTS, fZR_S = frac_T(ZS, RS, TS)

print("P window: T fraction =", fTP, "ZR fraction =", fZR_P)
print("S window: T fraction =", fTS, "ZR fraction =", fZR_S)

5.1.1 Visualizing Polarization Motion#

Particle motion plots show how the ground moves in different planes during the wave passage. This helps visualize:

  • P-wave: Should show rectilinear motion primarily in the Z-R plane (along the ray direction)

  • S-wave: Can show elliptical or rectilinear motion depending on SV/SH content

    • SV motion: appears in Z-R plane

    • SH motion: appears in the transverse (T) component

We’ll plot the particle motion in three planes: Z-R, Z-T, and R-T for both P and S windows.

# Create particle motion plots for P and S windows
fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# Time vectors for P and S windows
tP_window = times[iP0:iP1]
tS_window = times[iS0:iS1]

# Color by time to show motion direction
norm_P = plt.Normalize(vmin=tP_window.min(), vmax=tP_window.max())
norm_S = plt.Normalize(vmin=tS_window.min(), vmax=tS_window.max())
cmap = plt.cm.viridis

# ===== P-wave particle motion =====
# Z-R plane (primary plane for P-waves)
sc1 = axes[0, 0].scatter(RP, ZP, c=tP_window, cmap=cmap, s=20, alpha=0.6)
axes[0, 0].plot(RP, ZP, 'k-', alpha=0.3, linewidth=0.5)
axes[0, 0].plot(RP[0], ZP[0], 'go', markersize=8, label='Start')
axes[0, 0].plot(RP[-1], ZP[-1], 'ro', markersize=8, label='End')
axes[0, 0].set_xlabel('Radial (R)', fontsize=10)
axes[0, 0].set_ylabel('Vertical (Z)', fontsize=10)
axes[0, 0].set_title(f'P-wave: Z-R plane\n(Rectilinearity={rectP:.3f})', fontsize=11)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].legend(fontsize=8)
axes[0, 0].axis('equal')

# Z-T plane (should show minimal motion for pure P)
sc2 = axes[0, 1].scatter(TP, ZP, c=tP_window, cmap=cmap, s=20, alpha=0.6)
axes[0, 1].plot(TP, ZP, 'k-', alpha=0.3, linewidth=0.5)
axes[0, 1].plot(TP[0], ZP[0], 'go', markersize=8)
axes[0, 1].plot(TP[-1], ZP[-1], 'ro', markersize=8)
axes[0, 1].set_xlabel('Transverse (T)', fontsize=10)
axes[0, 1].set_ylabel('Vertical (Z)', fontsize=10)
axes[0, 1].set_title(f'P-wave: Z-T plane\n(T fraction={fTP:.3f})', fontsize=11)
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].axis('equal')

# R-T plane (horizontal plane)
sc3 = axes[0, 2].scatter(TP, RP, c=tP_window, cmap=cmap, s=20, alpha=0.6)
axes[0, 2].plot(TP, RP, 'k-', alpha=0.3, linewidth=0.5)
axes[0, 2].plot(TP[0], RP[0], 'go', markersize=8)
axes[0, 2].plot(TP[-1], RP[-1], 'ro', markersize=8)
axes[0, 2].set_xlabel('Transverse (T)', fontsize=10)
axes[0, 2].set_ylabel('Radial (R)', fontsize=10)
axes[0, 2].set_title('P-wave: R-T plane', fontsize=11)
axes[0, 2].grid(True, alpha=0.3)
axes[0, 2].axis('equal')

# ===== S-wave particle motion =====
# Z-R plane (SV motion appears here)
sc4 = axes[1, 0].scatter(RS, ZS, c=tS_window, cmap=cmap, s=20, alpha=0.6)
axes[1, 0].plot(RS, ZS, 'k-', alpha=0.3, linewidth=0.5)
axes[1, 0].plot(RS[0], ZS[0], 'go', markersize=8, label='Start')
axes[1, 0].plot(RS[-1], ZS[-1], 'ro', markersize=8, label='End')
axes[1, 0].set_xlabel('Radial (R)', fontsize=10)
axes[1, 0].set_ylabel('Vertical (Z)', fontsize=10)
axes[1, 0].set_title(f'S-wave: Z-R plane (SV)\n(Rectilinearity={rectS:.3f})', fontsize=11)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].legend(fontsize=8)
axes[1, 0].axis('equal')

# Z-T plane (SH motion shows on T)
sc5 = axes[1, 1].scatter(TS, ZS, c=tS_window, cmap=cmap, s=20, alpha=0.6)
axes[1, 1].plot(TS, ZS, 'k-', alpha=0.3, linewidth=0.5)
axes[1, 1].plot(TS[0], ZS[0], 'go', markersize=8)
axes[1, 1].plot(TS[-1], ZS[-1], 'ro', markersize=8)
axes[1, 1].set_xlabel('Transverse (T)', fontsize=10)
axes[1, 1].set_ylabel('Vertical (Z)', fontsize=10)
axes[1, 1].set_title(f'S-wave: Z-T plane (SH)\n(T fraction={fTS:.3f})', fontsize=11)
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].axis('equal')

# R-T plane (horizontal plane, shows both SV and SH)
sc6 = axes[1, 2].scatter(TS, RS, c=tS_window, cmap=cmap, s=20, alpha=0.6)
axes[1, 2].plot(TS, RS, 'k-', alpha=0.3, linewidth=0.5)
axes[1, 2].plot(TS[0], RS[0], 'go', markersize=8)
axes[1, 2].plot(TS[-1], RS[-1], 'ro', markersize=8)
axes[1, 2].set_xlabel('Transverse (T)', fontsize=10)
axes[1, 2].set_ylabel('Radial (R)', fontsize=10)
axes[1, 2].set_title('S-wave: R-T plane', fontsize=11)
axes[1, 2].grid(True, alpha=0.3)
axes[1, 2].axis('equal')

# Add colorbars
plt.colorbar(sc3, ax=axes[0, :], label='Time (s)', orientation='horizontal', pad=0.1, shrink=0.8)
plt.colorbar(sc6, ax=axes[1, :], label='Time (s)', orientation='horizontal', pad=0.1, shrink=0.8)

plt.tight_layout()
plt.show()

Interpretation Guide:

  • Rectilinear motion: Points align along a straight line (high rectilinearity)

  • Elliptical motion: Points trace an ellipse (common for S-waves with mixed SV/SH)

  • Scattered motion: Random pattern (low rectilinearity, noisy data)

For P-waves:

  • Z-R plane should show strong rectilinear motion along the ray path

  • Z-T plane should show minimal motion (little transverse energy)

  • Color progression from green (start) to red (end) shows the time evolution

For S-waves:

  • Z-R plane shows SV component (vertical-radial motion)

  • Z-T plane shows SH component (transverse-vertical motion)

  • R-T plane shows horizontal particle motion

  • S-waves often show more complex, elliptical patterns due to SV+SH mixing

5.2 Test 4 (concept): S window tends to be more transverse than P window#

if fTS <= fTP:
    print("NOTE: Transverse fraction not higher in S window. Diagnose:")
    print("- Are we looking at SV-dominant S (less T)?")
    print("- Is the back-azimuth correct? Station metadata correct?")
    print("- Is S window contaminated by later phases/coda?")
else:
    print("PASS (typical): S window is more transverse than P window.")

6. Linking back to Chapter 3 “core derivation” with a mini-derivation check#

This section is conceptual + computational and directly targets the multiple-choice learning outcomes:

Stress divergence as force density (Eq. 3.5–3.9)

Hooke’s law closure (Eq. 3.11–3.13)

Vector identity enabling P/S decoupling (Eq. 3.16–3.17)

Homogeneous vs source region (Eq. 3.9 vs 3.8)

6.1 Quick symbolic/numeric verification of the key vector identity#

Identity used in Chapter 3 (Eq. 3.16):

∇×∇×u=∇(∇⋅u)−∇2u.

We verify numerically on a random smooth vector field on a grid.

# Random smooth field on grid (reuse X,Y,Z from earlier, but smaller for speed)
dx = 1.0
x = np.arange(-8, 8+dx, dx)
y = np.arange(-8, 8+dx, dx)
z = np.arange(-8, 8+dx, dx)
X, Y, Z = np.meshgrid(x, y, z, indexing='ij')

u1 = np.sin(0.2*X)*np.cos(0.1*Y)
u2 = np.cos(0.15*Y)*np.sin(0.1*Z)
u3 = np.sin(0.1*X)*np.cos(0.12*Z)

# Compute div u
divu = np.gradient(u1, dx, axis=0) + np.gradient(u2, dx, axis=1) + np.gradient(u3, dx, axis=2)

# Compute grad(div u)
gdiv_x = np.gradient(divu, dx, axis=0)
gdiv_y = np.gradient(divu, dx, axis=1)
gdiv_z = np.gradient(divu, dx, axis=2)

# Compute Laplacian of u components
lap_u1 = (np.gradient(np.gradient(u1, dx, axis=0), dx, axis=0) +
          np.gradient(np.gradient(u1, dx, axis=1), dx, axis=1) +
          np.gradient(np.gradient(u1, dx, axis=2), dx, axis=2))
lap_u2 = (np.gradient(np.gradient(u2, dx, axis=0), dx, axis=0) +
          np.gradient(np.gradient(u2, dx, axis=1), dx, axis=1) +
          np.gradient(np.gradient(u2, dx, axis=2), dx, axis=2))
lap_u3 = (np.gradient(np.gradient(u3, dx, axis=0), dx, axis=0) +
          np.gradient(np.gradient(u3, dx, axis=1), dx, axis=1) +
          np.gradient(np.gradient(u3, dx, axis=2), dx, axis=2))

# Compute curl(curl(u))
# curl u:
curlu_x = np.gradient(u3, dx, axis=1) - np.gradient(u2, dx, axis=2)
curlu_y = np.gradient(u1, dx, axis=2) - np.gradient(u3, dx, axis=0)
curlu_z = np.gradient(u2, dx, axis=0) - np.gradient(u1, dx, axis=1)

# curl(curl u):
cc_x = np.gradient(curlu_z, dx, axis=1) - np.gradient(curlu_y, dx, axis=2)
cc_y = np.gradient(curlu_x, dx, axis=2) - np.gradient(curlu_z, dx, axis=0)
cc_z = np.gradient(curlu_y, dx, axis=0) - np.gradient(curlu_x, dx, axis=1)

# RHS = grad(div u) - Laplacian(u)
rhs_x = gdiv_x - lap_u1
rhs_y = gdiv_y - lap_u2
rhs_z = gdiv_z - lap_u3

# Compare norms
err = (np.mean(np.abs(cc_x - rhs_x)) +
       np.mean(np.abs(cc_y - rhs_y)) +
       np.mean(np.abs(cc_z - rhs_z)))

scale = (np.mean(np.abs(rhs_x)) + np.mean(np.abs(rhs_y)) + np.mean(np.abs(rhs_z)) + 1e-30)
rel_err = err/scale

print("Relative error (finite-diff check):", rel_err)
assert rel_err < 0.1, "Vector identity check failed beyond tolerance (try finer dx or smoother field)."
print("PASS: Numerical check of key vector identity (within FD tolerance).")
Contents