Elastic SDOF response spectra: TSL2PS • gmsp

# 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 — 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)
#>      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.42533 6.531097
#> 2:  0.02   0.1 493.6305 56.98416 6.596640
#> 3:  0.02   0.2 579.3027 58.97774 6.878137
#> 4:  0.05   0.0 487.7138 56.42533 6.531097
#> 5:  0.05   0.1 492.8057 56.98811 6.597085
#> 6:  0.05   0.2 571.3886 58.91904 6.839390

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")]
#>       OCID    Tn     ID           S
#>     <char> <num> <char>       <num>
#>  1:   D100   0.0    PSA 12.15984138
#>  2:   D100   0.1    PSA 16.32736580
#>  3:   D100   0.2    PSA 31.43589259
#>  4:   D100   0.0    PSV  0.64248700
#>  5:   D100   0.1    PSV  0.71579447
#>  6:   D100   0.2    PSV  1.36423235
#>  7:   D100   0.0     SD  0.03509437
#>  8:   D100   0.1     SD  0.03863945
#>  9:   D100   0.2     SD  0.07001303
#> 10:    D50   0.0    PSA 11.26617109
#> 11:    D50   0.1    PSA 14.12505433
#> 12:    D50   0.2    PSA 27.59419196
#> 13:    D50   0.0    PSV  0.59498360
#> 14:    D50   0.1    PSV  0.66637326
#> 15:    D50   0.2    PSV  1.24324408
#> 16:    D50   0.0     SD  0.03133454
#> 17:    D50   0.1     SD  0.03478889
#> 18:    D50   0.2     SD  0.05544259

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.