% Please input the data coordinates: x and y;
ss=size(x); N=ss(1); %This is the size of the input data, number of the sampling points.
Fs=2;  % Samping frequency，no smaller than 2!  
theta1=0; dx=0;dy=0;dxn=0;dyn=0;Y1=0;P1_1=0;P1_2=0; %Initialization   

%Compute turning angle per sampling step
for i=1:N-1
   dx(i)=x(i+1)-x(i);
   dy(i)=y(i+1)-y(i);
   dxn(i)=dx(i)/sqrt(dx(i)*dx(i)+dy(i)*dy(i));
   dyn(i)=dy(i)/sqrt(dx(i)*dx(i)+dy(i)*dy(i)); 
   dl(i)=sqrt(dx(i)*dx(i)+dy(i)*dy(i));  %Sampling distance
end
for j=1:N-2
 
    sint=(dxn(j)*dyn(j+1)-dxn(j+1)*dyn(j));%*180/pi;
    cost=(dxn(j)*dxn(j+1)+dyn(j)*dyn(j+1));
    if(cost<0 && sint>0)
        theta1(j)=acos(cost);%*180/pi;
    else if(cost<0 && sint<0)
        theta1(j)=-acos(cost);
        else theta1(j)=asin(sint);
        end
    end

    
end
Curvature=theta1*Fs;  % Curvature
min_radius_curvature=1.0/max(abs(Curvature));  %minimum radius of curvature

mC=(Curvature-mean(Curvature)).^4;
Kurtosis=mean(mC)/(var(Curvature)^2)-3.0;  %Excess kurtosis

%Fourier transform
M=N-2; m=1:M;
Y1=fft(theta1);
p1_2=abs(Y1); 
p1_2=p1_2(1:M/2+1); 
p1_1=p1_2.^2; 
p1_1(2:end-1)=2.0*p1_1(2:end-1);  % Frequency spectrum, See Wang & Imran, Palaeontology, 2024
f=Fs*(0:(M/2))/M;  

psquare=p1_1/(Fs*M);  % Spectral density
figure;
plot(f,psquare,'r-.'); 
P_un_smooth=sum(psquare(ceil(M/4+0.5):M/2+1))*(Fs/M);  % Unsmooth signal power (integrating the latter half of the spectral density)
P_all=sum(psquare(1:M/2+1))*(Fs/M);  % Total signal power (integrating the whole the spectral density)
Smoothness_index=1-P_un_smooth/P_all;   %Smoothness index
