function out = test_consrv 

global XXX

options = odeset('RelTol',1e-12,'AbsTol',1e-12,'BDF','on');

LA = 1 % length of air column 
LS = 5 % length of soil column
DA = 1
DS = 1e-2 
porosity = 0.1

interface = 60;
xhead = linspace(-LA,0,interface);
xsoil = linspace( 0,LS,2*interface);
dxhead = diff(xhead);
dxsoil = diff(xsoil);
new_grid_pt = min(dxhead(1),dxsoil(1));
x = union( union( xhead , xsoil ) , [new_grid_pt, -new_grid_pt] );
% makes grid uniform at interface

XXX = x;

C_0 = x<=0 ;
mass_0 =  C_0 * dxhead(1);
mass_0(1) = 1/2*mass_0(1);
mass_0(end) = 1/2*mass_0(interface);

[junk,C12] = ode15s(@mycdot,[0:100],C_0,options,DA,DS,0,0,porosity);

C_100 =  C12(end,:);

mass_100 = ((x<0) .* C_100)*dxhead(1) + ((x>=0) .* C_100)*dxsoil(end)*porosity ;

mass_100(1) = 1/2*mass_100(1);
mass_100(end) = 1/2*mass_100(end);

total_mass_0 = sum(mass_0);
total_mass_100 = sum(mass_100);

relative_loss = (total_mass_0 - total_mass_100)/total_mass_0;

out.m0 = mass_0;
out.m100 = mass_100;
out.tm0 = total_mass_0;
out.tm100 = total_mass_100;
out.loss = relative_loss
out.c0 = C_0;
out.c100 = C_100;
out.dxh = dxhead(1);
out.dxs = dxsoil(1);
out.idx = find(x==0);
out.air = x<=0;
out.soil = x>0

% keyboard

%----------------------------------------

function dcdt = mycdot(t,c,DA,DS,P,K,porosity);

% function dcdt = cdot(t,c,DA,DS,P,K,porosity);
%
% porosity * c' = (D*c_x)_x + P - K c
% with boundary conditions 
% DA*c_x(x=-Vh,t) = DS*c_x(x=Vs,t) = 0
% and P and K only active for x >0
% boundary interface @ x = 0 
% DA is diffusion in air
% DS is diffusion in soil

global XXX ;
x = XXX(:);
indx = find(x==0);

dx = diff(x);

aug_c = [c(1) ; c(:) ; c(end)];      % builds in boundary conditions
                                     % by repeated first and last entries
                                     % I make m(i+1) - m(i) = 0 -> m_x =0  

aug_x = [x(1)-dx(1) ; x(:) ; x(end)+dx(end)];
aug_xm = aug_x(1:end-1) + 1/2*diff(aug_x);

dxm = diff(aug_xm);

aug_d = (aug_x<=0) * DA + (aug_x>0)*DS ;  %the diffusion vector

dcdx = diff(aug_c)./diff(aug_x);

porosity_v = (x<0) + porosity*(x>=0);  % set por = soil_porosity at interface 

dcdt  = (diff( aug_d(1:end-1).*dcdx) ./ dxm + ...
      P.* (x>0) - ...
            K.* (x>0) .* aug_c(2:end-1)) ./ porosity_v ;

% keyboard
