MATLAB Answers

How can I visualise the side view of 3D matrix ?

1 view (last 30 days)
mohammed sait
mohammed sait on 22 Jul 2020
Answered: Swetha Polemoni on 6 Nov 2020
Hello,
I have a 3D matrix that represents the absorbance distribution of scattered laser light in a volume. The absorbance matrix is 3D A(i,j,k). The laser beam propagates in the z direction (corresponding to the k slice). I am able to see the planar distribution of xy direction at each distance z, but I need to visualise the overall propagation by viewing the z-y direction of the matrix. Please give me some suggestions as I am not proficient in MATLAB. The program below will simulate the scattered light distribution in xy at each z segment (please note that running this program takes couple of minutes)
clear all; clc; close all;
pi2 = 2*pi;
rad2deg = 180 /pi ;
deg2rad = pi /180;
ua = 0.114; us=20.7; ut= us+ ua; % optical properties of the medioum, [m^-1]
g = 0.9; % average cosine for scattering function
nw = 1.3; na = 1; % refraction index of water and air
length = 0.3 ; % length of the tank in m
width = .30; % in width of the tank
height = .30; % in height of the tank
Nphot = 1e6; % number of photons
dx=.005; % resolution in x m
dy=.005; % resolution in y m
dz=.005; % resolution in z m
Nx=width/dx; % number of indices in x
Ny=height/dy; % number of indices in y
Nz=length/dz; % number of indices in z
a = 0.005; % spot size of the laser in m
theta_rx =180 * deg2rad; % receiver zenith angle
Phi_tx = 0* deg2rad; % transmitter azimuth angle
Phi_rx = 0* deg2rad; % receiver azimuth angle
theta_div = 10 * deg2rad; % beam divergence angle
FOV = 180 * deg2rad; % field of view of the receiver
thres = 1e-4; % threshold value
m = 10;
ux_rx = sin(theta_rx)*cos(Phi_rx); uy_rx = sin(theta_rx)*sin(Phi_rx); uz_rx =cos(theta_rx); % direction cosine of receiver orientation
A = zeros(Nx,Ny,Nz); % matrix to store the absorbtion
% initiate
theta_tx = 0 * deg2rad; % in degrees
Dp = 0; % detector response, received photon counter
% tx rx initial direction
ux_tx = sin(theta_tx)*cos(Phi_tx); uy_tx = sin(theta_tx)*sin(Phi_tx); uz_tx =cos(theta_tx);
N=0;
while (N< Nphot)
% Photons inital angles
Phi_init = pi2 *rand; theta_init = acos( 1-rand*(1- cos(theta_div/2)))+ theta_tx;
%Phi_init = acos( 1-rand*(1- cos(theta_div/2)))+ Phi_tx; theta_init = acos( 1-rand*(1- cos(theta_div/2)))+ theta_tx;
% photon initial direction
ux = sin(theta_init)*cos(Phi_init); uy = sin(theta_init)*sin(Phi_init); uz =cos(theta_init);
% photon initial position
%x = width/2; y = 0; z =height/2;
x = 0; y = 0; z = 0;
W =1 ;
flag =1;
while (W ~= 0 && flag ==1)
% take a step
s = -log(rand)/ut;
% update the photon position (move)
x = x + s*ux; y = y + s*uy; z = z + s*uz;
%check if photon hit the boundary (ignore for now and terminate when photon is out of bound)
if x >= width/2 || y >= height/2 || z < 0
W = 0;
z = length +1;
end
% check if photon received
if norm([x,y,z]) >= length
if acos( [x,y,z].*[ux_rx,uy_rx,uz_rx] / norm([x,y,z])) <= FOV/2
Dp = Dp + W;
end
flag =0;
end
% Absorb
if(z < length && abs(y) < height/2 && abs(x) < width/2)
W = W*(1- ua/ut) ;
i = floor( (x+ width/2)/dx) + 1;
j = floor( (y + height/2)/dy) + 1;
k = floor( z /dz) + 1;
A(i,j,k) = A(i,j,k) + W;
end
% scatter
costheta = (1+g^2-((1-g^2)/1-g+2*g*rand)^2)*1/2*g;
sintheta = sqrt(1-costheta^2);
psi = pi2*rand;
if (psi < pi)
sinpsi = sqrt(1.0 - cos(psi)^2);
else
sinpsi = -sqrt(1.0 - cos(psi)^2);
end
if(abs(uz) > 0.99999)
uxx = sintheta * cos(psi);
uyy = sintheta * sinpsi;
uzz = sign(uz) * costheta;
else
temp = sqrt(1-uz^2);
uxx = sintheta*(ux*uz*cos(psi)-uy*sinpsi)/temp + ux*costheta;
uyy = sintheta*(uy*uz*cos(psi)+ux*sinpsi)/temp + uy*costheta;
uzz = -sintheta*cos(psi)*temp + uz*costheta;
end
ux = uxx;
uy = uyy;
uz = uzz;
% terminate roulete
if(W < thres)
if(rand > 1/m)
W = 0;
else
W = m * W;
end
end
end
N = N + 1;
end
A = A /(dz*dx*dy*Nphot*ua);
xx = -width/2:dx:width/2;
yy = -height/2:dy:height/2;
zz = 0:dz:length;
imagesc(zz,yy,A(:,:,60))
% colormap(blue);
colorbar;

  0 Comments

Sign in to comment.

Answers (1)

Swetha Polemoni
Swetha Polemoni on 6 Nov 2020
Hi.
It is my understanding that you want visualize side view of your 3d matrix after plotting. "View" command in matlab can be used to visualize 3d plots. Plotting in 3d can be done using "surf" command. Refer this example for better understanding.

  0 Comments

Sign in to comment.

Tags

Community Treasure Hunt

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

Start Hunting!