kill(all);

maxima:Construct Txy diagram Chloroform & hexane (= Mixture),version 1.0,part 2
peter.vlasschaert@gmail.com,30/08/2020

see website :' linspace(f),combspace '
------------------------------------------------------------------------------------------------------------------
maxima:Construct Txy diagram Chloroform & hexane (= Mixture),version 1.0,part 1
------------------------------------------------------------------------------------------------------------------
→ chloroform  (1)
→ hexane        (2)

load packages ,option
=================

load(lrats)$ ratprint:false$

physical parameters
================

total  pressure (atm)
----------------------------

P:1;  /* atm */

parameters  'Wilson Equation '
------------------------------------------

G12:1.2042 $ G21:0.39799 $

vapour pressure :(Antoine Equation : P(atm), T(kelvin))
---------------------------------------------------------------------------
ts : 'Kelvin'

p1sat(ts):= exp(9.33984-2696.79/(ts-46.14))$/*chloroform*/
p2sat(ts):= exp(9.20324-2697.55/(ts-48.78))$/*hexane    */

Estimate : boilling points
----------------------------------

Tb1:2696.79/(9.33984-log(1))+46.14;/*chloroform */
Tb2:2697.55/(9.20324-log(1))+48.78;/*hexane     */

341.88-273.15; /* temperature (unit °C) :hexane */

x-axis : range :  x ,y  mole fraction 
----------------------------------------------

p1:[0,0.002,0.004,0.006,0.008,0.01,0.015,0.02,0.04,0.06,0.08,0.1,0.12,0.14
 ,0.16,0.18,0.2,0.22,0.24,0.26,0.28,0.3,0.32,0.34,0.36,0.38,0.4,0.42,0.44,
 0.46,0.48,0.5,0.52,0.54,0.56,0.5800000000000001,0.6,0.62,0.64,0.66,0.68,
 0.7000000000000001,0.7200000000000001,0.74,0.76,0.78,0.8,0.8200000000000001,
 0.8400000000000001,0.86,0.88,0.9,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.985,
 0.99,0.995,1]$

p1a:copy(p1);

np:length(p1);

xp:p1$   /* liquid */
yp:p1$   /* vapour */
Tp:makelist(0,i,1,np)$   /*temperature(kelvin)*/

1e)  Tb2 > Tb1
============
set : values 

Tp[1]:Tb2$ Tp[np]:Tb1$

2e)  start value : simulation (dT = .01, T=Tb2)
===================================
set : values 

dT:0.001$  /* temperature step */
 /* kelvin */

load packages : "Newton Raphson", from Maxima,share.

load(mnewton)$

for i:2 thru np do(
      x1:xp[i],x2:1-x1,
/* Evaluate activity coefficients                */
term1:x1+x2*G12,term2:x2+x1*G21,
γ_1:exp(-log(term1)+ x2*(G12/term1-G21/term2)),
γ_2:exp(-log(term2)+ x1*(G21/term2-G12/term1)),
/* ********************************************* */        
/* find : f(T) = 0                               */   
    f:x1*γ_1*ev(p1sat(T))+x2*γ_2*ev(p2sat(T))-P ,     
    sol:mnewton([f],[T],[Tb2+(i-1)*dT]),
/* ********************************************* */          
/*  output →                                     */
    Tp[i]:rhs(flatten(sol)[1]),
    yp[i]:x1*γ_1*ev(p1sat(T),T=Tp[i])
/* ********************************************* */     
        )$   

use : flatten '   s :[[]]  = flatten remove the inner ones  s :[ ]

Tp;  /* temperature unit '= Kelvin ' */

yp; /* vapour mole fraction, chloroform */ 

p1a; /* liquid mole fraction, chloroform */

p2:lmin(Tp);  /* temperature for azeotrope in Kelvin ' chloroform , hexane '*/

p3:block(for i:1 thru np do(
if Tp[i] = p2 then return([i,yp[i]])
));

p3a:p3[1]; /* number in list 'p2' */
p3b:p3[2]; /* vapour mole fraction 'yp'*/

p1a[44];   /* liquid mole fraction 'p1a'*/

What we proved  x[i] ' = p1a '  =   y[i] ' = p3b ',  => azeotrope at this point 'chloroform,hexane'
see website
new: 25/08/2020:maxima : prove :result Azeotrope  x[i]=y[i], VLE (binary),version 1.0