%
%  w01wave.m -- (djm:  10 jan 2006)
%    - finite-difference (lax-wendroff)
%    - second-order interface?
%

%  wavespeeds
c0 = 1.0;  c1 = 0.5;

%  parameters
alpha = 0.5;  dt = 0.025;  N = 20;  M0 = 400;  M1 = 400;
dx0 = (c0*dt/alpha);   dx1 = (c1*dt/alpha);
x = [(-M0:-1)*dx0 0 (1:M1)*dx1];  xR = x(1);  xL = x(end);

%  theoretical coefficients
AR = (c1-c0)/(c1+c0);  AT = (2*c1)/(c1+c0);  
bot = min([AR 0])-0.1;  top = max([AT 1])+0.1;

%  compact support initial conditions
pc = -10;  pw = 2;
u0 = zeros(size(x));

u0 = (1+cos(pi*(x-pc      )/pw))/2.*(abs(x-pc      )<=pw);
u1 = (1+cos(pi*(x-pc-c0*dt)/pw))/2.*(abs(x-pc-c0*dt)<=pw);

figure(1);  plot(x,u0);  hold on;  plot(x,u0,'k.')
plot([0 0],[-1 2],'r--')
title('\bf scattering by a wavespeed discontinuity')
xlabel('\bf x-axis');  ylabel('\bf u(x,t)')
text(0.8*xR,1.05,['\bf time = 0.0'])
axis([xR xL bot top]);  pause

%  timestep loop
time = 0;

for loop = 1:40
    for j = 1:20
        time = time + dt;
        u = (alpha^2)*diff(u1,2) + 2*u1(2:end-1) - u0(2:end-1);
        u = [0 u 0];

        %  interface condition
        u(M0+1) = (4*u(M0+2)-u(M0+3))/(2*dx1) + (4*u(M0+0)-u(M0-1))/(2*dx0);
        u(M0+1) = u(M0+1)*(2/3)/((1/dx0) + (1/dx1));
        
        u0 = u1;  u1 = u;
    
    end
    
    hold off;  plot(x,u0);  hold on;  plot(x,u0,'k.')
    plot([0 0],[bot top],'r--')
    title('\bf scattering by a wavespeed discontinuity')
    xlabel('\bf x-axis');  ylabel('\bf u(x,t)')    
    text(0.8*xR,1.05,['\bf time = ' num2str(time,2)]) 
    axis([xR xL bot top]);  pause(0.25)
end

plot([xR xL xL xR],[AR AR AT AT],'r:')