% Wave Equation on a rectangle with damping and source
% utt=(uxx+uyy)-ut+f,0<x<10,0<y<10
% Four BC: u=0 (fixed BC)
% IC: u(x,y,0)=0, ut(x,y,0)=0
% Modified from Haroon Stephen, UNLV

clear all; clf;
Lx=10; Ly=10;
dx=0.1; dy=dx;
dt=.5*dx; %Required for stability of numerical method

%Number of intervals in x and y
nx=fix(Lx/dx);ny=fix(Ly/dy); %rounds to nearest integer
x=linspace(0,Lx,nx); y=linspace(0,Ly,ny); % rectangular grid

T=10; % end time

un=zeros(nx,ny); % initially all zero
unm1=un; unp1=un; % u at times t-1 and t+1

t=0;
while (t<T);
       un(:,[1 end])=0; %BC
       un([1 end],:)=0; %BC
       t=t+dt;
       unm1=un; %current value
       un=unp1; %previous value

       % Time-dependent Source f at center of rectangle
       un(50,50)=.06*sin(3*pi/2*t);  %Change position or put more than one
       % un(10,50)=.03*sin(3*pi/2*t); % Another source 
       for i=2:nx-1;
           for j=2:ny-1;
               unp1(i,j)=2*un(i,j)-unm1(i,j)+...
                   (dt/dx)^2*(un(i+1,j)+un(i,j+1)-4*un(i,j)+un(i-1,j)+un(i,j-1));
       %           -(dt)*(un(i,j)-unm1(i,j)); % Change to include damping
           end
       end
       clf;
       subplot(2,1,1);
       imagesc(x,y,un');
       colorbar;
       caxis([-.02 .02]);
       title(sprintf('t=%.2f', t));
       subplot(2,1,2);
       mesh(x,y,un');
       colorbar;
       caxis([-.02 .02]);
       axis([0,10,0,10,-.05,.05]) 
       pause(0.01);
end

    