clear all; clc; close all;
pi2 = 2*pi;
rad2deg = 180 /pi ;
deg2rad = pi /180;
ua = 0.114; us=20.7; ut= us+ ua;
g = 0.9;
nw = 1.3; na = 1;
length = 0.3 ;
width = .30;
height = .30;
Nphot = 1e6;
dx=.005;
dy=.005;
dz=.005;
Nx=width/dx;
Ny=height/dy;
Nz=length/dz;
a = 0.005;
theta_rx =180 * deg2rad;
Phi_tx = 0* deg2rad;
Phi_rx = 0* deg2rad;
theta_div = 10 * deg2rad;
FOV = 180 * deg2rad;
thres = 1e-4;
m = 10;
ux_rx = sin(theta_rx)*cos(Phi_rx); uy_rx = sin(theta_rx)*sin(Phi_rx); uz_rx =cos(theta_rx);
A = zeros(Nx,Ny,Nz);
theta_tx = 0 * deg2rad;
Dp = 0;
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)
Phi_init = pi2 *rand; theta_init = acos( 1-rand*(1- cos(theta_div/2)))+ theta_tx;
ux = sin(theta_init)*cos(Phi_init); uy = sin(theta_init)*sin(Phi_init); uz =cos(theta_init);
x = 0; y = 0; z = 0;
W =1 ;
flag =1;
while (W ~= 0 && flag ==1)
s = -log(rand)/ut;
x = x + s*ux; y = y + s*uy; z = z + s*uz;
if x >= width/2 || y >= height/2 || z < 0
W = 0;
z = length +1;
end
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
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
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;
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))
colorbar;
0 Comments
Sign in to comment.