# H2CO3 → H+ + HCO3-      4.5e-7
# HCO3- → H+ + CO3 2-     4.7e-11
# HSO4- → H+ + SO4 2-     1.0e-2
# H2O → H+ + OH-          1.0e-14

#variables: CO2, HCO3-, CO3 2-, H+, OH-, HSO4-, SO4 2-
#constants: H20

#"prod" denotes products, "reag" denotes reagents, "total" denotes the total reaction
#"var" denotes variable species; concentrations can change; "const" denotes constant species
#"rxn" denotes reactions matricies; define the parameters of the reactions
#"species" denotes actual concentration values.
#"terms" denotes terms in an equation
#"log" denotes the log base 10 of some matrix
# in all cases, each column is a species; and each row is a reaction.

#set up reactions.
varProdRxn = [0 1 0 1 0 0 0
              0 0 1 1 0 0 0
              0 0 0 1 0 0 1
              0 0 0 1 1 0 0];
varReagRxn = [1 0 0 0 0 0 0
              0 1 0 0 0 0 0
              0 0 0 0 0 1 0
              0 0 0 0 0 0 0];
constProdRxn = [0
                0
                0
                0];
constReagRxn = [0
                0
                0
                1];
prodRxn = [varProdRxn, constProdRxn];
reagRxn = [varReagRxn, constReagRxn];
totalRxn = prodRxn - reagRxn;
totalVarRxn = varProdRxn - varReagRxn;

k = [4.5e-7; 4.7e-11; 1.0e-2; 1.0e-14];
logk = log10(k);

nRxn = rows(varProdRxn);
nSpecies = columns(varProdRxn) + columns(constProdRxn);
rxnSize = size(nRxn, nSpecies);

#set up initial conditions
#note that no starting concentration can be 0, otherwise the program will get mad at you.
varSpecies = [1e-99 1e-99 1e-99 1e-99 1e-99 1 1e-99]; #CHANGE THIS
constSpecies = [1]; #SET/CALCULATE THIS
species = [varSpecies, constSpecies];

#simulation parameters
tolerance = 0.001; #maximum allowable error (percent/100) between k and q.

#slowly add carbonate to bisulfate
finalConcentration = 2;
step = 0.01;
steps = finalConcentration/step;
speciesArr = [];
t=step:step:finalConcentration;

#END OF SETUP. ACTUAL CALCULATIONS FOLLOW THIS LINE.

for i = 1:steps
  while true
    q = 10 .^ (totalRxn * transpose(log10(species)));

    error = (k - q) ./ k;
    if sum(abs(error) > tolerance) == 0 break end

    logProdSpecies = prodRxn * diag(log10(species));
    logProdCoeff = zeros(rxnSize);
    logReagSpecies = reagRxn * diag(log10(species));
    logReagCoeff = zeros(rxnSize);
    for rxnN = 1:nRxn
      for speciesN = 1:nSpecies
        logProdCoeff(rxnN, speciesN) = sum([logProdSpecies(rxnN,1:speciesN-1)...
                                         logProdSpecies(rxnN, speciesN+1:end)]);
        logReagCoeff(rxnN, speciesN) = sum([logReagSpecies(rxnN,1:speciesN-1)...
                                         logReagSpecies(rxnN, speciesN+1:end)]);
      end
    end

    prodCoeff = 10 .^ logProdCoeff;
    prodTerms = prodCoeff .* [varProdRxn, zeros(size(constProdRxn))];
    reagCoeff = 10 .^ logReagCoeff;
    reagTerms = reagCoeff .* [varReagRxn, zeros(size(constReagRxn))];

    linEq = zeros(nRxn, nRxn + 1);
    linEq(:, end) = (k - q) .* (10 .^ (reagRxn * transpose(log10(species))));
    for rxnN = 1:nRxn
      linEq(rxnN, 1:end-1) = prodTerms(rxnN,:) * transpose(totalRxn) - ...
                            (k(rxnN) .* (reagTerms(rxnN,:) * transpose(totalRxn)));
    end
    linEq = rref(linEq);

    dVarSpecies = transpose(linEq(:, end)) * totalVarRxn;
    if sum(varSpecies + dVarSpecies <= 0) > 0 #if the change will result in a negative concentration
      dVarSpecies = dVarSpecies ./ max(-dVarSpecies ./ varSpecies) + 1e-20;
      #scale the change to get as close to zero as possible
    end
    varSpecies = varSpecies + dVarSpecies;
    species = [varSpecies, constSpecies];
  end


  speciesArr = [speciesArr; species];

  varSpecies(3) += step;
  if(varSpecies(1) > 0.1) varSpecies(1) = 0.1; end
  species = [varSpecies, constSpecies];

end

co2 = speciesArr(:,1);
bicarb = speciesArr(:,2);
carb = speciesArr(:,3);
H = speciesArr(:,4);
pH = -log10(H);
OH = speciesArr(:,5);
bisulf = speciesArr(:,6);
sulf = speciesArr(:,7);