%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  code7c.m -- spectral wave w/ advective nonlinearity
%
%		u  + c u u  = d u    OR  d u
%		 t        x      xx         xxx
%
%  spectral, integrating-factor formulation:  U = F[u]; r=2,3
%
% 		(U exp(-d (ik)^r t))  = exp(-d (ik)^r t) (-0.5*c) (F[u^2]) 
%		                    t                                     x
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

global ik

choice = 1;
switch choice
  case {1}
	c  = 1.0;  d = 0.01;  r = 2;  par = [c,d,r];
	L  = 2*pi;  N = 256;
	dt = 0.005;  Nsteps = 200;  Ncycle = 5;
 	blurb = '\bf Burgers - diffusive regularization';
end

%  set grids: wavenumbers & coefficients
ik = i*(2*pi/L)*(-N/2:N/2-1)';  fact = exp(-d*dt*(ik.^r));
 
dx = L/N;  xg = (0:dx:L-dx)';  x = [xg ; L];

%  IV initialization & FFT
switch choice
  case {1}
    IVs  = sech(2*(xg-2)).^2;  pmin = -0.1;  pmax = 1.1;
  case {3}
    IVs  = exp(-2*(xg-2).^2);  pmin = -0.1;  pmax = 1.1;
  case {4}
    IVs  = sin(xg);  pmin = -1.5;  pmax = 1.5;
end

%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
%
%  characteristics plot for u_t + u u_x = 0
%

X0  = 0:1/12:7;
T0  = zeros(size(X0));
T1  =  ones(size(X0)) * 4*pi;
X1  = X0 + sech(2*(X0-2)).^2.*T1;

clf;  subplot(2,1,1);  hold on
axis([0 2*pi 0 5])

for k = 1:size(X0,2)
  plot([X0(k) X1(k)],[T0(k) T1(k)],'b')
end

title('\bf characteristics for advection nonlinearity')
xlabel('\bf x-axis');  ylabel('\bf t-axis')

subplot(2,1,2);  hold on;  axis([0 L pmin pmax]);
X   = -2*pi:pi/50:2*pi;
axis([0 2*pi -0.1 1.1])

for T = 1:1:4
  plot(X+T*sech(2*(X-2)).^2,sech(2*(X-2)).^2,'k:')
end
plot(X+sech(2*(X-2)).^2*5,sech(2*(X-2)).^2,'r:')

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

t = 0;  ug = IVs;  u = [ug ; ug(1)];
plot(x,u,'r');  plot(x,u,'k','markersize',7);

vg = fftshift(fft(ug))/N;

%  timestep loop

for k = 1:Ncycle
  for j = 1:Nsteps
	[vg,tj] = code7a_rk(vg,0,dt,par);

    vg = vg./fact;
    t  = t + dt;
  end
  
  ug = real(ifft(fftshift(vg))*N);
  plot(x,[ug ; ug(1)],'b');  pause(0.1)
end

plot(x,[ug ; ug(1)],'r');  
plot(x,[ug ; ug(1)],'k','markersize',7);
title(blurb);  xlabel('\bf x-axis');  ylabel('\bf u-axis')