Please help me to rewrite the following codes

%% Calculation of phi and f
 
    a11 =  ai(1);
    a12 = ((ai(1)*ai(2))^0.5)*(1-bin_PR(1,2));
    a13 = ((ai(1)*ai(3))^0.5).*(1-bin_PR(1,3));
    
    a21 = a12;
    a22 = ai(2);
    a23 = ((ai(2)*ai(3))^0.5)*(1-bin_PR(2,3));
    
    a31 = a13;
    a32 = a23;
    a33 = ai(3);
    
    sum_xa(1) = x(1)*a11 + x(2)*a12 +x(3)*a13;
    sum_xa(2) = x(1)*a21 + x(2)*a22 +x(3)*a23;
    sum_xa(3) = x(1)*a31 + x(2)*a32 +x(3)*a33;
  
    sum_ya(1) = y(1)*a11 + y(2)*a12 +y(3)*a13;
    sum_ya(2) = y(1)*a21 + y(2)*a22 +y(3)*a23;
    sum_ya(3) = y(1)*a31 + y(2)*a32 +y(3)*a33;
 
    phi_L=exp((bi./b_L).*(ZLt-1)-log(ZLt-B_L)-(A_L./(B_L.*(-2*sqrt(2)))).*(log((ZLt+(1-sqrt(2)).*B_L)/(ZLt+(1+sqrt(2)).*B_L)))*((2*sum_xa/a_L)-(bi./b_L)));
    
 
    phi_G=exp((bi./b_G).*(ZGt-1)-log(ZGt-B_G)-(A_G./(B_G.*(-2*sqrt(2)))).*(log((ZGt+(1-sqrt(2)).*B_G)/(ZGt+(1+sqrt(2)).*B_G)))*((2*sum_ya/a_G)-(bi./b_G)));
  
------------------------------------
-----------------------------------------
I want to write the above coding like this (below one):
for i=1:3
        for j=1:3
            el_L(i,j)=x(i).*x(j)*((ai(i).*ai(j)).^0.5).*(1-bin_PR(i,j));
            el_G(i,j)=y(i).*y(j)*((ai(i).*ai(j)).^0.5).*(1-bin_PR(i,j));
        end
    end
 
    a_L=sum(el_L(:));
    b_L=sum(x.*bi);

[madhan ravi]

provide all the values

[Pulok Deb]

clc
clear
format compact

T=input('Temperature(K):');        % Use T = 377.6
P=input('Pressure(MPa):');         % Use P = 27.58

C1  = 0.5;
C6  = 0.3;
C10 = 0.2;

Z =  [C1  C6  C10];
tc = [190.56  507.6   617.7];
pc = [4.599   3.025   2.11];
w =  [0.0115  0.3013  0.4923];
bin_PR =[0        0.0422    0.05;
         0.0422   0         0   ;
         0.05     0         0    ];
     
%%

K=(pc./P).*exp(5.37.*(1+w).*(1-(tc./T)));

err=1;   

while err>10^(-12)
    
    %% Calculation of nv
 
    nv=0.5;
    er=1;
    
    while er>10^(-12)
        
        f=(Z.*(K-1))./(1+nv.*(K-1));                   %f(nv) Equation
        f_prim=(-Z.*((K-1).^2))./(1+nv.*(K-1)).^2;     %f'(nv) Equation
        nv_new=nv;
        nv=nv-(sum(f)/sum(f_prim));
        er=(nv_new - nv)^2;  
    end
    
    x = Z./(1+(K-1).*nv);
    y = x.*K;

    %% Calculatiion of a,b,A,B in PR EOS

    m=0.3796+1.485.*w-0.1644.*(w.^2)+0.01667.*(w.^3);
    ac=(0.457235.*(0.008314.*tc).^2)./pc;
    alpha=(1+m.*(1-(T./tc).^0.5)).^2;
    ai=ac.*alpha;
    bi=0.077796.*0.008314.*tc./pc;
    
    for i=1:3
        for j=1:3
            el_L(i,j)=x(i).*x(j)*((ai(i).*ai(j)).^0.5).*(1-bin_PR(i,j));
            el_G(i,j)=y(i).*y(j)*((ai(i).*ai(j)).^0.5).*(1-bin_PR(i,j));
        end
    end

    a_L=sum(el_L(:));
    b_L=sum(x.*bi);
    A_L=(a_L.*P)./(0.008314.*T).^2;
    B_L=(b_L.*P)./(0.008314.*T);

    a_G=sum(el_G(:));
    b_G=sum(y.*bi);
    A_G=(a_G.*P)./(0.008314.*T).^2;
    B_G=(b_G.*P)./(0.008314.*T);

    %% Calculation of ZL and ZG

    a1L=-(1-B_L);
    a2L=(A_L-2*B_L-3*B_L^2);
    a3L=-((A_L*B_L)-(B_L^2)-(B_L^3));

    Eq_L=[1,a1L,a2L,a3L];
    Z_L=roots(Eq_L);
    ZL=Z_L(imag(Z_L)==0);
    Le=length(ZL);

    if Le==1
        ZLt=ZL;
    
    elseif Le==3
        if ZL(1)<=ZL(2) && ZL(1)<=ZL(3)
        ZLt=ZL(1);
        elseif ZL(2)<ZL(1)&& ZL(2)<=ZL(3)
        ZLt=ZL(2);
        else 
        ZLt=ZL(3);
        end
    end

    a1G=-(1-B_G);
    a2G=(A_G-2*B_G-3*B_G^2);
    a3G=-((A_G*B_G)-(B_G^2)-(B_G^3));

    Eq_G=[1,a1G,a2G,a3G];
    Z_G=roots(Eq_G);
    ZG=Z_G(imag(Z_G)==0);

    Le=length(ZG);

    if Le==1
        ZGt=ZG;
        
    elseif Le==3    
        if ZG(1)>=ZG(2) && ZG(1)>=ZG(3)
        ZGt=ZG(1);
        elseif ZG(2)>=ZG(1)&& ZG(2)>=ZG(3)
        ZGt=ZG(2);
        else 
        ZGt=ZG(3);
        end
    end

%% Calculation of phi and f

    a11 =  ai(1);
    a12 = ((ai(1)*ai(2))^0.5)*(1-bin_PR(1,2));
    a13 = ((ai(1)*ai(3))^0.5).*(1-bin_PR(1,3));
    
    a21 = a12;
    a22 = ai(2);
    a23 = ((ai(2)*ai(3))^0.5)*(1-bin_PR(2,3));
    
    a31 = a13;
    a32 = a23;
    a33 = ai(3);
    
    sum_xa(1) = x(1)*a11 + x(2)*a12 +x(3)*a13;
    sum_xa(2) = x(1)*a21 + x(2)*a22 +x(3)*a23;
    sum_xa(3) = x(1)*a31 + x(2)*a32 +x(3)*a33;
  
    sum_ya(1) = y(1)*a11 + y(2)*a12 +y(3)*a13;
    sum_ya(2) = y(1)*a21 + y(2)*a22 +y(3)*a23;
    sum_ya(3) = y(1)*a31 + y(2)*a32 +y(3)*a33;

    phi_L=exp((bi./b_L).*(ZLt-1)-log(ZLt-B_L)-(A_L./(B_L.*(-2*sqrt(2)))).*(log((ZLt+(1-sqrt(2)).*B_L)/(ZLt+(1+sqrt(2)).*B_L)))*((2*sum_xa/a_L)-(bi./b_L)));
    f_L=x.*P.*phi_L;

    phi_G=exp((bi./b_G).*(ZGt-1)-log(ZGt-B_G)-(A_G./(B_G.*(-2*sqrt(2)))).*(log((ZGt+(1-sqrt(2)).*B_G)/(ZGt+(1+sqrt(2)).*B_G)))*((2*sum_ya/a_G)-(bi./b_G)));
    f_G=y.*P.*phi_G;

%%
     nn=(1-(f_L./f_G)).^2;
     mm=(f_L./f_G);
     
     nn(isnan(nn))=0;
     mm(isnan(mm))=0;
     
     K=K.*mm;
     err=sum(nn(:));

end

disp('==============================================================================================')
disp('                                Flash Calculation')
disp('                                Peng-Robinson EOS')
disp('==============================================================================================')
Z
x
y

[Pulok Deb]

Above code is for 3 components. If I do calculation for 10 or 20 components, than how should I write the code at " Calculation of Phi and f " section. If I write like this, it will be large one. So I want to use loop and presented the coding simple with accurate result.

[madhan ravi]

i didn't get any error running your code

[Pulok Deb]

I know its right code. Here we use 3 components and use (i,j) from 1:3. But if I will do calculate 20 components, from 1:20.. It would be a large coding .. I want to make simple and smaller here in this part

[madhan ravi]

no you won't if you chane three to twenty you will get the results if and only if the other variables have the compatable sizes

[Pulok Deb]

so there is no way to write a simple loop or else for calculation of phi and f. right? but my instructor tell me its possible. I am confused now

Answer 1

[KSSV]

YOu need not to use a loop for your code. YOu can achieve it by simple matrix multiplication. Check the below code.

a11 =  ai(1);
a12 = ((ai(1)*ai(2))^0.5)*(1-bin_PR(1,2));
a13 = ((ai(1)*ai(3))^0.5).*(1-bin_PR(1,3));

a21 = a12;
a22 = ai(2);
a23 = ((ai(2)*ai(3))^0.5)*(1-bin_PR(2,3));

a31 = a13;
a32 = a23;
a33 = ai(3);
A = [a11 a12 a13 ; a21 a22 a23;a31 a32 a33] ;

b1 = x ;   % this should be a 3X1 array
b2 = y ;   % this should be a 3X1 array

sum_xa = A*b1 ;
sum_ya = A*b2 ;

phi_L=exp((bi./b_L).*(ZLt-1)-log(ZLt-B_L)-(A_L./(B_L.*(-2*sqrt(2)))).*(log((ZLt+(1-sqrt(2)).*B_L)/(ZLt+(1+sqrt(2)).*B_L)))*((2*sum_xa/a_L)-(bi./b_L)));
phi_G=exp((bi./b_G).*(ZGt-1)-log(ZGt-B_G)-(A_G./(B_G.*(-2*sqrt(2)))).*(log((ZGt+(1-sqrt(2)).*B_G)/(ZGt+(1+sqrt(2)).*B_G)))*((2*sum_ya/a_G)-(bi./b_G)));

[Pulok Deb]

I am looking for For loop &/or if-else loop. Thanks though

[KSSV]

a11 =  ai(1);
a12 = ((ai(1)*ai(2))^0.5)*(1-bin_PR(1,2));
a13 = ((ai(1)*ai(3))^0.5).*(1-bin_PR(1,3));
a21 = a12;
a22 = ai(2);
a23 = ((ai(2)*ai(3))^0.5)*(1-bin_PR(2,3));
a31 = a13;
a32 = a23;
a33 = ai(3);
A = [a11 a12 a13 ; a21 a22 a23;a31 a32 a33] ;

sum_xa = zeros(3,1) ; 
sum_ya = zeros(3,1) ;

for i = 1:3
    sum_xa(i) = A(i,1)*x(1)+A(i,2)*x(2)+A(i,3)*x(3) ;
    sum_ya(i) = A(i,1)*y(1)+A(i,2)*y(2)+A(i,3)*y(3) ;
end

phi_L=exp((bi./b_L).*(ZLt-1)-log(ZLt-B_L)-(A_L./(B_L.*(-2*sqrt(2)))).*(log((ZLt+(1-sqrt(2)).*B_L)/(ZLt+(1+sqrt(2)).*B_L)))*((2*sum_xa/a_L)-(bi./b_L)));
phi_G=exp((bi./b_G).*(ZGt-1)-log(ZGt-B_G)-(A_G./(B_G.*(-2*sqrt(2)))).*(log((ZGt+(1-sqrt(2)).*B_G)/(