Plotting inside for loop using if

Hello, I am trying to plot a formula in which if the input does not follow the constraint, calculation process change.I am trying to do this with if inside a for loop.My script is below.It works for single x values but It does not work for all the x values from 0 to 1 with increments of 0.01.Could anyone please fix the error I do here?

G=200;                 %mass flux[kg/m^2*s]
Pred=0.31052;          %reduced pressure
M=44.01;               %moleculer weight
q=10000;               %heat flux
gv=60.72775526;        %density of vapor phase
gl=1008.003926;        %density of liquid phase
g=9.81;                %gravitational acceleration
sig=7.53*10^-3;        %surface tension
vv=1.34741*10^-5;      %viscosity of vapor phase
vl=1.28*10^-4;         %viscosity of liquid phase
kl=0.128582692;        %conductivity of liqud phase
kv=0.0137;             %conductivity of vapor phase
hlv=236272.0191;       % enthalpy of vaporization
D=0.0061;              %inside diameter
Pr=2.2231;             % prandtl number for liquid flow
cpv=1387.633317;       %specific heat of vapor phase
for x=0.1:0.01:0.9;    %dry coefficient[-]
  Bo=q./(G.*hlv);
  Bd=(g.*(gl-gv).*D.^2)./sig;
  Rel=(G.*(1.-x).*D)./vl;
  x_critical=38.27.*Rel.^2.12.*(1000.*Bo).^1.64.*Bd.^-4.7;
  if x< x_critical 
    E=(1.+9.36.*10^3.*x.*Pr.*((gl./gv)-1)).^0.11;
    S=1./(1+0.00000162.*E.^0.69.*Rel.^1.11);
    hl=0.023.*Rel.^0.8.*Pr.^0.4.*(kl./D);
    hnb=55.*Pred.^0.12.*(-log10(Pred)).^-0.55.*M.^-0.5.*q.^0.67;
    htp_yoon=((S*hnb).^2+(E.*hl).^2).^0.5;
    plot(x,htp_yoon)
    hold on
  else
    Rel=(G.*(1-x).*D)./vl;
    Rev=(G.*x.*D./vv);
    hv=0.023.*Rev.^0.8.*(cpv.*vv./kv).^(0.4).*(kv./D);
    E=1+3000.*Bo.^0.86+1.12.*(x./(1-x)).^0.75.*(gl./gv).^0.41;
    hl=0.023.*Rel.^0.8.*Pr.^0.4.*(kl./D);
    Xtt=((1.-x)./x).^0.9.*(gv./gl).^0.5.*(vl./vv).^0.1;
    theta=2*pi*(36.23.*Rel.^3.47.*Bo.^4.84.*Bd.^-0.27.*(1./Xtt).^2.6);
    hwet=E.*hl;
    htp_yoon=((theta.*hv)+(2*pi-theta).*hwet)./(2*pi);
    plot(x,htp_yoon)
    hold on
  end
end

[dpb]

Doesn't appear to be anything untoward in the code; what, specifically, do you think is a problem?

Use the debugger and step through at the point in question...

[Jan]

@Berkay Cinar: Just a note: 9.36e3 is a cheap constant while 9.36.*10^3 is an expensive power operation and a multiplication.

Please mention, what the problem is. How do you see that "It does not work for all the x values"?

The code does not work for "all the x values from 0 to 1" because the loop runs over "0.1:0.01:0.9" only.

[dpb]

Good catch on the loop limits, Jan...didn't notice that.

[KSSV]

Try plot(x,htp_yoon,'.')

[dpb]

Oh, yeah, that too...he'd be better off building the whole output then plotting instead of inside the loop...

[Berkay Cinar]

After x=0.74 x is bigger than x critical so the calculation process should change.When I give single x values results are;

• x=0.1 5474
• x=0.2 5469
• x=0.3 5412
• x=0.4 5334
• x=0.5 5250
• x=0.6 5168
• x=0.7 5096
• x=0.74 3745
• x=0.8 3666
• x=0.9 3461but I need to get this result to a plot and when I plot it I think the calculation does not change when the x is higher than x critical.For every value of x I want to calculate htp_yoon according to x crit and record every result then plot everything.But I have no idea how to do that.Simple plotting at the end of the loop seems not to work.I tried your suggestion dpb.Thank you very much for your time but I dont understand why your script gives different result than mine. You got constant 5000 until 0.9 I have the results above.

[dpb]

Hmm....not sure where the code changed from what I copied prior time, but I also get that for the breakpoint now...

But, your plotting problem has to do with how you've structured the loop; by having only one point at a time, all you get is a dot on the screen that isn't visible because it isn't large enough. To see this, try changing your current plotting instruction to

plot(x,htp_yoon,'x','color','k')

and you'll see that before it did plot, you just couldn't seen the tiny little one-pixel dots.

But, the easier way is still to save up the computed array and plot in the end; you'll probably want those values for other purposes, anyway, won't you?

Anyway, the general idea is

1. make the range x values a vector and iterate over it in counted loop instead of using it as the loop value
2. create array for htp_yoon to match x
3. modify all references to variable x to use the subscripted value instead of the whole array reference
4. move the plotting to when you're otherwise done...

I fixed up the code in Answer to match; what I forgot before was that x was only a single value instead of the vector because you had used it for the iterator...

Answer 1

[dpb]

Reorganized your code to do as suggested above...

G=200;                 %mass flux[kg/m^2*s]
Pred=0.31052;          %reduced pressure
M=44.01;               %moleculer weight
q=10000;               %heat flux
gv=60.72775526;        %density of vapor phase
gl=1008.003926;        %density of liquid phase
g=9.81;                %gravitational acceleration
sig=7.53*10^-3;        %surface tension
vv=1.34741*10^-5;      %viscosity of vapor phase
vl=1.28*10^-4;         %viscosity of liquid phase
kl=0.128582692;        %conductivity of liqud phase
kv=0.0137;             %conductivity of vapor phase
hlv=236272.0191;       % enthalpy of vaporization
D=0.0061;              %inside diameter
Pr=2.2231;             % prandtl number for liquid flow
cpv=1387.633317;       %specific he
% make x values array 'cuz need 'em to plot with later...
x=0.1:0.01:0.9;        %dry coefficient[-]
L=length(x);              % find out how many there are
htp_yoon=zeros(size(x));  % preallocate for the output
for ix=1:L                % loop over thos values
  Bo=q./(G.*hlv);
  Bd=(g.*(gl-gv).*D.^2)./sig;
  Rel=(G.*(1.-x(ix)).*D)./vl;
  x_critical=38.27.*Rel.^2.12.*(1000.*Bo).^1.64.*Bd.^-4.7;
  if x(ix)<x_critical 
    E=(1.+9.36.*10^3.*x(ix).*Pr.*((gl./gv)-1)).^0.11;
    S=1./(1+0.00000162.*E.^0.69.*Rel.^1.11);
    hl=0.023.*Rel.^0.8.*Pr.^0.4.*(kl./D);
    hnb=55.*Pred.^0.12.*(-log10(Pred)).^-0.55.*M.^-0.5.*q.^0.67;
    htp_yoon(ix)=((S*hnb).^2+(E.*hl).^2).^0.5;
  else
    %Rel=(G.*(1-x).*D)./vl;  % NB:  Duplicated, not needed here
    Rev=(G.*x(ix).*D./vv);
    hv=0.023.*Rev.^0.8.*(cpv.*vv./kv).^(0.4).*(kv./D);
    E=1+3000.*Bo.^0.86+1.12.*(x(ix)./(1-x(ix))).^0.75.*(gl./gv).^0.41;
    hl=0.023.*Rel.^0.8.*Pr.^0.4.*(kl./D);
    Xtt=((1.-x(ix))./x(ix)).^0.9.*(gv./gl).^0.5.*(vl./vv).^0.1;
    theta=2*pi*(36.23.*Rel.^3.47.*Bo.^4.84.*Bd.^-0.27.*(1./Xtt).^2.6);
    hwet=E.*hl;
    htp_yoon(ix)=((theta.*hv)+(2*pi-theta).*hwet)./(2*pi);
  end
  %plot(x(ix),htp_yoon(ix),'x','color','k')  % same data variable
  %if ix==1,hold all, end  % either case co can plot here point wise if really want
end
plot(x,htp_yoon)  % now plot 'em all...

There's still a really big discontinuity at the breakpoint though...that expected?

[Berkay Cinar]

Thank you for your time dpb.That seems to do the trick.At the breakpoint it is supposed to be like that. Have a good day:)

[dpb]

OK, glad to help.

Now that have it; note that can refactor the above and eliminate the loop if you make the two cases into internal subroutines.

isCrit=(x>x_critical=38.27.*Rel.^2.12.*(1000.*Bo).^1.64.*Bd.^-4.7);  % logical vector
htp_yoon( isCrit)=htpCrit(xx(isCrit));
htp_yoon(~isCrit)=htpNonCrit(xx(~isCrit));

where the two functions are the code for the two if branches. You wrote them already to accept vector operations with the "dot" operators altho don't need for single elements, simply remove the subscript reference and pass the elements via the logical addressing.