arkottke / pyrotd Goto Github PK

import numpy as np
import matplotlib.pyplot as plt
import pyrotd

#%% system of units
m = 1.0
s = 1.0
kg = 1.0
N = kg*m/s/s

g = 9.80665*m/s/s

#%% read record
record = 'RSN8883_14383980_13849360.AT2'

record_content = open(record, mode='r').readlines()

# Flag indicating dt is found and that ground motion
# values should be read -- ASSUMES dt is on last line
# of header!!!
flag = 0

record_samples = []
for line in record_content:
    if line == '\n':
        # Blank line --> do nothing
        continue
    elif flag == 1:
        words = line.split()
        for word in words:
            record_samples.append(float(word))
    else:
        words = line.split()
        if words[0] == 'NPTS=':
            record_NPTS = words[1]
            record_DT = float(words[3])
            flag = 1

#%% compute spectral parameteres
f = np.logspace(-1, 2, 91)
T = 1/f
xi = 0.05
pyrotd_pSa = pyrotd.calc_spec_accels(record_DT, record_samples, f, xi, max_freq_ratio=5, osc_type='psa')

pyrotd_Sd = pyrotd.calc_spec_accels(record_DT, record_samples, f, xi, max_freq_ratio=5, osc_type='sd')

formulated_Sd = []
for index, psa in enumerate(pyrotd_pSa.spec_accel):
    formulated_Sd.append( psa*g * (T[index]**2) / (2 / np.pi)**2 )
formulated_Sd = np.asarray(formulated_Sd)

#%% plot
fig, axs= plt.subplots(layout='constrained')
fig.suptitle('Displacement spectra of ' + record)
axs.plot(T, pyrotd_Sd.spec_accel, label='Sd')
axs.plot(T, formulated_Sd, label=r'Sd = $pSa \cdot g \cdot \frac{T^2}{4\pi^2}$')
axs.set_xlabel('Period, s')
axs.set_ylabel('Spectral displacement')
axs.legend()
plt.show()

#%% plot
fig, axs= plt.subplots(layout='constrained')
fig.suptitle('Displacement spectra of ' + record)
axs.plot(T, pyrotd_Sd.spec_accel, label='$Sd$')
axs.plot(T, pyrotd_Sd.spec_accel*(g**3), label='$Sd \cdot g^3$')
axs.plot(T, formulated_Sd, label=r'$Sd = pSa \cdot g \cdot \frac{T^2}{4\pi^2}$')
axs.set_xlabel('Period, s')
axs.set_ylabel('Spectral displacement')
axs.legend()
plt.show()

import numpy as np
import pyrotd

#%% read record
record = 'RSN8883_14383980_13849360.AT2'

record_content = open(record, mode='r').readlines()

# Flag indicating dt is found and that ground motion
# values should be read -- ASSUMES dt is on last line
# of header!!!
flag = 0

record_samples = []
for line in record_content:
    if line == '\n':
        # Blank line --> do nothing
        continue
    elif flag == 1:
        words = line.split()
        for word in words:
            record_samples.append(float(word))
    else:
        words = line.split()
        if words[0] == 'NPTS=':
            record_NPTS = words[1]
            record_DT = float(words[3])
            flag = 1

#%% compute spectral acceleration
f = np.logspace(-1, 2, 91)
xi = 0.05
spec_accel = pyrotd.calc_spec_accels(record_DT, record_samples, f, xi, 5, 'sa').spec_accel