% hw4p234.m:
% Script for solving the traffic flow problem of HW4 
%  dq/dt + df/dx = 0, 0 < x < 1, f(q) = q*(1 - q)
%  q(x,0) = 0
%  q(0,t) = q0(t) = a sin^2(pi*t)
% for dx = 0.05, using:
%  A leapfrog-Asselin scheme (function hw4_leap)
%  A piecewise-linear FV scheme with downwind and MC slopes (function hw4_fv).
% We plot q(1,t) for all three methods vs. the exact soln for 0 < t < 3.

  a = 0.12;
  dx = 0.05;
  
  % Leapfrog parameters
  
  mumax_leap = 0.5;  % Max Courant number
  gam = 0.05;        % Asselin filter coeff.

  % FV parameters
  
  mumax_fv = 0.9;    % Max Courant number
  
  cmax = 1;          % Max wave speed
  T = 3;             % Length of integration
  L = 1;             % Domain size [0,L]
  nx = 1+round(L/dx);
  x = (0:nx-1)*dx;
  q  = zeros(1,nx); % Initial q(x,t=0)

  % Call leapfrog scheme
  
  dt_leap = mumax_leap*dx;
  t_leap = (0:dt_leap:T);
  q0_leap = a*sin(pi*t_leap).^2; % Define q(0,t) at the nt timesteps.
  q1_leap = hw4_leap(q,q0_leap,x,t_leap,gam);
  
  % Call FV scheme
  
  dt = mumax_fv*dx/cmax;
  t = (0:dt:T);
  q0 = a*sin(pi*t).^2; % Define q(0,t) at the nt timesteps.
  
  q1_LW = hw4_fv(q,q0,x,t,'L');
  q1_MC = hw4_fv(q,q0,x,t,'M');
  
  % Exact solution
  
  t0 = (-1:0.01:T); % Origin time of characteristic at x = 0
  qex = a*sin(pi*t0).^2; % qex(0,t0), t0>0.
  qex(t0<0) = 0;
  tex = t0 + 1./(1 - 2*qex); % qex conserved on dt/dx = 1/c(qex).

  % Make plot
  
  subplot(2,2,1)
  plot(tex,qex,'k-',t_leap,q1_leap,'r--',t,q1_LW,'g:',t,q1_MC,'b-.')
  xlabel('t')
  ylabel('q(1,t)')
  axis([0 3 -0.03 0.15])
  legend('Exact','Leap','LW','MC',2)

  qex_tleap = interp1(tex,qex,t_leap,'spline');
  maxerr_leap = norm(qex_tleap - q1_leap,inf)

  qex_t = interp1(tex,qex,t,'spline');
  maxerr_LW = norm(qex_t - q1_LW,inf)
  maxerr_MC = norm(qex_t - q1_MC,inf)