boiling_point (t0,x,AC,P)

function BoilingPoint = boiling_point (t0,x,AC,P)
% Finds boiling point for ideal liquids
%
% t0 - starting temperature
% x  - molar fractions
% AC - constants of Antoine equations
%       form:     N = number of compounds
%         [A1 B1 C1
%          A2 B2 C2
%          .  .  .
%          AN BN CN]
%
% P  - pressure

BoilingPointTolerance = 1e-6;

x = x(:);

if sum(x)~=1
    error('sum(x)~=1. function: boiling_point')
end

Ni   = length(x);   %number of compounds
p0   = zeros(1,Ni); %prealocation
t1   = t0;
iter = 0;     %counting iterations
MaxIter = 50; %maximum number of iterations

for i=1:Ni
    p0(i) = 10^(AC(i,1)-AC(i,2)/(AC(i,3)+t1));
end

criteria  = abs( P - p0*x );

while criteria > BoilingPointTolerance

    iter = iter +1;
    if iter == MaxIter
        error('MaxIter reached. function: boiling_point')
    end

    F = P - p0*x;

    %dFdt
    dFdt_i = zeros(1,Ni);
    for i=1:1:Ni
        dFdt_i(i) = x(i)*log(10)*AC(i,2)/((AC(i,3)+t1)^2)*(10^(AC(i,1)-AC(i,2)/(AC(i,3)+t1)));
    end
    dFdt = -sum(dFdt_i);

    t2 = t1 - F/dFdt;
    t1 = t2;

    p0 = zeros(1,Ni);
    for i=1:1:Ni
        p0(i)=10^(AC(i,1)-AC(i,2)/(AC(i,3)+t1));
    end

    criteria = abs(P-p0*x);
end

BoilingPoint = t1;