function [ q,BETA,z,phi ] = EOS_CAL(NC,PRESS,TEMP,phasetype,COMP,PROPS,BIC)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
%RK cubic EOS
RGAS=83.14;
SIGMA=1;EPSILON=0.0;OMEGA=0.08664;SI=0.42748;
a=zeros(NC,1);
b=zeros(NC,1);
ALPHA=zeros(NC,1);
TR=zeros(NC,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ic=1:NC
    TR(ic)=TEMP/PROPS(ic,1);
    ALPHA(ic)=TR(ic)^-0.5;
    a(ic)=SI*ALPHA(ic)*RGAS^2*PROPS(ic,1)^2/PROPS(ic,2);
    %a(ic) tabe yi ast az savabet va damaye mixture be damaye bohranie
    % ajzaye szande
    b(ic)=OMEGA*RGAS*PROPS(ic,1)/PROPS(ic,2);
    %b(ic) tabe yi ast savabetva damaye bohrani be feshare bohranie ajzaye
    %mixture
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
amix=0;
bmix=0;
qfc=zeros(NC,1);
for ic=1:NC
    bmix=bmix+COMP(ic)*b(ic);
    for jc=1:NC
        amix=amix+COMP(ic)*COMP(jc)*...
            (a(ic)*a(jc))^0.5*(1-BIC(ic,jc));
        qfc(ic)=qfc(ic)+COMP(jc)*a(jc)^0.5*(1-BIC(ic,jc));
        
    end
end
% vaed shodNE q & BETA
BETA=bmix*PRESS/(RGAS*TEMP);
q=amix/(bmix*RGAS*TEMP);
if(phasetype==0)
    %GAS
    z=0;
    znew=1;
    while(abs(znew-z)>10^-4)
        z=znew;
        znew=1+BETA-q*BETA*(z-BETA)/((z+EPSILON*BETA)*(z+SIGMA*BETA));
    end
else
    %LIQUID
    z=0;
    znew=BETA;
    while(abs(znew-z)>10^-4)
        z=znew;
        znew=BETA+(z+EPSILON*BETA)*(z+SIGMA*BETA)*(1+BETA-z)/(q*BETA);
    end
end
%vared shodane I & phi ba estefade az q & BETA & I
qbar=zeros(NC,1);
phi=zeros(NC,1);
I=log((z+BETA)/z);
for ic=1:NC
    qbar(ic)=q*((2*a(ic)^0.5*qfc(ic)/amix)-(b(ic)/bmix));
    phi(ic)=(b(ic)/bmix)*(z-1)-log(z-BETA)-qbar(ic)*I;
    phi(ic)=exp(phi(ic));
end

end