How can I visualise the side view of 3D matrix ?
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
Respuestas (1)
Ver también
Categorías
Más información sobre Photonics en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!