Kappa correction

src/pystrata/motion.py

@@ -23,15 +23,119 @@
 
 import enum
 import re
+import warnings
 
 import numpy as np
 import pyrvt
+import pykooh
 
 # Gravity in m/sec²
 from scipy.constants import g as GRAVITY
 
 _trapezoid = np.trapezoid
 
+# Integers and floats, including values without a leading digit (e.g., ".0100")
+# and Fortran style exponents (e.g., "1.0D-2").
+_RE_NUMBER = re.compile(r"[+-]?(?:\d+\.?\d*|\.\d+)(?:[eEdD][+-]?\d+)?")
+
+# Default frequency range for calculation kappa: [10,30]
+DEFAULT_KAPPA_FREQS = np.logspace(np.log10(10), np.log10(30), 100)
+
+def _compute_fourier_spectrum(time_step,
+                              accels,
+                              freqs = None,
+                              fa_length=None, 
+                              ko_bandwidth = None):
+    """Compute the Fourier Amplitude Spectrum of the time series."""
+
+    if fa_length is None:
+        # Use the next power of 2 for the length
+        n = 1
+        while n < accels.size:
+            n <<= 1
+    else:
+        n = fa_length
+
+    fft_freqs = np.fft.rfftfreq(n, d = time_step)
+
+    if freqs is None:
+        freqs = fft_freqs
+
+    if ko_bandwidth is None:
+        FAS = np.interp(freqs, 
+                        fft_freqs, 
+                        np.fft.rfft(accels, n))
+    else:
+        FAS = pykooh.smooth(freqs, 
+                            fft_freqs, 
+                            np.fft.rfft(accels, n),
+                            ko_bandwidth)
+
+    return freqs, FAS
+
+
+
+def _to_float(text):
+    """Convert a string to a float, permitting a Fortran style exponent."""
+    return float(text.replace("D", "E").replace("d", "e"))
+
+
+def _parse_at2_header(line):
+    """Parse the point count and time step from the header of an AT2 file.
+
+    Both of the PEER NGA layouts are supported::
+
+        4096    0.0100    NPTS, DT
+        NPTS=   5346, DT=   .0100 SEC,
+
+    as are variations that reverse the order of the two values, or that omit
+    the commas separating them.
+
+    Parameters
+    ----------
+    line: str
+        Fourth line of an AT2 file.
+
+    Returns
+    -------
+    npts: int
+        Number of points in the time series.
+    time_step: float
+        Time step of the time series [sec].
+    """
+    # Values that follow their label -- e.g., "NPTS= 5346" or "DT .0100". Each
+    # value is located by its own label, so their order does not matter.
+    found = {}
+    for key in ("NPTS", "DT"):
+        m = re.search(
+            r"\b" + key + r"\b\s*[=:]?\s*(" + _RE_NUMBER.pattern + ")",
+            line,
+            re.IGNORECASE,
+        )
+        if m:
+            found[key] = _to_float(m.group(1))
+
+    if len(found) < 2:
+        # Values that precede their labels -- e.g., "4096  0.0100  NPTS, DT".
+        values = [_to_float(v) for v in _RE_NUMBER.findall(line)]
+        if len(values) < 2:
+            raise ValueError(f"Unable to parse NPTS and DT from AT2 header: {line!r}")
+
+        values = values[:2]
+        upper = line.upper()
+        pos = {key: upper.find(key) for key in ("NPTS", "DT")}
+        if all(p >= 0 for p in pos.values()):
+            # Pair the values with the labels by order of appearance.
+            keys = sorted(pos, key=lambda key: pos[key])
+        else:
+            # Unlabeled, so rely on magnitude: the time step is the smaller of
+            # the two.
+            keys = ["DT", "NPTS"] if values[0] < values[1] else ["NPTS", "DT"]
+
+        found = dict(zip(keys, values))
+
+    return int(found["NPTS"]), found["DT"]
+
 
 class WaveField(enum.Enum):
     outcrop = 0
@@ -43,7 +147,7 @@ class Motion:
     def __init__(self, freqs=None):
         object.__init__(self)
 
-        self._freqs = np.array([] if freqs is None else freqs)
+        self._freqs = None if freqs is None else np.array(freqs)
         self._pga = None
         self._pgv = None
         self._arias_intensity = None
@@ -93,7 +197,7 @@ class TimeSeriesMotion(Motion):
     """Time-series motion for time series based site response analysis."""
 
     def __init__(
-        self, filename: str, description: str, time_step: float, accels, fa_length=None
+        self, filename: str, description: str, time_step: float, accels
     ):
         """Initialize the class from specified acceleration values.
 
@@ -120,8 +224,8 @@ def __init__(
         self._description = description
         self._time_step = time_step
         self._accels = np.asarray(accels)
-
-        self._calc_fourier_spectrum(fa_length)
+        self._kappa = None
+        self._fourier_amps = None
 
     @property
     def accels(self):
@@ -145,6 +249,7 @@ def times(self):
 
     @property
     def freqs(self):
+
         """Return the frequencies."""
         if self._freqs is None:
             self._calc_fourier_spectrum()
@@ -226,22 +331,43 @@ def calc_osc_accels(self, osc_freqs, osc_damping=0.05, tf=None):
         )
         return resp
 
-    def _calc_fourier_spectrum(self, fa_length=None):
+    def _calc_fourier_spectrum(self, 
+                               freqs = None, 
+                               fa_length = None, 
+                               ko_bandwidth = None):
         """Compute the Fourier Amplitude Spectrum of the time series."""
 
-        if fa_length is None:
-            # Use the next power of 2 for the length
-            n = 1
-            while n < self.accels.size:
-                n <<= 1
-        else:
-            n = fa_length
-
-        self._fourier_amps = np.fft.rfft(self._accels, n)
+        self._freqs, self._fourier_amps = _compute_fourier_spectrum(
+            self.time_step,
+            self._accels,
+            freqs = freqs,
+            fa_length= fa_length,
+            ko_bandwidth= ko_bandwidth
+        )
 
-        freq_step = 1.0 / (2 * self._time_step * (n / 2))
-        self._freqs = freq_step * np.arange(1 + n / 2)
+    @property
+    def kappa(self):
+
+        if self._kappa is None:
+            self._calc_kappa()
+
+        return self._kappa
+
+    def _calc_kappa(self,
+                    freqs_range = DEFAULT_KAPPA_FREQS, 
+                    fa_length=None, 
+                    ko_bandwidth = None):
+
+        _,fas = _compute_fourier_spectrum(
+            self.time_step,
+            self.accels,
+            freqs = freqs_range,
+            fa_length=fa_length,
+            ko_bandwidth=ko_bandwidth
+        )
 
+        self._kappa = -np.polyfit(freqs_range, np.log(abs(fas)),1)[0]/np.pi
+
     def _calc_sdof_tf(self, osc_freq, damping=0.05):
         """Compute the transfer function for a single-degree-of-freedom oscillator.
 
@@ -270,6 +396,15 @@ def _calc_sdof_tf(self, osc_freq, damping=0.05):
     def load_at2_file(cls, filename, scale=1.0):
         """Read an AT2 formatted time series.
 
+        The fourth line of the file provides the number of points and the time
+        step. Both of the PEER NGA layouts are read::
+
+            4096    0.0100    NPTS, DT
+            NPTS=   5346, DT=   .0100 SEC,
+
+        as are variations that reverse the order of the two values, or that
+        omit the commas separating them.
+
         Parameters
         ----------
         filename: str
@@ -281,11 +416,16 @@ def load_at2_file(cls, filename, scale=1.0):
             next(fp)
             description = next(fp).strip()
             next(fp)
-            parts = next(fp).split()
-            time_step = float(parts[1])
+            npts, time_step = _parse_at2_header(next(fp))
 
             accels = np.array([float(part) for line in fp for part in line.split()])
 
+        if accels.size != npts:
+            warnings.warn(
+                f"AT2 file '{filename}' specifies NPTS={npts}, but {accels.size} "
+                "accelerations were read."
+            )
+
         accels *= scale
         return cls(filename, description, time_step, accels)