function transfer = TRD_SF613g(inputs);

% K = inputs.K;    
% WFPS = inputs.WFPS;
% P_Units = inputs.P;

% function [x,t2,TM] = model2d(D,K)
%
%  NOVEMBER 11 2005 
%
% PHASE 0
% This code simulates methane profiles in gas/soil column
% A soil column is initially allowed to equilibrate for a specified time. 
% Both C12 and C13 methane profiles generated  
%
% PHASE I
% Tagged C13 methane is introduced to solution produced in PHASE 0
% This label is then measured at various points in time.
% Includes concentration of C12 and C13 over 
% the experimental window
%
% -------------------------------
%
% Numerical setup 
%
% PHASE 0 Experiment
%
% The initialization is modeled as 
% p(x) d/dt C(x,t) =  (D * C_x(x,t))_x + p(x) S - p(x) K * C(x,t)
% with boundary conditions
% C(0,t) = C12_atmos (known concentration Methane)
% and  D * d/dx C(Vh,t) = 0
% and initial profile
% C(x,0) = 0;
%
%
% PHASE I Experiment
%
% The second phase of experiment is modeled as
% p(x) d/dt C(x,t) =  (D * C_x(x,t))_x + p(x) S - p(x) K * C(x,t)
% C is calculated as moles/volume for each layer
% with no flux boundary conditions
% D * d/dx C(Vs,t) = 0  = D * d/dx C(Vh,t)`
% and initial profile
% C(x,0) = c12;
%
%
% D = diffusion coefficient
% P = Production rate 
% K = consumption rate
% Tmax = final time of simulation
% Vs = volume of headspace
% Vh = volume of soil
% C_atmos = concentration of methane in atmosphere
% p = porosity 
%
% Geometery  
% Vs < x < Vh,
% 0 < t < Tmax
% interface is grid pt nearest to zero

% declare x grid global
global XXX XMID XS XSMID

%----------------------------
% Begin parameters section
%----------------------------
%Incubation Physical parameters
JarVolume=1000; %mililiters of jar volume
SoilLength=10;  %length of soil core in cm
JarDiameter=4; %diameter of soil core in cm
WFPS=0.; %Percent of soil pores filled with water
Porosity=0.5; %Percent of soil volume that is pores

JarCrossSection=(JarDiameter/2)^2*pi; %area of soil core in cm^2
SoilArea = JarCrossSection;
SoilVolume= SoilArea*SoilLength; %Volume of soil core in mililiters (cm^3)
AirFilledPorosity=(1-WFPS)*Porosity;
VolumeSoilAir=SoilVolume*Porosity*(1-WFPS); %gives volume of soil pore space, ml
VolumeAirHeadspace= JarVolume-SoilVolume; %volume of air in the headspace, ml
LengthAirHeadspace=VolumeAirHeadspace/JarCrossSection; %change this from

AtmosphericMethanePPB=1800; %Concentration of Methane in the free atmosphere, parts per billion
HeadspaceMethane_nMoles=AtmosphericMethanePPB/22.4*.001*VolumeAirHeadspace;

Vh = VolumeAirHeadspace;   % Volume of headspace during the incubation
Vs = VolumeSoilAir;          % Volume of soil
Po = Porosity;           %porosity

% PDE parameters
%D = 1.2e+2;       % diffusion coeff
DA = 60*1.2e-2;    % diffusion in air, cm2/sec converted to cm2/min
DS = DA*0.9*(1-WFPS)^(2.3);   % diffusion in soil, cm2/sec converted to cm2/min, from Campbell (1985)
%DA = 0 ; DS = 0;
%P_Units=100; %Production in normal units of mg CH4-C/m2/day for the whole core
%P_Units=50 %Production in normal units of mg CH4-C/m2/day for the whole core
P_Units=0 %Production in normal units of mg CH4-C/m2/day for the whole core
PS = P_Units/12*1e6*22.4/10000*SoilArea/24/60/Vs*1000;% Production corrected for model units:ppb,cm,min 
%KS = 1;      % Consumption coeff
KS = 0.01;      % Consumption coeff
ConcSF6Init=100;  %Initial concentration of SF6 in the headspace

% Experimental parameters
Tinit = 100001;          % time allowed to equilibrate, minutes
TPreIncub=30;   %time for pre-incubation after addition of 13C-methane
Tmax = 300+TPreIncub;       % max simulation time, minutes
NumGasSamples=4;  % Number of gas samples taken from headspace
TGasSamples=round(linspace(TPreIncub,Tmax,NumGasSamples));  %time that gas samples are taken from the headspace

pC12_C13 = 0.0106;  % percent C13 in new production
p13Catm=0.01080;     % percent of 13C in atmospheric methane
Rpdb=0.01123;       % ratio of 13C to 12C in reference material
Fk = 0.990;         % fractionation during consumption against C13 vs. C12
FD = 0.999;         % fractionation during diffusion against C13 vs. C12
tagged = 5;      % multiplier for conc C13 methane in headspace tagged 

%-----------------
% Discretization parameters 
%Tstep = Tmax;                
%TstepEquil=1000;
%TstepIncub=Tmax+1;     % number of time samples in incubation, 1 per minute
Xstep =  2^7+1;  % min # cell in column
%-----------------

%----------------
% Plots
PLOTS = 1;        % 0 = off 1 = on

%---------------------------
% end parameters section
%---------------------------

% Build space grid 

% face location
x = linspace(-LengthAirHeadspace,SoilLength,Xstep);

% cell_center
x_mid = x(1:end-1) + diff(x)/2;
x = x(:); x_mid = x_mid(:);
[surface_val,surface_idx] = min(abs(x)); % first soil cell idx

xsoil_face = x(surface_idx:end);
xsoil_mid = x_mid(surface_idx:end);
face_idx = [1:Xstep]';
cell_idx = [1:Xstep-1]';

xe=x;
XXX = x;
XMID = x_mid;
XS = xsoil_face;
XSMID = xsoil_mid;

dx = diff(x);
e_porosity = 1 + 0*x_mid;
e_porosity(1:surface_idx-1) = 1;
e_porosity(surface_idx:end) = Porosity*(1-WFPS);
total_area = pi*(JarDiameter/2).^2; % cell
total_Vol = dx .* total_area; % cell
keyboard;
Vol = e_porosity(:) .* total_Vol; 
Vol_soil = Vol(surface_idx:end);

D = DA.*(x<=surface_val) + DS.*(x>surface_val);  % face
K = 0.*(x_mid<surface_val) + KS.*(x_mid>surface_val);  % cell
P = 0.*(x_mid<surface_val) + PS.*(x_mid>surface_val);  % cell
D = D(:); P= P(:); K = K(:);

% initial profile

c12conc = AtmosphericMethanePPB*(1-p13Catm);  % conc natural C12
c13conc = AtmosphericMethanePPB*(p13Catm) ;   % conc natural C13 

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

tspan0 = linspace(0,Tinit,10);
c120 = 0*xsoil_mid;
c130 = 0*xsoil_mid;
[t1,c12e] = ode15s(@cdotinit,tspan0,c120,options,DS,PS*(1-pC12_C13),KS,Vol_soil,c12conc);
[jk,c13e] = ode15s(@cdotinit,tspan0,c130,options,DS*FD,PS*(pC12_C13),KS*Fk,Vol_soil,c13conc);

% Final profiles
c12u= 0*x_mid;
c13u = 0*x_mid;
c12u(1:surface_idx-1) = c12conc;
c12u(surface_idx:end) = c12e(end,:)';  % equib conc natural 12C
c13u(1:surface_idx-1) = c13conc;        
c13u(surface_idx:end) = c13e(end,:)';  % equib conc natural 13C

c13t = (tagged*c13u) .* (x_mid<surface_val) ;          % conc added 13C
sf6 = ConcSF6Init.*(x_mid<surface_val); % sf6 profile
% set options if needed 
% type matlab>> help ode15s for solver options
options = odeset('RelTol',1e-9,'AbsTol',1e-10,'BDF','on');

tspan = [0:Tmax];

tic;

[junk,C12] = ode15s(@cdot,tspan,c12u,options,D,P*(1-pC12_C13),K,Vol);
[junk,C13U] = ode15s(@cdot,tspan,c13u,options,D.*FD,P*pC12_C13,K*Fk,Vol);
[t2,C13T] = ode15s(@cdot,tspan,c13t,options,D.*FD,P*0,K*Fk,Vol);
[t2, SF6] = ode15s(@cdot,tspan,sf6,options,D.*FD,P*0,K*0,Vol);
     
RunTime=toc
m_Vol = repmat(Vol(:)',length(t2),1);
z_check = sum(m_Vol.*SF6,2);
zz_check = (z_check - z_check(1))/z_check(1); % rel error in a
figure(10); plot(t2,zz_check);
TM=(C12+C13U+C13T);      %total methane in system
TMH=mean(TM(:,1:surface_idx-1),2); %Average concentration of methane in the headspace
HSFluxRatePPBMinute=polyfit(t2,TMH,1);
HSFluxRate=HSFluxRatePPBMinute(1)*60*24/22.4*12/1e6*VolumeAirHeadspace/1000*10000/SoilArea; %flux rate, mg/m2/day

TMS=mean(TM(:,surface_idx:end),2); %Average concentration of methane in the soil
SoilFluxRatePPBMinute=polyfit(t2,TMS,1);
SoilFluxRate=SoilFluxRatePPBMinute(1)*60*24/22.4*12/1e6*VolumeSoilAir/1000*10000/SoilArea; %flux rate, mg/m2/day
                                
TotalFluxRate=SoilFluxRate+HSFluxRate;

% changed from c13./c12;
Rall=c13u./c12u;
d13Call=((Rall-Rpdb)./Rpdb).*1000;

T13M=C13U+C13T;
RMH=mean(T13M(:,1:surface_idx-1),2)./mean(C12(:,1:surface_idx-1),2);
APMH=mean(T13M(:,1:surface_idx-1),2)./TMH;
d13CH=((RMH-Rpdb)./Rpdb).*1000;

GasSampled13C=d13CH(TGasSamples);
GasSampleConc=TMH(TGasSamples);
GasSampleAPMH=APMH(TGasSamples);
TGasSamples0= TGasSamples-TPreIncub;

SF6HS=mean(SF6(:,1:surface_idx-1),2);

%Calculate P and k as for 0-D model
options2 = optimset('TolX', 1e-8, 'TypicalX', [1e-10, 1e-10]);

t3=t2-TPreIncub;

TMHnM=TMH./22.4.*VolumeAirHeadspace./1000; %Number of nMoles CH4 in headspace
GasSamplesTMHnM=TMHnM(TGasSamples);  
APEMH=APMH-pC12_C13;   %Atom Percent Excess 13C-CH4 in headspace
LabelHSnM=  APEMH.*TMHnM;%Number of nMoles of label in headspace
GasSamplesLabelHSnM=LabelHSnM(TGasSamples)  %Number nM label in HS for Gas Samples only
TGasSamplesV=TGasSamples(:);

% Calculate the 0-D k and P parameters
k0 = -1 * polyfit(TGasSamplesV,log(GasSamplesLabelHSnM),1); %linear regression on labeled CH4
k0=k0(1);
k0 = k0(1)/Fk;  %k must be corrected for fractionation

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

[bestP0,min_cost0]=fminsearch(@cost0,[0],options2,transfer); %using linear regression for k and abs val for p

%%Variables to pass out of function
transfer.maxTMH=max(TMH);
transfer.bestP0_units=bestP0*12/1e6*60*24*10000/SoilArea;

%make some plots if switch is on
if PLOTS == 1
 close all;

  figure;
%  axis([min(x) max(x) min(c13u+c12u)*.9 max(c13u+c12u)*1.1]);
  subplot(3,2,1), plot(x_mid,[c13u+c12u,c12u]);
  title('Equilibrium profiles');
  xlabel('x');
  ylabel('Methane (ppb)');
  legend('C12+C13','C12',2);
  
  %axis([min(xe) max(xe) min(d13Call)*.9 max(d13Call)*1.1]);
%  subplot(3,2,2), plot(xe,d13Call);
%  title('Equilibrium d13C profile');
%  xlabel('x');
%  ylabel('d13C');
  
  subplot(3,2,3), plot(t2, TMH,TGasSamples,GasSampleConc,'ro');
  title('Methane concentration in headspace');
  xlabel('Time (minutes)');
  ylabel('Headspace Methane (ppb)');
  

  subplot(3,2,4), plot(t2,d13CH,TGasSamples,GasSampled13C,'ro');
  title('d13C of Methane in Headspace');
  xlabel('Time (minutes)');
  ylabel('d13C of Headspace CH4');
  
  subplot(3,2,5),  plot(TMH,d13CH,GasSampleConc,GasSampled13C,'ro')
  title('Headspace Methane d13C of vs. Concentration')
  xlabel('Concentration, ppb')
  ylabel('d13C of CH4')
  
  subplot(3,2,6), plot(t2,SF6HS);
  title('Headspace SF6')
  xlabel('time, min')
  ylabel('[SF6]')
  
% t = transfer.TGasSamples0;
% n0 = transfer.GasSamplesLabelHSnM(1);        %C13 orig
% m0 = transfer.GasSamplesTMHnM(1);
% alpha = transfer.Fk;     %fractionization constant
% pC12_C13=transfer.pC12_C13;
% 
% 
% M = P/k + (m0 - P/k).*exp(-k.*t); %model = capital
% N = n0.*exp(-k.*alpha.*t) ;
%   
  
  
  
  figure
  bestP0H=1.05*bestP0;
  bestP0L=0.95*bestP0;
  TMHnM_hat=bestP0/k0-(bestP0/k0-TMHnM(round(t3)==0))*exp(-k0.*t3); 
  TMHnM_hatL=bestP0L/k0-(bestP0L/k0-TMHnM(round(t3)==0))*exp(-k0.*t3);
  TMHnM_hatH=bestP0H/k0-(bestP0H/k0-TMHnM(round(t3)==0))*exp(-k0.*t3);
  plot(TGasSamples,GasSamplesTMHnM, 'bo',t2,TMHnM,t2,TMHnM_hat,'r:',t2,TMHnM_hatL,'g:',t2,TMHnM_hatH,'g:')
  xlabel('time')
  ylabel('CH4 in headspace, nM')
  legend('Gas Samples' , 'All data' , 'curve fit' , 'P+5%' , 'P-5%',0)
  
  
  %figure
  %LabelHSnM_hat = LabelHSnM(t3==0)
  
  figure
  APMH_hat=  (LabelHSnM(t3==0)*exp(-k0.* Fk .* t3))./TMHnM_hat+pC12_C13;  
  APMH_hatL=  (LabelHSnM(t3==0)*exp(-k0.* Fk .* t3))./TMHnM_hatL+pC12_C13;  
  APMH_hatH=  (LabelHSnM(t3==0)*exp(-k0.* Fk .* t3))./TMHnM_hatH+pC12_C13;  
  plot(TGasSamples,GasSampleAPMH, 'bo',t2,APMH,t2,APMH_hat,'r:',t2,APMH_hatH,'g:',t2,APMH_hatL,'g:')
  xlabel('time')
  ylabel('Atom Pecent 13C-CH4 in headspace')
  legend('Gas Samples' , 'All data' , 'curve fit',0)
  
  figure
  plot(log(GasSamplesTMHnM),100.*GasSampleAPMH,'bo'...
      ,log(TMHnM),100.*APMH,...
      log(TMHnM_hat),100.*APMH_hat,'r:'...
      ,log(TMHnM_hatL),100.*APMH_hatL,'g:'...
      ,log(TMHnM_hatH),100.*APMH_hatH,'g:')
  xlabel('log CH4 in headspace, nM')
  ylabel('Atom Pecent 13C-CH4 in headspace')
  legend('Gas Samples' , 'All data' , 'curve fit' , 'P+5%' , 'P-5%',0)
  
end


%----------------------------------------
function dcdt = cdotinit(t,c,D,P,K,vol,csurf);

%[jk,c13e] = ode15s(@cdotinit,tspan0,c130,options,DS,PS,KS,Vol,c13u);
% function dcdt = cdotinit(t,c,D,P,K,vol,csurf);
%
% solve equilibration problem, 
% which is limited to the soil column
% vol c' = (D*c_x)_x + vol*P - vol * K c
% with boundary conditions 
% c(x=0,t) = csurf and
% c_x(x=Vs,t) = 0
% boundary interface @ x = 0 
% D is diffusion
% P is a source term
% K is consumption term 
% vol = vol of cell


global XS XSMID;
x = XS(:);
xm = XSMID(:);
dx = diff(x);
dx_mid = diff(xm);

dx_mid = [dx_mid(1) ; dx_mid(:) ; dx_mid(end)];      % face quantity
aug_conc = [csurf ; c(:) ; c(end)];   % builds in boundary conditions
                                     % fix surface dirichlet cond csurf
                                     % and repeat last entries
                                     % to make c(i+1) - c(i) = 0 -> c_x =0  
% now find cdot

flux_face = D.*(diff(aug_conc)./dx_mid);          % face quantity
dcdt = 1./vol.*diff(flux_face)./dx - K.*c + P;    % back to cell quantity

dcdt = dcdt(:);
%----------------------------------------

function dcdt = cdot(t,c,D,P,K,vol);

%function dcdt = cdot(t,c,D,P,K,vol);

% vol c' = (D*c_x)_x + vol*P - vol * K c
% with boundary conditions 
% c_x(x=-Vh,t) = c_x(x=Vs,t) = 0
% and P and K only active for x >0
% boundary interface @ x = 0 
% D is diffusion
% P is a source term
% K is consumption term 
% vol = vol of cell

global XXX XMID;
x = XXX(:);
xm = XMID(:);
dx = diff(x);
dx_mid = diff(xm);

dx_mid = [dx_mid(1) ; dx_mid(:) ; dx_mid(end)];      % face quantity
aug_conc = [c(1) ; c(:) ; c(end)];   % builds in boundary conditions
                                     % by repeated first and last entries
                                     % I make c(i+1) - c(i) = 0 -> c_x =0  
% now find cdot

flux_face = D.*(diff(aug_conc)./dx_mid);             % face quantity
dcdt = 1./vol.*diff(flux_face)./dx - K.*c + P;    % back to cell quantity

dcdt = dcdt(:);
%----------------------------

function errlr = cost0(param,transfer); %linear regression error function for k & abs val for p


t3=transfer.t3;
P = param(1);
k = transfer.k0;
t = transfer.TGasSamples0;
n0 = transfer.LabelHSnM(t3==0);        %C13 orig
m0 = transfer.TMHnM(t3==0);
alpha = transfer.Fk;     %fractionization constant
pC12_C13=transfer.pC12_C13;


M = P/k + (m0 - P/k).*exp(-k.*t); %model = capital
N = n0.*exp(-k.*alpha.*t) ;

AP=N./M+pC12_C13;
AP=AP(:);

%AP =  N ./ (M + N) * 100;
M = M(:);
N = N(:);
m = transfer.GasSamplesTMHnM;      %lower case reflects measured data
n = transfer.GasSamplesLabelHSnM;  
m=m(:);
n=n(:);
ap=n./m+pC12_C13;

errlr = 20*sum(abs(M-m)./m) + 100*sum(abs(AP - ap)./ap);



%-----------------------------
% function l1error = cost1(transfer);
% 
% 
% P = param(1);
% k = param(2);
% t = extra.t;
% data.m = extra.m.data;   %C12 data
% m0 = extra.m.m0;        %C12 orig
% data.n = extra.n.data;   %C13 data
% data.n = extra.n.data;   %C13 data
% n0 = extra.n.n0;        %C13 orig
% alpha = extra.alpha;     %fractionization constant
% 
% M = P/k + (m0 - P/k).*exp(-k*t);
% N = n0*exp(-k*alpha*t) ;
% %AP =  N ./ (M + N) * 100;
% M = M(:);
% N = N(:);
% m = data.m(:);
% n = data.n(:);
% 
% l1error = 20*sum(abs(M-m)) + 100*sum(abs(N - n));
