function y=f(x)

global pH

%---------------
% Ca    -> x(1)
% SO4   -> x(2)
% CO3   -> x(3)
% H2CO3 -> x(4)
% HCO3  -> x(5)
% H     -> x(6)
% HSO4  -> x(7)
%---------------

Ksp1 = 10^-5.721;
Ksp2 = 10^-8.47;
K    = 10^1.47;
Ka1  = 10^-6.531;
Ka2  = 10^-10.35;
Ka3  = 10^-1.92;

Po2    = 10^-1.5;
SO4tot = 1e-4;

y = zeros(size(x));

y(1) = Ksp1 - x(1)*(x(2)+SO4tot);
y(2) = Ksp2 - x(1)*x(3);
y(3) = x(4) - Po2/K;
y(4) = Ka1*x(4) - x(5)*x(6);
y(5) = Ka2*x(5) - x(3)*x(6);
y(6) = Ka3*x(7) - x(6)*(x(2)+SO4tot);
y(7) = x(6) - 10^-pH;

