clear all;      close all;

global  R18w   R18OH

%__________________________________________________________________________
%              Input Initial (Experimental) Conditions

Data_Initial = 4;
            % Select 
            %    1: TS2 data & initial condition
            %    2: TS2-mnCA data & initial condition              
            %    3: TS2-CA2 data & initial condition
            %    4: TS2-3CA data & initial condition
            %    5: SS18 25C experiment data & initial condition 
Data_Entry;                                   
                            
                            
% Flag for enzymatic catalysis 
CA_catalysis = 0;   
      % Choose 0 for Michaelis-Menton model to account for CA catalysis
      % on CO2 hydration (based on Uchikawa & Zeebe 2012)
      
      %  or
      
      % Choose 1 to arbitrarily amplify k1f (CO2 hydration kinetic rate
      % constant (e.g., Chen et al.2018; Guo 2020) to fit the model to 
      % isotope data with no constraints on CA concentration such as 
      % biological systems. For this case, manually adjust the value of
      % "cat_factor" paramater, which can be found at Line 224 in
      % "Model_Constant.m"
                                     
% Temp. conversion from Celsius to Kelvin
Tk = Tc + 273.15;        


%__________________________________________________________________________
%                    Adjustable Parameters
KIF13_CO2H2O = 1.0176;    % 13C KIE for CO2 Hydration 
                          % (Yumol et al., 2020)  
KIF13_CO2OH = 1.0170;     % 13C KIE for CO2 Hydroxylation 
                          % (Christensen et al., 2020)

% Oxygen isotope KIFs  (Guo & Zhou, 2019, Table 3)
% 266 + 8 -> 2866 + H      
    F = ((1.5992E6/Tk^2)+(2.0211E3/Tk)-16.1);  % 8.7 @ 25C  
    KIF18_a1 = exp(-1*F/1000);                 % 0.9913 @ 25C
% 286 + 6 -> 2866 + H
    F = ((7.5198E5/Tk^2)-(4.050E3/Tk)+0.8);    % -4.3 @ 25C    
    KIF18_b1 = exp(-1*F/1000);                 % 1.0043 @ 25C
% 266 + '8 -> 2866
    F = ((-7.681E5/Tk^2)+(3.8436E4/Tk)-124.6); % -4.3257 @ 25C
    KIF18_a4 = exp(-1*F/1000);                 % 1.0043 @ 25C
% 286 + '6 -> 2866
    F = ((7.7752E5/Tk^2)-(1.8330E3/Tk)+0.9);   % 3.499 @ 25C
    KIF18_b4 = exp(-1*F/1000);                 % 0.9965 @ 25C    

% Clumped KIEs         (Guo, 2020, Table 2)
% 366 + 8 -> 3866 + H       0.999853 @ 25C
     KIE47_p1 = 1 + ( (4613.8393/Tk^2)-(4.0389/Tk)-0.185 )/1000;                                     
% 386 + 6 -> 3866 + H       0.999781 @ 25C        
     KIE47_s1 = 1 + ( (-5705.688/Tk^2)-(41.5925/Tk)-0.015 )/1000;     
% 366 + '8 -> 3866          0.999984 @ 25C
     KIE47_p4 = 1 + ( (-902.7635/Tk^2)+(157.1718/Tk)-0.533 )/1000;
% 386 + '6 -> 3866          0.999825 @ 25C
     KIE47_s4 = 1 + ( (-11771.2832/Tk^2)-(62.7060/Tk)+0.168 )/1000;   
   

%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



% Loading required constants
Model_Constants;

% d13C & d18O of CO2(aq)
d13C_co2d = ((d13C_SB + 1000) / Alp13_DCPd) - 1000;  % VPDB
  %d13C_co2d = d13C_SB - Eps13_DCPd;   % Sang's approach
d18O_co2d = ((d18O_dcp_t0 + 1000) / Alp18_DCPd) - 1000;  % VSMOW   

% 18O content in H2O & OH-.
R18w = ( (d18O_w / 1000)+1 ) * R18vsmow;
R18OH = Alp18_OHw * R18w;


% blank matrix into which initial condition will be loaded
C0 = zeros (1,8);
% giving numerical strings for the matrix for the initial condition
    co2d    = 1;        % [CO2(aq)]
    DCP     = 2;        % [DCP]  =[HCO3+CO3]
    co2d_13 = 3;        % 13R (13C/12C) of CO2(aq)
    DCP_13  = 4;        % 13R of DCP
    co2d_18 = 5;        % 18R (18O/16O) of CO2(aq)
    DCP_18  = 6;        % 18R of DCP
    co2d_47 = 7;        % 13C-18O clumps in CO2(aq)                                                   
    DCP_47  = 8;        % 13C-18O clumps in DCP

% initial condition/value for the variables    
C0(co2d) = (NaHCO3 * x1 * H) / K1;       
C0(DCP) = NaHCO3;                        
C0(co2d_13) = C0(co2d) * ((d13C_co2d / 1000) + 1) * R13vpdb;                           
C0(DCP_13) = NaHCO3 * ((d13C_SB / 1000) + 1) * R13vpdb; 
C0(co2d_18) = 2 * C0(co2d) * ((d18O_co2d / 1000) +1) * R18vsmow;               
C0(DCP_18) = 3 * NaHCO3 * ((d18O_1st / 1000) +1) * R18vpdb;                             
C0(co2d_47) = 2 * C0(co2d) * ((d13C_co2d / 1000) + 1) * R13vpdb * ...          
    ((d18O_co2d / 1000) +1) * R18vsmow * ((D47_co2 / 1000)+1);             
C0(DCP_47) = 3 * NaHCO3 * ((d13C_SB / 1000) + 1) * R13vpdb * ...                                 
    ((d18O_1st / 1000) +1) * R18vpdb * ((D47_DCP / 1000)+1);


TIME = [T0: 30:  600000];
options = odeset('MaxStep',40); 
         % lower the numner for faster equilibration ~100 for CAmn, CA2
         % 
[t,y] = ode15s(@IsoEquilibration, TIME, C0, options);

TIMEm = t/60;
co2d = y(:, co2d);
DCP = y(:, DCP);
co2d_13 = y(:, co2d_13);
DCP_13 = y(:, DCP_13);
co2d_18 = y(:, co2d_18);
DCP_18 = y(:, DCP_18);
co2d_47 = y(:, co2d_47);
DCP_47 = y(:, DCP_47);


d13C_DCP = ((DCP_13./DCP)./R13vpdb - 1).*1000;        
d18O_DCP_pdb = ((DCP_18./DCP)./R18vpdb/3 - 1).*1000;    % following JW
d18O_DCP_smow = ((DCP_18./DCP)./R18vsmow/3 - 1).*1000;  % following JW
      Eps18_DCP_w = ((d18O_DCP_smow+1000)./(d18O_w+1000) - 1).*1000;
Cap47_DCP = ((DCP_47.*DCP)./(DCP_13.*DCP_18) - 1).*1000;
Cap47_co2 = ((co2d_47.*co2d)./(co2d_13.*co2d_18) - 1).*1000;



figure(1)

subplot(2,1,1)
hold on; box on; set(gca,'FontWeight','bold');
plot(TIMEm, Eps18_DCP_w,'k-','LineWidth',1.5);
set(gca,'xscale','log');
xlabel('Time ( min )'); ylabel('^{18}\epsilon_{DCP-H2O} ( %^{_o} )');

subplot(2,1,2)
hold on; box on; set(gca,'FontWeight','bold');
plot(TIMEm, Cap47_DCP,'k-','LineWidth',1.5);
set(gca,'xscale','log');
xlabel('Time ( min )'); ylabel('\Delta_{47} of DCP ( %^{_o} )');



T = Data(:,2);
d18Opdb = Data(:,5);
   d18Osmow = 1.03091*(d18Opdb)+30.91;
Eps18_BaCO3_w = ((d18Osmow+1000)./(d18O_w+1000) - 1).*1000;
D47_BaCO3 = Data(:,9);
D47_error = Data(:,10);

figure(1)
subplot(2,1,1)
plot(T, Eps18_BaCO3_w,'ko','MarkerFaceColor','c');

subplot(2,1,2)
errorbar(T,D47_BaCO3,D47_error,'ko','MarkerFaceColor','c');