% Example m-file for number 3 from Problem Set 1.  
% Brian Hornbuckle, January 19, 2012.

clear;              % important to do this at the begining of every m-file

A_1 = 3;            % amplitude of wave y_1(x,t), cm
omega_1 = 20;       % radial frequency (angular velocity) of y_1, rad s^{-1}
beta_1 = 30;        % wavenumber (phase constant) of y_1, rad cm^{-1}

A_2 = -3;           % amplitude of wave y_2(x,t), cm
omega_2 = 20;       % radial frequency (angular velocity) of y_2, rad s^{-1}
beta_2 = 30;        % wavenumber (phase constant) of y_2, rad cm^{-1}

f_1 = omega_1./(2.*pi)          % frequency of y_1, Hz
lambda_1 = 2.*pi/beta_1         % wavelength of y_1, cm

f_2 = omega_2./(2.*pi)          % frequency of y_2, Hz
lambda_2 = 2.*pi/beta_2         % wavelength of y_2, cm

x_min = -1;  x_max = 1;             % limits of x, cm
x = linspace(x_min,x_max,10000);     % generates x vector,  cm

t = 0;              % time t, s

y_1 = A_1.*cos(omega_1.*t - beta_1.*x);      % wave y_1
y_2 = A_2.*cos(omega_2.*t + beta_2.*x);      % wave y_2


figure(1)

subplot(2,1,1)
plot(x,y_1,'b-');  grid on;
xlabel('x, cm');  ylabel('y_1(x,t=0), cm');
axis([-0.7,0.7,-5,5]);                          % zoom in and out using the axis command           

subplot(2,1,2)
plot(x,y_2,'r-');  grid on;
xlabel('x, cm');  ylabel('y_2(x,t=0), cm');
axis([-0.7,0.7,-5,5]);           

% -----------------------------------------------

t = 0:pi./100:2.*pi./50;        % time vector, s

y_1 = zeros(length(t),length(x));  y_2 = zeros(length(t),length(x));    % initialize

for n = 1:length(t);
    y_1(n,:) = A_1.*cos(omega_1.*t(n) - beta_1.*x);      % wave y_1(x,t)
    y_2(n,:) = A_2.*cos(omega_2.*t(n) + beta_2.*x);      % wave y_2(x,t)
end

figure(2)       % y_1 propagating

plot(x,y_1(1,:),'k-',x,y_1(2,:),'b-',x,y_1(3,:),'r-',x,y_1(4,:),'g-');
xlabel('x, cm');  ylabel('y_1(x,t), cm');  grid on;
axis([-0.4,0.4,-5,5]);           % zoom in and out using the axis command
legend(['t = ',num2str(t(1))],['t = ',num2str(t(2))],['t = ',num2str(t(3))],['t = ',num2str(t(4))]);