function [M,x,y,CCregular] = RegularGraph(N,k)
% N: Number of nodes
% k: Number of edges per node in initial regular graph
% Example: RegularGraph(10,5)

if nargin == 1
    k=N-1;
end

hold off

Angle=(2*pi)/N;

for v=0:N
    y(v+1)=sin(v*Angle);
    x(v+1)=cos(v*Angle);
end

plot(x([1:N]),y([1:N]),'or')
hold on

M=zeros(N,N);

for m=1:k/2
    M=M+diag(ones(N-m,1),m)+diag(ones(m,1),(N-m));
end

if rem(k,2)~=0
    M=M+diag(ones((N/2),1),(N/2));
end

for s=1:N-1
    for t=s+1:N
        if M(s,t)==1
            plot(x([s,t]),y([s,t]),'r')
            pause(0.05)
            hold on
        end
    end
end


for n=1:N
    if sum(M(n,:))+sum(M(N-n+1,:)) ~= k
        disp('Incorrect degree present!')
    end
end


InfOrNan=0;
for s=1:N
    C=M+M';
    for t=1:N
        if C(s,t) == 0 && (s~=t)
            C(:,t)=0;
            C(t,:)=0;
        end
    end
    C(s,:)=0;
    C(:,s)=0;
    c=0;
    for u=1:N
        c=c+sum(C(u,:));
    end
    e(s)=c/2;
    cc(s)=e(s)/sum(M(s,:));
    if isinf(cc(s))==1 || isnan(cc(s))==1
        cc(s)=0;
        InfOrNan=InfOrNan+1;
        % InfOrNan discounts the NaN's and Inf's
        % when calculating the mean below.
    end
end

CCregular=sum(cc)/(N-InfOrNan);

end