%Equilibrium calculation at given temperature, pressure and total
%composition in ideal mixture phases. The syntax is
%n = minGibb(n0,t,p,m0,A) where: 
%
%         n0 = initial mole number vector [mol]
%         t  = temperature [K]
%         p  = pressure [Pa]
%         m0 = standard chemical potential (ideal gas at T)
%         A  = formula matrix (atoms x substances)
%         n  = equilibrium mole number vector [mol]
%
%Copyright Tore Haug-Warberg 2003 (course TMT4140, NTNU)
%
 function[n] = minGibb(n0,t,p,m0,A)
%
%Parameters and physical quantities
 r       = 8.3143;                       % universal gas constant
 p0      = 1e5;                         % standard state pressure
 tol     = 1e-6;                          % convergence tolerance
 maxiter = 100;                    % maximum number of iterations
 maxstep = 0.9;                    % maximum relative step length
 epx     = 1e-6;                % arbitrarily small mole fraction
 [nc,ns] = size(A);               % number of atoms and substance
%
%Initialization
 n       = n0;
 i       = find(n==0);
 n(i)    = zeros(length(i),1) + sum(n)*epx;
 e       = ones(ns,1);
%
%Equilibrium calculation
 for iter=1:maxiter
   m    = m0 + r*t*log(p/p0*n/sum(n));
   D    = diag(n)/(r*t);
   lhs  = [A*D*A', A*n; (A*n)', 0];
   rhs  = [A*D*m-A*(n-n0); n'*m];
   x    = inv(lhs)*rhs;
   dn   =-D*(m-A'*x(1:nc)) + n*x(nc+1);
   step =-maxstep/min([-maxstep,min(dn./n)]);
   n    = n + step*dn;
   if max(abs(dn./n))<tol; return; end
 end
%
%Not converged
 n = NaN*n; return
