Changed Sea21 to g-s, the standard unit system of pyrvt.

src/pyrvt/motions.py

@@ -12,6 +12,7 @@
 
 import numpy as np
 import numpy.typing as npt
+from scipy.constants import g as gravity
 from scipy.interpolate import interp1d
 from scipy.stats import linregress
 
@@ -245,15 +246,12 @@ def calc_osc_accels(
             Peak pseudo-spectral acceleration of the oscillator
 
         """
-        if trans_func is not None:
-            tf = np.asarray(trans_func)
-        else:
-            tf = np.ones_like(self.freqs)
+        trans_func = 1 if trans_func is None else np.asarray(trans_func)
 
         resp = np.array(
             [
                 self.calc_peak(
-                    tf * calc_sdof_tf(self.freqs, of, osc_damping),
+                    trans_func * calc_sdof_tf(self.freqs, of, osc_damping),
                     osc_freq=of,
                     osc_damping=osc_damping,
                     site_tf=trans_func,
@@ -353,6 +351,7 @@ def __init__(
         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,
     ):
         """Initialize the motion.
@@ -383,6 +382,11 @@ def __init__(
         calc_kwds : dict, optional
             Keywords to be passed during the creation the peak calculator.
             These keywords are only required for some peak calculators.
+        freqs : array_like
+            frequencies for which the Fourier amplitude spectrum should be computed.
+            Defaults to `np.geomspace(0.05, 200, 512)`
+        disable_site_amp: bool, optional
+            if the crustal site amplification should be disable. Defaults to *False*.
 
         """
         super().__init__(peak_calculator=peak_calculator, calc_kwds=calc_kwds)
@@ -517,6 +521,10 @@ def __init__(
             * (self.stress_drop / self.seismic_moment) ** (1.0 / 3.0)
         )
 
+        # Combine the three components and convert from displacement to acceleration
+        self.calc_fourier_amps(freqs)
+        self._duration = self.calc_duration()
+
     def calc_duration(self) -> float:
         """Compute the duration by combination of source and path.
 
@@ -607,9 +615,10 @@ def calc_fourier_amps(self, freqs: npt.ArrayLike | None = None) -> np.ndarray:
         site_comp = 1 if self._disable_site_amp else (site_amp * site_dim)
 
         # Conversion factor to convert from dyne-cm into gravity-sec
-        conv = 1.0e-20 / 980.7
+        conv = 1.0e-20 / (100 * gravity)
         # Combine the three components and convert from displacement to
         # acceleration
+
         self._fourier_amps = (
             conv
             * (2.0 * np.pi * self._freqs) ** 2.0
@@ -789,8 +798,8 @@ def __init__(
 
         path_comp = geom_spread * anelastic_atten
 
-        # Convert to acceleration (cm/s)
-        conv = (2 * np.pi * self._freqs) ** 2
+        # Convert to acceleration (cm/s) and then into g-se
+        conv = (2 * np.pi * self._freqs) ** 2 / (gravity * 100)
 
         if disable_site_amp:
             site_tf = 1.0
@@ -852,15 +861,15 @@ def calc_depth_tor(mag: float, mechanism: str) -> float:
         """
         if mechanism == "RS":
             # Reverse and reverse-oblique faulting
-            depth_tor = 2.704 - 1.226 * max(mag - 5.849, 0)
+            fact = 2.704 - 1.226 * max(mag - 5.849, 0)
         else:
             # Combined strike-slip and normal faulting
-            depth_tor = 2.673 - 1.136 * max(mag - 4.970, 0)
+            fact = 2.673 - 1.136 * max(mag - 4.970, 0)
 
-        return max(depth_tor, 0) ** 2
+        return max(fact, 0) ** 2
 
     @staticmethod
-    def calc_duration(corner_freq, dist_ps):
+    def calc_duration(corner_freq: float, dist_ps: float) -> float:
         """Boore & Thomspson (2014) duration model."""
 
         # Source component. Equation 2