function [H_comp,D_comp,JRC_Turk,JRC_Lee] = frac_comp(p,radii)

%% 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 and Fractal Dimension using compass walking

%% --parameters-- %%
% p: array of profile, first column horizontal values, second column heights
% radii: array containing compass radii in mm, e.g. radii = [2 4 6 8 10];
% profile must be de-trended
% profile must have unique points
% profile must have sampling interval 1mm

%% calculate H and D
prof_length = max(p(:,1)); %length of profile in mm
resi = zeros(length(radii),2);
r = 0;
rchg = 0;
while r < max(radii)
    
    a = p(1,1); %center of first circle, x
    b = p(1,2); %center of first circle, y
    r = radii(rchg+1); %radius of the circle
    vs(1) = a; %initial guess for intersection
    vs(2) = b;
    circles = 0;
    start1 = 1; %line in v-vector for intersection search
    rchg = rchg+1;
    while vs(1)+r <= prof_length %many circles while
        
        circles = circles+1;
        a = vs(1);
        b = vs(2);
        %find the two points inside and outside of the circle: distance to origin
        x = zeros(1,2);
        y = zeros(1,2);
        for i = start1:length(p)
            
            if norm(p(i,:)-[a,b]) >= r
               x = [p(i-1,1),p(i,1)]';
               y = [p(i-1,2),p(i,2)]';
               start1 = i;
               break
            
            end
            
        end

        vs0 = [x,y]; %x1,y1;x2,y2: coordinates of the two points inside and outside

        %find intersection point of circle and profile
        vv = vs0(2,:)-vs0(1,:); %vector between inside and outside point
        vpos = vs0(1,:); %position vector
        phi1 = atand((y(1)-b)/(x(1)-a)); %angle from +hrztl first point
        if isnan(phi1)==1
            phi1=0;
        end
        phi2 = atand((y(2)-b)/(x(2)-a)); %angle from +hrztl second point
        phi = min(phi1,phi2):0.01:max(phi1,phi2); %testvector for angle
        t1 = (a+r*cosd(phi)-vpos(1))/vv(1); %x-coordinate
        t2 = (b+r*sind(phi)-vpos(2))/vv(2); %y-coordinate
        [~,I] = min(abs(t1-t2)); %I: Index of the min
        t = mean([t1(I),t2(I)]);
        vs = vpos+t*vv; %intersection point

        %rest f/r
        if norm(vs-[p(end,1),p(end,2)]) < r
            circles = circles+(norm(vs-[p(end,1),p(end,2)]))/r;
            break
        
        end

    end

    resi(rchg,1:2) = [r,circles]; %r and N+f/r
    vs=0;

end
    
X = [ones(length(log(resi(:,1))),1) log(resi(:,1))]; %log(r)
slope = X\log(resi(:,2)); %log(N+f/r)

% results
D_comp = -slope(2);
H_comp = 2+slope(2);

JRC_Lee = -0.87804+37.7844*((D_comp-1)/0.015)-16.9304*((D_comp-1)/0.015)^2;
JRC_Turk = -1138.6+1141.6*D_comp;

