%
%	Calculate the pH of a Boric Acid Solution using the Newton-Raphson Method
%
%	Variables are:  x1 = B(OH)3
%			        x2 = B(OH)4 -
%			        x3 = H3O +
%			        x4 = OH -
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear
%
% Input of the equilibrium constants and total concentrations
%
Ka = 7e-10;           % equilibrium constant for boric acid
Kw = 1e-14;            % equilibrium constant for water
CT = 5e-4;            % total concentration of boric acid
%
% Initial Guess is given for each species concentrations
%
X = [CT;0;1e-7;1e-7];
%
% Give the charge vector containing the charge of each species
%
Z = [0;-1;1;-1];
%
% Use the Newton-Raphson technique to determine the speciation: loop to minimize
% each function F
%
s =1;
while s > 1e-18
%
% Calculate the ionic strength
%
	I=1/2*sum(Z.^2.*X);
%
% Calculate the activity coefficients with the Guntelberg formula
%
    g0 = 1;                                 % for a neutral species
	g1 = 10^(-0.5*sqrt(I)/(1+sqrt(I)));     % for a 1-charge species
%
% Correct equilibrium constants for ideality deviation
%
	Kca = Ka*g0/g1^2;
	Kwc = Kw/g1^2;
%
% Calculation of the difference functions
%
	f1 = -X(2)+X(3)-X(4);                   % electroneutrality
	f2 = X(1)+X(2)-CT;                      % mass balance
	f3 = X(2)*X(3)-Kca*X(1);                % boric acid equilibrium
	f4 = X(3)*X(4)-Kwc;                     % water equilibrium

	F = [f1;f2;f3;f4];
%
% Construction of the Jacobien for the Newton-Raphson minimization
%
	J =[0 -1 1 -1;
	    1 1 0 0;
	    -Kca X(3) X(2) 0;
	    0 0 X(4) X(3)];
%
% Computation of the minimum d
%
	d = inv(J)*F;
	X = X-d;
%
% Test the convergence
%
	s = norm(F);
end
%
% Calculate the pH and check mass balance
%
disp ('pH =')
disp(-log10(X(3)))
disp('B(OH)3 = ')
disp(X(1))
disp('B(OH)4- = ')
disp(X(2))

Boric_tot = X(1)+X(2)           % Check mass balance