clear all;
clc;
close all;
Nsteps=30;
f=2000;
w=2*pi*f;
t=linspace(2*pi/Nsteps/w,2*pi/w,Nsteps);

rho1=1026;	% sea water  (see K&F p. 462)
c1=1500;
r1=rho1*c1;
k1=w/c1;

rho2=rho1/0.9*1.1;
c2=c1*0.9;	% slow bottom, figure 6.5
r2=r1*1.1;
k2=w/c2;

tetai=46.433/360*2*pi;	%input('enter incidence angle/n')
tetat=asin(c2/c1*sin(tetai));

k1x=k1*cos(tetai);
k1y=k1*sin(tetai);
k2x=k2*cos(tetat);
k2y=k2*sin(tetat);
ratio=cos(tetat)/cos(tetai);
R=(r2/r1-ratio)/(r2/r1+ratio);
T=1+R;
x=linspace(-1,1,200);
y=linspace(-1,1,200);

for mm=1:100
  	for nn=1:200
			pt(nn,mm)=exp(-i*(k1x*x(mm)+k1y*y(nn)))+R*exp(i*(k1x*x(mm)-k1y*y(nn)));
	end
end
   
for mm=101:200
  	for nn=1:200
			pt(nn,mm)=T*exp(-i*(k2x*x(mm)+k2y*y(nn)));
	end
end
   
pmag=sqrt(sqrt(abs(pt)));
pphase=angle(pt);
clim=max(max(pmag));
V=[-clim clim];							
pt=pmag.*exp(i*pphase);	

figure(1)
pt2=real(pt);
H=pcolor(x,y,pt2);
axis equal
axis([-1 1 -1 1])

M=moviein(Nsteps); 

for kk=1:Nsteps
   pt2=real(pt*exp(i*w*t(kk)));
   h=pcolor(x,y,pt2);
   axis equal
	xlabel('x');
   ylabel('y');
   title('Pressure')
   shading interp;
   caxis(V);
   axis([-1 1 -1 1])
   
   M(:,kk)=getframe;
% to store the images for web posting (bitmaps)
filnme=strcat('spher',int2str(kk));
eval(['print -dbitmap ' filnme ' -f']);
end

%movie(M,20,12);
