% Wave packet on an oblong domain. The wave speed is the same
% in the horizontal and vertical directions. The key to making this
% work properly is k2d, which defines the wavenumbers in such a way
% that the maximum kx and ky are equal.

clear
% Generate initial conditions
Nx = 32; % x-size of mesh
Ny = 64; % y-size of mesh
x = (1:Nx)' * ones(1,Ny) - Nx/2;
y = ones(Nx,1) * (1:Ny) - Ny/2;
% Wavepacket traveling along y:
psi0 = sin(y/3) .* exp(-y.^2/100);
% Time derivative is based on translation at c=1.
psidot0 = (1/3)*cos(y/3) .* exp(-y.^2/100) - (1/50) .* y .* psi0;
% pulse traveling along x:
psi1 = 0.1*exp(-x.^2/2);
% Time derivative is based on translation at c=1.
psidot1 = x .* psi1;
psi0 = psi0 + 2*psi1;
psidot0 = psidot0 + 2*psidot1;

% Set up the time axis
tfinal = Ny;
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
