function animson1d
clear all
close all
disp('Animation of sound waves')
cdef=1;
while 1
	c=input(sprintf('velocity c ? 9: quit, default: c= %0.3g ? ',cdef));
	if c==9
		break
	end
	if isempty(c)==1
		c=cdef;
	end	
	animwave(c)
end
%
% << animwave >> 
function animwave(c)
% Parameters to ajust
% ddx : distance between particules .1
dxx=.05;
% amp : amplitude of the particules oscillations .05
amp=.07;
% nlabel : printing label every nlabel dt  100 
nlabel=25
% l :  half length of the domain  1 
l=1;
% 
% Initialisation of the domain 
dx=.01;
x=-l:dx:l;
xx=-l:dxx:l;
zzero=0*xx;
k1=2*pi;
%c=2;
omega=k1*c;
pm=1; rhom=1; um=1;
% Time loop
figure
dt=.025;
tmax=1;
dtlabel=nlabel*dt;
nmax=round(tmax/dt)+1;
n=0;
for t=0:dt:tmax
	t_tmax=[t, tmax]
	s=sin(k1*x-omega*t);
	ss=sin(k1*xx-omega*t);
	p=pm*s;
	rho=rhom*s;
	u=um*s;
	uu=um*ss;
	% Premiere figure
	% pressure
	subplot(3,1,1)
	plot(x,0*u,'k-')
	hold on 
	plot(x,p,'m','LineWidth',2)
	grid on 
	if mod(t,dtlabel)==0
		title('pressure','FontSize',20,'Color','m')
		ylabel('p','FontSize',20,'Color','m')
	end
	hold off
	% density 
	subplot(3,1,2)
	plot(x,0*u,'k-')
	hold on 
	plot(x,rho,'g','LineWidth',2);
	grid on ; 
	if mod(t,dtlabel)==0
		title('density','FontSize',20,'Color','g')
		ylabel('\rho','FontSize',20,'Color','g')
	end
	hold off
	% velocity and particules
	subplot(3,1,3)
	zeta=amp*cos(-omega*t);
	xxp=xx+amp*cos(k1*xx-omega*t);
	xp=x+amp*cos(k1*x-omega*t);
	plot(x,0*u,'k-')
	hold on 
	plot(xp,u,'r','LineWidth',1);
	hold on
	for delh=-.8:.2:.8		
		plot(xxp,delh+zzero,'.','LineWidth',3)
		plot(zeta,delh,'o','LineWidth',1)
	end
	sca=.1;
	axis([-l l -1 1])
	grid on 
	if mod(t,dtlabel)==0
		title('velocity','FontSize',20,'Color','r')
		ylabel('u','FontSize',20,'Color','r')
		xlabel('x','FontSize',20)
		pause
	end
	hold off
	drawnow
	if n==0 
		M=moviein(nmax);
	end
	n=n+1
	M(:,n)=getframe; 
	% filew(gcf,n);
end
if n==nmax
	movie(M,2,4)
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')

