refactored calculation of intensity measures into measures.py

eqsig/__about__.py

@@ -1 +1 @@
-__version__ = "0.5.26"
+__version__ = "0.5.27"

eqsig/__init__.py

@@ -2,4 +2,5 @@
 from eqsig import design_spectra
 from eqsig.functions import *
 from eqsig.multiple import Cluster, combine_at_angle, compute_rotated
-from eqsig.measures import significant_duration, calculate_peak
+from eqsig.measures import calc_significant_duration, calculate_peak  # deprecated load
+from eqsig import measures

eqsig/measures.py

@@ -1,9 +1,12 @@
 import numpy as np
+import scipy.integrate
+from eqsig import duhamels
 
 
-def significant_duration(motion, dt, start=0.05, end=0.95):
+def calc_significant_duration(motion, dt, start=0.05, end=0.95):
     """
     Computes the significant duration using cumulative acceleration according to Trifunac and Brady (1975).
+
     :param motion: acceleration time series
     :param dt: time step
     :param start: threshold to start the duration
@@ -22,4 +25,159 @@ def significant_duration(motion, dt, start=0.05, end=0.95):
 
 def calculate_peak(motion):
     """Calculates the peak absolute response"""
-    return max(abs(min(motion)), max(motion))
+    return max(abs(min(motion)), max(motion))
+
+
+def calc_peak(motion):
+    """Calculates the peak absolute response"""
+    return max(abs(min(motion)), max(motion))
+
+
+def calc_sir(acc_sig):
+    """
+    Calculates the shaking intensity rate
+
+    ref:
+    Dashti, S., Bray, J. D., Pestana, J. S., Riemer, M., and Wilson, D., 2010. Centrifuge testing to
+    evaluate and mitigate liquefaction-induced building settlement mechanisms,
+    ASCE Journal of Geotechnical and Geoenv. Eng. 136, 918–929
+
+    :param acc_sig:
+    :return:
+    """
+    ai_total = acc_sig.arias_intensity
+    t5, t75 = calc_significant_duration(acc_sig.values, acc_sig.dt)
+    return 0.7 * ai_total / (t75 - t5)
+
+
+def _raw_calc_arias_intensity(acc, dt):
+    acc2 = acc ** 2
+    return np.pi / (2 * 9.81) * scipy.integrate.cumtrapz(acc2, dx=dt, initial=0)
+
+
+def calc_arias_intensity(acc_sig):
+    return _raw_calc_arias_intensity(acc_sig.values, acc_sig.dt)
+
+
+
+def calc_cav(acc_sig):
+    """
+    Calculates the Cumulative Absolute velocity
+
+    ref:
+    Electrical Power Research Institute. Standardization of the Cumulative
+    Absolute Velocity. 1991. EPRI TR-100082-1'2, Palo Alto, California.
+    """
+    abs_acc = np.abs(acc_sig.values)
+    return scipy.integrate.cumtrapz(abs_acc, dx=acc_sig.dt, initial=0)
+
+
+def calc_cav_dp(asig):
+    """
+    Calculates standardized cumulative absolute velocity
+
+    ref:
+    Campbell KW, Bozorgnia Y. Predictive equations for the horizontal component of standardized
+    cumulative absolute velocity as adapted for use in the shutdown of U.S. nuclear power plants.
+    Nucl Eng Des 2011;241:2558–69.
+
+    :param asig:
+    :return:
+    """
+    start = 0
+    pga_max = 0
+    cav_dp = 0
+    points_per_sec = (int(1 / asig.dt))
+    total_seconds = int(asig.time[-1])
+    cav_dp_1_series = []
+    acc_in_g = asig.values / 9.81
+
+    for i in range(0, total_seconds):
+
+        end = start + points_per_sec
+
+        interval_total_time = (start * asig.dt) + 1
+        interval_time = np.arange(start * asig.dt, interval_total_time, asig.dt)
+
+        acc_interval = []
+        for j in range(start, end + 1):
+            acc_interval.append(acc_in_g[j])
+
+        acc_interval = np.array(acc_interval)
+        abs_acc_interval = abs(acc_interval)
+
+        x_lower = start * asig.dt  # the lower limit of x
+        x_upper = end * asig.dt  # the upper limit of x
+        x_int = interval_time[np.where((x_lower <= interval_time) * (interval_time <= x_upper))]
+        y_int = np.abs(np.array(abs_acc_interval)[np.where((x_lower <= interval_time) * (interval_time <= x_upper))])
+        int_acc = scipy.integrate.trapz(y_int, x_int)
+
+        # calculation of pga (g)
+        pga = (max(abs_acc_interval))
+        if pga > pga_max:
+            pga_max = pga
+        if (pga - 0.025) < 0:
+            h = 0
+        if (pga - 0.025) >= 0:
+            h = 1
+
+        cav_dp = cav_dp + (h * int_acc)
+        cav_dp_1_series.append(cav_dp)
+        start = end
+    t1s = np.arange(total_seconds)
+    cav_dp_time_series = np.interp(asig.time, t1s, cav_dp_1_series)
+    return cav_dp_time_series
+
+
+def calc_isv(acc_sig):
+    """
+    Calculates the integral of the squared velocity
+
+     See Kokusho (2013)
+    :return:
+    """
+    return scipy.integrate.cumtrapz(acc_sig.velocity ** 2, dx=acc_sig.dt, initial=0)
+
+
+def cumulative_response_spectra(acc_signal, fun_name, periods=None, xi=None):
+
+    if periods is None:
+        periods = acc_signal.response_times
+    else:
+        periods = np.array(periods)
+    if xi is None:
+        xi = 0.05
+    resp_u, resp_v, resp_a = duhamels.response_series(acc_signal.values, acc_signal.dt, periods, xi)
+    if fun_name == "arias_intensity":
+        rs = _raw_calc_arias_intensity(resp_a, acc_signal.dt)
+    else:
+        raise ValueError
+    return rs
+
+
+def calc_max_velocity_period(asig):
+    from eqsig import AccSignal
+    periods = np.logspace(-1, 0.3, 100)
+    new_sig = AccSignal(asig.values, asig.dt)
+    new_sig.generate_response_spectrum(response_times=periods, xi=0.15)
+
+    max_index = np.argmax(new_sig.s_v)
+    max_period = periods[max_index]
+    return max_period
+
+
+def max_acceleration_period(asig):
+    from eqsig import AccSignal
+    periods = np.logspace(-1, 1, 100)
+    new_sig = AccSignal(asig.values, asig.dt)
+    new_sig.generate_response_spectrum(response_times=periods, xi=0)
+
+    max_index = np.argmax(new_sig.s_a)
+    max_period = periods[max_index]
+    return max_period
+
+
+def max_fa_period(asig):
+    max_index = np.argmax(asig.fa_spectrum)
+    max_period = 1. / asig.fa_frequencies[max_index]
+    return max_period

eqsig/single.py

@@ -8,6 +8,7 @@
 import eqsig.duhamels as dh
 import eqsig.displacements as sd
 import eqsig.measures as sm
+from eqsig import measures
 
 
 class Signal(object):
@@ -526,7 +527,7 @@ def generate_duration_stats(self):
             self.a_rms10 = -1.
 
         # Trifunac and Brady
-        self.sd_start, self.sd_end = sm.significant_duration(self.values, self.dt)
+        self.sd_start, self.sd_end = sm.calc_significant_duration(self.values, self.dt)
 
         self.t_595 = self.sd_end - self.sd_start
 
@@ -540,11 +541,9 @@ def generate_cumulative_stats(self):
         Absolute Velocity. 1991. EPRI TR-100082-1'2, Palo Alto, California.
         """
         # Arias intensity in m/s
-        acc2 = self.values ** 2
-        self.arias_intensity_series = np.pi / (2 * 9.81) * scipy.integrate.cumtrapz(acc2, dx=self.dt, initial=0)
+        self.arias_intensity_series = measures.calc_arias_intensity(self)
         self.arias_intensity = self.arias_intensity_series[-1]
-        abs_acc = np.abs(self.values)
-        self.cav_series = scipy.integrate.cumtrapz(abs_acc, dx=self.dt, initial=0)
+        self.cav_series = measures.calc_cav(self)
         self.cav = self.cav_series[-1]
 
     @property
@@ -571,18 +570,9 @@ def isv_series(self):
         if "isv_series" in self._cached_params:
             return self._cached_params["isv_series"]
         else:
-            return self._calculate_isv_series()
-
-    def _calculate_isv_series(self):
-        """
-        Calculates the integral of the squared velocity
-
-         See Kokusho (2013)
-        :return:
-        """
-        isv_series = scipy.integrate.cumtrapz(self.velocity ** 2, dx=self.dt, initial=0)
-        self._cached_params["isv_series"] = isv_series
-        return isv_series
+            isv_series = measures.calc_isv(self)
+            self._cached_params["isv_series"] = isv_series
+            return isv_series
 
     def generate_all_motion_stats(self):
         """

tests/conftest.py

@@ -1,9 +1,21 @@
 import os
 import sys
 
+import numpy as np
+
+from eqsig import AccSignal
+
 PACKAGE_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), '..'))
 
 sys.path.append(PACKAGE_DIR)
 TEST_DIR = os.path.dirname(__file__)
 OPENSEES_DATA_DIR = os.path.join(TEST_DIR, 'test_data/opensees_data/')
-TEST_DATA_DIR = os.path.join(TEST_DIR, 'unit_test_data/')
+TEST_DATA_DIR = os.path.join(TEST_DIR, 'unit_test_data/')
+
+
+def t_asig():
+    record_path = TEST_DATA_DIR
+    record_filename = 'test_motion_dt0p01.txt'
+    motion_step = 0.01
+    rec = np.loadtxt(record_path + record_filename)
+    return AccSignal(rec, motion_step)

tests/test_measures.py

@@ -0,0 +1,14 @@
+import numpy as np
+
+from eqsig import measures
+
+from tests import conftest
+
+
+def test_cavdp():
+    asig = conftest.t_asig()
+
+    cav_dp = measures.calc_cav_dp(asig)[-1]
+
+    # 0.810412819 from eqsig 0.5.26  tested with several motions
+    assert np.isclose(cav_dp, 0.81041282, rtol=0.001), cav_dp