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;