function[r, nnzRMM, nnzCMM, nnzRNM, nnzCNM]=getRadius(Data, nodesM, nodesN, coordM, coordN)

% getRadius: Computes the radiuses of the radial basis functions defined on 
% the interfaces.
% INPUT:
% Data:         Structure containing all the data parameters
% nodesM:       Set of interpolation nodes
% nodesN:       Set of evaluation nodes
% coordM:       Coordinate matrix for the interpolation nodes
% coordN:       Coordinate matrix for the evaluation nodes
% OUTPUT:
% r:            Radiuses of radial basis functions
% nnzRMM:       Number of off-diagonal nonzeros in each row of PhiMM
% nnzCMM:       Number of off-diagonal nonzeros in each column of PhiMM
% nnzRNM:       Number of nonzeros in each row of PhiNM
% nnzCNM:       Number of nonzeros in each column of PhiNM

M=length(nodesM);
N=length(nodesN);
r=zeros(M,1);

% Parameters
% gamma parameter: real number between ]0,1[
% Large values of gamma result in small supports
gamma=Data.Contact.gamma;
% Gamma parameter: real number between ]0,1[ with Gamma>gamma
% Used to ensure evaluation points are well within the
% support of radial basis functions (far enough from the boundary)
Gamma=Data.Contact.Gamma;
% Tolerance for contact detection
% If the tolerance is unknown, set it to some small value.
d0=Data.Contact.d0;
% Number of closest neighbors on current interface
c=Data.Contact.c;

% Number of off-diagonal nonzeros per row of PhiMM
nnzRMM=zeros(M,1);
% Number of off-diagonal nonzeros per column of PhiMM
nnzCMM=zeros(M,1);
% Number of nonzeros per row of PhiNM
nnzRNM=zeros(N,1);
% Number of nonzeros per column of PhiNM
nnzCNM=zeros(M,1);

maxS=inf;
f=0;
niter=0;
maxiter=10;

% Version 2
while maxS>f && niter<maxiter
    % Maximum number of supports to which an evaluation point may belong.
    % Condition for the matrix PhiMM to be diagonally dominant.
    f=floor(1./((1-gamma).^4.*(1+4*gamma)));
    
    for k=1:M
        point=coordM(k,:);
        neighbors=coordM;
        neighbors(k,:)=inf;
        distMM=vecnorm(neighbors-point,2,2);
        distMN=vecnorm(coordN-point,2,2);
        
        rMM=mink(distMM,c);
        rNM=sqrt(d0^2+0.25*rMM(end)^2);
        radius=max([rMM(end) rNM]);
        
        % Make sure the closest neighbor is at a distance >= gamma*radius.
        % If the point is too close (with respect to the radius), diagonal
        % dominance will be lost very quickly.
        if radius>rMM(1)/gamma
            radius=rMM(1)/gamma;
        end
        
        s1=distMM<radius;
        s2=distMN<Gamma*radius;
        
        nnzRMM(s1)=nnzRMM(s1)+1;
        nnzRNM(s2)=nnzRNM(s2)+1;
        nnzCMM(k)=sum(s1);
        nnzCNM(k)=sum(s2);
        
        r(k)=radius;
    end
    % Compute the maximum number of supports to which a point belongs
    maxS=max(nnzRMM);
    if maxS>f
        % Increase gamma (sequence converges to 1)
        gamma=0.5*(1+gamma);
        % Reinitialize counters
        nnzRMM=zeros(M,1);
        nnzCMM=zeros(M,1);
        nnzRNM=zeros(N,1);
        nnzCNM=zeros(M,1);
        niter=niter+1;
    end
end

if niter==maxiter && maxS>f
    warning('The maximum number of rounds for radius computations was reached. Try decreasing d0')
end
end