% Modes of the 2D Wave Equation in Polar Coordinates
% u_tt = k (r^2 u_rr+ ru_r + u_{theta theta})

% Initialize Data
%      a      = radius of membrane, 
%      T      = for Time Interval [0, T]
%      Nr, Nt = Number of Points in Space Grid, 
%      M      = Number of Time Steps
%      dr, dtheta = Delta r and Delta theta.
clear

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

% Root of Bessel function
Roots=besjroot(m,n,1e-3);
jmn=Roots(n);

% Spatial Grid
a=1;
Nr=30;
Ntheta=30;
dr=a/Nr;
dtheta=2*pi/Ntheta;
[r,theta]=meshgrid(0:dr:a,0:dtheta:2*pi);

% Times for one cylce
omega=c*jmn/a;
f=omega/2/pi;
T=2*pi/omega; 
M=50;
dt = T/M;

% Initialize u % Change to your initial condition
u=zeros(Nr+1,Ntheta+1);
u=besselj(m,jmn*r/a).*cos(m*theta); 
H=max(max(u));

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

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

movie(gcf, Mov,5,15)

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