%  waves6
%  Some elementary wave exercises to demonstrate the
%  effects of superposition, interferance, etc., between
%  two sinusoidal waves by way of a Matlab movie. 
%
%  Jim Price, January, 2000

clear
clear memory
close all

%  Set some default graphics things.
set(0,'DefaultTextFontSize',14)
set(0,'DefaultAxesFontSize',14)
set(0,'DefaultAxesLineWidth',1.6)
set(0,'DefaultLineLineWidth',1.4)


str2(1) =  {'Waves6:' };
str2(2) =  {'  '};
str2(3) =  {'Demonstrate the superposition of two sinusoidal '};
str2(4) =  {'propagating waves (the green and red lines). Select'};
str2(5) =  {'the frequency and wavenumber of the green wave by '};
str2(6) =  {'clicking the mouse on the dispersion diagram.'}; 
str2(7) =  {'The superposition of the waves is dashed blue, and '}; 
str2(8) =  {'the envelope is solid blue.  The phase speeds '};
str2(9) =  {'of the base waves are evident by the motion of  '};
str2(10) = {'green and red asterisks; the blue asterisk is the '};
str2(11) = {'carrier (or average).  The group moves at the speed'};
str2(12) = {'of the moving `g`. Hit any key to continue.'};
str2(13) = {''};
str2(14) = {'This code is public domain.'};
str2(15) = {''};
str2(16) = {'Jim Price, January, 2001.'};
str2(17) = {''};
hf3 = figure(10);
clf
set(hf3,'Position',[50 50 600 500])
set(gca,'Visible','off')
text(0, 0.50, str2,'FontSize', 12, 'HorizontalAlignment', 'left')
pause

figure(1)
clf reset
orient tall

x = -15:0.05:15;   %  the distance domain
dt = 0.1;          %  the time step
nstep = 200;       %  number of time steps to take
M = moviein(nstep);

%  define the wavelengths and wavenumbers of two pure waves
%  the base wave (the second pair is selected with the mouse)
k1 = 5.;
omega1 = 5.;

%  plot the dispersion diagram
subplot(2,1,1)
plot(k1, omega1, 'r*')
hold on; plot([0 10], [0 10], ':b')
axis([4 6 4 6])
axis('square')
grid
xlabel('wavenumber, k')
ylabel('angular freq., \omega')
title('Do you believe in group velocity?')

aaa = '  Use the mouse to select the (k, \omega) of the second wave.'

[k2, omega2] = ginput(1)
hold on
plot(k2, omega2, 'g*')


c1 = omega1/k1;
c2 = omega2/k2;
cgc = (omega2 + omega1)/(k2 + k1);
cg = (omega2 - omega1)/(k2 - k1);
deltg = 10; if cg <= 0; deltg = -10.; end


for j = 1:nstep
   
   t = j*dt;
   arg1 = k1*x - omega1*t;
   arg2 = k2*x - omega2*t;
   
   s1 = cos(arg1);
   s2 = cos(arg2);   
   sum = s1 + s2;

sg = 2*cos((k2-k1)*x/2 - (omega2 - omega1)*t/2);

   subplot(2,1,2)  
plot(x, sum, ':b', x, s1, 'r', x, s2, 'g', x, sg, 'b', x, -sg, 'b'); 

hold on
plot(t*c1-10, 0, '*r', t*c2-10, 0, '*g', t*cgc-10, 0, '*b')
hp = text(t*cg-deltg, 0, 'g');
hold off

 xlabel('x distance')
 ylabel('amplitude') 
 axis([-15 15 -2 2]);
   
 M(:,j) = getframe;
   
end

%  rerun the movie

movie(M,2,6);