Earthquake Focal Mechanisms — Exercises¶

EPS 130 — Introduction to Seismology

---

Setup¶

Run the cell below to load the plot_beachball function from the lecture.

In [ ]:

!pip install obspy cartopy -q

!pip install obspy cartopy -q

In [ ]:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D

def plot_beachball(strike, dip, rake, ax=None, title=None, show_compass=True):
    phi = np.radians(strike); delta = np.radians(dip); lam = np.radians(rake)
    n = np.array([-np.sin(delta)*np.sin(phi), np.sin(delta)*np.cos(phi), -np.cos(delta)])
    s = np.array([np.cos(lam)*np.cos(phi)+np.sin(lam)*np.cos(delta)*np.sin(phi),
                  np.cos(lam)*np.sin(phi)-np.sin(lam)*np.cos(delta)*np.cos(phi),
                 -np.sin(lam)*np.sin(delta)])
    M = np.outer(s, n) + np.outer(n, s)
    ng = 300; to = np.linspace(0, np.pi/2, ng); az = np.linspace(0, 2*np.pi, ng)
    TO, AZ = np.meshgrid(to, az)
    gn = np.sin(TO)*np.cos(AZ); ge = np.sin(TO)*np.sin(AZ); gd = np.cos(TO)
    amp = (M[0,0]*gn**2+M[1,1]*ge**2+M[2,2]*gd**2
           +2*M[0,1]*gn*ge+2*M[0,2]*gn*gd+2*M[1,2]*ge*gd)
    R = np.sqrt(2)*np.sin(TO/2); X = R*np.sin(AZ); Y = R*np.cos(AZ)
    if ax is None: fig, ax = plt.subplots(figsize=(5,5))
    ax.contourf(X,Y,amp,levels=[-1,0,1],colors=["white","#404040"])
    ax.contour(X,Y,amp,levels=[0],colors="k",linewidths=2)
    tc = np.linspace(0,2*np.pi,200)
    ax.plot(np.cos(tc),np.sin(tc),"k-",linewidth=2)
    if show_compass:
        for ang,lbl in [(0,"N"),(90,"E"),(180,"S"),(270,"W")]:
            ax.text(1.08*np.sin(np.radians(ang)),1.08*np.cos(np.radians(ang)),
                    lbl,fontsize=10,ha="center",va="center",fontweight="bold")
    if title is None: title=f"Strike={strike}, Dip={dip}, Rake={rake}"
    ax.set_title(title,fontsize=11)
    ax.set_xlim(-1.25,1.25); ax.set_ylim(-1.25,1.25)
    ax.set_aspect("equal"); ax.axis("off")
    return ax

print("Ready!")

import numpy as np import matplotlib.pyplot as plt from matplotlib.lines import Line2D def plot_beachball(strike, dip, rake, ax=None, title=None, show_compass=True): phi = np.radians(strike); delta = np.radians(dip); lam = np.radians(rake) n = np.array([-np.sin(delta)*np.sin(phi), np.sin(delta)*np.cos(phi), -np.cos(delta)]) s = np.array([np.cos(lam)*np.cos(phi)+np.sin(lam)*np.cos(delta)*np.sin(phi), np.cos(lam)*np.sin(phi)-np.sin(lam)*np.cos(delta)*np.cos(phi), -np.sin(lam)*np.sin(delta)]) M = np.outer(s, n) + np.outer(n, s) ng = 300; to = np.linspace(0, np.pi/2, ng); az = np.linspace(0, 2*np.pi, ng) TO, AZ = np.meshgrid(to, az) gn = np.sin(TO)*np.cos(AZ); ge = np.sin(TO)*np.sin(AZ); gd = np.cos(TO) amp = (M[0,0]*gn**2+M[1,1]*ge**2+M[2,2]*gd**2 +2*M[0,1]*gn*ge+2*M[0,2]*gn*gd+2*M[1,2]*ge*gd) R = np.sqrt(2)*np.sin(TO/2); X = R*np.sin(AZ); Y = R*np.cos(AZ) if ax is None: fig, ax = plt.subplots(figsize=(5,5)) ax.contourf(X,Y,amp,levels=[-1,0,1],colors=["white","#404040"]) ax.contour(X,Y,amp,levels=[0],colors="k",linewidths=2) tc = np.linspace(0,2*np.pi,200) ax.plot(np.cos(tc),np.sin(tc),"k-",linewidth=2) if show_compass: for ang,lbl in [(0,"N"),(90,"E"),(180,"S"),(270,"W")]: ax.text(1.08*np.sin(np.radians(ang)),1.08*np.cos(np.radians(ang)), lbl,fontsize=10,ha="center",va="center",fontweight="bold") if title is None: title=f"Strike={strike}, Dip={dip}, Rake={rake}" ax.set_title(title,fontsize=11) ax.set_xlim(-1.25,1.25); ax.set_ylim(-1.25,1.25) ax.set_aspect("equal"); ax.axis("off") return ax print("Ready!")

---

Exercise 1: Build a Beachball¶

Fill in the strike, dip, and rake to match each fault description. Recall:

• Rake: 0° = left-lateral, ±180° = right-lateral, −90° = normal, +90° = reverse
• Dip: 90° = vertical

1a. A vertical, right-lateral strike-slip fault striking East-West¶

In [ ]:

strike = ???
dip = ???
rake = ???

plot_beachball(strike, dip, rake)
plt.show()

strike = ??? dip = ??? rake = ??? plot_beachball(strike, dip, rake) plt.show()

1b. A normal fault striking North-South with a 45° dip¶

In [ ]:

strike = ???
dip = ???
rake = ???

plot_beachball(strike, dip, rake)
plt.show()

strike = ??? dip = ??? rake = ??? plot_beachball(strike, dip, rake) plt.show()

1c. The Hayward Fault (strike ≈ 325°, nearly vertical, right-lateral)¶

In [ ]:

strike = ???
dip = ???
rake = ???

plot_beachball(strike, dip, rake, title="Hayward Fault")
plt.show()

strike = ??? dip = ??? rake = ??? plot_beachball(strike, dip, rake, title="Hayward Fault") plt.show()

---

Exercise 2: Read Real Earthquake Beachballs¶

The four beachballs below are from real earthquakes. For each one, identify:

1. The fault type (strike-slip, normal, reverse/thrust, or oblique)
2. A tectonic setting where you would expect this type of faulting

Run the cell below to see the beachballs, then answer the questions.

Hint: If you get stuck, check the USGS event pages:

• A: nc75240492
• B: official20110311054624120_30
• C: usp000gvtu
• D: usp000g650

In [ ]:

fig, axes = plt.subplots(1, 4, figsize=(16, 4))

earthquakes = [
    (55,  85,  10, "Earthquake A\n(California, 2025)"),
    (193, 15,  63, "Earthquake B\n(Japan, 2011)"),
    (355, 47, -93, "Earthquake C\n(Italy, 2009)"),
    (325, 50,  90, "Earthquake D\n(China, 2008)"),
]

for ax, (s, d, r, title) in zip(axes, earthquakes):
    plot_beachball(s, d, r, ax=ax, title=title, show_compass=False)

plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 4, figsize=(16, 4)) earthquakes = [ (55, 85, 10, "Earthquake A\n(California, 2025)"), (193, 15, 63, "Earthquake B\n(Japan, 2011)"), (355, 47, -93, "Earthquake C\n(Italy, 2009)"), (325, 50, 90, "Earthquake D\n(China, 2008)"), ] for ax, (s, d, r, title) in zip(axes, earthquakes): plot_beachball(s, d, r, ax=ax, title=title, show_compass=False) plt.tight_layout() plt.show()

Earthquake A: What fault type? What tectonic setting?

Your answer here:

Earthquake B: What fault type? What tectonic setting?

Your answer here:

Earthquake C: What fault type? What tectonic setting?

Your answer here:

Earthquake D: What fault type? What tectonic setting?

Your answer here:

---

Exercise 3: The 2025 Berkeley Earthquake¶

Run the cells below to determine the focal mechanism from real P-wave polarity data — the same analysis as the lecture.

In [ ]:

from urllib.request import urlretrieve
import xml.etree.ElementTree as ET
from obspy.clients.fdsn import Client
from obspy import UTCDateTime
from obspy.geodetics import gps2dist_azimuth

# Earthquake parameters
origin_time = UTCDateTime("2025-09-22T09:56:12.52")
eq_lat, eq_lon, eq_depth = 37.863, -122.254, 7.6

# Download polarity data from USGS
url = ("https://earthquake.usgs.gov/product/phase-data/"
       "nc75240492/nc/1758546070900/quakeml.xml")
urlretrieve(url, "quakeml.xml")

tree = ET.parse("quakeml.xml")
root = tree.getroot()
ns = {"bed": "http://quakeml.org/xmlns/bed/1.2"}

polarity_picks = {}
for pick in root.findall(".//bed:pick", ns):
    pol_el = pick.find("bed:polarity", ns)
    if pol_el is None or pol_el.text not in ("positive", "negative"):
        continue
    wf = pick.find("bed:waveformID", ns)
    net = wf.get("networkCode")
    sta = wf.get("stationCode")
    polarity_picks[(net, sta)] = 1 if pol_el.text == "positive" else -1

# Get station coordinates
client = Client("NCEDC")
net_stations = {}
for net, sta in polarity_picks:
    net_stations.setdefault(net, []).append(sta)
station_coords = {}
for net, stations in net_stations.items():
    inv = client.get_stations(network=net, station=",".join(stations),
        starttime=UTCDateTime("2025-09-22"), endtime=UTCDateTime("2025-09-23"),
        level="station")
    for network in inv:
        for station in network:
            station_coords[(network.code, station.code)] = (
                station.latitude, station.longitude)

# Compute ray parameters
ray_params = []
for (net, sta), pol in polarity_picks.items():
    if (net, sta) not in station_coords:
        continue
    slat, slon = station_coords[(net, sta)]
    dist_m, az, baz = gps2dist_azimuth(eq_lat, eq_lon, slat, slon)
    dist_km = dist_m / 1000.0
    takeoff = np.degrees(np.arctan2(dist_km, eq_depth))
    ray_params.append({"polarity": pol, "takeoff": takeoff, "azimuth": baz})

print(f"Loaded {len(polarity_picks)} stations with polarity data")

from urllib.request import urlretrieve import xml.etree.ElementTree as ET from obspy.clients.fdsn import Client from obspy import UTCDateTime from obspy.geodetics import gps2dist_azimuth # Earthquake parameters origin_time = UTCDateTime("2025-09-22T09:56:12.52") eq_lat, eq_lon, eq_depth = 37.863, -122.254, 7.6 # Download polarity data from USGS url = ("https://earthquake.usgs.gov/product/phase-data/" "nc75240492/nc/1758546070900/quakeml.xml") urlretrieve(url, "quakeml.xml") tree = ET.parse("quakeml.xml") root = tree.getroot() ns = {"bed": "http://quakeml.org/xmlns/bed/1.2"} polarity_picks = {} for pick in root.findall(".//bed:pick", ns): pol_el = pick.find("bed:polarity", ns) if pol_el is None or pol_el.text not in ("positive", "negative"): continue wf = pick.find("bed:waveformID", ns) net = wf.get("networkCode") sta = wf.get("stationCode") polarity_picks[(net, sta)] = 1 if pol_el.text == "positive" else -1 # Get station coordinates client = Client("NCEDC") net_stations = {} for net, sta in polarity_picks: net_stations.setdefault(net, []).append(sta) station_coords = {} for net, stations in net_stations.items(): inv = client.get_stations(network=net, station=",".join(stations), starttime=UTCDateTime("2025-09-22"), endtime=UTCDateTime("2025-09-23"), level="station") for network in inv: for station in network: station_coords[(network.code, station.code)] = ( station.latitude, station.longitude) # Compute ray parameters ray_params = [] for (net, sta), pol in polarity_picks.items(): if (net, sta) not in station_coords: continue slat, slon = station_coords[(net, sta)] dist_m, az, baz = gps2dist_azimuth(eq_lat, eq_lon, slat, slon) dist_km = dist_m / 1000.0 takeoff = np.degrees(np.arctan2(dist_km, eq_depth)) ray_params.append({"polarity": pol, "takeoff": takeoff, "azimuth": baz}) print(f"Loaded {len(polarity_picks)} stations with polarity data")

In [ ]:

# Plot polarities on the focal sphere
def equal_area_project(takeoff_deg, azimuth_deg):
    r = np.sqrt(2) * np.sin(np.radians(takeoff_deg) / 2)
    return r * np.sin(np.radians(azimuth_deg)), r * np.cos(np.radians(azimuth_deg))

fig, ax = plt.subplots(figsize=(8, 8))
theta = np.linspace(0, 2*np.pi, 200)
ax.plot(np.cos(theta), np.sin(theta), "k-", linewidth=1.5)
for ang, lbl in [(0,"N"),(90,"E"),(180,"S"),(270,"W")]:
    ax.text(1.08*np.sin(np.radians(ang)), 1.08*np.cos(np.radians(ang)),
            lbl, fontsize=12, ha="center", va="center", fontweight="bold")
for r in ray_params:
    x, y = equal_area_project(r["takeoff"], r["azimuth"])
    if r["polarity"] == 1:
        ax.plot(x, y, "o", color="red", markersize=5,
                markeredgecolor="black", markeredgewidth=0.3)
    else:
        ax.plot(x, y, "o", color="white", markersize=5,
                markeredgecolor="blue", markeredgewidth=0.5)
ax.set_xlim(-1.25, 1.25); ax.set_ylim(-1.25, 1.25)
ax.set_aspect("equal"); ax.axis("off")
ax.set_title("P-Wave Polarities on the Focal Sphere", fontsize=13)
ax.legend(handles=[
    Line2D([0],[0], marker="o", color="w", markerfacecolor="red",
           markeredgecolor="black", markersize=8, label="Up (compression)"),
    Line2D([0],[0], marker="o", color="w", markerfacecolor="white",
           markeredgecolor="blue", markersize=8, label="Down (dilatation)"),
], loc="lower left", fontsize=11)
plt.tight_layout()
plt.show()

# Plot polarities on the focal sphere def equal_area_project(takeoff_deg, azimuth_deg): r = np.sqrt(2) * np.sin(np.radians(takeoff_deg) / 2) return r * np.sin(np.radians(azimuth_deg)), r * np.cos(np.radians(azimuth_deg)) fig, ax = plt.subplots(figsize=(8, 8)) theta = np.linspace(0, 2*np.pi, 200) ax.plot(np.cos(theta), np.sin(theta), "k-", linewidth=1.5) for ang, lbl in [(0,"N"),(90,"E"),(180,"S"),(270,"W")]: ax.text(1.08*np.sin(np.radians(ang)), 1.08*np.cos(np.radians(ang)), lbl, fontsize=12, ha="center", va="center", fontweight="bold") for r in ray_params: x, y = equal_area_project(r["takeoff"], r["azimuth"]) if r["polarity"] == 1: ax.plot(x, y, "o", color="red", markersize=5, markeredgecolor="black", markeredgewidth=0.3) else: ax.plot(x, y, "o", color="white", markersize=5, markeredgecolor="blue", markeredgewidth=0.5) ax.set_xlim(-1.25, 1.25); ax.set_ylim(-1.25, 1.25) ax.set_aspect("equal"); ax.axis("off") ax.set_title("P-Wave Polarities on the Focal Sphere", fontsize=13) ax.legend(handles=[ Line2D([0],[0], marker="o", color="w", markerfacecolor="red", markeredgecolor="black", markersize=8, label="Up (compression)"), Line2D([0],[0], marker="o", color="w", markerfacecolor="white", markeredgecolor="blue", markersize=8, label="Down (dilatation)"), ], loc="lower left", fontsize=11) plt.tight_layout() plt.show()

In [ ]:

# Grid search for best-fitting focal mechanism
def moment_tensor(strike_deg, dip_deg, rake_deg):
    phi=np.radians(strike_deg); delta=np.radians(dip_deg); lam=np.radians(rake_deg)
    n=np.array([-np.sin(delta)*np.sin(phi),np.sin(delta)*np.cos(phi),-np.cos(delta)])
    s=np.array([np.cos(lam)*np.cos(phi)+np.sin(lam)*np.cos(delta)*np.sin(phi),
                np.cos(lam)*np.sin(phi)-np.sin(lam)*np.cos(delta)*np.cos(phi),
               -np.sin(lam)*np.sin(delta)])
    return np.outer(s,n)+np.outer(n,s)

to_rad = np.radians([r["takeoff"] for r in ray_params])
az_rad = np.radians([r["azimuth"] for r in ray_params])
observed = np.array([r["polarity"] for r in ray_params])
gn = np.sin(to_rad)*np.cos(az_rad)
ge = np.sin(to_rad)*np.sin(az_rad)
gd = np.cos(to_rad)

best_score = 0
best_fm = (0, 0, 0)
for strike in range(0, 360, 5):
    for dip in range(5, 91, 5):
        for rake in range(-180, 181, 10):
            M = moment_tensor(strike, dip, rake)
            amp = (M[0,0]*gn**2+M[1,1]*ge**2+M[2,2]*gd**2
                   +2*M[0,1]*gn*ge+2*M[0,2]*gn*gd+2*M[1,2]*ge*gd)
            score = np.sum(np.sign(amp) == observed)
            if score > best_score:
                best_score = score
                best_fm = (strike, dip, rake)

print(f"Best-fitting focal mechanism:")
print(f"  Strike = {best_fm[0]},  Dip = {best_fm[1]},  Rake = {best_fm[2]}")
print(f"  Matching: {best_score}/{len(ray_params)} ({100*best_score/len(ray_params):.1f}%)")

# Grid search for best-fitting focal mechanism def moment_tensor(strike_deg, dip_deg, rake_deg): phi=np.radians(strike_deg); delta=np.radians(dip_deg); lam=np.radians(rake_deg) n=np.array([-np.sin(delta)*np.sin(phi),np.sin(delta)*np.cos(phi),-np.cos(delta)]) s=np.array([np.cos(lam)*np.cos(phi)+np.sin(lam)*np.cos(delta)*np.sin(phi), np.cos(lam)*np.sin(phi)-np.sin(lam)*np.cos(delta)*np.cos(phi), -np.sin(lam)*np.sin(delta)]) return np.outer(s,n)+np.outer(n,s) to_rad = np.radians([r["takeoff"] for r in ray_params]) az_rad = np.radians([r["azimuth"] for r in ray_params]) observed = np.array([r["polarity"] for r in ray_params]) gn = np.sin(to_rad)*np.cos(az_rad) ge = np.sin(to_rad)*np.sin(az_rad) gd = np.cos(to_rad) best_score = 0 best_fm = (0, 0, 0) for strike in range(0, 360, 5): for dip in range(5, 91, 5): for rake in range(-180, 181, 10): M = moment_tensor(strike, dip, rake) amp = (M[0,0]*gn**2+M[1,1]*ge**2+M[2,2]*gd**2 +2*M[0,1]*gn*ge+2*M[0,2]*gn*gd+2*M[1,2]*ge*gd) score = np.sum(np.sign(amp) == observed) if score > best_score: best_score = score best_fm = (strike, dip, rake) print(f"Best-fitting focal mechanism:") print(f" Strike = {best_fm[0]}, Dip = {best_fm[1]}, Rake = {best_fm[2]}") print(f" Matching: {best_score}/{len(ray_params)} ({100*best_score/len(ray_params):.1f}%)")

In [ ]:

# Compare our result with the USGS catalog
n_grid = 300
az_g = np.linspace(0, 360, n_grid); to_g = np.linspace(0, 90, n_grid//2)
AZ_g, TO_g = np.meshgrid(az_g, to_g)
R_g = np.sqrt(2)*np.sin(np.radians(TO_g)/2)
X_g, Y_g = R_g*np.sin(np.radians(AZ_g)), R_g*np.cos(np.radians(AZ_g))
gn_g = np.sin(np.radians(TO_g))*np.cos(np.radians(AZ_g))
ge_g = np.sin(np.radians(TO_g))*np.sin(np.radians(AZ_g))
gd_g = np.cos(np.radians(TO_g))

def radiation_on_grid(s, d, r):
    M = moment_tensor(s, d, r)
    return (M[0,0]*gn_g**2+M[1,1]*ge_g**2+M[2,2]*gd_g**2
            +2*M[0,1]*gn_g*ge_g+2*M[0,2]*gn_g*gd_g+2*M[1,2]*ge_g*gd_g)
def polarity_score(s, d, r):
    M = moment_tensor(s, d, r)
    a = (M[0,0]*gn**2+M[1,1]*ge**2+M[2,2]*gd**2
         +2*M[0,1]*gn*ge+2*M[0,2]*gn*gd+2*M[1,2]*ge*gd)
    return int(np.sum(np.sign(a)==observed))

fig, axes = plt.subplots(1, 2, figsize=(15, 7))
theta = np.linspace(0, 2*np.pi, 200)
for ax, ((s,d,r), label) in zip(axes, [
    (best_fm, "Our Result"), ((55,85,10), "USGS Catalog")]):
    AMP = radiation_on_grid(s, d, r)
    ax.contourf(X_g,Y_g,AMP,levels=[-1,0,1],colors=["white","lightcoral"],alpha=0.4)
    ax.contour(X_g,Y_g,AMP,levels=[0],colors="black",linewidths=2)
    ax.plot(np.cos(theta),np.sin(theta),"k-",linewidth=1.5)
    for ang,lbl in [(0,"N"),(90,"E"),(180,"S"),(270,"W")]:
        ax.text(1.08*np.sin(np.radians(ang)),1.08*np.cos(np.radians(ang)),
                lbl,fontsize=11,ha="center",va="center",fontweight="bold")
    for rp in ray_params:
        x,y = equal_area_project(rp["takeoff"],rp["azimuth"])
        if rp["polarity"]==1:
            ax.plot(x,y,"o",color="red",markersize=4,markeredgecolor="black",markeredgewidth=0.3)
        else:
            ax.plot(x,y,"o",color="white",markersize=4,markeredgecolor="blue",markeredgewidth=0.5)
    sc = polarity_score(s, d, r)
    ax.set_xlim(-1.25,1.25); ax.set_ylim(-1.25,1.25)
    ax.set_aspect("equal"); ax.axis("off")
    ax.set_title(f"{label}\nStrike={s}, Dip={d}, Rake={r}\n({sc}/{len(ray_params)} match)",fontsize=12)
plt.tight_layout()
plt.show()

# Compare our result with the USGS catalog n_grid = 300 az_g = np.linspace(0, 360, n_grid); to_g = np.linspace(0, 90, n_grid//2) AZ_g, TO_g = np.meshgrid(az_g, to_g) R_g = np.sqrt(2)*np.sin(np.radians(TO_g)/2) X_g, Y_g = R_g*np.sin(np.radians(AZ_g)), R_g*np.cos(np.radians(AZ_g)) gn_g = np.sin(np.radians(TO_g))*np.cos(np.radians(AZ_g)) ge_g = np.sin(np.radians(TO_g))*np.sin(np.radians(AZ_g)) gd_g = np.cos(np.radians(TO_g)) def radiation_on_grid(s, d, r): M = moment_tensor(s, d, r) return (M[0,0]*gn_g**2+M[1,1]*ge_g**2+M[2,2]*gd_g**2 +2*M[0,1]*gn_g*ge_g+2*M[0,2]*gn_g*gd_g+2*M[1,2]*ge_g*gd_g) def polarity_score(s, d, r): M = moment_tensor(s, d, r) a = (M[0,0]*gn**2+M[1,1]*ge**2+M[2,2]*gd**2 +2*M[0,1]*gn*ge+2*M[0,2]*gn*gd+2*M[1,2]*ge*gd) return int(np.sum(np.sign(a)==observed)) fig, axes = plt.subplots(1, 2, figsize=(15, 7)) theta = np.linspace(0, 2*np.pi, 200) for ax, ((s,d,r), label) in zip(axes, [ (best_fm, "Our Result"), ((55,85,10), "USGS Catalog")]): AMP = radiation_on_grid(s, d, r) ax.contourf(X_g,Y_g,AMP,levels=[-1,0,1],colors=["white","lightcoral"],alpha=0.4) ax.contour(X_g,Y_g,AMP,levels=[0],colors="black",linewidths=2) ax.plot(np.cos(theta),np.sin(theta),"k-",linewidth=1.5) for ang,lbl in [(0,"N"),(90,"E"),(180,"S"),(270,"W")]: ax.text(1.08*np.sin(np.radians(ang)),1.08*np.cos(np.radians(ang)), lbl,fontsize=11,ha="center",va="center",fontweight="bold") for rp in ray_params: x,y = equal_area_project(rp["takeoff"],rp["azimuth"]) if rp["polarity"]==1: ax.plot(x,y,"o",color="red",markersize=4,markeredgecolor="black",markeredgewidth=0.3) else: ax.plot(x,y,"o",color="white",markersize=4,markeredgecolor="blue",markeredgewidth=0.5) sc = polarity_score(s, d, r) ax.set_xlim(-1.25,1.25); ax.set_ylim(-1.25,1.25) ax.set_aspect("equal"); ax.axis("off") ax.set_title(f"{label}\nStrike={s}, Dip={d}, Rake={r}\n({sc}/{len(ray_params)} match)",fontsize=12) plt.tight_layout() plt.show()

Questions¶

The USGS reported two nodal planes for this earthquake:

	Strike	Dip	Rake
NP1	55°	85°	10°
NP2	324°	80°	174°

Question 3a: What type of faulting does the beachball show? Is this consistent with what you know about the Hayward Fault?

Your answer here:

Question 3b: The Hayward Fault strikes approximately 325° from North. Which nodal plane (NP1 or NP2) corresponds to the actual Hayward Fault? How do you know?

Your answer here:

Question 3c: From the P-wave polarity data alone (without knowing the Hayward Fault exists), could you determine which nodal plane is the real fault? What additional information would you need?

Your answer here: