clear; close all; clc;

%% Part a: Reproducing T-xy diagram CO2 + H2S
kij = 0;

P = 40;    % pressure system, atm
T = [280.15 281.15 282.15 283.15 288.15 293.15 298.15 303.15 308.15...
313.15 318.15 323.15];
X = [0.916 0.852 0.797 0.756 0.563 0.42 0.313 0.232 0.168 0.117 0.074...
0.035];
Yexp = fliplr([0.124 0.216 0.3 0.382 0.467 0.546 0.63 0.72 0.835 0.862...
0.899 0.942]);

options = optimoptions('fsolve', 'Algorithm', 'levenberg-marquardt',...
'Display', 'off');

TempCalc = 280:2:328;

MatrixSolutionA = zeros(length(TempCalc), 2);

z0 = [0.1 0.8];

for ii=1:length(TempCalc)
t = TempCalc(ii);
if ii > 1
   z0 = solucion;
end
fun = @(z) VLE(z, P, t, kij);
solucion = fsolve(fun, z0, options);
MatrixSolutionA(ii,:) = solucion;
end

compXA = MatrixSolutionA(:,1);
compYA = MatrixSolutionA(:,2);

%%
%
figure(1)
plot(X, T, 'ob')
hold on
grid on
plot(Yexp, T, '^k')
plot(compXA, TempCalc, 'LineWidth', 3)
plot(compYA, TempCalc, 'LineWidth', 3)
set(gca, 'FontSize', 18)
xlabel('Molar fraction, adim')
ylabel('Temperature, K')
title('Pressure 40 atm, kij = 0')
legend({'Liquiq - Exp','Vapor - Exp', 'Liquid - PR', 'Vapor - PR'},...
'Location','best')
print(gcf,'figTxy_P40_kij0.png','-dpng','-r300');

%% Part b: Reproducing T-xy diagram CO2 + H2S

kij = 0.097;
MatrixSolutionB = zeros(length(X), 2);
z0 = [0.1 0.8];
for ii=1:length(TempCalc)
t = TempCalc(ii);
if ii > 1
   z0 = solucion;
end
fun = @(z) VLE(z, P, t, kij);
solucion = fsolve(fun, z0, options);
MatrixSolutionB(ii,:) = solucion;
end

compXB = MatrixSolutionB(:,1);
compYB = MatrixSolutionB(:,2);

%%
%
figure(2)
plot(X, T, 'ob')
hold on
grid on
plot(Yexp, T, '^k')
plot(compXB, TempCalc, 'LineWidth', 3)
plot(compYB, TempCalc, 'LineWidth', 3)
set(gca, 'FontSize', 18)
xlabel('Molar fraction, adim')

Problem 1 (150 points): Manufacturing of parts from highly viscous...