Редактирование сценария изменения в MatLab

Question

Я пытаюсь дать/set 'rho_real' этому коду и получить P_Mpa как/набор выходов, однако исходный 'rho_real', похоже, не изменяется, поскольку я запускаю функция.

Любая помощь с благодарностью.

function P_Mpa = Density_vs_Pressure_Ethane (rho_real) 

% calculating Pressure vs Density curve using 

mi = [1.6069,1.6069]; 
sigmai = [3.5206,3.5206]; 1 
% epsilon = [2.873268218e-21,2.873268218e-21]; 
epsilon_ki = [191.42,191.42]; 
xi = [0.5,0.5]; % mole fraction of component i 
T = 298.15; % temperature of the system 
% k_bolt = 1.3806488e-23; % Boltzmann constant = 1.3806488 × 10-23 m2 kg s-2 K-1 

% Initial guess for density 

rho_real = 10000; 

rho = rho_real*6.022e-7; 

% Initial Kij 

Kij = 0; 

% Temperatuure-dependent segment diameter di of component i - matrix save 

di = zeros(1,2); 

giihs = zeros (1,2); 

rho_giihs = zeros (1,2); 


m_bar = 0; 

zi0=0; 
zi1=0; 
zi2=0; 
zi3=0; 
for i=1:2 
m_bar=m_bar+((xi(1,i))*(mi(1,i))); 


di(1,i) = (sigmai(1,i))*(1-(0.12*exp(-3*(epsilon_ki(1,i))/T))); 

% Calculating Zieeeeeeh for Hard Sphere compressibility (Rechcecked) 

zi0 = zi0+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^0)); 
zi1 = zi1+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^1)); 
zi2 = zi2+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^2)); 
zi3 = zi3+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^3)); 
end 

for i=1:2 

giihs (1,i)= (1/(1-zi3))+((((di(1,i))^2)/(2*(di(1,i))))*(3*zi2)/((1-zi3)^2))+((di(1,i)^2)/(2*(di(1,i)^2))^2)*((2*(zi2)^2)/((1-zi3)^3)); 

% Derevative of site to site radial distribution function of hard spheres 
% with respect to density 

rho_giihs (1,i) = ((zi3/((1-zi3)^2)))+(((((di(1,i))^2)/(2*(di(1,i)))))*... 
((((3*zi2)/((1-zi3)^2)))+((6*zi2*zi3)/((1-zi3)^3))))+... 
(((((di(1,i))^2)/(2*(di(1,i))))^2)*(((4*(zi2^2))/((1-zi3)^3))+... 
(((6*zi2^2)*zi3)/((1-zi3)^4)))); 
end 


eta = zi3; 

% Calculating the hard sphere compresibility factor (Rechecked) 

zhs = (zi3/(1-zi3)+((3*zi1*zi2)/((zi0*(1-zi3)^2)))+(((3*(zi2^3))-(zi3*(zi2^3)))/(zi0*((1-(zi3)^3))))); 

zhc=0; 

for i = 1:2 
% Calculating the hard chain compressibility factor (Rechecked) 

zhc= zhc+((xi(1,i)*(mi(1,i)-1)*((giihs(1,i))^-1)*((rho_giihs(1,i))))); 
end 
    zhc=m_bar*zhs-zhc; 

% a0 constants 

a0 = [0.9105631445,0.6361281449,2.6861347891,-26.547362491,97.759208784,-159.59154087,91.297774084]; % Checked and Correct Sadowski 2001 

% a1 constants 

a1 = [-0.3084016918,0.1860531159,-2.5030047259,21.419793629,-65.255885330,83.318680481,-33.746922930]; % Checked and Correct Sadowski 2001 

% a2 constants 

a2 = [-0.0906148351,0.4527842806,0.5962700728,-1.7241829131,-4.1302112531,13.776631870,-8.6728470368]; % Checked and Correct Sadowski 2001 

% b0 constants 

b0 = [0.7240946941,2.2382791861,-4.0025849485,-21.003576815,26.855641363,206.55133841,-355.60235612]; % Checked and Correct Sadowski 2001 

% b1 constants 

b1 = [-0.5755498075,0.6995095521,3.8925673390,-17.215471648,192.67226447,-161.82646165,-165.20769346]; % Checked and Correct Sadowski 2001 

% b2 constants 

b2 = [0.0976883116,-0.2557574982,-9.1558561530,20.642075974,-38.804430052,93.626774077,-29.666905585]; 

aim = zeros(1,7); 
bim = zeros(1,7); 

for i = 1:7 
aim(1,i)=a0(1,i)+(((m_bar-1)/m_bar)*a1(1,i))+(((m_bar-1)/m_bar)*((m_bar-2)/m_bar))*a2(1,i); % Checked and Correct Sadowski 2001 (rechecked) 
bim(1,i)=b0(1,i)+(((m_bar-1)/m_bar)*b1(1,i))+(((m_bar-1)/m_bar)*((m_bar-2)/m_bar))*b2(1,i); 


c1 = ((1 + (m_bar*(((8*eta)-(2*eta^2))/(1-eta)^4)+((1-m_bar)*(((20*eta)-(27*eta^2)+(12*eta^3)-(2*eta^4))/(((1-eta)*(2-eta))^2)))))^-1); 

c2 = ((-c1^2)*(m_bar*(((-4*eta^2)+(20*eta)+8)/((1-eta)^5))+((1-m_bar)*(((2*eta^3)+(12*eta^2)-(48*eta)+40)/(((1-eta)*(2-eta))^3))))); 

% Integrals of perturbation theory without/with respect to eta (rechecked) 
dI1 = 0; 
dI2 = 0; 
I1 = 0; 
I2 = 0; 
for j = 1:7 

I1 = I1+ aim(1,j)*eta^(j-1); 
I2 = I2+ bim(1,j)*eta^(j-1); 

dI1 = dI1 + (aim(1,j)*((j-1)+1)*(eta^(j-1))); 
dI2 = dI2 + (bim(1,j)*((j-1)+1)*(eta^(j-1))); 
end 
% Segment abriviations 

m2e = 0; 
m2e2 = 0; 
sigma_ij = zeros(1,2); 
epsilon_ij = zeros (1,2); 
for i = 1:2 
for j = 1:2 
    % Mixing rules for sigma and eta respectively. 
    sigma_ij(i,j)=0.5*(sigmai(1,i)+sigmai(1,j)); 
    epsilon_ij(i,j)=sqrt((epsilon_ki(1,i)*epsilon_ki(1,j))*(1-Kij)); 

    m2e=m2e+(xi(1,i)*xi(1,j)*mi(1,i)*mi(1,j)*((epsilon_ij(i,j)/(T)))*(sigma_ij(i,j)^3)); 
    m2e2=m2e2+(xi(1,i)*xi(1,j)*mi(1,i)*mi(1,j)*(((epsilon_ij(i,j)/(T)))^2)*(sigma_ij(i,j)^3)); 
end 

end 



% The dispersion contribution to the comprehensibility factor 

zdis = (-2*pi*rho*dI1*m2e)-(pi*rho*m_bar*((c1*dI2)+(c2*eta*I2))*m2e2); 

% Compresibility Factor (rechecked) 

z = 1 + zhc + zdis; 

zdis; 
zhc; 

P = z*8.314*T*rho/6.022e-7; % (rechecked) 
P_Mpa = P*1e-6; 

Answer 1

В строке ниже % Initial guess for density: вы установите значение rho_real в 10000. С этого момента, конечно, rho_real будет 10000 независимо от того, что входной аргумент был. Если вы удалите эту строку, ваша функция будет вести себя как ваш скрипт с переменной rho_real, определяемой входным аргументом.

У MATLAB есть способ предоставить значения по умолчанию для аргументов. Ваша функция может проверить, сколько аргументов было предоставлено. Если rho_real было опущено, то есть, если количество аргументов nargin составляет менее 1, только после этого вы устанавливаете rho_real.

if nargin < 1 
    rho_real = 10000 
end