% This is an attempt to generate a directed wave pulse.

clear
% Generate initial conditions
Nx = 32; % x-size of mesh
Ny = 128; % y-size of mesh
% Initial pulse
x = (1:Nx)' * ones(1,Ny) - Nx/2;
y = ones(Nx,1) * (1:Ny) - Ny/2;
psi0 = exp(-(x.^2 + 0.1*y.^2)/2);
% Time derivative is based on translation at c=1.
psidot0 = x .* psi0;

% Set up the time axis
tfinal = 2*Nx;
Nt = 2*tfinal
t = linspace(0, tfinal, Nt);

% Solve the problem
t0=cputime;
psi = waves2d(psi0, psidot0, t);
calculation_time_seconds = cputime - t0

% Graphics
unix('rm /tmp/waves*.png');   % delete existing movie frames from /tmp.
minpsi = min(min(min(psi)));  % for scaling the plot.
maxpsi = max(max(max(psi)));  % for scaling the plot.

colormap(hot(256));
t0=cputime;
for i=1:Nt
   image(255*(psi(:,:,i) - minpsi)/(maxpsi - minpsi)) % psi scaled 0-255
   axis('equal')
   filename = sprintf('/tmp/waves_image%03i.png',i);
   print(filename);
end
image_time_seconds = cputime - t0

colormap(cool(256));
t0=cputime;
for i=1:Nt
   mesh(psi(:,:,i))
   axis([1,max(Nx,Ny),1,max(Nx,Ny),-maxpsi,maxpsi])
   caxis([minpsi,maxpsi])
   shading('faceted')
   filename = sprintf('/tmp/waves_mesh%03i.png',i);
   print(filename);
end
mesh_time_seconds = cputime - t0
