%group velocity demo
%enter direction as:  
%       1.) waves traveling the same direction :     direction = 1
%       2.) waves traveling the opposite direction:  direction = -1
direction= -1;
Nsamples = 100;%Total number of displayed frames
Tlimit   = 10; %Number of cycles of HF
lambda1  = 500e-9;
lambda2  = 500e-9;
n1       = 1;  %index of refraction
n2       = 1;  %index of refraction


xlim     = 4e-6;%Spatial range
x        = linspace(-xlim,xlim,500);
c        = 3e8;
k1       = 2*pi/lambda1*n1;
k2       = 2*pi/lambda2*n2;
w1       = k1*c/n1;
w2       = k2*c/n2;
kbar     =(k1+k2)/2;
dk       =(k1-k2)/2;
wbar     =(w1+w2)/2;
dw       =(w1-w2)/2;
Tphase   = abs(2*pi/wbar);  %Period of HF

   tvec     = linspace(0,Tlimit*Tphase,Nsamples);

   for t=tvec

     u   = cos(k1*x-w1*t)+cos(k2*x-direction*w2*t);
     I   = abs(u).^2;
	  plot(x*1e6,u,x*1e6,abs(u).^2)
     axis([-xlim*1e6 xlim*1e6 -2 4])
     xlabel('microns')
	  title('combination of 2 waves')
     pause(0.1)
   end   %for