function surfaceWave(h,k,c,T,L)
% SurfaceWave.m
% Written by Jay Borthen (Dec 2010)
% Georgetown University, Mathematics and Statistics

% h is the mesh step size in both the x and y directions
% k is the time step size
% c is the wave celerity (velocity)
% T is the total number of time steps
% L is the length and width of the square membrane
% Example: SurfaceWave(1,1,0.7,300)

Lx=L;  %bx is the length of the membrane in the x direction
Ly=L;  %by is the length of the membrane in the y direction

%Check CFL condition
if c*(k/h)<=(1/sqrt(2))
sigma=c*(k/h);
else
    disp('The approxmation will not be stable!')
    return
end

%A=zeros(L);
%image(A);
%axis xy
%grid on
%[ry, rx]=getpts;
ry=20;
rx=20;

u=zeros(Lx,Ly,T);
for j=1:T-1
    for ii=1:Ly
        for i=1:Lx
            
            %Implement x-dimension boundary condition
            if i==1 | i==Lx
                u(i,ii,j)=0;
            else

                %Implement y-dimension boundary condition
                if ii==1 | ii==Ly
                    u(i,ii,j)=0;
                else

                    %Implement initial conditions
                    if j==1
                    u(i,ii,j)=0;
                    ut(i,ii,j)=(1/60*pi)*exp(-(1/2)*((((1+i*h)-rx)^2)+(((1+(ii)*h)-ry)^2)));
                    
                    %Finite difference
                    u(i,ii,j+1)=(sigma^2)*(u(i+1,ii,j)-4*u(i,ii,j)+u(i-1,ii,j)+u(i,ii+1,j)+u(i,ii-1,j))+2*u(i,ii,j)-u(i,ii,j+1)+2*k*ut(i,ii,j);
                    
                    else

                    u(i,ii,j+1)=(sigma^2)*(u(i+1,ii,j)-4*u(i,ii,j)+u(i-1,ii,j)+u(i,ii+1,j)+u(i,ii-1,j))+2*u(i,ii,j)-u(i,ii,j-1);
                    
                    end
                end
            end
        end
    end
    subplot(2,1,1)
    surf(u(:,:,j))
    title(['Homogeneous Wave Equation at ', num2str(j), ' timesteps of length ', num2str(k)])
    view(-37.5,70)
    colormap('winter')
    axis([0 Lx 0 Ly -0.07 0.07])
    subplot(2,1,2)
    contour(u(:,:,j))
    axis([0 Lx 0 Ly])
    pause(0.05)
end
fprintf('The origin of the wave occured at (%1.1f, %1.1f).\n',rx,ry)