assignment.m

% Definitions of values
P  = 760;
R = 1.9872;
A1 = 8.07131;
B1 = 1730.630;
C1 = 233.426;
A2 = 7.04115;
B2 = 1373.799;
C2 = 214.979;
alpha12 = 0.6092;
alpha21 = 0.6092;
A12 = 1999.3885;
A21 = 328.3587;

X = [];
T = [];
Y = [];

T1sat = B1 / (A1 - log10(P)) - C1;
T2sat = B2 / (A2 - log10(P)) - C2;

for x1=0:0.01:1
    X = [X x1];
    x2 = 1-x1;
    t = x1 * T1sat + x2 * T2sat;
    P1sat = power(10, A1 - B1 / (t + C1));

    er = 1;
    while er > 0.01
        P2sat = power(10, A2 - B2 / (t + C2));
        alpha = P1sat / P2sat;
        
        t12 = A12 / (R * (t + 273.15));
        t21 = A21 / (R * (t + 273.15));
        G12 = exp(-alpha12 * t12);
        G21 = exp(-alpha21 * t21);
        coeff1 = t21 * power(G21 / (x1 + x2*G21), 2) + t12*G12 / power(x2 + x1*G12, 2);
        coeff2 = t12 * power(G12 / (x2 + x1*G12), 2) + t21*G21 / power(x1 + x2*G21, 2);

        gamma1 = exp(coeff1 * x2 * x2);
        gamma2 = exp(coeff2 * x1 * x1);

        P1sat = P / (gamma1 * x1 + gamma2 * x2 / alpha);
        y = P1sat * gamma1 * x1 / P;
        t_prev = t;
        t = B1 / (A1 - log10(P1sat)) - C1;
        er = (t - t_prev) / t_prev * 100;
    end

    T = [T t];
    Y = [Y y];
end

T_exp = [115.5, 111.1, 109.3, 107.9, 107.7, 106.9, 105.6, 104.7, 103.7, 103.8, 101.7, 101.0, 99.7, 99.5, 97.5, 97.4, 96.2, 95.5, 95.4, 95.0, 94.7, 94.7, 94.6, 94.4, 94.5, 94.5, 94.8, 94.9, 94.9, 95.1, 95.3, 95.5, 96.0, 96.3, 96.6];
X_exp = [0.0, 0.0490, 0.0756, 0.0936, 0.1067, 0.1097, 0.1366, 0.1525, 0.1664, 0.1740, 0.1965, 0.2402, 0.2671, 0.2743, 0.3512, 0.3910, 0.4633, 0.5461, 0.5544, 0.6016, 0.6502, 0.6622, 0.6903, 0.7761, 0.8016, 0.8390, 0.8703, 0.8896, 0.8961, 0.9162, 0.9503, 0.9629, 0.9707, 0.9809, 0.9824];
Y_exp = [0.0, 0.1740, 0.2463, 0.2789, 0.3019, 0.3057, 0.3606, 0.3808, 0.4034, 0.4110, 0.4429, 0.4869, 0.5130, 0.5136, 0.5757, 0.6016, 0.6407, 0.6809, 0.6860, 0.7060, 0.7244, 0.7291, 0.7378, 0.7625, 0.7691, 0.7766, 0.7829, 0.7879, 0.7898, 0.7950, 0.8107, 0.8217, 0.8348, 0.8592, 0.8607];

 plot(T, X, T, Y, T_exp, X_exp, 'k*', T_exp, Y_exp,'b.');
 legend('X1', 'Y1', 'X1 Experimental', 'Y1 Experimental')

% plot(Y, X, Y_exp, X_exp, 'k.');
% legend('Calculated', 'Experimental');