MATLAB Answers

plotting a graph with different conditions using if statement

3 views (last 30 days)
Louisa Mezache
Louisa Mezache on 25 Oct 2020
Commented: dpb on 25 Oct 2020 at 22:11
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]')

  2 Comments

dpb
dpb on 25 Oct 2020
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
Louisa Mezache on 25 Oct 2020
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.

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 25 Oct 2020
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.
.

  2 Comments

Louisa Mezache
Louisa Mezache on 25 Oct 2020 at 20:35
Thank you, Star Strider! It looks like based on your and others' answers, I need to better understand vectorization. Thank you for your help.

Sign in to comment.

More Answers (2)

Anton Casas
Anton Casas on 25 Oct 2020
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

dpb
dpb on 25 Oct 2020
Edited: dpb on 25 Oct 2020
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...

  2 Comments

Louisa Mezache
Louisa Mezache on 25 Oct 2020 at 20:41
Thank you, dpb! I understand what you did there. That's very helpful.
dpb
dpb on 25 Oct 2020 at 22:11
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.

Sign in to comment.

Products


Release

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!