function animsurf2d
% Animation of surface waves
% O. Thual, 17/11/02
clear all
close all
disp('Animation of surface waves')
demo1
% << demo1 >> 
%
function demo1(h)
global Mla nla
% Parameters to ajust
nlamax=10;
delta=.1;
hdef=.4;
while 1
	leginput=sprintf('depth h [.1:.7] ? defaut :%0.3g, quit :9  ',hdef);
	h=input(leginput)
	if isempty(h)
		h=hdef
	end
	if h==9
		disp('Good Bye!')
		return
	end
	nla=0;
	close all
	animwave(h,nlamax,delta);
	nla
	movie(Mla,5,4)
end
%
% FUNCTIONS
% << animwave >>
function animwave(h,nlamax,delta)
global Mla nla
% ddx : distance between particules .1
ndxx=10; ndzz=10;
% nlabel : printing label every nlabel dt  100 
nlabel=100
% l :  half length of the domain  1 
l1=.5;  dl1=2*l1;
htop=h/5 ; htot=h+htop;
%axes('position',[.1 .1 .8 .8])
correc=1.2
posi=.1+correc*.7*htot/dl1;
positop=min([.8 posi]);
axes('position',[.1 .1 .8 positop])
% Initialisation of the domain
nptx=25; dx=l1/nptx; 
nptz=50; dz=h/nptz;
dxx=l1/ndxx; dzz=h/ndzz;
x=-l1:dx:l1;  z=-h:dz:0;
[X,Z]=meshgrid(x,z);
xx=-l1:dxx:l1; zz=-h:dzz:0;
[XX,ZZ]=meshgrid(xx,zz);
% others
zzerox=0*xx;
zmid=-h/2;
% frequency and wave number
omega=2*pi;  
T=2*pi/omega
k=2*pi;
pm=1;  um=1; 
% amp : amplitude of the particules oscillations .05
amp=.01*sinh(k*.25)/sinh(k*h)*h/.25;
% 
% Time loop
hold off
nmax=round(T/delta);
dt=T/nmax;
dtlabel=nlabel*dt;
for n=1:nmax
	t=n*dt
	% fine grid
	SI=sin(k*X-omega*t);
	CO=cos(k*X-omega*t);
	CH=cosh(k*(Z+h));
	SH=sinh(k*(Z+h));
	% particule grid
	SSI=sin(k*XX-omega*t);
	CCO=cos(k*XX-omega*t);
	CCH=cosh(k*(ZZ+h));
	SSH=sinh(k*(ZZ+h));
	% particule grid
	ssi=sin(k*xx-omega*t);
	cco=cos(k*xx-omega*t);
	cch=cosh(k*(zz+h));
	ssh=sinh(k*(zz+h));
	% all particules
	Xp=X-amp*CH.*CO;
	Zp=Z-amp*SH.*SI;
	% pressure
	p=pm*CH.*SI;
	contour(Xp,Zp,p,5);
	hold on 
	% plot horizontal straight line
	plot(xx,zzerox+zmid,'-k')
	plot(xx,zzerox,'-k')
    hold on	
	% all particules
	XXp=XX-amp*CCH.*CCO;
	ZZp=ZZ-amp*SSH.*SSI;
	UUp=-CCH.*SSI;
	WWp=SSH.*CCO;
	% special horizontal lines of particules
	xxhmid=xx-amp*cosh(k*(zmid+h))*cco;
	xxhtop=xx-amp*cosh(k*h)*cco;
	zzhmid=zmid-amp*sinh(k*(zmid+h))*ssi;
	zzhtop=-amp*sinh(k*h)*ssi;
	% special vertic lines of particules
	zzv=zz-amp*sinh(k*(zz+h))*sin(-omega*t);
	xxv=-amp*cosh(k*(zz+h))*cos(-omega*t);
	% plot all particules
	plot(XXp,ZZp,'.','Color','k','LineWidth',3)	
	quiver(XXp,ZZp,UUp,WWp)
	hold on
	% plot mid horizontal line of particules
	plot(xxhmid,zzhmid,'o','Color','b','LineWidth',2)
	hold on	
	% plot top horizontal line of particules
	plot(xxhtop,zzhtop,'o','Color','b','LineWidth',2)	
	hold on
	% plot  surface line
	plot(xxhtop,zzhtop,'-r','LineWidth',2)
	hold on
	% plot vertical line  of particules
	plot(xxv,zzv,'o','Color','b','LineWidth',2)
	% plot vertical line
	plot(0*zz,zz,'-k','LineWidth',1)
	hold on	
	% legend
	axis([-l1 l1 -h htop])
	leg=sprintf('pressure,   h= %0.3g',h);
	title(leg,'FontSize',20,'Color','m')
	ylabel('z','FontSize',20)
	xlabel('x','FontSize',20)
	hold off
	if nla==0 
		disp(' You can change the window size now !')
		pause
		Mla=moviein(nlamax);
	end
	nla=nla+1
	Mla(:,nla)=getframe; 
	% filew(gcf,n);
end
% 
% << filew >>
%
% writing image on  files 
function filew(gcf,n)
[a,map]=capture(gcf);
filename=sprintf('image%2.2d.jpg',n)
imwrite(a,map,filename,'jpg')

