function [nx ny nz] = normals2D(f,x,y)
%
% [nx ny nz] = normals2D(f,x,y)
%
% estimates surface normals of 2-d rough surface z=f(x,y) using central
% difference of a 3-point neighbourhood with the reciprocal of the squared 
% Euclidean distance as a weighting factor.
%
% Input:    f     - surface height
%           x     - surface point
%           y     - surface point
%
% Output:   nx    - x-components of surface normals
%           ny    - y-components of surface normals
%           nz    - z-components of surface normals
%
% Last updated: 2009-03-02 (David Bergström)
%

format long

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

ddx = zeros(nsp); % slopes in x-dir
ddy = zeros(nsp); % slopes in y-dir

% slopes calculation
for i=1:1:nsp; % i is x-component
    for j=1:1:nsp; % j is y-component
        if i == 1
            ddx(i,j) = (f(i+1,j) - f(i,j)) / (x(i+1) - x(i));
        elseif i == nsp
            ddx(i,j) = (f(i,j) - f(i-1,j)) / (x(i) - x(i-1));
        elseif i == 2 || i == nsp - 1
            ddx(i,j) = (f(i+1,j) - f(i-1,j))/(x(i+1) - x(i-1));
        elseif i == 3 || i == nsp - 2
            ddx(i,j) = 1/(1 + 1/4)*((f(i+1,j) - f(i-1,j))/(x(i+1) - x(i-1)) + (f(i+2,j) - f(i-2,j))/4/(x(i+2) - x(i-2)));
        else
            ddx(i,j) = 1/(1 + 1/4 + 1/9)*((f(i+1,j) - f(i-1,j))/(x(i+1) - x(i-1)) + (f(i+2,j) - f(i-2,j))/4/(x(i+2) - x(i-2)) + (f(i+3,j) - f(i-3,j))/9/(x(i+3) - x(i-3))); % y is columns
        end;
        
        if j == 1
            ddx(i,j) = (f(i,j+1) - f(i,j)) / (y(j+1) - y(i));
        elseif j == nsp
            ddx(i,j) = (f(i,j) - f(i,j-1)) / (y(j) - y(j-1));
        elseif j == 2 || j == nsp -1
            ddy(i,j) = (f(i,j+1) - f(i,j-1))/(y(j+1) - y(j-1));
        elseif j == 3 || j == nsp - 2
            ddy(i,j) = 1/(1 + 1/4)*((f(i,j+1) - f(i,j-1))/(y(j+1) - y(j-1)) + (f(i,j+2) - f(i,j-2))/4/(y(j+2) - y(j-2))); 
        else
            ddy(i,j) = 1/(1 + 1/4 + 1/9)*((f(i,j+1) - f(i,j-1))/(y(j+1) - y(j-1)) + (f(i,j+2) - f(i,j-2))/4/(y(j+2) - y(j-2)) + (f(i,j+3) - f(i,j-3))/9/(y(j+3) - y(j-3))); % x is rows
        end;
        
    end;
end;

nx = -ddx; ny = -ddy; nz = 1; % surface normal calculation
norm = sqrt(nx.^2+ny.^2+nz.^2); % norm of surface normal
nx = nx ./ norm; %
ny = ny ./ norm; % normalizations
nz = nz ./ norm; %
