%
%  w10oneway.m -- djm -- 06 mar 06
%    - one-way wave, method of lines
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global c dx

c = 1.0;  dt = 0.25;  T = 3;
Nx = 512;  Lx = 16;  dx = 2*Lx/Nx;

x = -Lx:dx:Lx-dx;

%  initial condition
u0 = zeros(size(x))';  u0(Nx/2+1) = 1.0;
%u0 = exp(-(x.^2));
%u0 = (x>-4).*(x<4);
u1 = u0;

figure(1);  clf;
plot(x,u0,'b.');  axis([-Lx Lx -0.2 0.2]);  pause

for time = 0:dt:T-dt
    [ti,u] = ode45(@w10oneway_ode,[time time+dt],u1);
    uE = u(end,:);
    plot(x,uE,'k.');  axis([-Lx Lx 1.05*min(uE) 1.05*max(uE)]);
    u1 = uE;  pause(0.1)
end
time = ti(end);
hold on
uE = u(end,:);
plot(x,uE,'b-');
plot(x,uE,'k.');  axis([-Lx Lx 1.05*min(uE) 1.05*max(uE)]);

%  airy plot
c1 = (2/c/time/dx^2)^(1/3);
uA = real((c1)*airy(c1*(x-c*time))) .* (x>0);
plot(x,dx*uA,'r-')
plot(x,dx*uA,'kx')
title('\bf airy dispersion')

figure(2);  clf;  hold on
plot(x,uE,'b-');
plot(x,uE,'ko');  axis([-Lx Lx 1.05*min(uE) 1.05*max(uE)]);

%  bessel solution
nn = round(x/dx);
ne = nn(1:2:end);  no = nn(2:2:end);

c2 = c*time/dx;

%  even terms
u_ev = real(besselj(abs(ne),c2));
u_od = real(besselj(abs(no),c2)) .* sign(no);

uB = [u_ev ; u_od];  uB = uB(:);
plot(x,uB,'rx');
title('\bf exact bessel solution')