function [acfx,clx,acfy,cly,lags] = acf2D(f,x,y,opt)
%
% [acfx,clx,acfy,cly,lags] = acf2D(f,x,y,opt)
%
% calculates the averaged autocovariance functions and correlation lengths
% of a 2-d surface f(x,y) in x and y directions.
%
% Input:    f    - surface heights
%           x    - surface points
%           y    - surface points
%           opt - optional parameter (type 'plot' for plotting the 
%                 normalized autocovariance functions)
%
% Output:   acfx - averaged autocovariance function in x direction
%           clx  - correlation length in x direction
%           acfy - averaged autocovariance function in y direction
%           cly  - correlation length in y direction
%           lags - lag length vector (useful for plotting the acfs)
%
% Last updated: 2010-07-29 (David Bergstrom)
%

format long

N = length(x); % number of sample points
lags = linspace(0,x(N)-x(1),N); % lag lengths

% autocovariance function calculation
acfx = zeros(1,N); acfy = zeros(1,N);
for i = 1:N
    cx = xcov(f(i,:),'coeff'); % the autocovariance function in x at row i
    acfx = acfx + cx(N:2*N-1)'; % right-sided version (cumulative sum)
    cy = xcov(f(:,i),'coeff'); % the autocovariance function in y at column i
    acfy = acfy + cy(N:2*N-1); % right-sided version (cumulative sum)
end;

acfx = acfx / N; acfy = acfy / N; % averaging of acfs

% correlation lengths calculation
% in x
k = 1;
while (acfx(k) > 1/exp(1))
    k = k + 1;
end;
clx = 1/2*(x(k-1)+x(k)-2*x(1)); % the correlation length in x
% in y
k = 1;
while (acfy(k) > 1/exp(1))
    k = k + 1;
end;
cly = 1/2*(y(k-1)+y(k)-2*y(1)); % the correlation length in y

% optional plotting
if nargin<4 || isempty(opt)
    return;
end;
if nargin==4 
    if ischar(opt)
        if strcmp(opt,'plot');
            subplot(1,2,1),plot(lags,acfx)
            xlabel('lag length (x)')
            ylabel('acf (normalized)')
            title('Plot of the average normalized acf in x-dir')
            subplot(1,2,2),plot(lags,acfy)
            xlabel('lag length (y)')
            ylabel('acf (normalized)')
            title('Plot of the average normalized acf in y-dir')
        else fprintf('%s is not a valid option. Type \''help acf2D\'' for further details.\n',opt); 
        end;
    else fprintf('Option must be a string. Type \''help acf2D\'' for further details.\n');
    end;
end;