kill(all)$

analytical method : binary distillation
peter.vlasschaert@gmail.com,05/04/2017
part2:stripping-section 'operating line'

alpha :constant relative volatility

assume(alpha >= 1);

p1:hyperbole       : equilbrium equation.
p2:linear equation : operating line.

p1:y=(alpha*x)/(1+(alpha-1)*x);
p2:y=m*x+b;

removes (p1,p2),solve for x

p3:ratsimp(rhs(p1)-rhs(p2))=0;
p4:(-1)*num(part(p3,1))=0;

p5:coeff(part(p4,1),x^2);
p6:coeff(part(p4,1),x);
p7:ev(part(p4,1),x=0);

solve : quadratic equation ,two roots.

p8:float(solve([p4], [x])[1]);
p9:float(solve([p4], [x])[2]);

m:0.501;
b:0.43;
alpha:5.5;

y=x ,45 degree line.'ref line'

p10:y=x;

p11:part(p1,2)$
p12:part(p10,2)$
p13:part(p2,2)$

plot2d([p11,p12,p13],
     [x,-.5,1.3],[y,-.5,1.3],
 [plot_format, gnuplot])$

x=zf =.4 : feed of tower ? yf value for zf

p14:zf=0.4;
p15:yf=m*part(p14,2)+b;

p141:part(p14,2);
p151:part(p15,2);

xb < zf : b=bottom 'tower'

p16:xb=0.2;
p17:yb=part(p16,2);

p161:part(p16,2);
p171:part(p17,2);

line (zf,yf)-(xb,yb)

p18:y=p151+((p171-p151)/(p161-p141))*(x-p141);

p181:part(p18,2);

plot2d([p11,p12,p13,p181],
     [x,-.5,1.3],[y,-.5,1.3],
 [plot_format, gnuplot])$

intersection'magenta line' and 'blue curve'
symbolic: y = m1*x+b1

p19:y=m1*x+b1;

p20:ratsimp(rhs(p1)-rhs(p19))=0;
p21:(-1)*num(part(p20,1))=0;

p22:coeff(part(p21,1),x^2);
p23:coeff(part(p21,1),x);
p24:ev(part(p21,1),x=0);

p25:float(solve([p21], [x])[1]);
p26:float(solve([p21], [x])[2]);

general equation : between two points (point1-point2)
point1 : (x1,y1) 
point2 : (x2,y2)
m1=(y2-y1)/(x2-x1)
y=y1+m1*(x-x1) => m1=slope
b1=y1-m1*x1    => y=m1*x+b1

p27:m1=(p151-p171)/(p141-p161);
p28:b1=p171-part(p27,2)*p161;

intersection :   'purple line' and 
'blue curve = equilibrium equation'
symbolic: y = m1*x+b1 : stripping 'operating line'
          

p271:part(p27,2);
p281:part(p28,2);

negative root: p25 (quadratic equation 'p21')

p29:ev(part(p25,2),p27,p28);

y value :'p29'

p291:y=p271*p29+p281;

stripping-line x<0 ,first point

positive root: p26 (quadratic equation 'p21')

p30:ev(part(p26,2),p27,p28);

y value :'p30'

p301:y=p271*p30+p281;

stripping-line x>0 ,second point