function [R,brdf,scat,shad,w] = goa2D(f,x,thd,n1,k1,n2,k2,nfrp,dstrb,br)
%
% [R,brdf,scat,shad,w] = goa2D(f,x,thd,n1,k1,n2,k2,nfrp,dstrb,br)
%
% calculates the directional-hemispherical reflectance for light
% incident from a dielectric medium with complex refractive index n1 + i*k1 on an
% opaque 1-d rough surface y=f(x) with complex refracttive index n2 + i*k2 with
% the angle of incidence thd (in degrees) using the geometric optics approximation (GOA)
% with nfrp first reflection points. Beam radius is set through br parameter
% while beam energy distribution is controlled via parameter dstrb ('uni' or 'gauss').
%
% Input f - surface height
% x - surface point
% thd - (global) angle of incidence in degrees
% n1 - refractive index of medium 1
% k1 - extinction coefficient of medium 1
% n2 - refractive index of medium 2
% k2 - extinction coefficient of medium 2
% nfrp - number of first reflection points (rays)
% dstrb - beam energy distribution (type 'uni' for uniform and 'gauss'
% for gaussian)
% br - beam radius
%
% Output: R - directional-hemispherical reflectance
% brdf - bidirectional reflectance distribution function resolved
% into first, second and higher order scattering
% scat - mean number of scattering events per incident ray
% shad - amount of shadowing (percentage)
% w - percentage of rays removed due to warning signals given
% during the simulation (indication of too few surface
% points)
% Last updated: 2010-09-30 (David Bergström)
%
tic
format long
nsp = length(x); % number of surface points
DX = (x(end)-x(1))/length(x); % sample distance
% check input
if nfrp > nsp
error('Too many rays! Unable to fit %g rays on surface with %g surface points.',nfrp,nsp);
end;
if br > (x(end) - x(1))
br = x(end) - x(1); % fit beam on surface
%error('Beam radius is too large! Unable to fit beam of radius %g on surface with length %g!',br,x(length(x)) - x(1));
end;
% initializations
theta = thd*pi/180; % AOI in radians
drcscat = zeros(3,181); % multiple scattering resolved differential reflection coefficient (1st order, 2nd, higher)
totalenergy = 0; % total incident energy on surface
bo = 0; % number of rays out-of-bounds
shp = zeros(1,nsp); % number of shadowed points
scatpts = zeros(1,nsp); % number of scattering points per incoming ray (1st refl. point)
w = zeros(1,nsp); % warning signals
% calculate slopes from central difference
ddx = conv([1 0 -1],f)/2/DX; ddx = [0 ddx(3:end-2) 0]; % (set slopes at border to 0)
% surface normals calculation
nx = -ddx; ny = 1;
norm = sqrt(nx.^2+ny.^2); % norm of surface normal
normals(:,1) = nx ./ norm; %
normals(:,2) = ny ./ norm; % normalizations
% mirrored surface
xM = x; fM = fliplr(f); normalsM = flipud(normals);
normalsM(:,1) = -normalsM(:,1); normalsM(:,2) = normalsM(:,2); % normals reversed
%%%%% Ray tracing algorithm starts here %%%%%
for frpi = 1:nfrp ; % loop over first reflection points using index
reflpoint = round(nsp/nfrp*frpi); % translate index to first reflection point
incomingray = [sin(theta) -cos(theta)]; % positive theta means from the left (going in positive x)
if shadow(x,f,reflpoint,theta) == 0 % Test for shadowing
scatpts(reflpoint) = 1; % scattering event found
incomingenergy = incidentenergy(incomingray,normals,x,ddx,reflpoint,dstrb,br); % calculation of energy of incident ray
incenergy = incomingenergy; % incidentenergy stored
totalenergy = totalenergy + incomingenergy; % add to total incident energy
% transform incident ray and energy to scattered ray and energy
[outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normals,reflpoint,n1,k1,n2,k2);
if outgoingenergy > incomingenergy
w(reflpoint) = 1; % report warning signal
outgoingenergy = incomingenergy;
outgoingray = incomingray;
end;
scattered = 0;
% check if new reflection point is available
nextreflpoint = findnewpoint(x,f,reflpoint,outgoingray);
if nextreflpoint == 0 && w(reflpoint) == 0 % no new point found
scattering_direction = scatdir(outgoingray);
if scattering_direction ~= -100 % definetely scattered
scattered = 1;
else % out-of-bounds
bo = bo + 1;
nextreflpoint = findmirrorpoint(x,f,xM,fM,reflpoint,outgoingray); % find new reflection point on mirrored surface
scatpts(reflpoint) = scatpts(reflpoint) + 1; % scattering event found on mirrored surface
incomingenergy = outgoingenergy; incomingray = outgoingray;
[outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normalsM,nextreflpoint,n1,k1,n2,k2); % transform incident ray and energy to scattered ray and energy
if outgoingenergy > incomingenergy
w(reflpoint) = 1; % report warning signal
outgoingenergy = incomingenergy;
outgoingray = incomingray;
end;
rayreturn = 0;
newpoint = findnewpoint(xM,fM,nextreflpoint,outgoingray);
while rayreturn == 0 && scattered == 0 && w(reflpoint) == 0% until ray either scattered or back to original surface
if newpoint == 0 % no new point found
scattering_direction = scatdir(outgoingray);
if incomingray(1)*outgoingray(1) < 0 % changing directions (i.e. returning to original surface)
rayreturn = 1;
nextreflpoint = findmirrorpoint(xM,fM,x,f,nextreflpoint,outgoingray);
if nextreflpoint == 0 % no new point found on original surface
scattering_direction = scatdir(outgoingray);
scattered = 1;
end;
else % finally scattered
scattered = 1;
end;
else
scatpts(reflpoint) = scatpts(reflpoint) + 1; % scattering event found on mirrored surface
incomingenergy = outgoingenergy; incomingray = outgoingray;
[outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normalsM,newpoint,n1,k1,n2,k2);
if outgoingenergy > incomingenergy
w(reflpoint) = 1; % report warning signal
outgoingenergy = incomingenergy;
outgoingray = incomingray;
end;
newpoint = findnewpoint(xM,fM,nextreflpoint,outgoingray);
end;
end;
end;
end;
while scattered == 0 && w(reflpoint) == 0
scatpts(reflpoint) = scatpts(reflpoint) + 1; % scattering event found
incomingenergy = outgoingenergy; incomingray = outgoingray;
% transform incident ray and energy to scattered ray and energy
[outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normals,nextreflpoint,n1,k1,n2,k2);
% error for grazing incident rays (decrease error by increasing number of
% surface points)
if outgoingenergy > incomingenergy
w(reflpoint) = 1; % report warning signal
outgoingenergy = incomingenergy;
outgoingray = incomingray;
end;
% check if new reflection point is available
np = findnewpoint(x,f,nextreflpoint,outgoingray);
if np == 0 && w(reflpoint) == 0 % no new point found
scattering_direction = scatdir(outgoingray);
if scattering_direction ~= -100 % definetely scattered
scattered = 1;
else % out-of-bounds (to mirrored surface)
bo = bo + 1;
nextreflpoint = findmirrorpoint(x,f,xM,fM,nextreflpoint,outgoingray);
scatpts(reflpoint) = scatpts(reflpoint) + 1; % scattering event found on mirrored surface
incomingenergy = outgoingenergy; incomingray = outgoingray;
[outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normalsM,nextreflpoint,n1,k1,n2,k2);
if outgoingenergy > incomingenergy
w(reflpoint) = 1; % report warning signal
outgoingenergy = incomingenergy;
outgoingray = incomingray;
end;
rayreturn = 0; originalray = incomingray;
newpoint = findnewpoint(xM,fM,nextreflpoint,outgoingray);
while rayreturn == 0 && scattered == 0 && w(reflpoint) == 0% until ray either scattered or back to original surface
if newpoint == 0 % no new point found
scattering_direction = scatdir(outgoingray);
if originalray(1)*outgoingray(1) < 0 % changing directions (i.e. returning to original surface)
rayreturn = 1;
nextreflpoint = findmirrorpoint(xM,fM,x,f,nextreflpoint,outgoingray);
if nextreflpoint == 0 % no new point found on original surface
R = -1; %
scat = -1; % error handling
shad = -1; %
return
end;
else % finally scattered
scattered = 1;
end;
else
scatpts(reflpoint) = scatpts(reflpoint) + 1; % scattering event found on mirrored surface
incomingenergy = outgoingenergy; incomingray = outgoingray;
[outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normalsM,newpoint,n1,k1,n2,k2);
if outgoingenergy > incomingenergy
w(reflpoint) = 1; % report warning signal
outgoingenergy = incomingenergy;
outgoingray = incomingray;
end;
newpoint = findnewpoint(xM,fM,newpoint,outgoingray);
end;
end;
end;
else
nextreflpoint = np;
end;
end;
if scattered == 1 && w(reflpoint) == 0% add to differential reflection coefficient
if scatpts(reflpoint) == 1
drcscat(1,round(scattering_direction*180/pi)+91) = drcscat(1,round(scattering_direction*180/pi)+91) + outgoingenergy; % first-order scattering
else
if scatpts(reflpoint) == 2
drcscat(2,round(scattering_direction*180/pi)+91) = drcscat(2,round(scattering_direction*180/pi)+91) + outgoingenergy; % second-order scattering
else drcscat(3,round(scattering_direction*180/pi)+91) = drcscat(3,round(scattering_direction*180/pi)+91) + outgoingenergy; % multiple-order scattering (> 2)
end;
end;
elseif w(reflpoint) == 1 % remove ray from calculation
totalenergy = totalenergy - incenergy;
scatpts(reflpoint) = 0;
end;
else shp(reflpoint) = shp(reflpoint) + 1; % first reflection point was shadowed
end;
end; % no more first reflection points (for)
%%%%%% RESULTS %%%%%%
drcscat = drcscat / totalenergy;
drc = sum(drcscat); % differential reflection coefficient
R = sum(drc); % directional-hemispherical reflectance
w = sum(nonzeros(w))/nfrp; % percentage of rays removed due to warnings
%brdf calculation
brdf(1,:) = drcscat(1,:)*180./cos((-90:90)*pi/180); % first order
brdf(2,:) = drcscat(2,:)*180./cos((-90:90)*pi/180); % second order
brdf(3,:) = drcscat(3,:)*180./cos((-90:90)*pi/180); % third and higher order
scat = mean(nonzeros(scatpts)); % mean number of scattering points per ray
shad = sum(nonzeros(shp))/(nfrp-w); % amount of shadowing (percentage)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shadowtesting function
function shadowed = shadow(x,f,reflpoint,theta)
endp = length(x);
shadowed = 0;
if theta > 0 % ray from the left of mean surface normal
j = reflpoint - 1;
while shadowed == 0 && j ~= 0
dx = x(reflpoint) - x(j); df = f(j) - f(reflpoint);
if atan(df/dx) > (pi/2 - theta)
shadowed = 1; % reflection point is shadowed
end;
j = j - 1;
end;
end;
if theta < 0 % ray from the right of mean surface normal
j = reflpoint + 1;
while shadowed == 0 && j ~= endp + 1
dx = x(j) - x(reflpoint); df = f(j) - f(reflpoint);
if atan(df/dx) > (pi/2 - (- theta))
shadowed = 1; % reflection point is shadowed
end;
j = j + 1;
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculation of energy for incident ray
function ie = incidentenergy(incomingray,normals,x,slopes,reflpoint,dstrb,br)
thetalocal = acos(dot(-incomingray,[normals(reflpoint,1) normals(reflpoint,2)])); % local AOI
segmentarea = ((x(length(x))-x(1))/length(x))./cos(atan(slopes(reflpoint)));
projectedarea = cos(thetalocal).*segmentarea; % area visible to ray
if strcmp(dstrb,'gauss') == 1
ie = projectedarea.*exp(-2.*x(reflpoint).^2/(br)^2); % gaussian intensity
elseif strcmp(dstrb,'uni') == 1
if x(reflpoint) >= -br/2 && x(reflpoint) <= br/2
ie = projectedarea; % unit intensity
else ie = 0;
end;
else fprintf('%s is not a valid distribution. Type \''help goa2D\'' for further details.\n',dstrb);
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% intensity reflection coefficients according to Fresnel's equations
function [Rs,Rp,R] = intrc(n1,n2,k2,theta)
p = sqrt(1/2*(sqrt((n2.^2-k2.^2-n1^2*sin(theta).^2).^2+4*n2.^2.*k2.^2)+(n2.^2-k2.^2-n1^2*sin(theta).^2)));
q = sqrt(1/2*(sqrt((n2.^2-k2.^2-n1^2*sin(theta).^2).^2+4*n2.^2.*k2.^2)-(n2.^2-k2.^2-n1^2*sin(theta).^2)));
Rs = ((n1*cos(theta)-p).^2+q.^2)./((n1*cos(theta)+p).^2+q.^2);
Rp = ((p-n1*sin(theta).*tan(theta)).^2+q.^2)./((p+n1*sin(theta).*tan(theta)).^2+q.^2).*Rs;
R = (Rs+Rp)./2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ray & energy scattering transformation (secondary reflection points)
function [outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normals,reflpoint,n1,k1,n2,k2)
thetalocal = acos(dot(-incomingray,[normals(reflpoint,1) normals(reflpoint,2)])); % local AOI
% Fresnel reflectance approximation
[~,~,R] = intrc(n1,n2,k2,thetalocal); outgoingenergy = incomingenergy*R;
% Snell's reflection law
outgoingray = incomingray - 2*[normals(reflpoint,1) normals(reflpoint,2)]*dot(incomingray,[normals(reflpoint,1) normals(reflpoint,2)]);
outgoingnormal = sqrt(dot(outgoingray,outgoingray));
outgoingray = outgoingray/outgoingnormal; % normalized
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find new reflection point on rough surface (or mirror surface)
function [newpoint] = findnewpoint(x,f,reflpoint,outgoingray)
newpoint = 0; endp = length(x);
if outgoingray(1) < 0 % left-moving
j = reflpoint - 1;
while newpoint == 0 && j ~= 0
Dx = x(reflpoint) - x(j); Df = f(j) - f(reflpoint);
if atan(outgoingray(2)/(-outgoingray(1))) < atan(Df/Dx)
newpoint = j;
end;
j = j - 1;
end;
elseif outgoingray(1) > 0 % right-moving
j = reflpoint + 1;
while newpoint == 0 && j ~= endp + 1
Dx = x(j) - x(reflpoint); Df = f(j) - f(reflpoint);
if atan(outgoingray(2)/outgoingray(1)) < atan(Df/Dx)
newpoint = j;
end;
j = j + 1;
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find new reflection point across surface/mirror surface boundary
function [newpoint] = findmirrorpoint(x,f,xM,fM,reflpoint,outgoingray)
newpoint = 0; endpx = length(x); endpxm = length(xM);
if outgoingray(1) > 0 % to right mirrored surface
j = 1;
while newpoint == 0 && j ~= endpxm
Dx = x(endpx) - x(reflpoint) + xM(j) - xM(1); Df = fM(j) - f(reflpoint);
if atan(outgoingray(2)/outgoingray(1)) < atan(Df/Dx)
newpoint = j;
end;
j = j + 1;
end;
else % to left mirrored surface
j = endpxm;
while newpoint == 0 && j ~= 0
Dx = x(reflpoint) - x(1) + xM(endpxm) - xM(j); Df = fM(j) - f(reflpoint);
if atan(outgoingray(2)/(-outgoingray(1))) < atan(Df/Dx)
newpoint = j;
end;
j = j - 1;
end;
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculates scattering direction (-90 < theta_s < 90)
function scd = scatdir(outgoingray)
if outgoingray(1) == 0 % ray going straight up
% finally scattered ray
scd = 0;
elseif outgoingray(1) > 0 % ray going right
if atan(outgoingray(2)/outgoingray(1)) < -10^(-17) % ray scattered outside edge of surface
scd = -100;
else
% finally scattered ray
scd = acos(dot(outgoingray,[0 1]));
end;
else % ray going left
if atan(outgoingray(2)/(-outgoingray(1))) < -10^(-17) % ray scattered outside edge of surface
scd = -100;
else
% finally scattered ray
scd = -acos(dot(outgoingray,[0 1]));
end;
end;