function retval = marchRHS(t, vCap)
% usage: retval = marchRHS(t, vCap)
% description: This function takes input of scalar time t,
%  current u and returns the rhs of the PDE
%          u_t = (D(x)*u_x)_x + P(u)_x + Q(u) + F(x,t)
% on the domain 0 <= x <= L, t > 0 with boundary conditions
% u(0,t) = B(t) and flags set nonzero if P, Q, B or F is to be used. 
% all x-values are defined in the global gxNodes.

% Method here is to convert to interval 2*pi, by change of variables u(x, t) = v(xi, t),
%  where xi  = (2*pi/L)*x, so that du/dx = (dv/dxi)*(2*pi/L) and the PDE becomes
%        v_t = (2*pi/L)^2*(D(x)*v_xi)_xi + (2*pi/L)*P(v)_xi + Q(v) + F((L/(2*pi))*xi,t)
% on the domain 0 <= x <= 2*pi, t > 0 with boundary condition u(0,t) = B(t).
 
global gL gxNodes gpFlag gqFlag gfFlag gbFlag;

% done with input, so calculate parameters of the problem
twoPiByL = 2*pi/gL;
N = length(vCap);
% set the exponents of Fourier approximation with DFT coefficients 
fExponents = [0:N/2, -(N/2-1):-1]';
%fExponents = [0:N/2-1 0 -N/2+1:-1]'; % slightly different if you use fft, ifft
% calculate DFT of v_xi
v_xiCap = i*fExponents.*vCap;
% calculate v, v^2, u_x, D*v, P(v), Q(v) and F(x,t) and transforms
% dealias vCap back to physical space with 2/3 rule before inverse transforming
uCap(ceil(N/3):(N - floor(N/3)+1)) = 0;
v = FFTransform(vCap, 1); % inverse FFT to compute v
%v = ifft(vCap); % if you need more speed, but be consistent 
% correct for left endpoint boundary condition, if requested
if (gbFlag)
  v = marchB(t);
end
v_xi = FFTransform(v_xiCap, 1); % inverse transform of v_xiCap
%v_xi = ifft(v_xiCap); % if you need more speed
% we won't need xi explicitly, do all evaluations using x
Dv_xi = marchD(gxNodes).*v_xi;
if (gpFlag || gqFlag)
  v2 = v.*v;
end
Dv_xiCap = FFTransform(Dv_xi, 0); % transform of Dv_xi
%Dv_xiCap = fft(Dv_xi); % if you need more speed
% ok, ready to compute rhs  -- start with FFT of (Dv_xi)_xi scaled
retval = twoPiByL*twoPiByL*i*fExponents.*Dv_xiCap;
% account for P(v), Q(v) and F(x,t) terms, if requested
if (gpFlag)
  retval = retval + i*fExponents.*FFTransform(marchP(v2, v, t), 0)*twoPiByL;
  %retval = retval + i*fExponents.*fft(marchP(v2, v, t))*twoPiByL; % if you need speed
end
if (gqFlag)
  retval = retval + FFTransform(marchQ(v2, v, t), 0);
  %retval = retval + fft(marchQ(v2, v, t));
end
if (gfFlag)
  retval = retval + FFTransform(marchF(gxNodes, t), 0);
  %retval = retval + fft(marchF(gxNodes, t)); % if you need speed
end
end    % marchRHS

function retval = marchP(v2, v, t)
% usage: retval = marchP(v2, v, t)
% description: This function accepts input vectors v2, v
% and scalar t, and computes polynomial P(u)  for
% the convection-diffusion eqn.  Coefficients in
% a*v2 + b*v may depend on t.

% default is constant value 0.0
%retval = 0.0*ones(size(v));
retval = -0.5*v2;
end    % marchP

function retval = marchQ(v2, v, t)
% usage: retval = marchQ(v2, v, t)
% description: This function accepts input vectors v2, v
% and scalar t, and computes polynomial Q(u)  for
% the convection-diffusion eqn.  Coefficients in
% a*v2 + b*v may depend on t.

% default is constant value 0.0
retval = 0.0*ones(size(v));
%retval = (v - v2);
end    % marchQ

function retval = marchB(t)
% usage: retval = marchB(t)
% description: This function accepts input scalar t and
% computes the left boundary condition for diffusion  eqn
% at time t.  marchB is assumed at both endpoints
% due to the assumption of periodicity in x of solution.

% default is constant value zero.
retval = 0.0;
end    % marchB

function retval = marchD(x)
% usage: retval = marchD(x)
% description: This function accepts input vector x
% and computes diffusion coefficient for the convection-
% diffusion eqn.  march(0)  should match the initial
% conditions. We assume IC and soln are periodic
% of period L that is specified in the main problem.

% default is constant value 1.0
retval = 1.0*ones(size(x));
%retval = 0.5*ones(size(x));
%retval = 0.25*ones(size(x));
end   % marchD

function retval = marchF(x,t)
% usage: retval = marchFf(x,t)
% description: This function accepts input vector x, scalar
% time t, and computes a time and space dependent
% forcing term for the convection-diffusion problem.

% default is constant value zero.
retval = zeros(size(x));
%retval = -0.5*exp(-t/2)*(4*marchD(x).*cos(2*x) + sin(x).*sin(x));
%retval = sin(x).*(1 + 0.25*t + t*t*cos(x));
%retval = cos(x).*(1 + 0.25*t - t*t*sin(x) );
%retval = 0.5*(t + x.*(4*pi*(1-t) + x.*((8*pi*pi*t*t + 2*t - 2) + x.*(- 8*pi*t*t + 2*t*t*x))));
%retval =  0.5*(-2*pi*t + x.*(3*t + x.*(4*pi*(1-t) + x.*(2*(t-1) + x.*(8*pi*pi*t*t + x.*(-8*pi*t*t + 2*t*t*x))))));
end    % marchF

