function [Rp,lambda,Sdi,SF] = stat_selec(p)
% calculation of 
%   Rp: roughness profile index 
%   lambda: ultimate slope of profile 
%   Sdi: standard deviation of the slope angles
%   SF: structure function
%
%   [Rp,lambda,Sdi,SF] = stat_selec(P)
%
%   P: profile array (first column x, second column y)
%   !!! profile must be de-trended
%   !!! profile must have constant intervals of x
%
%
% Please cite: 
% Marsch K & Fernandez-Steeger (2021) 
% https://doi.org/10.1007/s00603-020-02348-0
%
%
% License:
% "This work is licensed under the Creative Commons Attribution 4.0
% International (CC BY 4.0) License. To view a copy of this license, visit:
% http://creativecommons.org/licenses/by/4.0/"

%% general
ex = diff(p(:,1));
ey = diff(p(:,2));
L = p(end,1)-p(1,1);

%% calculate Rp
Le = sqrt(ex.^2+ey.^2);
Rp = sum(Le)/L;

%% calculate lamda
lambda = (max(p(:,2))-min(p(:,2)))/L;

%% calculate Sdi
i_e = ey./ex;
i_avg = mean(i_e);
Sdi = atand(sqrt(mean((i_e-i_avg).^2)));

%% calculate SF
samp_int = abs(p(1,1)-p(2,1));
SF = 1/L*sum(ey.^2*samp_int);
