function retval = convectRHS(t, u)
% usage: retval = convectRHS(t, u)
% 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 u(x,0) = uIC(x),
% with flags set nonzero if P, Q or F is to be used. 

% 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/initial conditions
%         u(0,t) = B(t) and v(xi,0) = uInit((2*pi/L)*x).

% programmer: t. shores
% last update: 3/1/16
 % octave and matlab compatible
 
global gL gxNodes;

% user set choice of using extra terms or not here
pFlag = false;    % true if P(u) is nonzero, user must define convectP(x,t) below
qFlag = false;    % true if Q(u) is nonzero, user must define convectQ(x,t) below
fFlag = true;      % true if a forcing term is used
bFlag = true;     % true if boundary conditions are to be met

% done with input, so calculate parameters of the problem
v = u;		% if we ignore arguments, values of both PDEs are the same
vCap = FFTransform(v, 0); %  initial FFT value of v, transforms wrt xi
twoPiByL = 2*pi/gL;
N = length(v);
% set the exponents of Fourier approximation with DFT coefficients 
fExponents = [0:N/2, -(N/2-1):-1]';
% 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
v_xi = FFTransform(v_xiCap, 1); % inverse transform of v_xiCap
% we won't need xi explicitly, do all evaluations using x
Dv_xi = convectD(gxNodes).*v_xi;
if (pFlag || qFlag)
  v2 = v.*v;
end
Dv_xiCap = FFTransform(Dv_xi, 0); % transform of Dv_xi
% 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 (pFlag)
  retval = retval + i*fExponents.*FFTransform(convectP(v2, v, t), 0)*twoPiByL;
end
if (qFlag)
  retval = retval + FFTransform(convectQ(v2, v, t), 0);
end
if (fFlag)
 retval = retval + FFTransform(convectF(gxNodes, t),0);
end
retval = FFTransform(retval, 1); % inverse transform of transformed rhs
% correct for left endpoint boundary condition, if requested
if (bFlag)
  retval(1) = convectB(t);
end

end

function retval = convectP(v2, v, t)
% usage: retval = convectP(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.

% Octave and Matlab compatible file
% default is constant value 0.0
retval = 0.0*ones(size(v));
%retval = -0.5*v2;
end

function retval = convectQ(v2, v, t)
% usage: retval = convectQ(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.

% Octave and Matlab compatible file
% default is constant value 0.0
retval = 0.0*ones(size(v));
%retval = (v - v2);
end

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

% Octave and Matlab compatible file
% default is constant value zero.
retval = 0.0;
end

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

% Octave and Matlab compatible file
% 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

function retval = convectF(x,t)
% usage: retval = convectFf(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.

% Octave and Matlab compatible file
% default is constant value zero.
retval = zeros(size(x));
%retval = -0.5*exp(-t/2)*(4*convectD(x).*cos(2*x) + sin(x).*sin(x));
%retval = -exp(-t/2)/2;
%retval = zeros(size(x));
end
