Z-toolz

Z-toolz è una semplice libreria per scilab che sto scrivendo per poter risolvere diverse tipologie di esercizi,per caricarla in scilab è sufficiente copincollare questo codice in un file .sci e dare “load” da scilab.C’è da precisare che non è esente da errori,e che non può essere presa come riferimento per progetti”seri”,cmq è in continuo aggiornamento e perfezionamento.

//FUNZIONI DI SERVIZIO

//POSITIVE_ex
//(prende le componenti positive di un vettore e le mette in
//un vettore di pari dimensione le cui altre componenti saranno nulle
//ordina poi in maniera decrescente le componenti)
function [vp]=p_ex(vett)
s=size(vett);
N=s(1);
vp=zeros(N);
for j=1:3
if(vett(j)>0)
vp(j)=vett(j);
end;
end;
vp=sort(vp);
endfunction

//POSITIVE (restituisce true se tutte le componenti del vettore sono positive)
function [test]=positive(v)
sv=size(v);
s=sv(1);
test=%T;
for i=1:s
if(v(i)<0)
test=%F;
end;
end;
endfunction

//CHECK_FALSE (se almeno una componente è “FALSE” restituisce “FALSE”)

function [c]=check_f(v)
s=size(v);
c=%t;
for i=1:s(1)
if(v(i)==%f)
c=%f;
break;
end;
end
endfunction

//FUNZIONE DERIVATA (di una funzione “funz”)

function [z]=drv(x,funz,eps)
z=(funz(x+eps)-funz(x-eps))/(2*eps);
endfunction

//FUNZIONE GRADIENTE

function [g]=grad(x,funz,N,eps)

for i=1:N
h=zeros(N,1);
h(i)=eps;
g(i)=(funz(x+h)-funz(x-h))/(2*eps);
end

endfunction

//FUNZIONE HESSIANO

function [he]=hes(x,funz,N,eps)
he=zeros(N,N);
//saturo k,k
for i=1:N
h=zeros(N,1);
h(i)=eps;
he(i,i)=(funz(x+h)-2*funz(x)+funz(x-h))/eps^2;
end;
//

for i=1:N
for j=1:N
h=zeros(N,1);
h1=zeros(N,1);
h(i)=eps;
h1(j)=eps;
f1=funz(x-h-h1);
f2=funz(x+h-h1);
f3=funz(x+h+h1);
f4=funz(x-h+h1);
he(i,j)=1000*(f1+f3-f4-f2)/(4*eps);
end
end;

endfunction

//FUNZIONE ARMIJO
function [stp]=armijo(x,g,d,funz,delta,gamm,c)

stp=abs(g’*d)/norm(d)^2;
k=0;
while (funz(x+stp*d)>funz(x)+gamm*stp*g’*d)

stp=delta*stp;
if (k==200)
break;
end;
k=k+1;
end;

endfunction

//METODO DEL GRADIENTE CON PASSO Armijo
function [min_g]=m_grad(x0,funz,delta,gamm,c)
s=size(x0);
N=s(2);
x=x0′;

g=grad(x,funz,N,0.00001);
d=-g/norm(g);
k=0;
while (norm(g)>0.01)

d=-g/norm(g);
alpha=armijo(x,g,d,funz,delta,gamm,c);
printf(”\n passo\n %f”,alpha);
printf(”\n”);
x=x+alpha*d;
g=grad(x,funz,N,0.0001);
printf(”\n iterata\n %d”,k);
printf(”\n”);
if (k==1000)
break;
end;
k=k+1;

end;
min_g=x;
endfunction

//METODO DI BISEZIONE

function [r]=bisez(funz,a0,b0,eps)
a=a0;
b=b0;
while (abs(funz(a)-funz(b))>eps)

r=(a+b)/2;
value=funz(a)*funz(r)
if (value>0)
a=r;
end;
if value<0
b=r;
end;

end;
x=[a0:0.1:b0];
y=funz(x);
plot2d(x,y);
plot2d(x,y+100);

endfunction

//EULERO (eq.differenziali del tipo y’=f(x,y(x)))

function eul(funz,x0,N,a,b)

//N:numero di intervalli
//a,b:intervallo di integrazione

x=zeros(N+1);
y=zeros(N+1);
x(1)=x0;
y(1)=funz(x0);
h=abs(b-a)/N;

for i=2:N+1
x(i)=x(i-1)+h;
y(i)=y(i-1)+funz(x(i-1),y(i-1))*h;
end;
plot2d(x,y);
endfunction

//CARDANO
//a è il vettore dei coefficienti del polinomio cubico
function [x]=card(a)
x=zeros(3);
q=(a(1)^2-3*a(2))/9;
r=(2*a(1)^3-9*a(2)*a(1)+27*a(3))/54;
if (q^3-r^2<0)
x(1)=-sign(r)*(((r^2-q^3)^0.5+abs(r))^(1/3)+(q)/((r^2-q^3)^0.5+abs(r))^(1/3))-a(1)/3;
end;
if (q^3-r^2>0)
d=acos(r/q^(3/2));
x(1)=-2*q^0.5*cos(d/3)-a(1)/3;
x(2)=-2*q^0.5*cos((d+2*%pi)/3)-a(1)/3;
x(3)=-2*q^0.5*cos((d+4*%pi)/3)-a(1)/3;
end;
endfunction

//INTEGRAZIONE NUMERICA PER RETTANGOLI (di un set di dati)
function [I]=int_r(x,y)
s=size(x);
I=0;
for i=1:s(2)-1
temp(i)=(y(i)+y(i+1))/2*(x(i+1)-x(i));
I=I+temp(i);
end;
endfunction

//DERIVATA NUMERICA di un set di dati
function [d]=td(x,y)
s=size(x);
for i=2:s(2)-1
d(i-1)=(y(i+1)-y(i-1))/(x(i+1)-x(i-1));
end;
plot(x,y,d);
endfunction

///////////////////////////////////////////////////////////////////////////////////////////

//TENSIONE DI VAPORE DA VAN DER WAALS
function [td]=ps_vdw(T,Tc,Pc,plt)
//dichiarazione di costanti (eps è la tolleranza per il check)
eps=0.00000001;
R=8.31;
delta=0;
m=0.01;
k=0;
//calcolo a e b dai dati critici:
a=27*R^2*Tc^2/(64*Pc);
b=R*Tc/(8*Pc);
//dichiarazione della funzione vdw:
function [p]=vdw(v,T,a,b)
p=(8.31*T/(v-b))-a/v^2;
endfunction
//vettore coefficienti del polinomio per vmax vmin:
c_1=zeros(3);
c_1(1)=-2*a/(R*T);
c_1(2)=4*a*b/(R*T);
c_1(3)=-2*a*b^2/(R*T);
//vettore dei coefficienti del polinomio per vv e vl (dipende da pa,ps o pb)
function [c_2]=cp(T,p,a,b)
R=8.31;
c_2=zeros(3);
c_2(1)=-R*T/p-b;
c_2(2)=a/p;
c_2(3)=-a*b/p;
endfunction
//funzione phi
function [f]=phi(T,p,vl,vv,a,b)
R=8.31;
f=p*(vv-vl)-R*T*log((vv-b)/(vl-b))-a*(1/vv-1/vl);
endfunction
//calcolo limiti
vm=card(c_1);
vm=sort(vm);
vmin=vm(2);
vmax=vm(1);
pmin=vdw(vmin,T,a,b);
pmax=vdw(vmax,T,a,b);
//intervallo
delta=m*(pmax-pmin);
pa=pmax-delta;
while(pa<0)
m=m/2;
delta=m*(pmax-pmin);
pa=pmax-delta;
end;
pb=0+delta;
//affinamento pb
printf(”affino pb \n”);
cb=cp(T,pb,a,b);
v_b=card(cb);
sb=size(v_b);
while(sb(1)==1|positive(v_b)==%F)
printf(”*”);
pb=pb+delta;
if(pb>pa-2*delta)
break;
end;
cb=cp(T,pb,a,b);
v_b=card(cb);
sb=size(v_b);
end;
printf(”\n”);
//inizilizzazione
ca=cp(T,pa,a,b);
v_a=card(ca);
va_l=min(v_a);
va_v=max(v_a);
phi_a=phi(T,pa,va_l,va_v,a,b);
vb_l=min(v_b);
vb_v=max(v_b);
phi_b=phi(T,pb,vb_l,vb_v,a,b);
value=phi_a*phi_b;
printf(”inizializzazione: \n”);
printf(”pa=%f,pb=%f \n”,pa,pb);
printf(”phi(pa)*phi(pb)=%f \n”,value);
//bisezione
printf(”bisezione: \n”);
while(abs(phi_a-phi_b)>eps)
ps=(pa+pb)/2;
c=cp(T,ps,a,b);
v_s=card(c);
vs_l=min(v_s);
vs_v=max(v_s);
phi_s=phi(T,ps,vs_l,vs_v,a,b);
//check
value=phi_a*phi_s;
printf(”value =%f \n”,value);
if (value>0)
pa=ps;
c=cp(T,pa,a,b);
v_a=card(c);
va_l=min(v_a);
va_v=max(v_a);
phi_a=phi(T,pa,va_l,va_v,a,b);
end;
if (value<0)
pb=ps;
c=cp(T,pb,a,b);
v_b=card(c);
vb_l=min(v_b);
vb_v=max(v_b);
phi_b=phi(T,pb,vb_l,vb_v,a,b);
end;

end;
printf(”value =%f \n”,value);
printf(”\n”);
//plot
if(plt~=0)
h=(vmax*1.5-vmin*0.8)/100;
z=[vmin*0.8:h:vmax*1.5]‘;
sz=size(z);
s_z=sz(1);
for i=1:s_z
w1(i)=vdw(z(i),T,a,b);
w2(i)=pmax;
w3(i)=pmin;
w4(i)=ps;
end;
plot(z,[w1 w2 w3 w4]);
end;
lambda=T*8.31*log((vs_v-b)/(vs_l-b));
td=[ps vs_l vs_v a b lambda];
endfunction

//FUNZIONE PLOTTA PS vs T
function ps_vdw_T(Tc,Pc,T1,T2)
plt=0;
N=30;
h=(T2-T1)/N;
x=[T1:h:T2]‘;
sx=size(x);
for i=1:sx(1)
printf(”T=%f \n”,x(i));
p=ps_vdw(x(i),Tc,Pc,plt);
y(i)=p(1)/101300;
end;
plot(x-273,y,’r.->’);
endfunction

//FUNZIONE PLOTTA lambda vs T
function ps_vdw_T_l(Tc,Pc,T1,T2)
plt=0;
N=30;
h=(T2-T1)/N;
x=[T1:h:T2]‘;
sx=size(x);
for i=1:sx(1)
printf(”T=%f \n”,x(i));
p=ps_vdw(x(i),Tc,Pc,plt);
y(i)=p(1);
a(i)=p(4);
b(i)=p(5);
vs_l(i)=p(2);
vs_v(i)=p(3);
lambda(i)=p(6);
end;
plot(x,lambda,’b.->’);
endfunction

///////////////////////////////////////////////////////////////////////////////////////////

//PENG-ROBINSON
function [p]=pr_crit(v,T,Tc,Pc,w)
R=8.31;
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k*(1-(T/Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2)*alpha/Pc;
b=0.07780*R*Tc/Pc;
p=R*T/(v-b)-a/(v*(v+b)+b*(v-b));
endfunction

//PENG-ROBINSON
function [p]=pr(v,T,a,b,w)
R=8.31;
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k*(1-(T/Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2)*alpha/Pc;
b=0.07780*R*Tc/Pc;
p=R*T/(v-b)-a/(v*(v+b)+b*(v-b));
endfunction

//CALCOLO VOLUME CON PENG_ROBINSON
function [v]=vcalc(T,p,Tc,Pc,w)
R=8.31;
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k*(1-(T/Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2)*alpha/Pc;
b=0.07780*R*Tc/Pc;
a_1=zeros(3);
a_1(1)=b-R*T/p;
a_1(2)=-2*b*R*T/p+a/p-3*b^2;
a_1(3)=(b^2*R*T/p)-(a*b/p)+b^3;
v=card(a_1);
endfunction

//TENSIONE DI VAPORE CON PENG-ROBINSON
function [ps]=ps_pr(T,Tc,Pc,w)
R=8.31;
m=0.001;
eps=0.01;
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k*(1-(T/Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2)*alpha/Pc;
b=0.07780*R*Tc/Pc;
//coefficienti polinomio quarto grado (vmax,vmin)
a_2=zeros(5);
a_2(1)=b^4-2*a*b^3/(R*T);
a_2(2)=2*a*b^2/(R*T)-4*b^3;
a_2(3)=2*a*b/(R*T)+2*b^2;
a_2(4)=4*b-2*a/(R*T);
a_2(5)=1;
v_4=poly(a_2,”v”,”coeff”);
r_v_4=roots(v_4);
r_r_v_4=real(r_v_4);
vmin=r_r_v_4(3);
vmax=r_r_v_4(4);
pmin=pr(vmin,T,a,b,w);
pmax=pr(vmax,T,a,b,w);
//funzione phi
function [f]=phi(T,p,vl,vv,a,b)
R=8.31;
f=p*(vv-vl)-(sqrt(2)*a*log(b*(1-sqrt(2))+vl)/(4*b)-sqrt(2)*a*log(b*(1-sqrt(2))+vv)/(4*b)-R*T*log(vl-b)+R*T*log(vv-b)-sqrt(2)*a*log(b*(sqrt(2)+1)+vl)/(4*b)+sqrt(2)*a*log(b*(sqrt(2)+1)+vv)/(4*b));
endfunction
//intervallo
delta=m*(pmax-pmin);
pa=pmax-delta;
while(pa<0)
m=m/2;
delta=m*(pmax-pmin);
pa=pmax-delta;
end;
pb=0+delta;
//affinamento pb
v_b=vcalc(T,pb,Tc,Pc,w);
sb=size(v_b);
while(sb(1)==1|positive(v_b)==%F)
pb=pb+delta;
if(pb>pa-2*delta)
break;
end;
v_b=vcalc(T,pb,Tc,Pc,w);
sb=size(v_b);
end;
//inizializzazione
v_a=vcalc(T,pa,Tc,Pc,w);
va_l=min(v_a);
va_v=max(v_a);
phi_a=phi(T,pa,va_l,va_v,a,b);
v_b=vcalc(T,pb,Tc,Pc,w);
vb_l=min(v_b);
vb_v=max(v_b);
phi_b=phi(T,pb,vb_l,vb_v,a,b);
value=phi_a*phi_b;
//bisezione
counter=0;
while(abs(phi_a-phi_b)>eps)
counter=counter+1;
ps=(pa+pb)/2;
v_s=vcalc(T,ps,Tc,Pc,w);
vs_l=min(v_s);
vs_v=max(v_s);
phi_s=phi(T,ps,vs_l,vs_v,a,b);
//check
value=phi_a*phi_s;
if (value>0)
pa=ps;
v_a=vcalc(T,pa,Tc,Pc,w);
va_l=min(v_a);
va_v=max(v_a);
phi_a=phi(T,pa,va_l,va_v,a,b);
end;
if (value<0)
pb=ps;
v_b=vcalc(T,pb,Tc,Pc,w);
vb_l=min(v_b);
vb_v=max(v_b);
phi_b=phi(T,pb,vb_l,vb_v,a,b);
end;
if(counter==1000)
break;
end;
end;//(endwhile)

//plot
//h=(vmax*1.5-vmin*0.75)/100;
//z=[vmin*0.75:h:vmax*1.5]‘;
//sz=size(z);
//s_z=sz(1);
// for i=1:s_z
// w1(i)=pr(z(i),T,a,b);
// w2(i)=pr(vmax,T,a,b);
// w3(i)=pr(vmin,T,a,b);
// w4(i)=ps;
// end;
// plot(z,[w1,w2,w3,w4]);

endfunction

//FUNZIONE PLOTTA PS vs T (peng_robinson)
function ps_pr_T(Tc,Pc,w,T1,T2)
plt=0;
N=30;
h=(T2-T1)/N;
x=[T1:h:T2]‘;
sx=size(x);
for i=1:sx(1)
printf(”\n T=%f \n”,x(i));
p=ps_pr(x(i),Tc,Pc,w);
y(i)=p(1)/101300;
end;
plot(x-273,y,’r.->’);
endfunction

//COEFF.DI FUGACITÀ DI UN COMPONENTE PURO CON PENG-ROBINSON (liquido)
function [c]=c_liq(p,Tc,Pc,w)
R=8.31;
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k.*(1-(T./Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2).*alpha./Pc;
b=0.07780*R.*Tc./Pc;
v=vcalc(T,p,Tc,Pc,w);
vl=min(v);
zl=p*vl/(R*T);
B=b*p/(R*T);
A=a*p/(R*T)^2;
l_phi=(zl-1)-log(zl-B)-A/(2*sqrt(2)*B)*log((zl+(1+sqrt(2))*B)/(zl+(1-sqrt(2))*B));
c=exp(l_phi);
endfunction;

//COEFF.DI FUGACITÀ DI UN COMPONENTE PURO CON PENG-ROBINSON (vapore)
function [c]=c_vap(p,Tc,Pc,w)
R=8.31;
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k.*(1-(T./Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2).*alpha./Pc;
b=0.07780*R.*Tc./Pc;
v=vcalc(T,p,Tc,Pc,w);
vv=max(v);
zv=p*vv/(R*T);
B=b*p/(R*T);
A=a*p/(R*T)^2;
l_phi=(zv-1)-log(zv-B)-A/(2*sqrt(2)*B)*log((zv+(1+sqrt(2))*B)/(zv+(1-sqrt(2))*B));
c=exp(l_phi);
endfunction;

//COEFFICIENTI DI FUGACITA’ IN MISCELA CON PENG-ROBINSON EOS (liquido)
function [phi]=phi_liq(Tc,Pc,w,P,T,zm)
c_1=size(Tc);
c=c_1(2);
phi=zeros(c);
ac=zeros(c,c);
bc=zeros(c,c);
amix=0;
bmix=0;
//calcolo costanti:
R=8.31;
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k.*(1-(T./Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2).*alpha./Pc;
b=0.07780*R.*Tc./Pc;
for i=1:c
for j=1:c
ac(i,j)=sqrt(a(i)*a(j));
bc(i,j)=(b(i)+b(j))/2;
if(i==j)
ac(i,j)=a(i);
bc(i,j)=b(i);
end;
end;
end;

for i=1:c
for j=1:c
amix=amix+zm(i)*zm(j)*ac(i,j);
bmix=bmix+zm(i)*zm(j)*bc(i,j);
end;
end;
//calcolo volume liquido
coef=zeros(3);
coef(1)=bmix-R*T/P;
coef(2)=-2*bmix*R*T/P+amix/P-3*bmix^2;
coef(3)=(bmix^2*R*T/P)-(amix*bmix/P)+bmix^3;
v=card(coef);
vl=min(v);
//fattore di compressibilità della miscela
z=P*vl/(R*T);
//calcolo degli addendi per ln(phi)
ln_phi=zeros(c);
s=zeros(8);
for i=1:c
temp=0;
s(1)=b(i)/bmix*(z-1);
s(2)=log(z-bmix*P/(R*T));
s(3)=amix/(2*sqrt(2)*bmix*R*T);
for j=1:c
temp=temp+2*(zm(j)*ac(i,j));
end;
s(4)=temp;
s(5)=s(4)/amix-b(i)/bmix;
s(6)=z+(sqrt(2)+1)*bmix*P/(R*T);
s(7)=z-(sqrt(2)-1)*bmix*P/(R*T);
s(8)=s(6)/s(7);
ln_phi(i)=s(1)-s(2)-s(3)*s(5)*log(s(8));
phi(i)=exp(ln_phi(i));
end;
endfunction

//COEFFICIENTI DI FUGACITA’ IN MISCELA CON PENG-ROBINSON EOS (vapore)
function [phi]=phi_vap(Tc,Pc,w,P,T,zm)
R=8.31;
c_1=size(Tc);
c=c_1(2);
phi=zeros(c);
ac=zeros(c,c);
bc=zeros(c,c);
amix=0;
bmix=0;
//calcolo costanti:
k=0.37464+1.54226*w-0.26992*w^2;
alpha=(1+k.*(1-(T./Tc)^0.5))^2;
a=0.45724*R^2*(Tc^2).*alpha./Pc;
b=0.07780*R.*Tc./Pc;
for i=1:c
for j=1:c
ac(i,j)=sqrt(a(i)*a(j));
bc(i,j)=(b(i)+b(j))/2;
if(i==j)
ac(i,j)=a(i);
bc(i,j)=b(i);
end;
end;
end;

for i=1:c
for j=1:c
amix=amix+zm(i)*zm(j)*ac(i,j);
bmix=bmix+zm(i)*zm(j)*bc(i,j);
end;
end;
//calcolo volume liquido
coef=zeros(3);
coef(1)=bmix-R*T/P;
coef(2)=-2*bmix*R*T/P+amix/P-3*bmix^2;
coef(3)=(bmix^2*R*T/P)-(amix*bmix/P)+bmix^3;
v=card(coef);
vv=max(v);
//fattore di compressibilità della miscela
z=P*vv/(R*T);
//calcolo degli addendi per ln(phi)
ln_phi=zeros(c);
s=zeros(8);
for i=1:c
temp=0;
s(1)=b(i)/bmix*(z-1);
s(2)=log(z-bmix*P/(R*T));
s(3)=amix/(2*sqrt(2)*bmix*R*T);
for j=1:c
temp=temp+2*(zm(j)*ac(i,j));
end;
s(4)=temp;
s(5)=s(4)/amix-b(i)/bmix;
s(6)=z+(sqrt(2)+1)*bmix*P/(R*T);
s(7)=z-(sqrt(2)-1)*bmix*P/(R*T);
s(8)=s(6)/s(7);
ln_phi(i)=s(1)-s(2)-s(3)*s(5)*log(s(8));
phi(i)=exp(ln_phi(i));
end;
endfunction
//////////////////////////////////////////////
//COEFFICIENTI DI ATTIVITÀ DA PENG-ROBINSON
function [g]=gamma_pr(Tc,Pc,w,P,T,x)
s=size(x);
c=s(2);
g=zeros(c);
for i=1:c
g(i)=phi_liq(Tc(i),Pc(i),w(i),P,T,x(i))/(x(i)*c_liq(P,Tc(i),Pc(i),w(i)));
end;

endfunction;

//ATTIVITÀ DA PENG-ROBINSON
function [g]=act_pr(Tc,Pc,w,P,T,x)
s=size(x);
c=s(2);
g=zeros(c);
for i=1:c
g(i)=phi_liq(Tc(i),Pc(i),w(i),P,T,x(i))/(c_liq(P,Tc(i),Pc(i),w(i)));
end;

endfunction;

//PLOT COEFFICIENTI DI ATTIVITÀ
function plot_act(Tc,Pc,w,P,T)
n=30;
x=zeros(2,n-1);
for j=1:n-1
x(1,j)=(j)/(n);
x(2,j)=1-x(1,j);
end;

gamm=zeros(2,n-1);
x_1=x’;
perc=0;
winId=waitbar(’wait’);
realtimeinit(0);
for i=1:n-1
gamm_=gamma_pr(Tc,Pc,w,P,T,x_1(i,:));
gamm(1,i)=gamm_(1);
gamm(2,i)=gamm_(2);
realtime(0.1*i);
waitbar(i/(n-1),winId);
end;
winclose(winId);
plot(x(1,:),log((gamm(1,:))),’b.->’);
plot(x(1,:),log((gamm(2,:))),’r.->’);
pause
endfunction

//CALCOLO DEL FLASH peng_robinson (da correggere)
function [zeq]=flash_pr(P,T,zm,Tc,Pc,w)
//costanti
s=size(zm);
c=s(2);
pb=0;
pvap=zeros(c);
eps=0.0001;
//calcolo approx di P inizio ebollizione
for i=1:c
pvap(i)=ps_pr(T,Tc(i),Pc(i),w(i));
pb=pb+zm(i)*pvap(i);
end;
//calcolo k di primo tentativo
k=pvap/pb;
//definisco la funzione phi
function [f]=psi(alpha_1)
f=0;
for i=1:c
f=f+zm(i)*(1-k(i))/(alpha_1+k(i)*(1-alpha_1));
end
endfunction
check=%f;
counter=0;
while(check==%f)
//calcolo alpha con bisezione
printf(”biseziono \n”);
alpha_2=bisez(psi,0,1,0.00001);
if(alpha_2(2)==%f)
printf(”radice non trovata” \n);
break;
end;
alpha=alpha_2(1);
//calcolo x e y
x=zeros(c);
y=zeros(c);
for i=1:c
x(i)=zm(i)/(alpha+k(i)*(1-alpha));
y(i)=k(i)*zm(i)/(alpha+k(i)*(1-alpha));
end;
//calcolo coefficienti fugacità
phi_l=phi_liq(Tc,Pc,w,P,T,zm);
phi_v=phi_vap(Tc,Pc,w,P,T,zm);
//ricalcolo k
k=phi_l./phi_v;
//ricalcolo x,y
for i=1:c
x_1(i)=zm(i)/(alpha+k(i)*(1-alpha));
y_1(i)=k(i)*zm(i)/(alpha+k(i)*(1-alpha));
end;
deltax=abs(x-x_1);
deltay=abs(y-y_1);
check_1=deltax<eps;
check_2=deltay<eps;
check_3=check_1&check_2;
check=check_f(check_3);
counter=counter+1;
if(counter==1000)
printf(”raggiunto il numero max di iterate” \n);
break;
end;

end;//endwhile
zeq=[x y];
endfunction

//CALCOLO DEL FLASH PER PIÙ TEMPERATURE PRESSIONE COSTANTE
function [z]=flash_pr_T(P,T_v,zm,Tc,Pc,w)
s_z=size(zm);
c=s_z(2);
s_t=size(T_v);
for i=1:s_t(2)
z_T=flash_pr(P,T_v(i),zm,Tc,Pc,w)
x(:,i)=z_T(:,1);
y(:,i)=z_T(:,2);
end
z=[x;y];
plot(x(1,:),y(1,:),’b.->’);
pause
endfunction

//CALCOLO PRESSIONE INIZIO EBOLLIZIONE
function [pb]=pib_pr(T,x,Tc,Pc,w)
//costanti
s=size(x);
c=s(2);
pb=0;
pvap=zeros(c);
eps=0.0001;
//calcolo approx di P inizio ebollizione
for i=1:c
pvap(i)=ps_pr(T,Tc(i),Pc(i),w(i));
pb=pb+x(i)*pvap(i);
end;
//calcolo k di primo tentativo
k=pvap/pb;
//calcolo y
for i=1:c
y(i)=k(i)*x(i);
end;
//è possibile che le yi siano maggiori di 1 allora correggo
scarto=(sum(y)-1)/sum(y);
y=y*(1-scarto);

check_est=%f;
check_int=%f;
count_int=0;
count_est=0;
while(check_est==%f)
phi_l=phi_liq(Tc,Pc,w,pb,T,x);
// printf(”est \n”);
while(check_int==%f)
phi_v=phi_vap(Tc,Pc,w,pb,T,y’);
//calcolo y_1
//printf(”int \n”);
y_1=(y.*phi_l)./phi_v;
//è possibile che le y_1i siano maggiori di 1 allora correggo
scarto=(sum(y_1)-1)/sum(y_1);
y_1=y_1*(1-scarto);
deltay=abs(y-y_1);
check_int=check_f(deltay<eps);
y=y_1;
count_int=count_int+1;
if(count_int==100)
break;
end;
end;//while interno
check_est=check_f(abs(sum(y_1)-1)<eps);
pb=pb*sum(y_1);
count_est=count_est+1;
if(count_est==100)
break;
end;
end;//while esterno

endfunction

//CALCOLO PRESSIONE INIZIO CONDENSAZIONE
function [pd]=pic_pr(T,y,Tc,Pc,w)
//costanti
s=size(y);
c=s(2);
pb=0;
pvap=zeros(c);
eps=0.0001;
//calcolo approx di P inizio ebollizione
pd_1=0;
for i=1:c
pvap(i)=ps_pr(T,Tc(i),Pc(i),w(i));
pd_1=pd_1+y(i)/pvap(i);
end;
pd=1/pd_1;
//calcolo k di primo tentativo
k=pvap/pd;
//calcolo y
for i=1:c
x(i)=y(i)/k(i);
end;
//è possibile che le xi siano maggiori di 1 allora correggo
scarto=(sum(x)-1)/sum(x);
x=x*(1-scarto);
check_est=%f;
check_int=%f;
count_int=0;
count_est=0;
while(check_est==%f)
phi_l=phi_liq(Tc,Pc,w,pd,T,x’);
// printf(”est \n”);
while(check_int==%f)
phi_v=phi_vap(Tc,Pc,w,pd,T,y);
//calcolo y_1
//printf(”int \n”);
x_1=(x.*phi_l)./phi_v;
//è possibile che le x_1i siano maggiori di 1 allora correggo
scarto=(sum(x_1)-1)/sum(x_1);
x_1=x_1*(1-scarto);
deltay=abs(x-x_1);
check_int=check_f(deltay<eps);
x=x_1;
count_int=count_int+1;
if(count_int==100)
break;
end;
end;//while interno
check_est=check_f(abs(sum(x_1)-1)<eps);
pd=pd*sum(x_1);
count_est=count_est+1;
if(count_est==100)
break;
end;
end;//while esterno

endfunction

//GRAFICO ISOTERMO MISCELA BINARIA
function [pb]=pib_pr_x(T,Tc,Pc,w)
//x_v è matrice di vettori composizione
n=30;
x=zeros(2,n);
for j=1:n+1
x(1,j)=(j-1)/n;
x(2,j)=1-x(1,j);
end;
pb=zeros(1,n+1);
pd=zeros(1,n+1);
x_1=x’;
perc=0;
winId=waitbar(’wait’);
realtimeinit(0);
for i=1:n+1
pb(i)=pib_pr(T,x_1(i,:),Tc,Pc,w);
pd(i)=pic_pr(T,x_1(i,:),Tc,Pc,w);
realtime(0.1*i);
waitbar(i/(n+1),winId);
end;
winclose(winId);
plot(x(1,:),pd/101300,’b.->’);
plot(x(1,:),pb/101300,’r.->’);
endfunction

//FUNZIONI PROVA
deff(”z=objf(x)”,”z=2*(x(1)+10)^2+3*x(2)^2+x(2)”);
deff(”z=objf1(x)”,”z=x(1)^3+x(1)*x(2)-x(1)*x(2)^2″);
deff(”z=objf2(x)”,”z=x(1)^2+x(2)^2″);
deff(”y=objf3(x)”,”y=2*x^3+x^2+x-3+2″);
deff(”y1=objf4(x,y)”,”y1=5*x”);
deff(”y2=objf5(x,y)”,”y2=2*x^2+y”);