Reflection & Transmission Coefficients: Impedance Contrast vs Incidence Angle#
Colab note: This notebook is designed to run on Google Colab. The first code cell installs dependencies.
Learning Objectives
Use Snell’s law (constant \(p\)) to relate incidence angle to slowness in each layer.
Compute SH reflection/transmission coefficients \(S^{\downarrow}S^{\uparrow}\) and \(S^{\downarrow}S^{\downarrow}\) as functions of \(p\) (or \(\theta_1\)).
Predict and identify polarity reversals, critical angles, and postcritical (complex) coefficients.
Distinguish amplitude ratios (“raw” coefficients) from energy partitioning (energy-normalized coefficients).
Prerequisites:
Complex numbers
Basic trigonometry
Plotting with Matplotlib
Reference: Shearer (Chapter 6.3), Reflection and transmission coefficients
Notebook Outline:
This notebook explores how seismic waves reflect and transmit at an interface between two media, focusing on how impedance contrast and incidence angle control the behavior of the coefficients.
# 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
# Install missing packages if in Colab/User environment
# Minimal dependencies for this notebook
required_packages = ["numpy", "matplotlib"]
missing_packages = []
for package in required_packages:
try:
__import__(package)
except ImportError:
missing_packages.append(package)
if missing_packages:
print(f"Installing missing packages: {missing_packages}")
try:
import subprocess
subprocess.check_call([sys.executable, "-m", "pip", "install"] + missing_packages)
print("Installation complete.")
except Exception as e:
print(f"Failed to install packages: {e}")
else:
print("All required packages are installed.")
import numpy as np
import matplotlib.pyplot as plt
1. Utilities: ray parameter \(p\) and vertical slowness \(\eta\)#
Following the chapter:
Horizontal slowness (ray parameter) is constant across the interface: $\(p = u\,\sin\theta = \frac{\sin\theta}{c}\)$
Vertical slowness: $\(\eta = \sqrt{u^2 - p^2} = \sqrt{\frac{1}{c^2} - p^2}\)$
For postcritical incidence, \(\eta\) becomes imaginary; we keep complex arithmetic so this happens naturally.
def ray_parameter(theta1_rad, c1):
"""Horizontal slowness p (s/km) given incidence angle theta1 and velocity c1 (km/s)."""
return np.sin(theta1_rad) / c1
def vertical_slowness(p, c):
"""Vertical slowness eta = sqrt(1/c^2 - p^2). Complex for postcritical."""
return np.sqrt((1.0 / c**2) - p**2 + 0j)
def cos_theta_from_p(p, c):
"""Define cos(theta) := c * eta (works pre- and postcritical)."""
return c * vertical_slowness(p, c)
def theta_from_p(p, c):
"""Incidence angle theta (from vertical). Returns complex if p>1/c."""
# sin(theta)=p*c; for postcritical this is >1 so theta becomes complex; keep real part for plotting.
return np.arcsin(np.clip(np.real(p*c), -1, 1))
2. Vertical incidence: impedance contrast (quick calibration)#
From the chapter (vertical incidence):
For SH (and SV) $\(S^{\downarrow}S^{\uparrow}_{\rm vert} = \frac{\rho_1\beta_1 - \rho_2\beta_2}{\rho_1\beta_1 + \rho_2\beta_2}, \qquad S^{\downarrow}S^{\downarrow}_{\rm vert} = \frac{2\rho_1\beta_1}{\rho_1\beta_1 + \rho_2\beta_2}.\)$
For P (note the sign convention difference in the book): $\(P^{\downarrow}P^{\uparrow}_{\rm vert} = -\frac{\rho_1\alpha_1 - \rho_2\alpha_2}{\rho_1\alpha_1 + \rho_2\alpha_2}, \qquad P^{\downarrow}P^{\downarrow}_{\rm vert} = \frac{2\rho_1\alpha_1}{\rho_1\alpha_1 + \rho_2\alpha_2}.\)$
Here \(Z=\rho c\) is the impedance. The chapter emphasizes that a 10% impedance increase corresponds to about a 5% reflection coefficient magnitude at vertical incidence.
def R_T_vertical_SH(rho1, beta1, rho2, beta2):
Z1, Z2 = rho1*beta1, rho2*beta2
R = (Z1 - Z2) / (Z1 + Z2) # Shearer eq. 6.48
T = (2*Z1) / (Z1 + Z2) # Shearer eq. 6.49 (note: transmitted amplitude in layer 2)
return R, T
def R_T_vertical_P(rho1, alpha1, rho2, alpha2):
Z1, Z2 = rho1*alpha1, rho2*alpha2
R = -(Z1 - Z2) / (Z1 + Z2) # Shearer eq. 6.50
T = (2*Z1) / (Z1 + Z2) # Shearer eq. 6.51
return R, T
# Base model (students can edit)
rho1, beta1 = 2.6, 3.2
rho2, beta2 = 2.9, 3.9
Rsh0, Tsh0 = R_T_vertical_SH(rho1, beta1, rho2, beta2)
print(f"Vertical SH: R = {Rsh0:.3f}, T = {Tsh0:.3f} (incident amplitude = 1)")
2.1 Exercise 1 (prediction-first): sweep impedance contrast#
Horizontal axis: \(\Delta Z / \bar{Z}\) (relative to mean impedance)
Vertical axis: reflection coefficient
Mark where you expect polarity reversals.
def impedance(rho, c):
return rho*c
# Sweep by scaling lower-layer impedance while holding upper fixed.
# (This is pedagogical: in reality rho and velocity covary; we'll vary Z2 directly.)
Z1 = impedance(rho1, beta1)
contrast = np.linspace(-0.6, 0.6, 301) # ±60% relative to mean, wide enough to see behavior
Zbar = Z1 # use Z1 as reference for convenience
Z2 = Z1 * (1 + contrast)
R = (Z1 - Z2)/(Z1 + Z2) # SH vertical (same algebra as eq. 6.48)
T = (2*Z1)/(Z1 + Z2) # eq. 6.49
plt.figure()
plt.plot(contrast*100, np.real(R), label="R (SH, vertical)")
plt.axhline(0, linewidth=1)
plt.axvline(0, linewidth=1)
plt.xlabel("Impedance contrast ΔZ/Z1 (%)")
plt.ylabel("Reflection coefficient")
plt.title("Vertical incidence: reflection vs impedance contrast")
plt.legend()
plt.show()
plt.figure()
plt.plot(contrast*100, np.real(T), label="T (SH, vertical)")
plt.axhline(1, linewidth=1)
plt.axvline(0, linewidth=1)
plt.xlabel("Impedance contrast ΔZ/Z1 (%)")
plt.ylabel("Transmission amplitude (raw)")
plt.title("Vertical incidence: transmission vs impedance contrast")
plt.legend()
plt.show()
3. SH reflection/transmission vs incidence angle#
references: (Shearer 6.3.1)
For an incident downgoing SH wave in layer 1, Shearer derives:
Angles are measured from vertical and are linked by constant \(p\) (Snell’s law).
def SH_coefficients_raw(theta1_deg, rho1, beta1, rho2, beta2):
"""Raw (amplitude) SH reflection/transmission for incident downgoing SH in layer 1."""
theta1 = np.deg2rad(theta1_deg)
p = ray_parameter(theta1, beta1) # p = sin(theta1)/beta1
c1, c2 = beta1, beta2
cos1 = cos_theta_from_p(p, c1)
cos2 = cos_theta_from_p(p, c2)
num = rho1*c1*cos1 - rho2*c2*cos2
den = rho1*c1*cos1 + rho2*c2*cos2
R = num/den
T = (2*rho1*c1*cos1)/den
return p, R, T
def SH_coefficients_energy_norm(theta1_deg, rho1, beta1, rho2, beta2):
"""Energy-normalized SH coefficients (Shearer 6.56-6.57)."""
p, Rraw, Traw = SH_coefficients_raw(theta1_deg, rho1, beta1, rho2, beta2)
c1, c2 = beta1, beta2
cos1 = cos_theta_from_p(p, c1)
cos2 = cos_theta_from_p(p, c2)
# Reflection norm equals raw for same-wave-type reflection (Shearer 6.56)
Rnorm = Rraw
# Transmission norm includes the sqrt factor (Shearer 6.57)
Tnorm = Traw * np.sqrt((rho2*c2*cos2)/(rho1*c1*cos1))
return p, Rnorm, Tnorm, Rraw, Traw
3.1 Exercise 2#
We reproduce the chapter’s key qualitative behaviors:
Near-vertical: coefficients ~ vertical-incidence values.
At some intermediate angle, \(|S^{\downarrow}S^{\uparrow}|\) may pass through zero (no reflected pulse).
Near the critical angle (when \(\theta_2\to 90^\circ\)), the transmitted ray becomes horizontal; coefficients can become large in amplitude.
Postcritical: coefficients become complex, and the reflected pulse is distorted by a phase advance.
# Choose a model. (This example is similar in spirit to Shearer’s Moho example, but editable.)
rho1, beta1 = 2.9, 3.9
rho2, beta2 = 3.38, 4.49
angles = np.linspace(0, 80, 401)
Rraw = np.zeros_like(angles, dtype=complex)
Traw = np.zeros_like(angles, dtype=complex)
Rnorm = np.zeros_like(angles, dtype=complex)
Tnorm = np.zeros_like(angles, dtype=complex)
pvals = np.zeros_like(angles, dtype=float)
for i, th in enumerate(angles):
p, Rn, Tn, Rr, Tr = SH_coefficients_energy_norm(th, rho1, beta1, rho2, beta2)
pvals[i] = np.real(p)
Rraw[i], Traw[i] = Rr, Tr
Rnorm[i], Tnorm[i] = Rn, Tn
# Critical angle for transmission into layer 2: p = 1/beta2 -> sin(theta1c)/beta1 = 1/beta2
# => sin(theta1c) = beta1/beta2 (only if beta1<beta2)
if beta1 < beta2:
theta1c = np.rad2deg(np.arcsin(beta1/beta2))
else:
theta1c = np.nan
plt.figure()
plt.plot(angles, np.abs(Rraw), label="|S↓S↑| (raw)")
plt.plot(angles, np.abs(Traw), label="|S↓S↓| (raw)")
if np.isfinite(theta1c):
plt.axvline(theta1c, linestyle="--", label=f"critical θ1 ≈ {theta1c:.1f}°")
plt.xlabel("Incidence angle θ1 (deg, from vertical)")
plt.ylabel("Magnitude")
plt.title("SH raw coefficients vs incidence angle")
plt.legend()
plt.show()
plt.figure()
plt.plot(angles, np.angle(Rraw, deg=True), label="phase(S↓S↑) [deg]")
plt.plot(angles, np.angle(Traw, deg=True), label="phase(S↓S↓) [deg]")
if np.isfinite(theta1c):
plt.axvline(theta1c, linestyle="--", label="critical θ1")
plt.xlabel("Incidence angle θ1 (deg, from vertical)")
plt.ylabel("Phase (deg)")
plt.title("SH phase vs incidence angle (complex postcritical behavior)")
plt.legend()
plt.show()
3.2 Exercise 3: impedance contrast vs incidence angle#
We now turn this into a toy parameter study with two knobs:
Impedance ratio \(Z_2/Z_1\) where \(Z=\rho\beta\) (SH impedance).
Incidence angle \(\theta_1\) (from vertical) in layer 1.
We will visualize \(S^{\downarrow}S^{\uparrow}\) as a heatmap in \((\theta_1, Z_2/Z_1)\) space.
Prediction prompt: Where do you expect:
\(S^{\downarrow}S^{\uparrow}=0\)? (no reflection)
\(|S^{\downarrow}S^{\uparrow}|=1\)? (total reflection; typically postcritical)
regions of strong phase rotation (complex coefficients)?
# Fix beta1 and beta2 ratio so critical behavior exists; vary rho2 to tune impedance.
rho1, beta1 = 2.9, 3.9
beta2 = 4.49
# We'll vary Z2/Z1 by varying rho2 while holding beta2 fixed (toy knob).
Zratio = np.linspace(0.6, 1.6, 201) # Z2/Z1
angles = np.linspace(0, 80, 321)
Rmag = np.zeros((len(Zratio), len(angles)))
Rph = np.zeros((len(Zratio), len(angles)))
Z1 = rho1*beta1
for i, zr in enumerate(Zratio):
Z2 = zr*Z1
rho2 = Z2/beta2
for j, th in enumerate(angles):
_, R, _ = SH_coefficients_raw(th, rho1, beta1, rho2, beta2)
Rmag[i, j] = np.abs(R)
Rph[i, j] = np.angle(R, deg=True)
plt.figure()
plt.imshow(Rmag, aspect="auto", origin="lower",
extent=[angles.min(), angles.max(), Zratio.min(), Zratio.max()])
plt.colorbar(label="|S↓S↑|")
plt.xlabel("Incidence angle θ1 (deg)")
plt.ylabel("Impedance ratio Z2/Z1")
plt.title("Magnitude of SH reflection coefficient vs angle and impedance contrast")
plt.show()
plt.figure()
plt.imshow(Rph, aspect="auto", origin="lower",
extent=[angles.min(), angles.max(), Zratio.min(), Zratio.max()])
plt.colorbar(label="phase(S↓S↑) [deg]")
plt.xlabel("Incidence angle θ1 (deg)")
plt.ylabel("Impedance ratio Z2/Z1")
plt.title("Phase of SH reflection coefficient vs angle and impedance contrast")
plt.show()
4. Energy partitioning#
The chapter highlights that amplitude is not conserved, but energy is.
Energy flux density on the interface for a wave of amplitude \(A\) is proportional to: $\(\tilde{E}_{\rm flux} \propto \rho\,c\,A^2\,\cos\theta.\)$ So we define an energy-normalized coefficient (Shearer 6.55) so that the sum of squared coefficients equals 1.
4.1 Exercise 4: verify energy conservation (precritical)#
Pick an incidence angle below critical and verify: $\(|R_{\rm norm}|^2 + |T_{\rm norm}|^2 \approx 1.\)$
Then repeat above critical and verify that transmission energy goes to zero and \(|R|\to 1\) (total internal reflection).
def energy_check(theta1_deg, rho1, beta1, rho2, beta2):
p, Rn, Tn, Rr, Tr = SH_coefficients_energy_norm(theta1_deg, rho1, beta1, rho2, beta2)
return {
"theta1_deg": theta1_deg,
"p": float(np.real(p)),
"Rraw": Rr,
"Traw": Tr,
"Rnorm": Rn,
"Tnorm": Tn,
"energy_sum": np.abs(Rn)**2 + np.abs(Tn)**2
}
rho1, beta1 = 2.9, 3.9
rho2, beta2 = 3.38, 4.49
for th in [10, 30, 50, 65]:
out = energy_check(th, rho1, beta1, rho2, beta2)
print(f"θ1={out['theta1_deg']:>4.0f}°, energy sum |R|^2+|T|^2 = {out['energy_sum']:.4f}")
print(f" Rraw={out['Rraw']:.3f}, Traw={out['Traw']:.3f} (complex allowed)")
5. Synthesis questions#
Impedance-only vs angle-dependent: At vertical incidence, coefficients depend only on \(Z=\rho c\). Why does angle dependence introduce \(\cos\theta\) factors? What physical quantity is being projected?
Polarity (sign) conventions: In the chapter, why does \(P^{\downarrow}P^{\uparrow}_{\rm vert}\) have the opposite sign convention to \(S^{\downarrow}S^{\uparrow}_{\rm vert}\)?
Zero reflection angle: For the model you chose, find (numerically) the angle where \(S^{\downarrow}S^{\uparrow}=0\). What combination of \(\rho\beta\cos\theta\) balances makes that happen?
Postcritical distortion: When \(S^{\downarrow}S^{\uparrow}\) is complex, how does multiplying a spectrum by \(e^{i\phi}\) distort a time-domain pulse? Relate to the chapter’s discussion of phase advances and Hilbert transforms.