Elastic SDOF response spectra: TSL2PS

Input and output

psw <- PSL2PSW(ps)
psl_again <- PSW2PSL(psw)
head(psl_again)
#>      OCID    Tn     ID         S
#>    <char> <num> <char>     <num>
#> 1:     H1   0.0    PSA  487.7138
#> 2:     H1   0.1    PSA  492.8057
#> 3:     H1   0.2    PSA  571.3886
#> 4:     H1   0.5    PSA 1021.4305
#> 5:     H1   1.0    PSA 3026.4257
#> 6:     H1   2.0    PSA  316.6346

ps_xi <- TSL2PS(tsl[OCID == "H1"], xi = c(0.02, 0.05),
                Tn = c(0.1, 0.2), output = "PSW")
head(ps_xi)
#>       xi    Tn   PSA.H1     PSV.H1       SD.H1
#>    <num> <num>    <num>      <num>       <num>
#> 1:  0.02   0.0 487.7138 56.4253259 6.531096805
#> 2:  0.02   0.1 493.6305  0.9069311 0.001670949
#> 3:  0.02   0.2 579.3027  1.8773197 0.006969010
#> 4:  0.05   0.0 487.7138 56.4253259 6.531096805
#> 5:  0.05   0.1 492.8057  0.9069939 0.001671061
#> 6:  0.05   0.2 571.3886  1.8754513 0.006929751

SDOF model

For each period \(T_n\) the code computes

\[\omega_n = \frac{2\pi}{T_n}, \qquad C = 2\,\xi\,\omega_n, \qquad K = \omega_n^2.\]

A 2D linear state-space system is integrated:

\[\dot{\mathbf{y}} = \mathbf{A}\,\mathbf{y} + \mathbf{B}\,u(t), \qquad \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -K & -C \end{bmatrix}, \qquad \mathbf{B} = \begin{bmatrix} 0 \\ 1 \end{bmatrix}.\]

The input \(u(t)\) is the series value s[k]. The sign convention (e.g., \(u(t) = \pm a_g(t)\)) is not enforced by the code; since the final outputs use absolute maxima, the sign does not affect PSA / PSV / SD.

Discretisation via matrix exponential

For a time step \(\Delta t\), assuming piecewise-constant input within the step, the exact update is

\[\mathbf{y}_k = e^{\mathbf{A}\Delta t}\,\mathbf{y}_{k-1} + \left(e^{\mathbf{A}\Delta t} - \mathbf{I}\right) \mathbf{A}^{-1}\,\mathbf{B}\,u_k.\]

Implementation:

• Ae = expm::expm(A * dt)
• AeB = (Ae - I) %*% pracma::pinv(A) %*% B

The state is then updated by looping over \(k = 2 \ldots N\).

How PSA / PSV / SD are built

In the single-signal SDOF worker, after simulation the code takes the absolute maximum of displacement (first state):

\[SD = \max_t |y_1(t)|.\]

It then reports

\[PSV = \omega_n\,SD, \qquad PSA = \omega_n^2\,SD.\]

For canonical TSL input, TSL2PS() runs the worker on the matching source series and keeps the matching spectral ID: PSA from AT, PSV from VT, and SD from DT. The function prepends \(T_n = 0\) with the corresponding peak value (PGA, PGV, or PGD).

D50 and D100 horizontal spectra

For canonical TSL input, TSL2PS(D50 = TRUE) can add the median rotated horizontal component as an ordinary derived OCID named D50. TSL2PS(D100 = TRUE) can add the maximum rotated horizontal component as OCID = "D100". The input must contain OCID = "H1" and OCID = "H2" for each source ID (AT, VT, and DT). The vertical component UP remains an ordinary component and is not used to compute D50 or D100.

D50 and D100 follow the same spectral-ID contract as the component spectra:

• PSA.D50 is computed from rotated AT;
• PSV.D50 is computed from rotated VT;
• SD.D50 is computed from rotated DT.
• PSA.D100, PSV.D100, and SD.D100 use the same source-ID mapping.

D100 is computed independently for each period and spectral ID. The angle that maximizes PSA can differ from the angle that maximizes PSV or SD.

t_d50 <- seq(0, 1, by = 0.01)
at_wide <- data.table(
  t = t_d50,
  H1 = 10 * sin(2 * pi * 3 * t_d50),
  H2 = 7 * cos(2 * pi * 4 * t_d50),
  UP = 3 * sin(2 * pi * 2 * t_d50)
)

tsl <- AT2TS(at_wide, units.source = "mm", isRaw = FALSE,
             output = "TSL", audit = FALSE)
rot_psl <- TSL2PS(tsl, Tn = c(0.1, 0.2),
                  output = "PSL", D50 = TRUE, D100 = TRUE, nTheta = 12L)
rot_psl[OCID %chin% c("D50", "D100")]
#>       OCID    Tn     ID            S
#>     <char> <num> <char>        <num>
#>  1:   D100   0.0    PSA 1.215984e+01
#>  2:   D100   0.1    PSA 1.632737e+01
#>  3:   D100   0.2    PSA 3.143589e+01
#>  4:   D100   0.0    PSV 6.424870e-01
#>  5:   D100   0.1    PSV 1.139222e-02
#>  6:   D100   0.2    PSV 4.342486e-02
#>  7:   D100   0.0     SD 3.509437e-02
#>  8:   D100   0.1     SD 9.787487e-06
#>  9:   D100   0.2     SD 7.093803e-05
#> 10:    D50   0.0    PSA 1.126617e+01
#> 11:    D50   0.1    PSA 1.412505e+01
#> 12:    D50   0.2    PSA 2.759419e+01
#> 13:    D50   0.0    PSV 5.949836e-01
#> 14:    D50   0.1    PSV 1.060566e-02
#> 15:    D50   0.2    PSV 3.957369e-02
#> 16:    D50   0.0     SD 3.133454e-02
#> 17:    D50   0.1     SD 8.812128e-06
#> 18:    D50   0.2     SD 5.617509e-05

Notes (limitations and corner cases)

• If the user provides Tn that includes 0, the call errors because TSL2PS() adds the Tn = 0 peak-value anchor internally.
• Complexity is nested over period and sample: long records and dense period grids can be slow.
• Time regularization is handled inside the worker; the series is interpolated to a rationalised dt.
• TSL2PS() derives grouping metadata from the TSL schema. It does not expose BY, COL.s, COL.t, or COL.ID.

References

• Nigam, N. C., & Jennings, P. C. (1969). Calculation of Response Spectra from Strong-Motion Earthquake Records. Bulletin of the Seismological Society of America, 59(2), 909–922.
• Chopra, A. K. (2017). Dynamics of Structures: Theory and Applications to Earthquake Engineering (5th ed.). Pearson.