function [H_rmscor,sdh1p,sdh1mm,D_rmscor] = frac_rmscor(p,perp)

%% 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/"

%% please cite 
% --> Kristofer Marsch (2020)
% --> DOI for this data available at https://orcid.org/0000-0002-1857-506X


%% calculating Hurst Exponent, sdh1p, sdh1mm and Fractal Dimension using the standard deviation of the correlation funtion
% based on Renard et al. (2006) and Candela et al. (2009)
% based on Stigsson & Mas Ivars (2019)

%% --parameters-- %%
% p: array of profile, first column horizontal values, second column heights
% perp: percent of profile length cutoff as integer 10%-->0.1
% profile must be de-trended
% profile must be equally spaced, e.g. constant sampling interval

%% calculate H and D
x = p(:,1)-min(p(:,1)); %move to first quadrant
y = p(:,2)-min(p(:,2));

numcut = floor(perp*length(x)); %Malinverno (1990): steps should smaller 0.2*N
steps = linspace(1,numcut,numcut);

count = 1;
for i = steps
    dy = y(1:1:end-i)-y(i+1:1:end);
    dx = mean(abs(x(1:1:end-i)-x(i+1:1:end)));
    result(count,1) = i;
    result(count,2) = std(dy,1); %standard deviation with 1/n factor
    result(count,3) = dx; %samling length in same unit as dy
    count = count+1;
end

X2 = [ones(length(log10(result(:,3))),1) log10(result(:,3))];
c = X2\log10(result(:,2));

H_rmscor = c(2);
D_rmscor = 2-H_rmscor;

%% calculate sdh1p and sdh1mm
sdh1p = 10^c(1);

sampstep = mean(diff(x));
sdh1mm = sdh1p/(sampstep/1)^H_rmscor;
