REQUIRE "atoms.a4l";

(* problem 10.3, pp 359 VanNess, 7th Edition *)

MODEL all;

T IS_A temperature;
P, P1S, P2S IS_A pressure;

A IS_A factor;
gamma1, gamma2 IS_A positive_factor;
A1, B1, C1 IS_A positive_factor;
A2, B2, C2 IS_A positive_factor;      
x1, x2, y1, y2 IS_A fraction;

alpha IS_A positive_factor;

vp1: P1S/1000.0{Pa} = exp(A1 - B1/(T/1{K} - C1));
vp2: P2S/1000.0{Pa} = exp(A2 - B2/(T/1{K} - C2));

sumx: x1 + x2 = 1.0;
sumy: y1 + y2 = 1.0;

vle1: x1*gamma1*P1S = y1*P;
vle2: x2*gamma2*P2S = y2*P;

activity1: gamma1 = exp(A*x2*x2);
activity2: gamma2 = exp(A*x1*x1);

eqA: A = 2.771 - 0.00523*T/1.0{K};

eqforalpha: alpha*x1*y2 = y1*x2;

METHODS
        METHOD clear;
        T.fixed := FALSE;
        P.fixed := FALSE;
        P1S.fixed := FALSE;
        P2S.fixed := FALSE;
        A1.fixed := FALSE;
        B1.fixed := FALSE;
        C1.fixed := FALSE;
        A2.fixed := FALSE;
        B2.fixed := FALSE;
        C2.fixed := FALSE;
        x1.fixed := FALSE;
        x2.fixed := FALSE;
        y1.fixed := FALSE;
        y2.fixed := FALSE;
        A.fixed := FALSE;
        gamma1.fixed := FALSE;
        gamma2.fixed := FALSE;
        alpha.fixed := FALSE;
        END clear;

        METHOD antoine;
        A1.fixed := TRUE;
        A1 := 16.59158;
        B1.fixed := TRUE;
        B1 := 3643.31;
        C1.fixed := TRUE;
        C1 := 33.424;
        A2.fixed := TRUE;
        A2 := 14.25326;
        B2.fixed := TRUE;
        B2 := 2665.54;
        C2.fixed := TRUE;
        C2 := 53.424;
        END antoine;

        METHOD parta;
        T.fixed := TRUE;
        T := 318.15 {K};
        x1.fixed := TRUE;
        x1 := 0.25;
        END parta;

        METHOD partb;
        T.fixed := TRUE;
        T := 318.5 {K};
        y1.fixed := TRUE;
        y1 := 0.6;
        END partb;

        METHOD partc;
        P.fixed := TRUE;
        P := 1.0{atm};
        x1.fixed := TRUE;
        x1 := 0.85;
        END partc;
 
        METHOD partd;
        P.fixed := TRUE;
        P := 1.0{atm};
        y1.fixed := TRUE;
        y1 := 0.4;
        END partd;

        METHOD parte;
        alpha.fixed := TRUE;
        alpha := 1.0;
        T.fixed := TRUE;
        T := 318.15 {K};
        END parte;

        METHOD values;
        x1 := 0.5;
        y1 := 0.5;
        T := 320 {K};
        P := 1.0 {atm};
        gamma1 := 1.0;
        gamma2 := 1.0;
        END values;

        METHOD doparta;
        RUN clear;
        RUN antoine;
        RUN parta;
        END doparta;

        METHOD dopartb;
        RUN clear;
        RUN antoine;
        RUN partb;
        END dopartb;

        METHOD dopartc;
        RUN clear;
        RUN antoine;
        RUN partc;
        END dopartc;

        METHOD dopartd;
        RUN clear;
        RUN antoine;
        RUN partd;
        END dopartd;

        METHOD doparte;
        RUN clear;
        RUN antoine;
        RUN parte;
        END doparte;
END all;