function EqnBIC = EqnBIC_TransCD_Ex12MovingGauss
%% EqnBIC: Transient CD: A moving Gaussian "plume" w/ lin.vel. & no-source 
% Example 1.2 on [0,1]^2;  Allowing constant velocity or no-diffusion
% James Liu, ColoState; 2012/07--2017/02 

% These parameters can be adjusted 
T = 0.35;  % Time period [0,T] 
D = 1e-3;  % Diffusion coefficient 
a = 1;  a1 = 0.1;  b = 1;  b1 = 0;  % Velocity v=[a+a1*x,b+b1*y] 
xc = 0.25;  yc = 0.25;  % Center of the initial "plume"
sigma22 = 0.01;  % 2*sigma^2: Size of the initial "plume" 

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

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

% The right-hand function in the PDE: No source 
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 = fxnc(pt,t);
end

% Neumann boundary condition for the concentration 
% --None 

% Known initial concentration 
function c0 = fxnc0(pt)
  x = pt(:,1);  y = pt(:,2);
  tmp = (x-xc).^2 + (y-yc).^2;
  c0 = exp(-tmp/sigma22);
end

% Known exact concentration 
function c = fxnc(pt,t)
  x = pt(:,1);  y = pt(:,2);
  if (a1~=0) 
    xstar = (1/a1) * (exp(-a1*t)*(a+a1*x) - a);
  else 
    xstar = x - a * t;
  end
  if (b1~=0) 
    ystar = (1/b1) * (exp(-b1*t)*(b+b1*y) - b);
  else 
    ystar = y - b * t;
  end
  tmp = (xstar-xc).^2 + (ystar-yc).^2;
  c = sqrt(sigma22/(sigma22+4*D*t)) * exp(-tmp/(sigma22+4*D*t));
end

% Known exact characteristics (backward, true for forward also) 
% from head (x,y,t1) backward to foot (xstar,ystar,t0) 
function qt = fxnExactChara(pt,t1,t0)
  x = pt(:,1);  y = pt(:,2);  dt = t1 - t0;
  if (a1~=0) 
    xstar = (1/a1) * (exp(-a1*dt)*(a+a1*x) - a);
  else 
    xstar = x - a * dt;
  end
  if (b1~=0) 
    ystar = (1/b1) * (exp(-b1*dt)*(b+b1*y) - b);
  else 
    ystar = y - b * dt;
  end
  qt = [xstar,ystar];
end

end