function [R,brdf,scat,shad,w] = goa3D(f,x,y,thd,phd,n1,k1,n2,k2,nfrp,dstrb,br)
%
% [R,brdf,scat,shad,w] = goa3D(f,x,y,thd,phd,n1,k1,n2,k2,nfrp,dstrb,br)
%
% calculates the directional-hemispherical reflectance for light
% incident from a medium with complex refractive index n1 + i*k1 on a 2-d 
% rough surface f(x,y) with index n2 + i*k2 with the angles of incidence 
% thd (polar) and phd (azimuth) (both in degrees) using the 
% geometric optics approximation (GOA) with nfrp^2 first reflection points. 
% Beam radius is set through br parameter while beam energy distribution is 
% controlled via parameter dstrb ('gauss' or 'uni').
%
% Input:    f     - surface height
%           x     - surface point
%           y     - 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) along square
%                   side
%           dstrb - beam energy distribution (type 'uni' for uniform and 'gauss'
%                   for gaussian)
%           br    - beam radius
%
% Output:   R    - directional-hemispherical reflectance
%           brdf - bidirectional reflection 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)
%

format long
nsp = length(x); % number of surface points along edge

% 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
end;

% initializations
theta = thd*pi/180; % Polar (zenith) AOI in radians
phi = phd*pi/180; % Azimuth AOI in radians

endpx = length(x); endpy = length(y); % endpoints

brdf = zeros(3,91,360); % multiple scattering resolved brdf (1st order)

totalenergy = 0; % total incident energy on surface
reflectedenergy = 0;

shp = zeros(nsp,nsp); % number of shadowed points
scp = zeros(nsp,nsp); % number of scattering points per incoming ray (for every 1st refl. point)
w = zeros(nsp,nsp); % number of warnings due to invalid AOI

DY = (y(end)-y(1))/length(y); % sample distance
slask = conv2(f(:,:),[1 0 -1],'same')/2/DY;
ddy = zeros(nsp);
ddy(2:end-1,2:end-1) = slask(2:end-1,2:end-1);

DX = (x(end)-x(1))/length(x); % sample distance
slask = conv2(f(:,:),[1 0 -1]','same')/2/DX;
ddx = zeros(nsp);
ddx(2:end-1,2:end-1) = slask(2:end-1,2:end-1);

normals_y = -ddy; normals_x = -ddx; normals_z = ones(nsp,nsp);
norm_of_normal = sqrt(normals_x.^2+normals_y.^2+normals_z.^2); % norm of surface normal
normals_x = normals_x ./ norm_of_normal; %
normals_y = normals_y ./ norm_of_normal; % normalizations
normals_z = normals_z ./ norm_of_normal; %

XX = x; YY = y; ZZ = f; normals_xx = normals_x; normals_yy = normals_y; normals_zz = normals_z;

%%%%% Ray tracing algorithm starts here %%%%%
for nyi = 1:nfrp; % loop over y
    ny = round(nsp/nfrp*nyi); % translate index to first reflection point
    for nxi = 1:nfrp; % loop over x
        nx = round(nsp/nfrp*nxi); % translate index to first reflection point

        % reset surface
        x = XX; y = YY; f = ZZ; normals_x = normals_xx; normals_y = normals_yy; normals_z = normals_zz;  
        
        reflection_point = [x(nx) y(ny) f(nx,ny)];

        incomingray = -[sin(theta)*cos(phi) sin(theta)*sin(phi) cos(theta)];
       
        if shadow(x,y,f,nx,ny,theta,phi) == 0 % Test for shadowing
            
            scp(nx,ny) = 1; % scattering event found
            w(nx,ny) = 0;
            
            incomingenergy = incidentenergy(incomingray,normals_x,normals_y,normals_z,x,y,ddx,ddy,nx,ny,dstrb,br); % calculation of energy of incident ray
            incomingenergy_orig = incomingenergy;
                      
            % transform incident ray and energy to scattered ray and energy
            [outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normals_x,normals_y,normals_z,nx,ny,n1,k1,n2,k2);
            
            scattered = 0;
            
            if outgoingenergy > incomingenergy
                w(nx,ny) = 1;
                outgoingenergy = incomingenergy;
                outgoingray = incomingray;
                scp(nx,ny) = 0;
                scattered = 1;
            else totalenergy = totalenergy + incomingenergy; % add to total incident energy
            end;
            
            if scattered == 0 % not a weird behaving ray

                % check if new reflection point is available
                [nextreflpoint,oobpoint] = findnewpoint(x,y,f,nx,ny,outgoingray);

                if nextreflpoint(1) == 0 && nextreflpoint(2) == 0  % no new relection point found

                    [scdir_theta, scdir_phi] = scatdir(outgoingray);

                    if scdir_theta <= pi/2 % definetely scattered
                        scattered = 1;
                    else % out-of-bounds
                    end;

                end;
            
            end;
    
            
            while scattered == 0
                
                if oobpoint(1) == 0 && oobpoint(2) == 0 % ray not out-of-bounds
                    
                    reflection_point = [x(nextreflpoint(1)) y(nextreflpoint(2)) f(nextreflpoint(1),nextreflpoint(2))];
                    
                    scp(nx,ny) = scp(nx,ny) + 1; % new scattering event found

                    incomingenergy = outgoingenergy; incomingray = outgoingray;

                    % transform incident ray and energy to scattered ray and
                    % energy
                    [outgoingray,outgoingenergy] = rayscat(incomingray,incomingenergy,normals_x,normals_y,normals_z,nextreflpoint(1),nextreflpoint(2),n1,k1,n2,k2);

                    % error for grazing incident rays (decrease error by increasing number of
                    % surface points)
                    if outgoingenergy > incomingenergy
                        w(nx,ny) = 1;
                        outgoingenergy = incomingenergy;
                        outgoingray = incomingray;
                        scp(nx,ny) = 0; 
                        totalenergy = totalenergy - incomingenergy_orig;
                        scattered = 1;
                    end;

                    % check if new reflection point is available
                    [np,oobpoint] = findnewpoint(x,y,f,nextreflpoint(1),nextreflpoint(2),outgoingray);

                else
                    
                    %bo(rp) = bo(rp) + 1;
                    
                    switch oobpoint(1)
                        case endpx
                            switch oobpoint(2)
                                case 1     % +x -y
                                    
                                    % mirrored surface for beams out-of-bounds to lower right and upper left

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,rot90(f.',2),oobpoint(1),oobpoint(2),outgoingray,8,reflection_point);
                                    
                                    % new surface
                                    f = rot90(f.',2); normals_x = -rot90(normals_x.',2); normals_y = -rot90(normals_y.',2); normals_z = rot90(normals_z.',2);

                                case endpy % +x +y
                                    
                                    % mirrored surface for beams out-of-bounds to upper right and lower left

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,f.',oobpoint(1),oobpoint(2),outgoingray,5,reflection_point);
                                    
                                    % new surface
                                    f = f.'; normals_x = -normals_x.'; normals_y = -normals_y.'; normals_z = normals_z.';

                                otherwise  % +x
                                    
                                    % mirrored surface for beams out-of-bounds in x-direction

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,flipud(f),oobpoint(1),oobpoint(2),outgoingray,1,reflection_point);

                                    % new surface
                                    f = flipud(f); normals_x = -flipud(normals_x); normals_y = flipud(normals_y); normals_z = flipud(normals_z);

                            end
                        case 1
                            switch oobpoint(2)
                                case 1     % -x -y
                                    
                                    % mirrored surface for beams out-of-bounds to upper right and lower left

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,f.',oobpoint(1),oobpoint(2),outgoingray,7,reflection_point);
                                    
                                    % new 
                                    f = f.'; normals_x = -normals_x.'; normals_y = -normals_y.'; normals_z = normals_z.';

                                case endpy % -x +y
                                    
                                    % mirrored surface for beams out-of-bounds to lower right and upper left

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,rot90(f.',2),oobpoint(1),oobpoint(2),outgoingray,6,reflection_point);
                                    
                                    % new surface
                                    f = rot90(f.',2); normals_x = -rot90(normals_x.',2); normals_y = -rot90(normals_y.',2); normals_z = rot90(normals_z.',2);

                                otherwise  % -x
                                    
                                    % mirrored surface for beams out-of-bounds in x-direction

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,flipud(f),oobpoint(1),oobpoint(2),outgoingray,3,reflection_point);

                                    % new surface
                                    f = flipud(f); normals_x = flipud(-normals_x); normals_y = flipud(normals_y); normals_z = flipud(normals_z);

                            end
                        otherwise
                            switch oobpoint(2)
                                case 1     % -y
                                    
                                    % mirrored surface for beams out-of-bounds in y-direction

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,fliplr(f),oobpoint(1),oobpoint(2),outgoingray,4,reflection_point);
                                    
                                    % new surface
                                    f = fliplr(f); normals_x = fliplr(normals_x); normals_y = -fliplr(normals_y); normals_z = fliplr(normals_z);

                                case endpy % +y
                                    
                                    % mirrored surface for beams out-of-bounds in y-direction

                                    [np,oobpoint,ro] = findmirrorpoint(x,y,fliplr(f),oobpoint(1),oobpoint(2),outgoingray,2,reflection_point);

                                    % new surface
                                    f = fliplr(f); normals_x = fliplr(normals_x); normals_y = -fliplr(normals_y); normals_z = fliplr(normals_z);

                                otherwise
                                    error('Not out-of-bound!');
                            end
                           
                    end
                    
                    if (np(1) == 0 && np(2) == 0) && oobpoint(1) ~= 0 && oobpoint(2) ~= 0
                        reflection_point = ro; % update ray origin in new coordinate system
                    end;

                
                end;
                

                if (np(1) == 0 && np(2) == 0) % no new point found and ray not out-of-bounds

                    [scdir_theta, scdir_phi] = scatdir(outgoingray);

                    if scdir_theta <= pi/2 % definetely scattered
                        scattered = 1;
                    else % out-of-bounds
                    end;

                else
                    nextreflpoint = np; % continue tracing from found reflection point
                end;

            end; % while

            if scattered == 1 && w(nx,ny) == 0 % add to differential reflection coefficient
                
                reflectedenergy = reflectedenergy + outgoingenergy;
                
                scdt = round(scdir_theta*180/pi);
                scdp = round(scdir_phi*180/pi);
                
                if scdp == 0
                    scdp = 360;
                end;
                
                % solid angle calculation
                if scdt+1 == 1
                    solidangle = pi/180*(1-cos(pi/360));
                elseif scdt+1 == 91
                    solidangle = pi/180*(cos(pi/2-pi/360));
                else solidangle = pi/180*(cos((scdt)*pi/180-pi/360) - cos((scdt)*pi/180+pi/360));
                end;
                    
                if scp(nx,ny) == 1
                    brdf(1,scdt+1,scdp) = brdf(1,scdt+1,scdp)+ outgoingenergy / solidangle / cos(scdir_theta); % first-order scattering
                else
                    if scp(nx,ny) == 2
                        brdf(2,scdt+1,scdp) = brdf(2,scdt+1,scdp) + outgoingenergy / solidangle / cos(scdir_theta); % second-order scattering
                    else brdf(3,scdt+1,scdp) = brdf(3,scdt+1,scdp) + outgoingenergy / solidangle / cos(scdir_theta); % multiple-order scattering (> 2)
                    end;
                end;

            end;

            
        else shp(nx,ny) = shp(nx,ny) + 1;
        
        end;
        
    end;

end; % no more first reflection points


%%%%%% RESULTS %%%%%%

brdf = brdf / totalenergy;

R = reflectedenergy / totalenergy; % directional-hemispherical reflectance

w = sum(nonzeros(w))/nfrp^2; % percentage of rays removed due to warnings
scat = mean(nonzeros(scp)); % mean number of scattering points per ray
shad = sum(nonzeros(shp))/(nfrp^2-w);  % amount of shadowing (percentage)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shadowtesting function

function shadowed = shadow(x,y,f,nx,ny,theta,phi)

shadowed = 0; 
dx = x(2) - x(1); dy = y(2) - y(1);
endpx = length(x); endpy = length(y); % endpoints

if (phi >= 7*pi/4) || (phi <= pi/4) % trace along positive x
    
    i = nx + 1; j = ny;
    while shadowed == 0 && i ~= (endpx + 1) && j ~= (endpy + 1) && j ~= 0
        
        delta_x = x(i) - x(nx); % diff change in x
        diff_y = delta_x*tan(phi); % diff change in y
        j = ny + round(diff_y/dy); % increment in y
        
        if j == (endpy + 1)
            return
        end;
        
        if j == 0
            return
        end;
        
        delta_y = y(j) - y(ny); % point diff change in y
        delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)
        
        if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                shadowed = 1;
        end;
        
        i = i + 1;
    end;
    
else if (phi >= 3*pi/4) && (phi <= 5*pi/4) % trace along negative x
    
        i = nx - 1; j = ny;
        while shadowed == 0 && i ~= 0 && j ~= (endpy + 1) && j ~= 0
            
            delta_x = x(i) - x(nx); % diff change in x
            diff_y = delta_x*tan(phi); % diff change in y
            j = ny + round(diff_y/dy); % increment in y
            
            if j == (endpy + 1)
                return
            end;
            
            if j == 0
                return
            end;
            
            delta_y = y(j) - y(ny); % point diff change in y
            delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)

            if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                shadowed = 1;
            end;
            
            i = i - 1;
        end;

    else if (phi > pi/4) && (phi < 3*pi/4) % trace along positive y
            
            i = nx; j = ny + 1;
            while shadowed == 0 && j ~= endpy + 1 && i ~= endpx + 1 && i ~= 0
                
                delta_y = y(j) - y(ny); % diff change in y 
                diff_x = delta_y*tan(pi/2-phi); % diff change in x
                i = nx + round(diff_x/dx); % increment in x

                if i == (endpx + 1)
                    return
                end;
                
                if i == 0
                    return
                end;
                
                delta_x = x(i) - x(nx); % diff change in x
                delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)

                if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                    shadowed = 1;
                end;

                j = j + 1;
            end;
            
        else % trace along negative y
            
            i = nx; j = ny - 1;
            while shadowed == 0 && j ~= 0 && i ~= endpx + 1 && i ~= 0

                delta_y = y(j) - y(ny); % diff change in y 
                diff_x = delta_y*(tan(3*pi/2-phi)); % diff change in x
                i = nx + round(diff_x/dx); % increment in x
                
                if i == (endpx + 1)
                    return
                end;
                
                if i == 0
                    return
                end;
                
                delta_x = x(i) - x(nx); % diff change in x
                delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)

                if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                    shadowed = 1;
                end;
                
                j = j - 1;
                
            end;

        end;
        
    end;
    
end;

    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate energy for incident ray

function ie = incidentenergy(incomingray,normals_x,normals_y,normals_z,x,y,ddx,ddy,rpx,rpy,dstrb,br)

thetalocal = acos(dot(-incomingray,[normals_x(rpx,rpy) normals_y(rpx,rpy) normals_z(rpx,rpy)])); % local AOI

segmentarea = (x(length(x))-x(1))/length(x)*(y(length(y))-y(1))/length(y)/cos(atan(ddx(rpx,rpy)))/cos(atan(ddy(rpx,rpy))); % area of tangent surface element
projectedarea = cos(thetalocal)*segmentarea; % area visible to incident ray
if strcmp(dstrb,'gauss') == 1
    ie = projectedarea.*exp(-2.*(x(rpx).^2+y(rpy).^2)/(br)^2); % gaussian intensity
elseif strcmp(dstrb,'uni') == 1
    if x(rpx) >= -br/2 && x(rpx) <= br/2 && y(rpy) >= -br/2 && y(rpy) <= br/2
        ie = projectedarea; % unit intensity
    else ie = 0;
    end;
else fprintf('%s is not a valid distribution. Type \''help goa3D\'' 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_x,normals_y,normals_z,rpx,rpy,n1,k1,n2,k2)

thetalocal = acos(dot(-incomingray,[normals_x(rpx,rpy) normals_y(rpx,rpy) normals_z(rpx,rpy)])); % local AOI

% Fresnel reflectance approximation
[~,~,R] = intrc(n1,n2,k2,thetalocal); 
outgoingenergy = incomingenergy*R;

% Snell's reflection law
outgoingray = incomingray - 2*(dot(incomingray,[normals_x(rpx,rpy) normals_y(rpx,rpy) normals_z(rpx,rpy)])*[normals_x(rpx,rpy) normals_y(rpx,rpy) normals_z(rpx,rpy)]);
outgoingnorm = sqrt(dot(outgoingray,outgoingray));
outgoingray = outgoingray / outgoingnorm; % normalized


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find new reflection point on rough surface

function [newpoint,oobpoint] = findnewpoint(x,y,f,nx,ny,outgoingray)

newpoint = [0 0]; oobpoint = [0 0];
dx = x(2) - x(1); dy = y(2) - y(1);
endpx = length(x); endpy = length(y); % endpoints

[theta,phi] = scatdir(outgoingray);

if theta ~= 0 % (& phi ~= 0) % ray pointed straight up (definetely scattered)
    
    if (phi >= 7*pi/4) || (phi <= pi/4) % trace along positive x

        i = nx + 1; j = ny;

        while (newpoint(1) == 0 && newpoint(2) == 0) && i ~= (endpx + 1) && j ~= (endpy + 1) && j ~= 0

            delta_x = x(i) - x(nx); % point diff change in x
            diff_y = delta_x*tan(phi); % ray diff change in y
            j = ny + round(diff_y/dy); % point increment in y (best appr. to ray)

            if j >= (endpy + 1)
                oobpoint = [i endpy];
                % out-of-bound in positive y
                return
            end;
            
            if j <= 0
                oobpoint = [i 1];
                % out-of-bound in negative y
                return
            end;
            
            delta_y = y(j) - y(ny); % point diff change in y
            delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)

            if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                newpoint = [i j];
            end;

            i = i + 1;
                      
        end;
        
        if i == endpx + 1
            oobpoint = [endpx j];
        end;

    elseif (phi >= 3*pi/4) && (phi <= 5*pi/4) % trace along negative x

        i = nx - 1; j = ny;
        
        while (newpoint(1) == 0 && newpoint(2) == 0) && i ~= 0 && j ~= (endpy + 1) && j ~= 0

            delta_x = x(i) - x(nx); % point diff change in x
            diff_y = delta_x*tan(phi); % ray diff change in y
            j = ny + round(diff_y/dy); % point increment in y (best appr. to ray)
            
            if j >= (endpy + 1)
                oobpoint = [i endpy];
                % out-of-bound in positive y
                return
            end;
            
            if j <= 0
                oobpoint = [i 1];
                % out-of-bound in negative y
                return
            end;
            
            delta_y = y(j) - y(ny); % point diff change in y
            delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)
            
            if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                newpoint = [i j];
            end;

            i = i - 1;
            
            if i == 0
                oobpoint = [1 j];
            end;
            
        end;

    elseif (phi > pi/4) && (phi < 3*pi/4) % trace along positive y
        
        i = nx; j = ny + 1;
        while (newpoint(1) == 0 && newpoint(2) == 0) && j ~= endpy + 1 && i ~= endpx + 1 && i ~= 0

            delta_y = y(j) - y(ny); % point diff change in y
            diff_x = delta_y*tan(pi/2-phi); % diff change in x
            i = nx + round(diff_x/dx); % increment in y
            
            if i >= (endpx + 1)
                oobpoint = [endpx j];
                % out-of-bound in positive x
                return
            end;
            
            if i <= 0
                oobpoint = [1 j];
                % out-of-bound in negative x
                return
            end;
            
            delta_x = x(i) - x(nx); % point diff change in x
            delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)

            if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                newpoint = [i j];
            end;

            j = j + 1;
            
            if j == endpy + 1
                oobpoint = [i endpy];
            end;

        end;

    else % trace along negative y
        
        i = nx; j = ny - 1;
        while (newpoint(1) == 0 && newpoint(2) == 0) && j ~= 0 && i ~= endpx + 1 && i ~= 0

            delta_y = y(j) - y(ny); % point diff change in y
            diff_x = delta_y*(tan(3*pi/2-phi)); % diff change in x
            i = nx + round(diff_x/dx); % increment in y
            
            if i >= (endpx + 1)
                oobpoint = [endpx j];
                % out-of-bound in positive x
                return
            end;
            
            if i <= 0
                oobpoint = [1 j];
                % out-of-bound in negative x
                return
            end;

            delta_x = x(i) - x(nx); % point diff change in x
            delta_z = f(i,j) - f(nx,ny); % point diff change in z (height)

            if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                newpoint = [i j];
            end;

            j = j - 1;
            
            if j == 0
                oobpoint = [i 1];
            end;
            
        end;

    end;
    
end;
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find new reflection point across surface/mirror surface boundary

function [newpoint,oobpoint,RO] = findmirrorpoint(x,y,f,oobx,ooby,outgoingray,oobdir,rayorigin)

newpoint = [0 0]; oobpoint = [0 0];
dx = x(2) - x(1); dy = y(2) - y(1);
endpx = length(x); endpy = length(y); % endpoints

[theta,phi] = scatdir(outgoingray);

if oobdir == 1

    startpoint = [1 ooby];
    RO = [-(endpx - 1)*dx 0 0] + rayorigin; % ray origin in mirror system

elseif oobdir == 2

    startpoint = [oobx 1];
    RO = [0 -(endpy - 1)*dy 0] + rayorigin; % ray origin in mirror system

elseif oobdir == 3

    startpoint = [endpx ooby];
    RO = [(endpx - 1)*dx 0 0] + rayorigin; % ray origin in mirror system

elseif oobdir == 4
    startpoint = [oobx endpy];
    RO = [0 (endpy - 1)*dy 0] + rayorigin; % ray origin in mirror system

elseif oobdir == 5

    startpoint = [1 1];
    RO = [-(endpx - 1)*dx -(endpy - 1)*dy 0] + rayorigin; % ray origin in mirror system

elseif oobdir == 6

    startpoint = [endpx 1];
    RO = [(endpx - 1)*dx -(endpy - 1)*dy 0] + rayorigin; % ray origin in mirror system

elseif oobdir == 7

    startpoint = [endpx endpy];
    RO = [(endpx - 1)*dx (endpy - 1)*dy 0] + rayorigin; % ray origin in mirror system

else

    startpoint = [1 endpy];
    RO = [-(endpx - 1)*dx (endpy - 1)*dy 0] + rayorigin; % ray origin in mirror system

end;

if (phi >= 7*pi/4) || (phi <= pi/4) % trace along positive x

    i = startpoint(1) + 1; j = startpoint(2);

    while (newpoint(1) == 0 && newpoint(2) == 0) && i ~= (endpx + 1) && j ~= (endpy + 1) && j ~= 0

        delta_x = x(i) - RO(1); % point diff change in x
        diff_y = delta_x*tan(phi); % ray diff change in y
        j = startpoint(2) + round(diff_y/dy) - round((y(ooby) - rayorigin(2))/dy); % point increment in y

        if j >= (endpy + 1)
            oobpoint = [i endpy];
            % out-of-bound in positive y
            return
        end;
        if j <= 0
            oobpoint = [i 1];
            % out-of-bound in negative y
            return
        end;

        delta_y = y(j) - RO(2); % point diff change in y
        delta_z = f(i,j) - rayorigin(3); % point diff change in z (height)
        
        if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
            newpoint = [i j];
        end;

        i = i + 1;

    end;

    if i == endpx + 1
        delta_x = (x(i-1) + dx - RO(1)); % diff change in x
        diff_y = delta_x*tan(phi); % diff change in y
        j = startpoint(2) + round(diff_y/dy) - round((y(ooby) - rayorigin(2))/dy); % increment in y
        
        if j >= (endpy + 1)
            oobpoint = [endpx endpy];
            % out-of-bound in positive x and y
            return

        elseif j <= 0
            oobpoint = [endpx 1];
            % out-of-bound in negative y
            return
        else
            oobpoint = [endpx j];
        end;
    end;


elseif (phi >= 3*pi/4) && (phi <= 5*pi/4) % trace along negative x

    i = startpoint(1) - 1; j = startpoint(2);

    while (newpoint(1) == 0 && newpoint(2) == 0) && i ~= 0 && j ~= (endpy + 1) && j ~= 0

        delta_x = x(i) - RO(1); % diff change in x
        diff_y = delta_x*tan(phi); % diff change in y
        j = startpoint(2) + round(diff_y/dy) - round((y(ooby) - rayorigin(2))/dy); % increment in y

        if j >= (endpy + 1)
            oobpoint = [i endpy];
            % out-of-bound in positive y
            return
        end;

        if j <= 0
            oobpoint = [i 1];
            % out-of-bound in negative y
            return
        end;
        
        delta_y = y(j) - RO(2);
        delta_z = f(i,j) - rayorigin(3); % point diff change in z (height)

        if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
            newpoint = [i j];
        end;

        i = i - 1;

        if i == 0
            delta_x = x(i+1) - dx - RO(1); % diff change in x
            diff_y = delta_x*tan(phi); % diff change in y
            j = startpoint(2) + round(diff_y/dy) - round((y(ooby) - rayorigin(2))/dy); % increment in y
            
            if j >= (endpy + 1)
                oobpoint = [1 endpy];
                % out-of-bound in positive x and y
                return

            elseif j <= 0
                oobpoint = [1 1];
                % out-of-bound in negative y
                return

            else
                oobpoint = [1 j];
            end;
        end;

    end;

elseif (phi > pi/4) && (phi < 3*pi/4) % trace along positive y

    i = startpoint(1); j = startpoint(2) + 1;
    
    while (newpoint(1) == 0 && newpoint(2) == 0) && j ~= endpy + 1 && j ~= endpx + 1 && j ~= 0

        delta_y = y(j) - RO(2); % point diff change in y
        diff_x = delta_y*tan(pi/2-phi); % ray diff change in x
        i = startpoint(1) + round(diff_x/dx) - round((x(oobx) - rayorigin(1))/dx); % point increment in y
        
        if i >= (endpx + 1)
            oobpoint = [endpx j];
            % out-of-bound in positive x
            return
        end;

        if i <= 0
            oobpoint = [1 j];
            % out-of-bound in negative x
            return
        end;
        
        delta_x = x(i) - RO(1); % point diff change in x
        delta_z = f(i,j) - rayorigin(3); % point diff change in z (height)

        if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
            newpoint = [i j];
        end;

        j = j + 1;

        if j == endpy + 1
            delta_y = y(j-1) + dy - RO(2); % diff change in y
            diff_x = delta_y*tan(pi/2-phi); % diff change in x
            i = startpoint(1) + round(diff_x/dx) - round((x(oobx) - rayorigin(1))/dx); % increment in y
            
            if i >= (endpx + 1)
                oobpoint = [endpx endpy];
                % out-of-bound in positive x and y
                return

            elseif i <= 0
                oobpoint = [1 endpy];
                % out-of-bound in negative y
                return

            else
                oobpoint = [i endpy];
            end;
        end;

    end;

else % trace along negative y
    
    i = startpoint(1); j = startpoint(2) - 1;
    
    while (newpoint(1) == 0 && newpoint(2) == 0) && j ~= 0 && i ~= endpx + 1 && i ~= 0
        
        delta_y = y(j) - RO(2); % diff change in y
        diff_x = delta_y*(tan(3*pi/2-phi)); % diff change in x
        i = startpoint(1) + round(diff_x/dx) - round((x(oobx) - rayorigin(1))/dx); % increment in y
        
        if i >= (endpx + 1)
            oobpoint = [endpx j];
            % out-of-bound in positive x
            return
        end;
        if i <= 0
            oobpoint = [1 j];
            % out-of-bound in negative x
            return
        end;
        
        delta_x = x(i) - RO(1); % point diff change in x
        delta_z = f(i,j) - rayorigin(3); % point diff change in z (height)

        if theta >= pi/2 - atan(delta_z/sqrt(delta_x^2 + delta_y^2))
                newpoint = [i j];
        end;

        j = j - 1;

        if j == 0
            delta_y = y(j+1) - dy - RO(2); % diff change in y
            diff_x = delta_y*(tan(3*pi/2-phi)); % diff change in x
            i = startpoint(1) + round(diff_x/dx) - round((x(oobx) - rayorigin(1))/dx); % increment in y
            if i >= (endpx + 1)
                oobpoint = [endpx 1];
                % out-of-bound in positive x and y
                return

            elseif i <= 0
                oobpoint = [1 1];
                % out-of-bound in negative y
                return

            else
                oobpoint = [i 1];
            end;
        end;

    end;

end;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculates scattering direction (0 <= sctheta <= pi/2, 0 <= scphi <= 2*pi)

function [sctheta scphi] = scatdir(outgoingray)

%sctheta = acos(outgoingray(3));
sctheta = acos(dot(outgoingray,[0 0 1]));

if outgoingray(1) < 0
    scphi = pi + atan(outgoingray(2)/outgoingray(1));
elseif outgoingray(1) > 0
    if outgoingray(2) >= 0
        scphi = atan(outgoingray(2)/outgoingray(1));
    else scphi = 2*pi + atan(outgoingray(2)/outgoingray(1));
    end;
elseif outgoingray(1) == 0
    if outgoingray(2) > 0
        scphi = pi/2;
    elseif outgoingray(2) < 0
        scphi = 3*pi/2;
    else scphi = 0;
    end;
end;
