% Leapfrog with the fourth-order Laplacian, otherwise identical to
% models/leapfrog.m. The only change is `lap4` in place of `lap2`.
%
% What it buys is grid dispersion. A discrete Laplacian gets the wave speed
% slightly wrong at wavelengths that the grid barely resolves, and the error
% grows with distance travelled: a pulse made of many wavelengths slowly comes
% apart, its short components lagging behind, leaving a trailing ripple that
% is entirely numerical. The 5-point stencil's phase error goes as (k*h)^2,
% the 9-point one's as (k*h)^4, so at the eight or ten points per wavelength
% where a run is cheap enough to be interactive, the difference is easy to
% see: raise the frequency until models/leapfrog.m visibly trails, then switch
% to this one.
%
% It is not free. The stencil reads twice as many neighbours, and its
% stability limit is tighter (c*dt/h <= sqrt(3/8) rather than 1/sqrt(2)), so
% the app also takes a slightly smaller timestep. Roughly: a little more work
% per step, in exchange for needing a much coarser grid at the same accuracy.

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)
  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)));

  sd = sig * dt;
  lap = lap4(p);
  pn = (2*p - (1 - sd) .* pm + (c*dt).^2 .* lap + (dt*dt) * s) ./ (1 + sd);

  pold = p;
  tn = t + dt;
end