% TERMO.M %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute vapour composition using Wilson's model and extended Raoult's law

function [y,dT] = termo(x,T,flag,Nt)

global termodata

A = termodata.A;
V = termodata.V;
W = termodata.W;
R = termodata.R;
P = termodata.P;

Nc = 4;

x = x';
x = [
	x
	1 - sum(x)
	];

Ps = zeros(1,Nc);
ratio = zeros(1,Nc);
LAMBDA = zeros(Nc,Nc);
LNGAMMA = zeros(1,Nc);

switch flag
    case 1
        y = zeros(Nt,Nc);
        dT = zeros(Nt,1);
        
        for i = 1:Nt
            % i: Stage
            for j = 1:Nc
                % j: Component
                
                % Compute vapour pressure using Antoine equation
                Ps(j) = 10^(A(j,1) - A(j,2) / (T(i) - 273.15 + A(j,3)));
                % Wilson model
                for k = 1:Nc
                    % k: Component
                    if j == k
                        LAMBDA(j,k) = 1;
                    else
                        LAMBDA(j,k) = (V(k) / V(j)) * exp(-W(j,k) / (R * T(i)));
                    end
                end
            end
            
            for j = 1:Nc
            	% j: Component
            	for k = 1:Nc
            		% k: Component
            		ratio(k) = LAMBDA(k,j) / (x(:,i)' * LAMBDA(k,:)');
                end
            	LNGAMMA(j) = -log(x(:,i)' * LAMBDA(j,:)') + 1 - x(:,i)' * ratio';
            end
            GAMMA = exp(LNGAMMA)';
            % Using VLE formula to compute vapour mole fraction
            for j = 1:Nc
            	% j: Component
            	y(i,j) = Ps(j) * x(j,i) * GAMMA(j) / (x(:,i)' .* Ps * GAMMA);
            end
            % Update law for temperature
            dT(i) = 1 * (P - x(:,i)' .* Ps * GAMMA ./ 760);
        end
        y = y(:,1:(Nc - 1));
    case 3
    	y = zeros(1,Nc);
    	
        for i = 1:Nc
        	% i: Component
        	Ps(i) = 10^(A(i,1) - A(i,2) / (T(1) - 273.15 + A(i,3)));
        end
        for i = 1:Nc
        	% i: Component
        	for j = 1:Nc
        		% j: Component
        		if (i == j)
        			LAMBDA(i,j) = 1;
        		else
        			LAMBDA(i,j) = (V(j) / V(i)) * exp(-W(i,j) / (R * T(1)));
        		end
        	end
        end
        for i = 1:Nc
        	% i: Component
        	for j = 1:Nc
        		% j: Component
        		ratio(j) = LAMBDA(j,i) / (x' * LAMBDA(j,:)');
            end
        	LNGAMMA(i) = -log(x' * LAMBDA(i,:)') + 1 - x'*ratio';
        end
        GAMMA = exp(LNGAMMA)';
        for i = 1:Nc
        	% i: Component
        	y(i) = Ps(i) * x(i) * GAMMA(i) / (x' .* Ps * GAMMA);
        end
        y = y(1:(Nc - 1));
        % Update law for temperature
        dT = 1 * (P - x' .* Ps * GAMMA./760);
end
