Elastic SDOF response spectra: TSL2PS
In gmsp: Ground Motion Signal Processing

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)

# Minimal executable example
library(gmsp)
library(data.table)
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)

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 — domain-normalised acceleration response from AT,
PSV — domain-normalised velocity response from VT,
SD — domain-normalised displacement response from DT.

These are three independent domain operators. They are not three classical pseudo-ordinates derived from one oscillator driven only by acceleration.

Input and output

Input

.x: canonical TSL data.table with t, s, ID, OCID, plus optional metadata.
xi: damping ratio ($0 \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)

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)

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

Let $d \in {A,V,D}$ identify the acceleration, velocity, or displacement domain and let $x_d(t)$ be the matching AT, VT, or DT input. The worker solves an independent state equation in each domain,

$$ \ddot z_d + 2\xi\omega_n\dot z_d + \omega_n^2 z_d = x_d(t), $$

whose transfer function is

$$ H_d(p) = \frac{1}{p^2 + 2\xi\omega_n p + \omega_n^2}. $$

gmsp uses the domain-normalised operator

$$ G_d(p) = \omega_n^2 H_d(p), \qquad S_d(T_n) = \omega_n^2 \max_t |z_d(t)|. $$

The multiplier is therefore $\omega_n^2$ in all three branches. If $x_d$ has units $U_d$, the first state has units $U_d\,\mathrm{s}^2$ and $S_d$ retains $U_d$. Consequently PSA has acceleration units, PSV has velocity units, and SD has displacement units. Because $z_A$, $z_V$, and $z_D$ are responses to different source histories, identities such as $PSV = \omega_n SD$ do not apply across these columns.

The function prepends $T_n = 0$ with the corresponding unfiltered peak input value: PGA from AT, PGV from VT, and PGD from DT.

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. The $\omega_n^2$ domain normalisation is applied to every rotated response before the D50 quantile or D100 maximum is evaluated.

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")]

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.