function [u, v, acc] = plate_sdof(omega, sigma, p, dt)
% PLATE_SDOF  Exact response of a damped oscillator to piecewise-linear forcing.
%
%   [u, v, acc] = plate_sdof(omega, sigma, p, dt) integrates
%
%       u'' + 2 sigma u' + omega^2 u = p(t)
%
%   from rest, where p is sampled every dt seconds and varies linearly between
%   samples.  Returns displacement u, velocity v and acceleration acc at the
%   sample times.  The update is the exact solution over each interval
%   (Nigam and Jennings, 1969; Chopra, Dynamics of Structures, sec. 5.2), so
%   there is no time-step stability limit and no phase error: the only
%   approximation is the linear interpolation of p between samples.
%
%   Requires sigma < omega (underdamped).

if sigma >= omega
    error('plate_sdof: the oscillator must be underdamped (sigma = %g >= omega = %g)', sigma, omega);
end
wd = sqrt(omega^2 - sigma^2);          % damped natural frequency
e = exp(-sigma * dt);
s = sin(wd * dt);
c = cos(wd * dt);
w2 = omega^2;
zr = sigma / wd;                        % zeta / sqrt(1 - zeta^2)

% Coefficients of  u(n+1) = A u(n) + B v(n) + C p(n) + D p(n+1)
%                  v(n+1) = A2 u(n) + B2 v(n) + C2 p(n) + D2 p(n+1)
A = e * (zr * s + c);
B = e * s / wd;
C = (2 * sigma / (w2 * dt) + e * (((w2 - 2 * sigma^2) / (w2 * wd * dt) - zr) * s - (1 + 2 * sigma / (w2 * dt)) * c)) / w2;
D = (1 - 2 * sigma / (w2 * dt) + e * ((2 * sigma^2 - w2) / (w2 * wd * dt) * s + 2 * sigma / (w2 * dt) * c)) / w2;
A2 = -e * w2 / wd * s;
B2 = e * (c - zr * s);
C2 = (-1 / dt + e * ((w2 / wd + sigma / (dt * wd)) * s + c / dt)) / w2;
D2 = (1 - e * (zr * s + c)) / (w2 * dt);

n = numel(p);
u = zeros(size(p));
v = zeros(size(p));
for k = 1:n - 1
    u(k + 1) = A * u(k) + B * v(k) + C * p(k) + D * p(k + 1);
    v(k + 1) = A2 * u(k) + B2 * v(k) + C2 * p(k) + D2 * p(k + 1);
end
acc = p - 2 * sigma * v - w2 * u;
end