Skip to contents

Elastic SDOF response spectra: TSL2PS

spectra.Rmd
# Minimal executable example
library(gmsp)
library(data.table)
#> 
#> Attaching package: 'data.table'
#> The following object is masked from 'package:base':
#> 
#>     %notin%
t_vec <- seq(0, 20, by = 0.01)
dt_acc <- data.table(
  t = t_vec,
  H1 = 500 * sin(2 * pi * 1 * t_vec) * exp(-0.1 * t_vec),
  H2 = 300 * cos(2 * pi * 1.5 * t_vec) * exp(-0.1 * t_vec),
  UP = 100 * sin(2 * pi * 0.8 * t_vec) * exp(-0.1 * t_vec)
)
tsl <- AT2TS(dt_acc, units.source = "mm", isRaw = FALSE,
             output = "TSL", audit = FALSE)
ps <- TSL2PS(tsl[OCID == "H1"], xi = 0.05,
             Tn = c(0.1, 0.2, 0.5, 1.0, 2.0),
             output = "PSL")
head(ps)
#>      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

TSL2PS() computes elastic single-degree-of-freedom (SDOF) response spectra for canonical TSL input. Source time-series ID values map to spectral ID values: AT -> PSA, VT -> PSV, and DT -> SD.

The reported spectral IDs are:

  • PSA — pseudo-spectral acceleration,
  • PSV — pseudo-spectral velocity,
  • SD — spectral displacement.

Classic references: Nigam & Jennings (1969), Chopra (2017).

Input and output

Input

  • .x: canonical TSL data.table with t, s, ID, OCID, plus optional metadata.
  • xi: damping ratio (0≤ξ≤10 \le \xi \le 1), default 0.05.
  • Tn: period vector in seconds. If NULL, a default grid is used. Do not include 0; the function prepends the Tn = 0 peak-value anchor internally.

Grouping metadata is derived from the TSL schema: all columns except t, s, ID, and OCID are metadata keys. OCID remains the component/channel key.

Output

output = "PSL" returns the canonical long table. output = "PSW" returns the wide projection:

  • PSW: wide projection with columns PSA.<OCID>, PSV.<OCID>, SD.<OCID>.
  • PSL: canonical long table with columns Tn, ID ∈ {PSA, PSV, SD}, OCID, S, plus grouping keys.

TSL2PS() is the public spectra interface. It consumes canonical TSL; it does not expose BY, COL.s, COL.t, or COL.ID.

PSL2PSW() and PSW2PSL() expose the same long/wide projection for existing spectra tables:

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

When xi is a vector, TSL2PS() runs the scalar spectra path once per damping ratio and adds xi as metadata. Scalar xi keeps the historical schema without an xi column.

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 TnT_n the code computes

ωn=2πTn,C=2ξωn,K=ωn2.\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:

𝐲̇=𝐀𝐲+𝐁u(t),𝐀=[01−K−C],𝐁=[01].\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)u(t) is the series value s[k]. The sign convention (e.g., u(t)=±ag(t)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 Δt\Delta t, assuming piecewise-constant input within the step, the exact update is

𝐲k=e𝐀Δt𝐲k−1+(e𝐀Δt−𝐈)𝐀−1𝐁uk.\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…Nk = 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=maxt|y1(t)|.SD = \max_t |y_1(t)|.

It then reports

PSV=ωnSD,PSA=ωn2SD.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 Tn=0T_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.