% Modes of the 2D Wave Equation in Cartesian Coordinates
% u_tt = k (u_xx + u_yy)

% Initialize Data
%      L, W   = Length and Width of membrane, 
%      T      = for Time Interval [0, T]
%      Nx, Ny = Number of Points in Space Grid, 
%      M      = Number of Time Steps
%      dx, dy = Delta x and Delta y.

clear

% Select Mode (n,m = 1,2,3, ...)
nx=3;
ny=2;
c=1;

% Spatial Grid
L=1;
W=2;
Nx=50;
Ny=50;
dx=L/Nx;
dy=W/Ny;
[x,y]=meshgrid(0:dx:L,0:dy:W);

% Time steps for one cycle
omega=c*sqrt((nx*pi/L)^2+(ny*pi/W)^2); 
f=omega/2/pi;
T=2*pi/omega;
M=50;
dt = T/M;

% Initialize u % Change to your initial condition
u=zeros(Nx+1,Ny+1);
u=u+sin(nx*pi*x/L).*sin(ny*pi*y/W); 
H=max(max(u));

% Plot initial condition
surf(x,y,u,'FaceColor','red','EdgeColor','none')
camlight left;
lighting phong
xlabel('x')
ylabel('y')
title(['Frequency = ' num2str(f) 'Hz. Solution at t = 0'])
axis([0,L,0,W,-H,H])
frame=1;
Mov(frame)=getframe(gcf); 
%pause(0.5) 

% Time evolution
for j=1:M
    t=j*dt;  
    u=sin(nx*pi*x/L).*sin(ny*pi*y/W)*cos(omega*t);   
   % Plot 3D solution every 10 time steps
   %if mod(m,10)==0
       frame=frame+1;
       surf(x,y,u,'FaceColor','red','EdgeColor','none')
       camlight left;
       lighting phong
       xlabel('x')
       ylabel('y')
       title(['Frequency = ' num2str(f) 'Hz. Solution at t = ' num2str(t)])
       axis([0,L,0,W,-H,H])
       Mov(frame)=getframe(gcf); 
       %pause (0.5)
   %end
end

movie(gcf, Mov,5,15)
%movie(gcf, Mov,5)

% Extra - show or create a movie
%    movie(gcf, Mov)         % plays movie
%    movie(gcf, Mov,1,2)     % plays at 2 fps 
%    movie(gcf, Mov,10,5)    % repeats 10 times at 5 fps
% Create movie 
%    movie2avi(Mov, 'heat2d.avi', 'compression', 'None');