function transfer = current  

% This is a test for conservation
% in solving the system
% p(x) u_t = (D(x) u_x)_x - p(x) K u + p(x) P
% with no flux BC's

TESTS = 1 

LH = 2;
LS = 5;
total_num = 2^6+1;
%total_num = 2^4+1;

% uniform grid (non-uniform causes trouble)
x = linspace(-LH,LS,total_num);                          % face location (n) 
x_mid = x(1:end-1) + diff(x)/2;                          % cell center (n-1)

x = x(:);
x_mid = x_mid(:);

face_idx = [1:length(x)]';
cell_idx = [1:length(x)-1]';
[interface_val,interface_idx] = min(abs(x));% interface_idx=first soil cell_idx

DA = 1;                                             % face quantity
DS = 0.25;                                          % face quantity
D = DA.*(x<=interface_val) + DS.*(x>interface_val); % face quantity
D = D(:); 

Porosity = .6;                                      % cell quantity
WFPS = 0.25;                                        % cell quantity
PS = 0; PA = 0;                                  % cell quantity
KS = 0; KA = 0;                                     % cell quantity
P = PA.*(x_mid<interface_val) + PS.*(x_mid>interface_val); % cell quantity
K = KA.*(x_mid<interface_val) + KS.*(x_mid>interface_val); % cell quantity
dx = diff(x);                                       % cell quantity (length)
dx = dx(:);
P = P(:); K = K(:);

[interface_val,interface_idx] = min(abs(x));   % interface_idx = first soil cell_idx
area = pi*(1-(WFPS*Porosity*(cell_idx>=interface_idx)));   % cell quantity
area = area(:);                                   % cell quantity   
vol = dx .* area;                                 % cell quantity

m0 = 1*(cell_idx<interface_idx);                  % mass
c0 = m0 ./ vol;                                   % concentration
m0 = m0(:);
c0 = c0(:);

%options = [];
%options = odeset('BDF','on');
options = odeset('RelTol',1e-6,'AbsTol',1e-8,'BDF','on');

ti = 0;
tf = 50;

tic;
% mass calculation
[t1,m] = ode15s(@cdot1,[ti tf],m0,options,D,area,vol,x,x_mid);
toc

tic;
% concentration calculation
[t2,c] = ode15s(@cdot2,[ti tf],c0,options,D,vol,x,x_mid);
toc

tic;
% concentration calculation including production (P) and consumption (K)
[t3,c_full] = ode15s(@cdot3,[ti tf],c0,options,D,P,K,vol,x,x_mid);
toc


% TESTS
if TESTS == 1
   c00 = rand.*(x_mid>interface_val);
   alpha = rand;
   j = []; jj = [];
   [j,jj] = ode15s(@cdot3,[ti tf],c00,options,0,P+alpha*(x_mid<interface_val),0,vol,x,x_mid);
   p = polyfit(j,jj(:,interface_idx-1),1);
   p(1) - alpha 
   p(2) -  P(interface_idx-1)
   pause(1);

   polyfit(j,jj(:,interface_idx),1);
   p(1) - alpha 
   p(2) - P(interface_idx)
   pause(1);
   
   polyfit(j,jj(:,interface_idx+1),1);
   p(1) - alpha 
   p(2) - P(interface_idx+1)

   c00 = rand.*(x_mid>interface_val);
   beta = randn/2;
   j = []; jj = [];
   [j,jj] = ode15s(@cdot3,[ti tf],c00,options,0,0,-beta*(x_mid>interface_val),vol,x,x_mid);
   p = polyfit(j,log(jj(:,interface_idx-1)),1);
   p(1) - K(interface_idx-1)
   exp(p(2)) - c00(interface_idx-1)
   pause(1);
   p = polyfit(j,log(jj(:,interface_idx)),1);
   p(1) - beta
   exp(p(2)) - c00(interface_idx)
   pause(1);
   p = polyfit(j,log(jj(:,interface_idx+1)),1);
   (p(1) - beta) ./ beta 
   (exp(p(2)) - c00(interface_idx+1)) ./ c00(interface_idx+1) 
   pause(1);

end

[M,N] = size(m);
vol_m = repmat(vol(:)',M,1);
conc = m ./ vol_m;

[M,N] = size(c);
vol_m = repmat(vol(:)',M,1);
mass = c .* vol_m;

[M,N] = size(c_full);
vol_m_full = repmat(vol(:)',M,1);
mass_full = c_full .* vol_m_full;

% write all variables to 'transfer' pointer array
vars = who;
for i = 1:length(vars)
    cname = ['transfer.',vars{i}];
    eval([cname,'=',vars{i},';']);
end


flag = input('Plot?');
if flag == 1
   figure;
   [xxx,ttt] = meshgrid(x_mid,t1);
   mesh(xxx,ttt,conc);
   ylabel('time');
   xlabel('space');
   title('Concentration (from mass) profile over space and time');

   figure;
   [xxx,ttt] = meshgrid(x_mid,t2);
   mesh(xxx,ttt,c);
   ylabel('time');
   xlabel('space');
   title('Concentration (2) profile over space and time');

   figure;
   [xxx,ttt] = meshgrid(x_mid,t3);
   mesh(xxx,ttt,c_full);
   ylabel('time');
   xlabel('space');
   title('Concentration (3) profile over space and time');

   figure;
   total_m0 = sum(m0); 
   subplot(3,1,1);
   plot(t1,sum(m,2) - total_m0);
   xlabel('t');
   ylabel('mass-m_0');
   title('(mass-m_0) vs. time');
   subplot(3,1,2);
   plot(t2,sum(mass,2)-total_m0);
   xlabel('t');
   ylabel('mass-m_0');
   title('(conc*vol) = (mass-m_0) vs. time');
   subplot(3,1,3);
   plot(t3,sum(mass_full,2)-total_m0);
   xlabel('t');
   ylabel('mass-m_0');
   title('(conc_{full}*vol) = (mass-m_0) vs. time');

end
