function [nx ny] = normals1D(f,x)
%
% [nx ny] = normals1D(f,x)
%
% Estimates surface normals of 1-d rough surface y=f(x) sing central
% difference of a 3-point neighbourhood with the reciprocal of the squared 
% Euclidean distance as a weighting factor.
%
% Input:    f - surface heights
%           x - surface points
%
% Output:   nx    - x-components of surface normals
%           ny    - y-components of surface normals
%
% Last updated: 2009-03-02 (David Bergström)
%

format long

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

ddx = zeros(1,nsp); % slopes in x-dir

% slopes calculation
for i=1:1:nsp;
    if i == 1
        ddx(i) = (f(i+1) - f(i)) / (x(i+1) - x(i));
    elseif i == nsp
        ddx(i) = (f(i) - f(i-1)) / (x(i) - x(i-1));
    elseif i == 2 || i == nsp - 1
        ddx(i) = (f(i+1) - f(i-1)) / (x(i+1) - x(i-1));
    elseif i == 3 || i == nsp - 2
        ddx(i) = 1/(1 + 1/4)*((f(i+1) - f(i-1))/(x(i+1) - x(i-1)) + (f(i+2) - f(i-2))/4/(x(i+2) - x(i-2)));
    else
        ddx(i) = 1/(1 + 1/4 + 1/9)*((f(i+1) - f(i-1))/(x(i+1) - x(i-1)) + (f(i+2) - f(i-2))/4/(x(i+2) - x(i-2)) + (f(i+3) - f(i-3))/9/(x(i+3) - x(i-3)));
    end;
end;

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




