function animint2d
% Animation of internal waves
% O. Thual, 17/11/02
clear all
close all
disp('Animation of internal waves')
demodef=1
while 1
	% 1: demo1, 2:demo2
	demo=input('1:demo1, 2:demo2, 9:quit ?' );
	if demo==1 
		demo1
	end
	if demo==2
		demo2
	end
	if isempty(demo)==1
		demo1
	end
	if demo==9
		return
	end
end
% << demo1 >> 
%
function demo1
global Mla nla
% Parameters to ajust
nlamax=20; 
% theta : angle of the wave vecteur
nla=0;
lambdadef=.7;
lambda=input('lambda ? 9 to quit   ')
if isempty(lambda)==1
	lambda=lambdadef
end
if lambda==9
	return
end
T=1/lambda;
delta=T/nlamax;
animwave(lambda,nlamax,delta);
nla
movie(Mla,3,4)
%
% << demo2 >> 
%
function demo2
global Mla nla
% Parameters to ajust
nlamax=20; 
% delta : time step
delta=.25
% theta : angle of the wave vecteur
nla=0; 
for lambda=1:-.3:.3
	animwave(lambda,nlamax,delta);
	nla
end
newplot
movie(Mla,3,3)
%
% FUNCTIONS
% << animwave >>
function animwave(lambda,nlamax,delta)
global Mla nla
theta=acos(lambda); sthe=sin(theta) ; cthe=cos(theta);
% ddx : distance between particules .1
dxx=.05; dzz=.05;
% amp : amplitude of the particules oscillations .05
amp=.07;
% nlabel : printing label every nlabel dt  100 
nlabel=100
% l :  half length of the domain  1 
l=.5; 
% Initialisation of the domain 
dx=.01; dz=dx;
x=-l:dx:l;  z=x;
[X,Z]=meshgrid(x,z);
xx=-l:dxx:l; zz=xx;
[XXa,ZZa]=meshgrid(xx,zz);
XX=XXa*cthe-ZZa*sthe;
ZZ=XXa*sthe+ZZa*cthe;
zzero=0*xx;
omega=2*pi*cthe;  
T=2*pi/omega
k1=2*pi; k3=k1*tan(theta)
pm=1; rhom=1; um=1;
% 
% Time loop
hold off
nmax=round(T/delta);
dt=T/nmax;
dtlabel=nlabel*dt;
% line to draw // to k
xxk=zz*cthe;
zzk=zz*sthe;
% line to draw orthogonal to k
xxe=-zz*sthe;
zze=zz*cthe;
for n=1:nmax
	t=n*dt
	S=sin(k1*X+k3*Z-omega*t);
	SS=sin(k1*XX+k3*ZZ-omega*t);
	C=cos(k1*X+k3*Z-omega*t);
	CC=cos(k1*XX+k3*ZZ-omega*t);
	p=pm*C;
	%	pfix=pm*cos(k1*X+k3*Z);;
	%contour(X,Z,p,[-1:.2:1]);
	contour(X,Z,p,5);
	% particules
	XXp=XX+amp*SS*sthe;
	ZZp=ZZ-amp*SS*cthe;
	% line of particules // to k
	ssv=sin(k1*xxk+k3*zzk-omega*t);
	ccv=cos(k1*xxk+k3*zzk-omega*t);
	xxv=xxk+amp*ssv*sthe;
	zzv=zzk-amp*ssv*cthe;
	%
	hold on
	% plot all particules
	plot(XXp,ZZp,'.','Color','k','LineWidth',2)	
	% plot straight lines
	plot(xxk,zzk,'-k')	
	plot(xxe,zze,'-k','LineWidth',1)
	% plot line of particules // to k
	plot(xxv,zzv,'o','Color','b','LineWidth',2)	
	% legend
	leg=sprintf('pressure,   lambda= %0.3g',lambda);
	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);
		% M=moviein(nmax);
	end
	% M(:,n)=getframe;
	nla=nla+1
	Mla(:,nla)=getframe; 
	%filew(gcf,n);
end
%	movie(M,2,4)
% 
% << 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')
