function out=IndiciAssOrd(N,mododicalcolo)
% Funzione per il calcolo degli indici di associazione per variabili
% ordinali

% Se il numero di argomenti di input è minore di 2 allora si utilizza il
% modo di calcolo 1 per trovare C e D
if nargin<2
    mododicalcolo=1;
end

% Estrai il numero di righe e di colonne delle tabella di input N
[I,J] = size(N);

n = sum(sum(N)); %sample size

% Trovo le coppie discordanti e concordanti tramite una delle due formule
    C=0;
    D=0;
    for i=1:I
        for j=1:J
            if i<I && j<J
                xsel=N(i+1:I,j+1:J);
                C=C+N(i,j)*sum(xsel(:));
            end
            if j>1 && i<I
                xsel=N(i+1:I,1:j-1);
                D=D+N(i,j)*sum(xsel(:));
            end
        end
    end
    

% totpairs= total number of pairs
totpairs=0.5*n*(n-1);
% nidot = row sums
nidot=sum(N,2);
% ndotj = columns sums
ndotj=sum(N,1);

% Gamma index
gam = (C - D)/(C + D);
% Tau-a index
taua = (C - D)/totpairs;
% Taub-b index
% Tx = number of pairs tied on X
Tx=sum(nidot.*(nidot-1)/2);
% Ty = number of pairs tied on Y
Ty=sum(ndotj.*(ndotj-1)/2);
taub=(C - D)/sqrt( (totpairs-Tx)*(totpairs-Ty) );
% Tau-c index
m=min(I,J);
tauc=m*2*(C-D)/( n^2 *(m-1));
% Somers index
wr=n^2-sum(nidot.^2);
som = 2*(C-D)/wr;

out=struct;
out.gam=gam;
out.taua=taua;
out.taub=taub;
out.tauc=tauc;
out.som=som;

end