% Senior Honors II course: Problem Set #2
%
% Example solutions by Alex Barnett 2/10/04. Cut and paste into
% your command window, chunk by chunk. Or, just run the whole script.

% -------------------- Question 1 -----------------------
clear;   % clear all variables

% choose range of alpha values which can be used for either plot...
alpha = logspace(-3, 3, 100);  % alpha is now a vector of values

% Rewrite a,b formulae from lecture 3 in terms of alpha: fill a,b vectors
a = (1-alpha)./(1+alpha);   % note ./ (element-by-element division), not /
b = 2./(1+alpha);
% since a,b real don't need complex conjugate:
R = a.^2;                   % note .^2 (each element is squared).
T = alpha .* b.^2;

figure; plot(alpha, [a; b; R; T], '-'); % plotting multiple curves use ;
legend('a', 'b', 'R', 'T'); title('HW #2, 1a');
xlabel('\alpha \equiv k''/k'); % note double-quote to get ' inside quotes
ylabel('coefficients');
axis([0 5 -1 2]);

figure; loglog(alpha, [R; T], '-');
legend('R', 'T');  title('HW #2, 1b');
xlabel('\alpha \equiv k''/k'); ylabel('power coefficients');

% -------------------- Question 2 -----------------------
clear;
alpha = 5;  % alpha now is just a scalar (Matlab calls 1-by-1 object)
a = (1-alpha)/(1+alpha);   % compute coefficients a,b
b = 2/(1+alpha);
c = 1;                % given
kl = 1;               % left side
kr = kl*alpha;        % right side
omega = c*kl;         % omega = c.k using left side

xl = -10:.1:0;        % x values on left side (including 0)
xr = 0:.1:10;         % x values on right side (including 0)

% time-independent complex wave, different vectors of values for
% each side...
ul = exp(i*kl*xl) + a*exp(-i*kl*xl);  % i=sqrt(-1) if not defined
ur = b*exp(i*kr*xr);

% t=0 so the e^iwt factors are all 1, can ignore.
figure; plot(xl, real(ul), 'b-', xr, real(ur), 'b-'); % both b- linetype
title('HW #2, 2a');
xlabel('x'); ylabel('u(x,t=0)');

t = pi/2;
% include time evolution factor
ult = ul.*exp(-i*omega*t);
urt = ur.*exp(-i*omega*t);
figure; plot(xl, real(ult), 'b-', xr, real(urt), 'b-'); % both b- linetype
title('HW #2, 2b');
xlabel('x'); ylabel('u(x,t=\pi/2)');
text(-8, 1.8, 'Notice how the amplitude ratio appears very different at this time');

% animation... reuse above code with t increasing in a loop...
figure;
for t=0:0.1:100
  ult = ul.*exp(-i*omega*t);
  urt = ur.*exp(-i*omega*t);
  plot(xl, real(ult), 'b-', xr, real(urt), 'b-');
  axis([-10 10 -2 2]);   % need to fix range otherwise will jump around
  title('HW #2, 2c');
  xlabel('x'); ylabel('u(x,t)'); text(1,1, sprintf('t=%g', t));
  drawnow;
  pause(.05);
end
