function EqnBIC = EqnBIC_TransCD_ExConstVelPtInfsn
%% EqnBIC: Trans.conv.-diff.: Const.vel., Point infusion, No source 
% Example on [0,1]^2 
% Adopted from JohnNovo_JCP_2012 
% James Liu, ColoState; 2012/07--2017/02

% These parameters can be adjusted
T = 0.5;  % Time period [0,T]  % T=1 or T>1 are more interesting 
D = 1e-6;  % Diffusion coefficient 
a = 1;  b = 0;  % Velocity vector v=[a,b]
xc = 0;  yc = 0.5;  % Position of the constant point infusion 

EqnBIC = struct('T',T, ...
  'fxnv',@fxnv, 'D',D, 'fxnf',@fxnf, ...
  'fxncD',@fxncD, 'fxncN',@fxncN, 'fxnc0',@fxnc0, ...
  'fxnc',@fxnc, ...
  'fxnBkwdChara',@fxnExactChara);

% Known velocity field 
function v = fxnv(pt) 
  x = pt(:,1);  % y = pt(:,2);
  v = [a*ones(size(x)), b*ones(size(x))];
end

% Known no-source: The right-hand side function in the conv.-diff. eqn. 
function f = fxnf(pt,t)
  x = pt(:,1);  % y = pt(:,2);
  f = zeros(size(x));
end

% Dirichlet boundary condition for the concentration 
function cD = fxncD(pt,t) 
  x = pt(:,1);  y = pt(:,2);
  cD = zeros(size(x));
  k = find(x==xc & y==yc);  cD(k) = 1;
end

%% JL20151228: TO BE FINISHED 
% Neumann boundary condition for the concentration 
% Homo. Neumann cond. on the right bndry. x=1 

% Known initial concentration
function c0 = fxnc0(pt)
  x = pt(:,1);  y = pt(:,2);
  c0 = zeros(size(x));
  k = find(x==xc & y==yc);  c0(k) = 1;
end

% No known exact concentration or just set to 0 
function c = fxnc(pt,t)
  x = pt(:,1);  % y = pt(:,2);
  c = zeros(size(x));
end

% Known exact characteristics (backward, true for forward also)
% from head (x,y,t1) backward to foot (xstar,ystar,tnew,t0)
function qt = fxnExactChara(pt,t1,t0)
  x = pt(:,1);  y = pt(:,2);  dt = t1 - t0;
  xstar = x - a * dt;
  ystar = y - b * dt;
  qt = [xstar,ystar];
end

end