How to run a for loop with two functions (nested for loop?)

I am trying to perform a series of iterative calculations given two inputs. First I want to calculate for three pressures. Second, I want to calculate over a range of values from 0:0.01:1.0. I am using two functions: one from the XSteam.m package and one user-defined function. I am getting an error with the second function but I don't really know how to resolve this. First is the input in the script.

P = [50.0, 100.0, 150.0];     % Pressure (bars)
alpha = 0.0:0.01:1.0;         % Void fraction from 0 to 1 in increments of 0.01
dh = 0.01;                    % Pipe diameter in meters (1 cm = 0.01 m)

for i = 1:numel(P);
    sigma(i) = XSteam('st_p',P(i));
    rhof(i) = XSteam('rhoL_P',P(i));
    rhov(i) = XSteam('rhoV_P',P(i));
    for n = 1:numel(alpha);
        a(n) = calc_a( alpha, sigma(i), rhof(i), rhov(i), dh );
    end
end

The values with XSteam are correctly generated if I put them in their own loop but I am unable to get the second function working when running through a for loop. Here is the calc_a function

function [ a ] = calc_a( alpha, sigma, rhof, rhov, dh )
g = 9.81;
db = sqrt(sigma /(g*(rhof - rhov)));
if alpha < 0.5            
    a = (6*alpha)/ db
elseif alpha >0.8 
    a = (4*sqrt(alpha))/ dh
else  
    omega = (alpha - 0.5)/(0.8 - 0.5)
    a = ((6*alpha)/ db)*(1-omega) + ((4*alpha)/ dh)*omega
end
end

Is there a better way to do this or what am I missing? Thanks for the help.

[SALAH ALRABEEI]

I think u might have a mistake in alpha in the 2nd loop. If Iam correct you should send one value of apha at foe each iteration. So it is alpha(i) not alpha

Answer 1

[DGM]

Pretty sure that you're just intending to process alpha one element at a time, so just do that. Of course, you're overwriting a each time the outer loop cycles, so you'll want to fix that too.

P = [50.0, 100.0, 150.0];     % Pressure (bars)
alpha = 0.0:0.01:1.0;         % Void fraction from 0 to 1 in increments of 0.01
dh = 0.01;                    % Pipe diameter in meters (1 cm = 0.01 m)

for m = 1:numel(P)
    sigma(m) = XSteam('st_p',P(m));
    rhof(m) = XSteam('rhoL_P',P(m));
    rhov(m) = XSteam('rhoV_P',P(m));
    for n = 1:numel(alpha)
        a(m,n) = calc_a( alpha(n), sigma(m), rhof(m), rhov(m), dh );
    end
end
% user-defined function is unchanged

The output is now a matrix, with one row for each value of P. If you leverage implicit array expansion available since R2016b, you can avoid the inner loop.

P = [50.0, 100.0, 150.0];     % Pressure (bars)
alpha = 0.0:0.01:1.0;         % Void fraction from 0 to 1 in increments of 0.01
dh = 0.01;                    % Pipe diameter in meters (1 cm = 0.01 m)

M = numel(P);
sigma = zeros(M,1);
rhof = zeros(M,1);
rhov = zeros(M,1);
for m = 1:M
    % i guess xsteam doesn't handle vector inputs
    sigma(m) = XSteam('st_p',P(m));
    rhof(m) = XSteam('rhoL_P',P(m));
    rhov(m) = XSteam('rhoV_P',P(m));
end

% this is everything that the calc_a() used to do
% but now it's vectorized
g = 9.81;
db = sqrt(sigma ./(g*(rhof - rhov)));
omega = (alpha - 0.5)/(0.8 - 0.5);
m1 = alpha<0.5;
m2 = alpha>0.8;

% this should be identical to a
A = ((6*alpha)./ db).*m1 + ((4*sqrt(alpha))/ dh).*m2 ...
    + (((6*alpha)./ db).*(1-omega) + ((4*alpha)/ dh).*omega).*(~m1 & ~m2);

Whichever method you use is up to you.

[Kevin Burg]

Thanks. I actually got to your first solution while messing around with it some more but the second method is interesting since I wasn't aware that could be done.

[DGM]

All that's really necessary for the second method to work is that sigma, etc be oriented orthogonal to alpha. Then when elementwise arithmetic is performed on them, the result is a matrix. The legacy way of doing the same is to use bsxfun().