Undefined function or variable C

I've been given a function by my professor that calculates PT flash for non -ideal systems. I'm trying to use it to simulate an experiment. But continously get error.

I think it may be the way i've set my input up.

The sub-functions are given below:

function [PD,xDew]=PDew_Raoult(Ps,y,T)

% compute dew pressure at defined T for a mixture with a user given vapour
% pressure using the Raoult Law
c=length(y);
for i=1:c
    Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));

PD=1/sum(y./Ps);
xDew=y*PD./Ps;

end

function[PB,Ybub]=Pbubble_Raoult(Ps,x)
%this is a function that calculates the bubble point of a fluid with given
%vapour pressure and composition 
 Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)))
i=length(x);
PB=Ps.*x;
Ybub=Ps.*x/PB;
end

function [ alphaV ] = RACHFORDRICE_BISECT( K,Z )
%This function applies the bisection method to the Rachford-Rice equation
%for alphaV in (0,1).
%Input/output are on molar basis. 

a=0;
b=1;
%psi_0=sum(z.*(K-1))
%psi_1=sum((z.*(K-1))./K)

while b-a>0.000001
    alphaV=(a+b)/2;
    psi_alphaV=sum((Z.*(K-1))./((1-alphaV)+alphaV*K));
    psi_b=sum((Z.*(K-1))./((1-b)+b*K));
    if psi_alphaV*psi_b<0
        a=alphaV;
    else
        b=alphaV;
    end
end

function[PB,Ybub]=PB_VLE_Wilson(Z,C,MW,rhoL,T,BIP)
% this is a function that calculates the bubble point using the wilson
% activity coeefficient model.
c=length(Z);
for i=1:c
    Psi=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
    [gamma,~]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z);
    PB=sum(Psi.*Z.*gamma);
    Ybub=Psi.*gamma/PB;
end

function [PD,xDew]=PD_VLE_wilson (C,MW,rhoL,BIP,T,Z)
% this is a function that calculates the dewpoint of a mixture with the
% wilson correlation for generating activity coefficients for non ideal 
%liquid phase
c=length(Z);
for i=1:c
     Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
end 
% initial guess of Pdew using PD_VLE_ID(C,T,y)
[PD,xDew]=PDew_Raoult(Ps,Z,T); %initial guess of PD(1) and calculation of x1

[gamma,~]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z);% calculation of gamma at the initial guess of x1
iter=0;
while  max(abs((PD*Z-Ps.*x.*gamma)./(PD*Z)))>0.001 &&iter<1000  %isofugacity conditions and control on the  number of iterations
    PD=1/sum(Z./Ps.*gamma);
    xDew=(PD*Z)./(Ps.*gamma);
    [gamma,~]=ACTIVITY_WILSON(MW,rhoL,BIP,T,x);% initial guess of gamma
    iter =iter +1;
end

if iter<1000
    PD=1/sum(Z./Ps.*gamma);
    xDew=(PD*Z)./(Ps.*gamma);
else
    PD=0;
    xDew=0;
    disp("Dew point not found");
end
end

function [gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z)
%This program calculates the activity coefficients (gamma) and the 
%activities (a) of each component of a mixture of c components using the 
%Wilson model. 
%INPUT PARAMETERS: MW: vector (1xc) reporting the molecular weights of the 
%c components; rhoL, vector (1xc) reporting the liquid density of the c 
%pure components at temperature T; BIP is a matrix cxc reporting the energy   
%interaction parameters (BIP(i,j)=lambda_ij-lambda_ii, J/mol). The energy
%interaction parameters are temperature independent; T: temperature of the 
%system; x vector (1xc) reporting the mole fractions of the components of 
%the mixture.   
%OUTPUT PARAMETERS: gamma: vector 1xc reporting the activity 
%coefficients of the components of the mixture; a: vector 1xc reporting the 
%activities of the components of the mixture. 
%Unless otherwise stated, all input/output parameters are expressed
%according to MKS.

R=8.314;
c=length(Z);

%Molar volumes of the pure liquid components composing the mixture
Vl=1./((rhoL*1000)./MW);

%Lambda terms (dimensionless) of the Wilson formula
for i=1:c
    for j=1:c
        Lambda(i,j)=(Vl(j)/Vl(i))*exp(-BIP(i,j)/(R*T));
    end
end

for i=1:c
    for j=1:c
        A=sum(Z.*Lambda(j,:));
        C(j)=Z(j)*Lambda(j,i)/A;
    end
    lngamma(i)=1-log(sum(Z.*Lambda(i,:)))-sum(C);
    gamma(i)=exp(lngamma(i));
    a(i)=gamma(i)*Z(i);
end

end
%% the main function to be ran is given as


function[Xeq, Yeq, alphaV, Fl, Fv]=PTFLASH_VLE_Wilson(C,MW,rhoL,BIP,P,T,Z)
%this is a function that calculates PT flash for non -ideal systems where K
%cannot be explicit since it depends on X which is also a variable
c=length(Z);
for i=1:c
  
    Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));
end

% validating the possibility of having VLE

[PB,y]=PB_VLE_Wilson(Z,C,MW,rhoL,T,BIP);

[PD,x]=PD_VLE_wilson (C,MW,rhoL,BIP,T,Z);

% checking if there is VLE
if  P<=PD           %% y here represents the mole fraction of componenets
    Xeq=0;
    Yeq=y;
    alphaV=1;
    Fl=0;
    Fv=P*Z;    
elseif P>=PB           %% x here represents the mole fraction of components
    Xeq=x;
    Yeq=0;
    alphaV=0;
    [gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Xeq);
    Fl=Ps.*Xeq.*gamma;
    Fv=0;
else
    Xeq=(x+y)/2 ;      % first guess on x
    [gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Xeq);   %first guess on gamma
    K=(Ps.*gamma)/P;  % first guess on K
     Psi_0=sum(z.*(K-1));
        Psi_1=sum(z.*(K-1)./K);
        if Psi_0*Psi_1>=0    % bad initial guess
    Xeq=0;
    Yeq=0;
    alphaV=0;
    Fl=0;
    Fv=0;
    disp("bad initial guess")
        else       %supposedly good initial guess
            Fl=zeros(1,c);     % trick to enter the while cycle at the begininng
            Fv=Fl+1;
            iter=0;            % initialization of iteration count
            while max(abs((Fv-Fl)./Fv))>0.00001 && iter<10000 && Psi_0*Psi_1<0
                [ alphaV ] = RACHFORDRICE_BISECT( K,Z );
     xeq=Z./((1-alphaV)+alphaV*K);
            Yeq=K.*xeq;
           [gamma,a]=ACTIVITY_WILSON(MW,rhoL,BIP,T,Z);   %Gamma new guess
            K=(Ps.*gamma)/P;
            psi_0=sum(z.*(K-1));
            psi_1=sum(z.*(K-1)./K);
            Fv=P*Yeq;
            Fl=Ps.Xeq.*gamma;
            iter=iter+1;
            end
        end
end
end

[Rik]

There are many issues with the code you posted.

• There are some extra or missing end keywords, so we can't run this code at all
• You didn't provide reasonable inputs for PDew_Raoult
• There are still several mlint warnings
• The indentation suggests that some code was intended to be in a loop
• You assing a value to c (lower case), but then also use C (upper case) without defining it

And this is already from the first few lines. That last one will already cause the error you got.

Answer 1

[Guillaume]

If this exactly the code you've been given without you having made any modification:

function [PD,xDew]=PDew_Raoult(Ps,y,T)
% compute dew pressure at defined T for a mixture with a user given vapour
% pressure using the Raoult Law
c=length(y);
for i=1:c
    Ps=exp((C(i,1))+(C(i,2)/T)+(C(i,3)*log(T))+(C(i,4)*T^C(i,5)));

then it was never working and you need to go back to your professor to ask for code that actually work.

C is indeed undefined (and clearly not meant to be the same as the lowercase c). Until you know what C is meant to be, there's nothing you can do to fix it. I would suspect that C is meant to be y with y a Nx5 matrix (with N>5 otherwise, you'll hit a bug in the code!)

Also completely absurd is that the function takes Ps as an input and then immediately discard that input to replace it with a calculated value.

Also, completely absurd is tha the function calculates values in a loop and at each step of the loop, completely throws away the values calculated in the previous iteration. In the end, it only keeps the values calculated in the last, the previous steps having been a waste of time.

So, again, that code never worked, even if it had worked, it had bugs, made no sense and was poorly written.