s = 340; %array size (spacial resolution of the simulation)
l = 440;

%create arrays
p = zeros(l,s);  %past
c = zeros(l,s);  %current
f = zeros(l,s);  %future
r = zeros(l,s);

speed = 10;  			%propagation speed
dt = 0.0004; 			%timestep (small!!) (time resolution of the simulation)
dx = 0.01; 				%distance between elements
rr = speed * dt / dx; 
r(1:l, 1:s)   = speed * dt / dx;
r(1:l, 320:s)  = 0.001*speed * dt / dx;
r(1:l, 1:20)   = 0.001*speed * dt / dx;
stopTime = 0.6;        	%stop time

%for dynamic boundary
n = 300;

%initial loop
%for i = s/2 - s/16 : s/2 + s/16
	%initial conditions
%	c(i,1:s) = 0.5 * sin(2*pi/s*[1:s]) * sin(i *2*pi/s);
	%c(i, s/2 - s/16 : s/2 + s/16) = - 2 * cos(0.5 * 2 * pi / (s/8) * [s/2 - s/16 : s/2 + s/16]) * cos(0.5 * 2 * pi / (s/8) * i);
%	p(i,1:s) = c(i,1:s);
%end

%initial condition 2
p(l/2 + l/4 -1 : l/2 + l/4 + 1 , 1:2) = 0;%1.5 * sin(t*0.0*n);
c(l/2 + l/4 -1 : l/2 + l/4 + 1 , 1:2) = 0;%1.5 * sin(t*dt*n);


%main loop
for t = 0 : dt : stopTime
	%wave equation
  %for i=2:1:s-1
    %for j=2:1:s-1
      %f(i,j) = 2*c(i,j) - p(i,j) + r(i,j)*( c(i-1,j) + c(i+1,j) + c(i,j-1) + c(i,j+1) - 4*c(i,j));
      f(2:l-1, 2:s-1) = 2 * c(2:l-1,2:s-1) - p(2:l-1,2:s-1) + r(2:l-1,2:s-1).*( c(1:l-2,2:s-1) + c(3:l,2:s-1) + c(2:l-1,1:s-2) + c(2:l-1,3:s) - 4*c(2:l-1,2:s-1) );
	  %end
  %end
	%f(s/2-s/8 : s/2+s/8 , 1:2) = 1.5 * sin(t*n);

	%body force (dynamic boundary)
	%f(l/2 + l/4-1 : l/2+l/4+1 , 1:2) = 1.5 * sin(t*n);
	%f(s/2 - s/4-1 : s/2-s/4+1 , 1:2) = 1.5 * sin(t*n);
  %f(1:2, s/2 - s/4-1 : s/2 - s/4+1) = 1.5 * sin(t*n);
  f(1:2, s/2 - 1 : s/2 + 1) = 1.5 * sin(2*t*n);
	%f(2:s-1,1:2) = 1;
	
	%transparent boundaries (no reflection)
	f(2:l-1,1) = (2 * c(2:l-1,1) + (rr-1) * p(2:l-1,1) + 2.*r(2:l-1,1).*r(2:l-1,1).*(c(2:l-1,2)   + c(3:l,1) + c(1:l-2,1) - 3 * c(2:l-1,1)))/(1+rr); % Y:1
	f(2:l-1,s) = (2 * c(2:l-1,s) + (rr-1) * p(2:l-1,s) + 2.*r(2:l-1,s).*r(2:l-1,s).*(c(2:l-1,s-1) + c(3:l,s) + c(1:l-2,s) - 3 * c(2:l-1,s)))/(1+rr); % Y:s
	f(1,2:s-1) = (2 * c(1,2:s-1) + (rr-1) * p(1,2:s-1) + 2.*r(1,2:s-1).*r(1,2:s-1).*(c(2,2:s-1)   + c(1,3:s) + c(1,1:s-2) - 3 * c(1,2:s-1)))/(1+rr); % X:1
	f(l,2:s-1) = (2 * c(l,2:s-1) + (rr-1) * p(l,2:s-1) + 2.*r(l,2:s-1).*r(l,2:s-1).*(c(l-1,2:s-1) + c(l,3:s) + c(l,1:s-2) - 3 * c(l,2:s-1)))/(1+rr); % Y:s
	
	
	%special corners (transparent boundary)
	f(1,1) = ( 2 * c(1,1) + (rr-1) * p(1,1) + 2*rr*rr*( c(2,1)   + c(1,2)   - 2*c(1,1))) / (1+rr);  % X:1 ; Y:1
  f(l,s) = ( 2 * c(l,s) + (rr-1) * p(l,s) + 2*rr*rr*( c(l-1,s) + c(l,s-1) - 2*c(l,s))) / (1+rr);   % X:s ; Y:s
  f(1,s) = ( 2 * c(1,s) + (rr-1) * p(1,s) + 2*rr*rr*( c(2,s)   + c(1,s-1) - 2*c(1,s))) / (1+rr);   % X:1 ; Y:s
  f(l,1) = ( 2 * c(l,1) + (rr-1) * p(l,1) + 2*rr*rr*( c(l-1,1) + c(l,2)   - 2*c(l,1))) / (1+rr);   % X:s ; Y:1
  
	%update values
	p(1:l,1:s) = c(1:l,1:s);
	c(1:l,1:s) = f(1:l,1:s);
	
	%plot
	if mod( t / dt , 4 ) == 0 

    %h=imagesc((1:1:l),(1:1:s)',c',[-1,1]); 
    %daspect  ([1 1 1]);
    %set(gca,'FontSize',20);
    %drawnow
    
		b = surf(c); 				%surface plot
    colormap(bone (64));
    light;
		set(b, 'edgecolor','none') %turn mesh edges off
		axis([0 l 0 s -2 2])  		 %axis scale
		caxis([-0.5 0.5])         	   %color axis scale
		pause ( .03 )
	end
	
end

hold off