function [rhoPort, tauPort, alphaPort, tauDirect] = simulPort(rho, xi, aP, aS, N)
% Input parameters:
% rho --- surface reflectivity
% xi --- port aspect ratio
% aP --- port radius
% aS --- sphere radius
% N --- number of rays
% Output:
% rhoPort --- port reflectance
% tauPort --- port transmittance
% alphaPort --- port absorptance
% tauDirect --- port transmittance (straight-through)
MAXITER = 1e4; % Just in case, to avoid getting stuck in an endless loop
hP = 2 * xi * aP; % Port height
% Bottom and top face coordinates. Sphere is centred at O.
zBot = realsqrt(aS^2-aP^2);
zTop = zBot + hP;
% Trivial thin-port case
if hP == 0
rhoPort = 0;
alphaPort = 0;
tauPort = 1;
tauDirect = 1;
return
end
% Initializing ray counters
nbEscape = 0;
nbReturn = 0;
nbAbsorbed = 0;
nbDirect = 0;
% Loop over rays
for nRay = 1:N
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find initial ray entering a cylindrical port %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Take random point on spherical cap inside port, this will be a ray
% position cos-distributed between cos theta_P and 1
theta = acos(1-(1 - cos(asin(aP/aS)))*rand);
phi = 2 * pi * rand;
xP = aS * sin(theta) * cos(phi);
yP = aS * sin(theta) * sin(phi);
zP = aS * cos(theta);
% Take random point A on sphere surface (outside cap)
thetaA = acos(-1+(1 + cos(asin(aP/aS)))*rand);
phiA = 2 * pi * rand;
xA = aS * sin(thetaA) * cos(phiA);
yA = aS * sin(thetaA) * sin(phiA);
zA = aS * cos(thetaA);
% Direction vector
dAP = sqrt((xP - xA)^2+(yP - yA)^2+(zP - zA)^2);
uX = (xP - xA) / dAP;
uY = (yP - yA) / dAP;
uZ = (zP - zA) / dAP;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main loop over reflections inside the port %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nIter = 0;
bAlive = true;
while (nIter < MAXITER) && bAlive
% Find intersection point on cylindrical port (assuming we are
% inside the finite cylinder)
% First intersect with infinite long cylinder
ca = uX^2 + uY^2;
cb = uX * xP + uY * yP;
cc = xP^2 + yP^2 - aP^2;
t1 = (-cb + sqrt(cb^2 - ca * cc)) / ca;
zP = zP + uZ * t1; % Coordinate of intersect (will be new ray position)
if (zP >= zTop); % Escape through top face (detector)
nbEscape = nbEscape + 1;
if nIter == 0
nbDirect = nbDirect + 1;
end
bAlive = false;
elseif (zP <= zBot); % Returns through bottom face (sphere)
nbReturn = nbReturn + 1;
bAlive = false;
else % Intersection must be on side of finite cylinder
% Test first if ray is absorbed by wall
if (rand > rho) % absorbed rays
nbAbsorbed = nbAbsorbed + 1;
bAlive = false;
else % If not, implement random Lambertian reflection
% Intersection point (update ray position)
xP = xP + uX * t1;
yP = yP + uY * t1;
% Inward normal for Lambertian reflection
nx = -xP / aP;
ny = -yP / aP;
% nz = 0;
% Choose point randomly on unit sphere
rnphi = 2 * pi * rand; % 0