Skip to content

Add option to define custom parameter for Seismological Source Theory RVT Motion #53

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jimxia7 wants to merge 3 commits into
base: main
Choose a base branch
from
Open

Add option to define custom parameter for Seismological Source Theory RVT Motion #53

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0a03d61
Updated test_calc_compatible_spectra
jimxia7 Apr 3, 2026
1487465
Added option for custom parameter for Source Theory RVT Motion
jimxia7 Apr 3, 2026
2a372d3
Format/lint/test
jimxia7 Apr 3, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/conf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@
from datetime import date

# Add the src directory to the path
sys.path.insert(0, os.path.abspath('../src'))
sys.path.insert(0, os.path.abspath("../src"))

import pyrvt

Expand Down
229 changes: 141 additions & 88 deletions src/pyrvt/motions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,21 +109,35 @@ def calc_sdof_tf(
)


def calc_stress_drop(magnitude: float) -> float:
"""Stress drop using Atkinson & Boore (2011) model.
def calc_stress_drop(magnitude: float, method: str = "Atkinson&Boore2011") -> float:
"""Calculating Stress drops.

Parameters
----------
magnitude : float
Moment magnitude of the stress drop.
Moment magnitude of the stress drop.

Method: String
Atkinson&Boore2011: Atkinson & Boore (2011) model
Stafford2022: Stafford et al. (2022) model

Returns
-------
stress_drop : float
Stress drop (bars).

"""
return 10 ** (3.45 - 0.2 * max(magnitude, 5.0))

if method == "Atkinson&Boore2011":
return 10 ** (3.45 - 0.2 * max(magnitude, 5.0))
elif method == "Stafford2022":
return np.exp(2.296 + 0.4624 * np.min([magnitude - 5.0, 0.0])) * 10
elif method is None:
return 10 ** (3.45 - 0.2 * max(magnitude, 5.0))
else:
raise NotImplementedError(
f"Method {method} not implemented for stress drop calculation."
)


def calc_geometric_spreading(
Expand Down Expand Up @@ -478,12 +492,12 @@ def __init__(
magnitude: float,
distance: float,
region: str,
stress_drop: float | None = None,
depth: float | None = 8,
peak_calculator: str | peak_calculators.Calculator | None = None,
calc_kwds: dict | None = None,
freqs: npt.ArrayLike | None = None,
disable_site_amp: bool = False,
**kwargs,
) -> None:
"""Initialize the motion.

Expand All @@ -496,11 +510,6 @@ def __init__(
region : str
Region for the parameters. Either 'cena' for Central and Eastern North
America, or 'wna' for Western North America.
stress_drop : float, optional
Stress drop of the event (bars). If `None`, then the default value
is used. For `region` is 'cena', the default value is computed by
the model, while for `region` is 'wna' the
default value is 100 bars.
depth : float, optional
Hypocenter depth (km). The `depth` is combined with the `distance`
to compute the hypocentral distance.
Expand All @@ -519,6 +528,25 @@ def __init__(
disable_site_amp: bool, optional
if the crustal site amplification should be disable. Defaults to *False*.

**kwargs: optional arguments
The following parameters can be specified:
- `shear_velocity`: shear velocity (km/s).
- `path_atten_coeff`: path attenuation coefficient.
Quality Factor = path_atten_coeff * f^path_atten_power
- `path_atten_power`: path attenuation power.
Quality Factor = path_atten_coeff * f^path_atten_power
- `density`: density (g/cm³).
- `site_atten`: site attenuation (s).
- `geometric_spreading`: geometric spreading parameters.
- 'stress_drop': float, optional
Stress drop of the event (bars). If `None`, then the default value
is used. For `region` is 'cena', the default value is computed by
the model, while for `region` is 'wna' the
default value is 100 bars.
- 'stress_drop_method': define how the stress drop is computed
useless if stress_drop is specified.
- 'site_amp': site amplification function.

"""
super().__init__(peak_calculator=peak_calculator, calc_kwds=calc_kwds)

Expand All @@ -529,111 +557,136 @@ def __init__(

if self.region == "wna":
# Default parameters for the WUS from Campbell (2003)
self.shear_velocity = 3.5
self.path_atten_coeff = 180.0
self.path_atten_power = 0.45
self.density = 2.8
self.site_atten = 0.04
self.shear_velocity = kwargs.get("shear_velocity", 3.5)
self.path_atten_coeff = kwargs.get("path_atten_coeff", 180.0)
self.path_atten_power = kwargs.get("path_atten_power", 0.45)
self.density = kwargs.get("density", 2.8)
self.site_atten = kwargs.get("site_atten", 0.04)

self.geometric_spreading = kwargs.get(
"geometric_spreading", [(1, 40), (0.5, None)]
)

self.geometric_spreading = [(1, 40), (0.5, None)]
stress_drop = kwargs.get("stress_drop", None)
stress_drop_method = kwargs.get("stress_drop_method", None)

if stress_drop:
self.stress_drop = stress_drop
elif stress_drop_method:
self.stress_drop = calc_stress_drop(
magnitude, method=stress_drop_method
)
else:
self.stress_drop = 100.0

# Crustal amplification from Campbell (2003) using the
# log-frequency and the amplification based on a quarter-wave
# length approximation
self.site_amp = make_interp_spline(
np.log(
self.site_amp = kwargs.get(
"site_amp",
make_interp_spline(
np.log(
[
0.01,
0.09,
0.16,
0.51,
0.84,
1.25,
2.26,
3.17,
6.05,
16.60,
61.20,
100.00,
]
),
[
0.01,
0.09,
0.16,
0.51,
0.84,
1.25,
2.26,
3.17,
6.05,
16.60,
61.20,
100.00,
]
1.00,
1.10,
1.18,
1.42,
1.58,
1.74,
2.06,
2.25,
2.58,
3.13,
4.00,
4.40,
],
k=1,
),
[
1.00,
1.10,
1.18,
1.42,
1.58,
1.74,
2.06,
2.25,
2.58,
3.13,
4.00,
4.40,
],
k=1,
)

elif self.region == "cena":
# Default parameters for the CEUS from Campbell (2003)
self.shear_velocity = 3.6
self.density = 2.8
self.path_atten_coeff = 680.0
self.path_atten_power = 0.36
self.site_atten = 0.006
self.shear_velocity = kwargs.get("shear_velocity", 3.6)
self.path_atten_coeff = kwargs.get("path_atten_coeff", 680.0)
self.path_atten_power = kwargs.get("path_atten_power", 0.36)
self.density = kwargs.get("density", 2.8)
self.site_atten = kwargs.get("site_atten", 0.006)

self.geometric_spreading = kwargs.get(
"geometric_spreading", [(1, 70), (0, 130), (0.5, None)]
)

self.geometric_spreading = [(1, 70), (0, 130), (0.5, None)]
stress_drop = kwargs.get("stress_drop", None)
stress_drop_method = kwargs.get("stress_drop_method", None)

if stress_drop:
self.stress_drop = stress_drop
elif stress_drop_method:
self.stress_drop = calc_stress_drop(
magnitude, method=stress_drop_method
)
else:
self.stress_drop = calc_stress_drop(magnitude)

# Crustal amplification from Campbell (2003) using the
# log-frequency and the amplification based on a quarter-wave
# length approximation
self.site_amp = make_interp_spline(
np.log(
self.site_amp = kwargs.get(
"site_amp",
make_interp_spline(
np.log(
[
0.01,
0.10,
0.20,
0.30,
0.50,
0.90,
1.25,
1.80,
3.00,
5.30,
8.00,
14.00,
30.00,
60.00,
100.00,
]
),
[
0.01,
0.10,
0.20,
0.30,
0.50,
0.90,
1.25,
1.80,
3.00,
5.30,
8.00,
14.00,
30.00,
60.00,
100.00,
]
1.00,
1.02,
1.03,
1.05,
1.07,
1.09,
1.11,
1.12,
1.13,
1.14,
1.15,
1.15,
1.15,
1.15,
1.15,
],
k=1,
),
[
1.00,
1.02,
1.03,
1.05,
1.07,
1.09,
1.11,
1.12,
1.13,
1.14,
1.15,
1.15,
1.15,
1.15,
1.15,
],
k=1,
)

else:
Expand Down
5 changes: 3 additions & 2 deletions src/pyrvt/peak_calculators.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,6 @@
from scipy.interpolate import LinearNDInterpolator
from scipy.signal import argrelmax


# From: https://github.com/scipy/scipy/issues/4831#issuecomment-258501648
# Set the signature of the future types function
# The second argument must be a pointer or it won't work right,
Expand Down Expand Up @@ -1311,7 +1310,9 @@ def _calc_peak_factor(
"""
m0, m2 = sspectrum.moments(0, 2)

tE = 4 * np.trapezoid(np.square(sspectrum.squared_fa), x=sspectrum.freqs) / m0**2
tE = (
4 * np.trapezoid(np.square(sspectrum.squared_fa), x=sspectrum.freqs) / m0**2
)

# Central frequency
freq_cent = np.sqrt(m2 / m0) / (2 * np.pi)
Expand Down
4 changes: 2 additions & 2 deletions tests/test_tools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -280,8 +280,8 @@ def test_calc_compatible_spectra(method):
method, periods, events[:1], 0.05
)

# Test that the fit is within 2% of the target
assert_allclose(events[0]["psa"], events_mod[0]["psa_calc"], rtol=0.05)
# Test that the fit is within 5% of the target
assert_allclose(events[0]["psa"], events_mod[0]["psa_calc"], rtol=0.05, atol=1e-4)


@pytest.mark.parametrize("fixed_spacing", [True, False])
Expand Down