function EqnBIC = EqnBIC_TransCD_Ex23RotCylin
%% EqnBIC: Trans.CD: A clockwise rotating cylinder "plume" w/ no source
% Example 2.3 on [0,1]^2;  Rot. around dom. center;  Allowing no-diffusion
% See BochevPeterGao_CMAME_2013, DengGin_NLAA_2014 
% James Liu, ColoState; 2012/07--2017/02 

% These parameters can be adjusted 
T = 2*pi;  % Time period [0,T], 2*pi for a complete rotation 
D = 1e-5;  % Diffusion coefficient 
xc = 0.25;  yc = 0.25;  % Center of the initial "plume" 
r0 = 0.2;  % Size/Radius 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 = [y-0.5, -(x-0.5)];
end

% The right-hand function in the PDE 
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));
end

% Neumann boundary condition for the concentration 
% --None 

% Known initial concentration 
function c0 = fxnc0(pt)
  x = pt(:,1);  y = pt(:,2);
  r = sqrt((x-xc).^2 + (y-yc).^2);
  c0 = zeros(size(x));
  k = find(r<=r0);  c0(k) = 1; 
end

% ??? Known exact concentration 
function c = fxnc(pt,t)
  x = pt(:,1);  y = pt(:,2);
  X = 2*(x-0.5);  Y = 2*(y-0.5);
  Xstar = cos(t)*X - sin(t)*Y;
  Ystar = sin(t)*X + cos(t)*Y;
  xstar = 0.5 * (Xstar+1);
  ystar = 0.5 * (Ystar+1);
  tmp = (xstar-xc).^2 + (ystar-yc).^2;
  % c = sqrt(sigma22/(sigma22+4*D*t)) * exp(-tmp/(sigma22+4*D*t));
  c = 0;
end

end