function N = nprfun(T, V, P, y)
% N = NPRFUN(T, V, P, y)
%
% Beregner totalt antall mol for en gassblanding av O2, He og N2 med temperatur
% T [°K], volum V [m3], trykk P [Pa] og molfraksjoner y ved en Peng-Robinson
% tilstandsligning.
%
% Skrevet av Olaf Trygve Berglihn <olafb@pvv.org> 5. mai 2005
% Beskyttet av GPL v2 lisens <http://www.gnu.org/licenses/gpl.txt>

% gasskonstant
R = 8.31;

% Peng-Robinson parametre
% komponent-rekkefølge: O2, He, N2
Tc = [-118.38; -267.96; -146.96] + 273.15; % °K
Pc = [5080.02e3; 226.97e3; 3394.37e3];     % Pa
Vc = [0.07320; 0.05730; 0.09000];          % m3/kmol
w = [0.01900; -0.39003; 0.04000];
d = [0,        0,      -0.01200
     0,        0,       0.06850 
     -0.01200, 0.06850, 0];

% N = sum(n);
% y = n/N;
Tr = T./Tc;
yij = y*y';

% Peng-Robinson konstanter
b=0.07780*R*Tc./Pc;
k=0.37464 + 1.54226*w - 0.26992*w.^2;
alpha = (1 + k.*(1-sqrt(Tr))).^2;
ac = 0.45724*R^2*Tc.^2./Pc;
a = ac.*alpha;
aij = (1-d).*(a*a').^(0.5);

% Mikseregler
am = sum(sum(triu(yij).*triu(aij),1));
bm = sum(y.*b);

% løser ligning for N med Newtons algoritme
% bruker ideell gass som startverdi
N = P*V/(R*T);
reltol = 1e-5;
relres = reltol*10;
maxit = 100;
it = 0;

while (relres > reltol) & (it < maxit)
	f = P - ( N*R*T/(V-N*bm) - N^2*am/(V*(V+N*bm) + N*bm*(V-N*bm)));
	df = -( R*T*(V-N*bm) - N*R*T*(V-N*bm)*(-bm) ) / (V-N*bm)^2 ...
		+ ( 2*N*am*(V*(V+N*bm)+N*bm*(V-N*bm)) ...
		-   N^2*am*(V*bm+bm*(V-N*bm)+N*bm*(-bm)) ) ...
		/ (V*(V+N*bm) + N*bm*(V-N*bm))^2;
    Ni = N - f/df;
    relres = abs(Ni-N)/N;
    N = Ni; 
    it = it + 1;
end

