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

global ik


choice = 4;
switch choice
  case {4}
	c  = 1.0;   d = -0.00064;  r = 4;  par = [c,d,r];
	L  = 2*pi;  N = 128;
	dt = 0.01;  Nsteps = 50;  Ncycle = 8;
 	blurb = '\bf 4^{th}-order hyper-diffusive regularization';
  case {3}
	c  = 1.0;   d = -0.00008;  r = 4;  par = [c,d,r];
	L  = 2*pi;  N = 256;
	dt = 0.005;  Nsteps = 100;  Ncycle = 8;
 	blurb = '\bf 4^{th}-order hyper-diffusive regularization';
  case {1}
	c  = 1.0;   d = -0.00001;  r = 4;  par = [c,d,r];
	L  = 2*pi;  N = 512;
	dt = 0.0025;  Nsteps = 200;  Ncycle = 8;
 	blurb = '\bf 4^{th}-order hyper-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 {5}
    IVs  = exp(-2*(xg-pi).^2);  pmin = -1.1;  pmax = 1.1;
  case {1,2,3,4}
    IVs  = sin(xg);  pmin = -1.5;  pmax = 1.5;
end

t = 0;  ug = IVs;  u = [ug ; ug(1)];
%clf;  hold on;  axis([0 L pmin pmax]);
hold on;  axis([0 L pmin pmax]);
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');  plot(x,[ug ; ug(1)]);  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');  hold off

deriv = (ug(N/2+2)-ug(N/2))/(2/N)