% Second-order leapfrog for the 2D acoustic wave equation.
%
%   p_tt + 2*sig*p_t = c^2 * lap(p) + s(x, y, t)
%
% Pressure only, at constant density, so the medium is the sound speed c
% (metres per second) and the absorption sig (inverse seconds) that the scene
% defines. Everything here is SI: x, y in metres, t and dt in seconds, f in
% hertz — none of it is declared as such anywhere, it simply follows from x,
% y and t being metres and seconds, which is the scene's and the app's doing,
% not this file's. Centring both the second time
% derivative and the damping term on step n,
%
%   (p+ - 2p + p-)/dt^2 + sig*(p+ - p-)/dt = c^2*lap(p) + s
%
% and solving for p+ gives the update below. It is explicit: the only thing
% that couples neighbouring points is `lap2`, the 5-point Laplacian, which the
% host supplies as a single GPU dispatch. Everything else here is element-wise
% and compiles to one kernel per line.
%
% Stability wants c*dt/h <= 1/sqrt(2); the app sets dt from the fastest speed
% anywhere in the scene, so a fast scatterer slows the whole run down.
%
% Time is carried as a field rather than a number, because a batch of steps is
% one replay of a fixed sequence of GPU operations and nothing the host writes
% between frames can change inside it. `tn = t + dt` makes the model keep its
% own clock, which is what lets the source term be correct however many steps
% are batched into a submit.

function [p, pm, t] = init(npts)
  p  = zeros(npts, 1);
  pm = zeros(npts, 1);
  t  = zeros(npts, 1);
end

function [pn, pold, tn] = step(p, pm, t, x, y, c, sig, dt, f, t0, tw, cw, x0, y0, w, point)
  % The source. `cw` blends between a Gaussian pulse (0) and a wave that
  % turns on smoothly and stays on (1); `point` blends between a line source
  % spanning the grid in y, whose far field is a plane wave, and a point
  % source at (x0, y0).
  %
  % The om^2 is only a choice of units, not a physical amplitude — this is a
  % body force of arbitrary strength, not a source with a real acoustic power
  % rating. Such a force drives a response that falls off as 1/om^2 — two
  % time integrations — so without the factor the field would shrink tenfold
  % every time the frequency slider tripled. The equation is linear; scaling
  % the source scales the answer and nothing else.
  u = (t - t0) / tw;
  env = (1 - cw) * exp(-u .* u) + cw * (0.5 * (1 + tanh(u)));
  gx = ((x - x0) / w) .^ 2;
  gy = point * (((y - y0) / w) .^ 2);
  om = 2*pi*f;
  s = (om*om) * (env .* sin(om*(t - t0)) .* exp(-(gx + gy)));

  % One step. The damping is what the absorbing layer acts through: sig is
  % zero over the interior, so there p+ is the plain leapfrog update.
  sd = sig * dt;
  lap = lap2(p);
  pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);

  % This step's field becomes the next step's history. A line on its own, so
  % it plans as a copy rather than being folded into the update above.
  pold = p;
  tn = t + dt;
end