plotting a graph with different conditions using if statement

Hello,

I'm attempting to plot a graph based on different functions of x. I've attempted an if else statement for the different ranges of x, but nothing shows up. I'm not sure why or how else I should try to plot this. Below is my code. Please let me know if there's anything I should clarify.

Thank you!

Best,

Louisa

%clear
%clc

mu_a1 = 0.1; %absorption coefficient for 1D geometry in mm^-1
sigma = 1.0; %backscattering coefficient in mm-1
D = 5.0; %thickness of tissue in mm
E = 1; % incident laser irradiance in W/mm^2
m = (mu_a1 + sigma)/sigma;
b = sqrt(m^2 - 1);

x = 0;
C = 1;
if (x >= 0) && (x <= C)
    x = x + 0.1;
    R = sinh(b*sigma*x) / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %reflectance
    T = m / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %transmittance
    phi = 2*E/(1-R^2);
elseif C < x < (C + D)
    x = x + 0.1;
    R = sinh(b*sigma*x) / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %reflectance
    T = m / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %transmittance
    phi = (E/(1-R)) + (E*T/(1-R));
elseif  x >= (C + D)
    x = x + 0.1;
    R = sinh(b*sigma*x) / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %reflectance
    T = m / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %transmittance
    phi = E*T/(1-R);
end
plot(x,phi,'b')
xlabel('X [mm]')
ylabel('Fluence Rate [W/mm^2]')

[dpb]

There's nothing to repeat the if...end block so all that happens is you calculate the x=0; condition, increment x, but then exit.

You would need to iterate over the desired range of x and save the computed x,phi values in arrays to plot.

[Louisa Mezache]

Thank you for your response. Does this mean my 'x = x + 0.1' is in the wrong position? I'm not sure I know how to same the computed x,phi values in an array. Should this be done within each condition of the if statement? I'm new to MATLAB, so I appreciate your help.

Answer 1

[Star Strider]

Your code can be made significantly more efficient by using vectorisation and ‘logical indesing’:

mu_a1 = 0.1; %absorption coefficient for 1D geometry in mm^-1
sigma = 1.0; %backscattering coefficient in mm-1
D = 5.0; %thickness of tissue in mm
E = 1; % incident laser irradiance in W/mm^2
m = (mu_a1 + sigma)/sigma;
b = sqrt(m^2 - 1);
C = 1;

x = linspace(0, 1, 100);

R = sinh(b*sigma*x) ./ (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %reflectance
T = m ./ (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %transmittance

phifcn = @(x) (2*E./(1-R.^2)).*((x >= 0) & (x <= C)) + ((E./(1-R)) + (E*T./(1-R))).*((C < x) & (x < (C + D))) + (E*T./(1-R)).*(x >= (C + D));

figure
plot(x,phifcn(x),'b')
xlabel('X [mm]')
ylabel('Fluence Rate [W/mm^2]')
grid

producing:

See: Array vs. Matrix Operations for the vectorisation expalanation, and Matrix Indexing to understand how logical indexing works. The if block is not necessary.

.

[Louisa Mezache]

Thank you, Star Strider! It looks like based on your and others' answers, I need to better understand vectorization. Thank you for your help.

[Star Strider]

As always, my pleasure!

Answer 2

[Anton Casas]

First of all, you need vectors on x and phi to plot a normal graph. The way you are using your variables, x and phi are numbers, so you are only plotting a point.

Even efore that, as dpb said, the main problem in your code is that there are not any loops, so the evaluation of the variables only occur once.

Here you have 2 options:

1. either you hold on and plot point by point, or
2. you save x and phi values on vectors and plot them at the end.

IMO, the second one is much cleaner, but requires a little bit more of modification on your code. Judging from your code (apologies for that), you do not have much understanding of vectors, so that could be slightly more difficult for you.

The first option would look somehow like this (I didn't execute it so I can't assure it will work 100%):

for x = 0:0.1:10 % This way the plot will end in x = 10
    C = 1;
    if (x >= 0) && (x <= C)
        % x = x + 0.1; Updates of x are now done by the loop
        R = sinh(b*sigma*x) / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %reflectance
        T = m / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %transmittance
        phi = 2*E/(1-R^2);
    elseif C < x < (C + D)
        % x = x + 0.1;
        R = sinh(b*sigma*x) / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %reflectance
        T = m / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %transmittance
        phi = (E/(1-R)) + (E*T/(1-R));
    elseif  x >= (C + D)
        % x = x + 0.1;
        R = sinh(b*sigma*x) / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %reflectance
        T = m / (m*sinh(b*sigma*x) + b*cosh(b*sigma*x)); %transmittance
        phi = E*T/(1-R);
    end
    plot(x,phi,'b+') % You may want to change + to other symbol
    hold on % This will prevent your plot from clearing
    xlabel('X [mm]')
    ylabel('Fluence Rate [W/mm^2]')
end

[Louisa Mezache]

Thank you, Anton! This was helpful.

Answer 3

[dpb]

What's the desired range for X?

Read the "Getting Started" section on array indexing and operations...

In short, in MATLAB if you write something like

xStart=0; xEnd=10;
x=linspace(xStart,xEnd);  % generates vector of 100 elements between xStart,xEnd, inclusive
R = sinh(b*sigma*x)./(m*sinh(b*sigma*x) + b*cosh(b*sigma*x));   %reflectance
T = m ./ (m*sinh(b*sigma*x) + b*cosh(b*sigma*x));               %transmittance

you'll have two vectors R, T of the same length as x over those ranges. MUCH simpler than writing the loop.

Then, there's the conditional on the phi parameter to deal with...that's a little more effort, but not a lot--use logical addressing.

isLE_C=(x <= C);                        % the first conditional logical addressing vector
phi(isLE_C) = 2*E./(1-R(isLE_C).^2);    % compute those -- NB: the "dot" operators
isBT_C_CD = (x>C) & (x<(C+D));          % ML requires two logical expressions
phi(isBT_C_CD) = E./(1-R(isBT_C_CD)) + E(isBT_C_CD).*T(isBT_C_CD)./(1-R(isBT_C_CD));
isGT_CD = (x >= (C + D));
phi(isGT_CD) = E(isGT_CD).*T(isGT_CD)./(1-R(isGT_CD));

Then you have vectors x,phi to plot...

NB: Air code; may be some typos; may have missed a "dot" somewhere, but is the idea...

[Louisa Mezache]

Thank you, dpb! I understand what you did there. That's very helpful.

[dpb]

Star S's "trick" of using the summation of the pieces shortens the overall code but is the same result -- the above is explicit to illustrate the manipulation by section.