Added responsespectra source code files

src/functions.py

@@ -0,0 +1,87 @@
+# -*- coding: utf-8 -*-
+"""
+- Function to compute the response spectra of time series using the Duhamel integral technique.
+
+Inputs:
+   - data: acceleration data in the time domain
+   - delta: Sampling rate of the time-series (in Hz)
+   - T: Output period range in second, Example (if delta>=20 Hz): T = np.concatenate((np.arange(.1, 1, .01), np.arange(1, 20, .1)))
+   - xi: Damping factor (Standard: 5% -> 0.05)
+   - Resp_type: Response type, choose between:
+                - 'SA'  : Acceleration Spectra
+                - 'PSA' : Pseudo-acceleration Spectra
+                - 'SV'  : Velocity Spectra
+                - 'PSV' : Pseudo-velocity Spectra
+                - 'SD'  : Displacement Spectra
+
+ Output:
+     - Response spectra in the unit specified by 'Resp_type'
+"""
+
+#  -----------------
+#  Copyright © 2023 ACK Cyfronet AGH, Poland.
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+#
+#  This work was partially funded by DT-GEO Project.
+#  -----------------
+
+import numpy as np
+from numba import jit
+
+
+@jit(nopython=True)
+def rs_func(data, dt, T, xi, Resp_type):
+
+    n = len(data)    
+    n_T = len(T)
+    
+    # dt = 1 / delta
+    w = 2 * np.pi / T
+
+    mass = 1  #  constant mass (=1)
+    c = 2 * xi * w * mass
+    wd = w * np.sqrt(1 - xi**2)
+    p1 = -mass * data   
+    
+    # predefine output matrices
+    S = np.zeros(n_T)
+    for j in np.arange(n_T):
+        # Duhamel time domain matrix form
+        I0 = 1/w[j]**2*(1-np.exp(-xi*w[j]*dt)*(xi*w[j]/wd[j]*np.sin(wd[j]*dt)+np.cos(wd[j]*dt)))
+        J0 = 1/w[j]**2*(xi*w[j]+np.exp(-xi*w[j]*dt)*(-xi*w[j]*np.cos(wd[j]*dt)+wd[j]*np.sin(wd[j]*dt)))
+        
+        AA = [[np.exp(-xi*w[j]*dt)*(np.cos(wd[j]*dt)+xi*w[j]/wd[j]*np.sin(wd[j]*dt)) , np.exp(-xi*w[j]*dt)*np.sin(wd[j]*dt)/wd[j] ] , 
+               [-w[j]**2*np.exp(-xi*w[j]*dt)*np.sin(wd[j]*dt)/wd[j] ,np.exp(-xi*w[j]*dt)*(np.cos(wd[j]*dt)-xi*w[j]/wd[j]*np.sin(wd[j]*dt)) ]]
+        BB = [[I0*(1+xi/w[j]/dt)+J0/w[j]**2/dt-1/w[j]**2 , -xi/w[j]/dt*I0-J0/w[j]**2/dt+1/w[j]**2 ] ,
+            [J0-(xi*w[j]+1/dt)*I0, I0/dt] ]
+
+        u1 = np.zeros(n)
+        udre1 = np.zeros(n)
+        
+        for xx in range(1,n,1) :
+    
+            u1[xx] = AA[0][0]*u1[xx-1] + AA[0][1]*udre1[xx-1] + BB[0][0]*p1[xx-1] + BB[0][1]*p1[xx]
+            udre1[xx] = AA[1][0]*u1[xx-1]+AA[1][1]*udre1[xx-1] + BB[1][0]*p1[xx-1]+BB[1][1]*p1[xx]
+
+        if Resp_type == "SA":
+            S[j] = np.max(np.abs(-(w[j] ** 2 * u1 + c[j] * udre1)))
+        elif Resp_type == "PSA":
+            S[j] = np.max(np.abs(u1)) * w[j] ** 2
+        elif Resp_type == "SV":
+            S[j] = np.max(np.abs(udre1))
+        elif Resp_type == "PSV":
+            S[j] = np.max(np.abs(u1)) * w[j]
+        elif Resp_type == "SD":
+            S[j] = np.max(np.abs(u1))
+    return S

src/responsespectra.py

@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+
+#  -----------------
+#  Copyright © 2023 ACK Cyfronet AGH, Poland.
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+#
+#  This work was partially funded by DT-GEO Project.
+#  -----------------
+
+import numpy as np
+from obspy import read
+from functions import rs_func
+
+
+def main(miniseed_file, Damping = 0.05):
+    """
+    Python application that reads the acceleration miniseed file and calculates the response spectra.
+    Arguments:
+        miniseed_file: path to input file of type 'MiniSEED in accleration unit'
+        Damping: Damping factor. 0.05 means 5% damping ratio
+    Returns:
+        CSV files of the response spectra values.
+    """
+    Resp_types = ['SA','PSA','SV','PSV','SD']
+    T = np.concatenate(
+        (
+            np.arange(0.05, 0.1, 0.005), #0.05s to 0.1s by 0.005s intervals
+            np.arange(0.1, 0.5, 0.01), #0.1s to 0.5s by 0.01s intervals
+            np.arange(0.5, 1, 0.02),
+            np.arange(1, 5, 0.1),
+            np.arange(5, 10.5, 0.5),
+        )
+    )
+    stream_raw = read(miniseed_file)
+    for resp in Resp_types:
+        stream = stream_raw.copy() #start with fresh stream
+        csv_header = ["Period"] #start of header for csv output file
+        Sfin = np.zeros((len(stream), len(T)))
+        for i, st in enumerate(stream):
+            st.data = st.data * 100 #convert from m/s/s to cm/s/s
+            dt = round(st.stats.delta, 3)  # Get dt sampling rate in seconds
+            tr_data = st.data # extract time series
+            Sfin[i] = rs_func(tr_data, dt, T, Damping, Resp_type=resp)
+            net = st.stats.network
+            stn = st.stats.station
+            loc = st.stats.location
+            chan = st.stats.channel
+            nslc = net + '.' + stn + '.' + loc + '.' + chan
+            csv_header.append(nslc + " " + resp)
+        csv_array = np.vstack((T,Sfin))
+        csv_array = csv_array.T          
+        np.savetxt('spectra_'+resp+'.csv', csv_array, delimiter=",", header=",".join(csv_header))

src/responsespectra_wrapper.py

@@ -0,0 +1,38 @@
+# -*- coding: utf-8 -*-
+
+#  -----------------
+#  Copyright © 2023 ACK Cyfronet AGH, Poland.
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+#
+#  This work was partially funded by DT-GEO Project.
+#  -----------------
+
+import sys
+import argparse
+from responsespectra import main as responsespectra
+
+
+def main(argv):
+    parser = argparse.ArgumentParser()
+    parser.add_argument("miniseed_file", help="Path to input file of type 'MiniSEED in accleration unit'")
+    parser.add_argument("--Damping", help="Damping factor. 0.05 means 5% damping ratio",
+                        type=float, default=0.05, required=False)
+
+    args = parser.parse_args()
+    responsespectra(**vars(args))
+    return
+
+
+if __name__ == '__main__':
+    main(sys.argv)