3 views (last 30 days)

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]')

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.

.

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:

- either you hold on and plot point by point, or
- 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
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...

dpb
on 25 Oct 2020 at 22:11

Opportunities for recent engineering grads.

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

Start Hunting!
## 2 Comments

## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/625633-plotting-a-graph-with-different-conditions-using-if-statement#comment_1085188

⋮## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/625633-plotting-a-graph-with-different-conditions-using-if-statement#comment_1085188

## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/625633-plotting-a-graph-with-different-conditions-using-if-statement#comment_1085213

⋮## Direct link to this comment

https://uk.mathworks.com/matlabcentral/answers/625633-plotting-a-graph-with-different-conditions-using-if-statement#comment_1085213

Sign in to comment.