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,400,100)

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;
% close(gcf)

ry=0.45*L;
rx=0.62*L;

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
                    
                    %if i<Lx/2+4 && i>Lx/2-4 && ii<Ly/2+4 && ii>Ly/2-4
                    %    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)=exp(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)+3.1*k*ut(i,ii,j);
                            
                        else
                            
                            %Undamped: 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)); 
                            u(i,ii,j+1)= exp(-j/(25*T)).*((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
        
    end
    
    surf(-u(:,:,j))
    %text(3,12,-0.09,'Unity Analytics LLC','Color',[0 0 0],'FontWeight','bold')
	%hold on
	%text(1,12,-0.12,'(www.unityanalytics.net)','Color',[0 0 0])
    hold on
    colormap('winter')
    axis([0 L 0 L -0.15 0.15])
    %Z=0.02.*ones(L,L);
    %Z(:,L/2+3)=-0.04;
    %Z(:,L/2-3)=-0.04;
    %Z(L/2-3,:)=-0.04;
    %Z(L/2+3,:)=-0.04;
    %mesh((L/2-3):(L/2+3),(L/2-3):(L/2+3),Z((L/2-3):(L/2+3),(L/2-3):(L/2+3)))
   

    hold off
    set(gcf,'Renderer','zbuffer')

    set(gcf,'Color',[1 1 1]);
%figure('Color',[0.8 0.8 0.8]);
if j>50
im=frame2im(getframe);
imB = imresize(im,[144 210]);
imwrite(imB,strcat('C:/Users/Jay/Desktop/Wave Screenshots/Wave', num2str(j),'.png'))
end
%M(j)= getframe;
    pause(0.05)
    
    
%      if j==10 || j==20 || j==50 || j==100 || j==200 || j==300 || j==500
%      	figure
%          subplot(2,1,1)
%          surf(-u(:,:,j))
%          hold on
%          %Z=0.04.*ones(L,L);
%          %Z(:,L/2+3)=-0.04;
%          %Z(:,L/2-3)=-0.04;
%          %Z(L/2-3,:)=-0.04;
%          %Z(L/2+3,:)=-0.04;
%          %mesh((L/2-3):(L/2+3),(L/2-3):(L/2+3),Z((L/2-3):(L/2+3),(L/2-3):(L/2+3))) 
%          hold off
%          title(['Homogeneous Damped Wave Equation at Timestep ', num2str(j)])
%          view(-37.5,70)
%          colormap('winter')
%          axis([0 Lx 0 Ly -0.1 0.1])
%          subplot(2,1,2)
%          contour(-u(:,:,j))
%          hold on
%          grid on
%          axis([0 Lx 0 Ly])
%      end
end
%movie(M)
%movie2avi(M,'SimpleWaveSim.avi','compression','Cinepak')
%fprintf('The origin of the wave occured at (%1.1f, %1.1f).\n',rx,ry)