% Gas membrane separation calculation
%
Ph=[15 25 35 45 55];
xf=.229;
ypv=[.342 .394 .415 .431 .427;.367 .411 .456 .472 .476;.374 .436 .478 .510 .518];
xrv=[.115 .101 .081 .075 .060;.178 .153 .135 .125 .111;.200 .189 .180 .174 .161];
Kr=2;a=6;yi=0.2;theta=0.5;
dy=.01;
f1='yi*((1-xr)*r-(1-yi))-a*(1-yi)*(xr*r-yi)';
f2='xr*(1-theta)+yp*theta-xf';
dt1='((xf*r-yp)*(xr*r-yi)*(xf*r-yp+xr*r-yi)/2)^(1/3)';
f3='yp*Kr*theta-(1-theta)*a*eval(dt1)';
dt2='(((1-xf)*r-(1-yp))*((1-xr)*r-(1-yi))*((1-xf)*r-(1-yp)+(1-xr)*r-(1-yi))/2)^(1/3)';
f4='(1-yp)*Kr*theta-(1-theta)*eval(dt2)';
%
df1dk=0;df1dt=0;
df2dk=0;df2da=0;df2dy=0;
df4da=0;
for ni=1:3;
for nj=1:5;
    xr=xrv(ni,nj);yp=ypv(ni,nj);r=(Ph(nj)+14.7)/14.7;
%
x=[Kr a yi theta];
disp('Guessed values')
fprintf('KR = %8.3f, alfa = %8.3f, yi = %8.3f, theta = %8.3f\n',Kr,a,yi,theta)
for i=1:20
%
f=[eval(f1) eval(f2) eval(f3) eval(f4)];
df1da=(1-yi)*(xr*r-yi);
df1dy=(1-xr)*r-(1-yi)+yi+(xr*r-yi)+a*(1-yi);
df2dt=-xr+yp;
dlm1=eval(dt1);dlm2=eval(dt2);
df3dk=yp*theta;
df3da=-(1-theta)*dlm1;
df3dt=yp*Kr+a*dlm1;
df4dk=(1-yp)*theta;
df4dt=(1-yp)*Kr+dlm2;
yi=yi+dy;
df3dy=(eval(f3)-f(3))/dy;
df4dy=(eval(f4)-f(4))/dy;
yi=yi-dy;
Jt=[df1dk df1da df1dy df1dt;
    df2dk df2da df2dy df2dt;
    df3dk df3da df3dy df3dt;
    df4dk df4da df4dy df4dt;];
dx=Jt\f';
x=x-dx';
Kr=x(1);a=x(2);yi=x(3);theta=x(4);
a_exp=yp*(1-xr)/xr/(1-yp);
if max(abs(dx))<.001, break, end
end
fprintf('For P/p = %8.3f,Calculated values from %g iterations\n',r,i)
fprintf('KR = %8.3f, alfa = %8.3f, alfa exp = %8.3f, yi = %8.3f, theta = %8.3f\n',Kr,a,a_exp,yi,theta)
disp(' ')
end
end