% doppler_sirene.m
%
% P. Flandrin

N = 100;
t = 1:N;
M = 100;
ph = linspace(0,2*pi,M);
a = 1.8;% 1.8 .75
r = a*t;
p = 5;
Q = 40;
for k = 1:Q
	eval(['r' num2str(k) ' = [zeros(1,k*p) a*(1:N-k*p)];'])
end

figure(1)
clf

pause

for k = 1:N
	plot(-N+r(k)*cos(ph),r(k)*sin(ph))
	hold on
	plot(-N,0,'o')
	for m = 1:Q
		eval(['plot(-N-p+1*p+r' num2str(m) '(k)*cos(ph),r' num2str(m) '(k)*sin(ph))'])
	end
	plot([-N/2 -N+1-1],[-3*N/4 0],'LineWidth',2)
	plot(-N/2,-3*N/4,'.','MarkerSize',25)
	plot([-3*N/2 N/2],[0 0])
	axis('square')
	axis([-3*N/2 N/2 -N N])
	set(gca,'XTick',[],'YTick',[])
	pause(.5)
	hold off
end

pause

figure(1)
clf

for k = 1:N
	plot(-N+r(k)*cos(ph),r(k)*sin(ph))
	hold on
	plot(-N+k-1,0,'o')
	for m = 1:Q
		eval(['plot(-N+m*p+r' num2str(m) '(k)*cos(ph),r' num2str(m) '(k)*sin(ph))'])
	end
	plot([-N/2 -N+k-1],[-3*N/4 0],'LineWidth',2)
	plot(-N/2,-3*N/4,'.','MarkerSize',25)
	plot([-3*N/2 N/2],[0 0])
	axis('square')
	axis([-3*N/2 N/2 -N N])
	set(gca,'XTick',[],'YTick',[])
	title(['Doppler effect (V/c = ',num2str(1/a,2),')'],'FontSize',16)
	pause(.5)
	hold off
end
