Slip Distribution Code

slip.py

import numpy as np
import matplotlib.pyplot as plt


# Parameters
def von_karman_transfer(field, L, dx, H=0.75):
    # Get field dimensions
    ny, nx = field.shape
    kx = np.fft.fftfreq(nx, d=dx) * 2 * np.pi  # Wavenumbers in x
    ky = np.fft.fftfreq(ny, d=dx) * 2 * np.pi  # Wavenumbers in y
    kx, ky = np.meshgrid(kx, ky, indexing="xy")

    # Radial wavenumber

    # Von Kármán transfer function
    return (
        kx,
        ky,
        np.sqrt(L.prod() / (1 + (L[0] * kx) ** 2 + (L[1] * ky) ** 2) ** (H + 1)),
    )


def von_karman_filter_2d(field, L, dx, H=0.75):
    """
    Apply a von Kármán filter to a 2D field.

    Parameters:
    - field: 2D input array (e.g., spatial field).
    - L: Correlation length (outer scale).
    - dx: Sampling interval (assumed equal in both dimensions).

    Returns:
    - Filtered 2D field.
    """

    # Fourier transform of the field
    field_ft = np.fft.fft2(field)

    # Apply the filter
    filtered_ft = field_ft * von_karman_transfer(field, L, dx, H)

    # Inverse Fourier transform
    filtered_field = np.fft.ifft2(filtered_ft).real

    return filtered_field


def genslip_filter(geometry, field, mw, dx):
    a = 10 ** (np.array([0.5 * mw - 1.7, 0.333 * mw - 0.7]))
    kx, ky, h_k = von_karman_transfer(field, a, dx)
    ny, nx = field.shape
    c = np.array([nx * dx / 2, ny * dx / 2])

    f = np.reciprocal(1 + c[0] ** 2 * kx**2 + c[1] ** 2 * ky**2)
    field_ft = np.fft.fft2(field)
    dc = field_ft[0, 0]
    theta = np.random.uniform(-np.pi, np.pi, size=field.shape)
    field_st = dc / np.sqrt(a.prod()) * h_k * np.exp(1j * theta)
    spatial_slip = np.fft.ifft2(field_ft * f + field_st * (1 - f)).real
    mean_slip = spatial_slip.mean()
    slip_stddev = np.std(spatial_slip)
    target_stddev = 0.85 * mean_slip
    spatial_slip = mean_slip + (spatial_slip - mean_slip) * (
        target_stddev / slip_stddev
    )
    spatial_slip -= spatial_slip.min()

    area = np.prod(field.shape) * (dx * 1000) ** 2  # m^2
    moment = 3.3e10 * area * (spatial_slip / 100).mean()  # mu * m^2 * m

    print(moment)
    target_moment = 10 ** ((3 / 2) * (mw + 6.03333))
    print(target_moment)
    spatial_slip *= target_moment / moment
    return spatial_slip


m = 6.9
taper = 0.05

length = 400

width = 200

taper_width = 10


# Create the initial slip array

slip = np.ones((width, length))


# Create taper weights for the bottom

bottom_taper = np.ones(width)

bottom_taper[-taper_width:] = np.linspace(1, taper, taper_width)


# Apply bottom taper

slip *= bottom_taper[:, np.newaxis]


# Create taper weights for the left and right

left_right_taper = np.ones(length)

left_right_taper[:taper_width] = np.linspace(taper, 1, taper_width)  # Left taper

left_right_taper[-taper_width:] = np.linspace(1, taper, taper_width)  # Right taper


# Apply left and right tapers

slip *= left_right_taper[np.newaxis, :]

genslip = genslip_filter(slip, m, 0.1)
hypocentre = np.array([0.5, 0.5])
vs = 3.5  # km/s


print(genslip.mean())
print(np.std(genslip))


plt.figure(figsize=(10, 5))
plt.imshow(genslip, cmap="hot", origin="lower", aspect="auto")
plt.colorbar(label="Slip Value")
plt.title("Heatmap of Tapered Slip Array")
plt.xlabel("Length")
plt.ylabel("Width")
plt.savefig("test.png")