% Usage: [u, v] = pset3prob2(x2, T, b, c, lambda, sigma2, stfreq, stcen, sxloc)
%
% Solves the leap-frog wave equation for periodic boundaries with a source:
%        exp(-(t-stcen)^2) * sin(stfreq * (t-stcen))
% located at x=sxloc.  (For the problem set, stcen = stfreq = 5, and sxloc = 10.)
%
% x2 is the grid of x points in the unit cell, at *twice* the resolution.
% That is, it includes *both* the u and v staggered-grid points.
% (For the problem set, x2 = [0:0.05:19.95].)
%
% T is the time to run for.  We start with initial conditions u = v = 0.
% b and c are the constants in the u and v equations, and lambda = dt/dx.
%
% sigma2 is the PML periodicity sigma(x) at the x2 points.
%
% Returns u and v at time T.  u is at x2(1:2:end) and v is at x2(2:2:end).

function [u, v] = pset3prob2(x2, T, b, c, lambda, sigma2, stfreq, stcen, sxloc)
  xu = x2(1:2:end);
  xv = x2(2:2:end);
  u = zeros(1,length(xu));
  v = zeros(1,length(xv));
  % note: we must have length(xu) == length(xv)

  sigu = sigma2(1:2:end);
  sigv = sigma2(2:2:end);

%  h = semilogy(xu, abs(u));
%  axis([min(x2) max(x2) 1e-12 1]);
  h = plot(xu, u);
  axis([min(x2) max(x2) -1 1]);
  set(h,'EraseMode', 'xor');
  dx = (max(xu) - min(xu)) / (length(xu) - 1);
  dt = dx * lambda;
  t = 0;
  
  sx = (abs(xu - sxloc) == min(abs(xu - sxloc))) / dx;
  sx0 = find(sx ~= 0); sx0 = sx0(2:end); sx(sx0) = 0;

  while (t < T)
     t = t + dt;
    stval = exp(-(t-stcen)^2) * sin(stfreq*(t-stcen));

    u = (u .* (1 - dx*sigu/2) + lambda * b * (v - [v(end),v(1:end-1)]) + dt * sx * stval) ./ (1 + dx*sigu/2);
    v = (v .* (1 - dx*sigv/2) + lambda * c * ([u(2:end),u(1)] - u)) ./ (1 + dx*sigv/2);

    set(h, 'XData', xu, 'YData', u);
%    set(h, 'XData', xu, 'YData', abs(u));
    drawnow;
  end