% Numerically calculate soln. to advection eqn. dq/dt + dq/dx = 0 on
% domain 0<x<1 with periodic BCs using spectral method with RK4
% time-differencing scheme.

  N = 32; % Number of modes
  nu = 0.5; % Courant number 
  L = 1;  % Domain size
  x = L*(0:(N-1))/N;  %  x-gridpoints [1xN]
  dx = L/N;
  M = [0:(N/2-1) (-N/2):(-1)];  
  k = 2*pi*M/L;  % Wavenumbers [1xN].  

  q0 = sign(0.25-abs(x-0.4));  % Initial condition
  tf = 1; % Final time
  dt = nu*dx; % Timestep
  nt = round(tf/dt); % Number of timesteps to take

  qhat = fft(q0);  % Initial Fourier expansion coeffs [1xN].
  for it = 1:nt

    % March forward dqhat/dt = -Shat using RK4
    % where Shat is DFT of S(q) = dq/dx

    d1 = -dt*1i*k.*qhat; 
    d2 = -dt*1i*k.*(qhat + 0.5*d1);
    d3 = -dt*1i*k.*(qhat + 0.5*d2);
    d4 = -dt*1i*k.*(qhat + d3);
    qhat = qhat + (d1 + 2*d2 + 2*d3 + d4)/6; % New qhat
  end
  q = real(ifft(qhat)); % Numerical q at tf [1xN] (should be real, but take
			% real part for plotting.)

%----------Make plots
% Plot 1: Initial condition showing Gibbs ringing
% Plot 2: Power spectrum of ICs
% Plot 3: Numerical solution at t=1
% Plot 4: Numerical phase speed

  Nf = 256;   % Number of plotting points
  xf = (0:(Nf-1))*L/Nf;
  q0f = sign(0.25-abs(xf-0.4));  % Exact IC on plotting grid

  subplot(2,2,1)
  qif = interpft(q0,Nf);  % Spectrally truncated IC on plotting grid
  plot([xf L],[q0f q0f(1)],'b-',[x L],[q0 q0(1)],'r+',...
       [xf L],[qif qif(1)],'r-') 
  xlabel('x')
  ylabel('q_0')
  axis([0 1 -1.5 1.5])
  title(['Spectral ringing of a square wave IC, N = ',num2str(N)])
  legend('Exact','Gridpoints','Truncated')

  subplot(2,2,2)
  q0hat = fft(q0);
  cfit = 350;
  Mpos = M(M>0);
  semilogy(M,abs(q0hat).^2,'r+',Mpos,cfit*Mpos.^(-2),'r--')
  xlabel('K_n/(2\pi)')
  ylabel('|qhat_n|^2')
  title('Power spectrum of q_0 with k^{-2} fit curve')

  subplot(2,2,3)
  qf = interpft(q,Nf);  % Numerical solution on plotting grid
  plot([xf L],[q0f q0f(1)],'b-',[x L],[q q(1)],'r+',[xf L],[qf qf(1)],'r-') 
  xlabel('x')
  ylabel('q')
  axis([0 1 -1.5 1.5])
  title(['Spectral method, RK4, dt/dx = ',num2str(nu),', t = ' num2str(tf)])
  legend('Exact','Gridpoints','Continuous')

  subplot(2,2,4)
  phasespeed = 1 - angle(qhat./q0hat)./(k*tf);
  plot(M,phasespeed,'+')
  xlabel('K_n/(2\pi)')
  ylabel('Numerical phase speed')  
  title('Phase speed error vs. wavenumber')