function [Mc,V,DP,GDe,CDc]=Masses(rho,coord,nm,Nsm,Ntm,Num,GyroFlag)

Mc=zeros(24);

xn=coord(:,1);
yn=coord(:,2);
zn=coord(:,3);

V=0;
mindetJ=+inf;

if GyroFlag >=1
    
GDe=zeros(24);
CDc=zeros(24);
% Calculations necessary for Coriolis matrix
% centers of mass position calculation
pos=0.5;

Bpos=[-pos,-pos,-pos;
    pos,-pos,-pos;
    pos,pos,-pos;
    -pos,pos,-pos;
    -pos,-pos,pos;
    pos,-pos,pos;
    pos,pos,pos;
    -pos,pos,pos];

XX=zeros(8);
YY=zeros(8);

for tt=1:8
    
    
    s=Bpos(tt,1);
    t=Bpos(tt,2);
    u=Bpos(tt,3);
    
    N1 = (1-s)*(1-t)*(1-u)/8;
    N2 = (1+s)*(1-t)*(1-u)/8;
    N3 = (1+s)*(1+t)*(1-u)/8;
    N4 = (1-s)*(1+t)*(1-u)/8;
    N5 = (1-s)*(1-t)*(1+u)/8;
    N6 = (1+s)*(1-t)*(1+u)/8;
    N7 = (1+s)*(1+t)*(1+u)/8;
    N8 = (1-s)*(1+t)*(1+u)/8;
    
    nt=[N1,N2,N3,N4,N5,N6,N7,N8];
    Bt=nt*coord;
    
    for rr=1:8
        
        
        
        s=Bpos(rr,1);
        t=Bpos(rr,2);
        u=Bpos(rr,3);
        
        N1 = (1-s)*(1-t)*(1-u)/8;
        N2 = (1+s)*(1-t)*(1-u)/8;
        N3 = (1+s)*(1+t)*(1-u)/8;
        N4 = (1-s)*(1+t)*(1-u)/8;
        N5 = (1-s)*(1-t)*(1+u)/8;
        N6 = (1+s)*(1-t)*(1+u)/8;
        N7 = (1+s)*(1+t)*(1+u)/8;
        N8 = (1-s)*(1+t)*(1+u)/8;
        
        nr=[N1,N2,N3,N4,N5,N6,N7,N8];
        Br=nr*coord;
        
        XX(tt,rr)=(Bt(1)+Br(1))/2;
        YY(tt,rr)=(Bt(2)+Br(2))/2;
        
        
    end
end
elseif GyroFlag < 1 % Gyroscopic matrices not requested

    GDe=0; XX=0; YY=0;
    CDc=0;
else
        XX=0; YY=0;        GDe=0;
    CDc=0;
end
%% inertial matrices computation

n=zeros(3,24);
Ns=zeros(1,8);
Nt=zeros(1,8);
Nu=zeros(1,8);

for ii=1:8
        
    % shape function matrix
    n=nm(:,:,ii);
    
    % shape functions derivatives
    Ns=Nsm(ii,:);
    Nt=Ntm(ii,:);
    Nu=Num(ii,:);
    
    xs=Ns*xn;
    
    xu=Nu*xn;
    
    xt=Nt*xn;
    
    ys=Ns*yn;
    
    yu=Nu*yn;
    
    yt=Nt*yn;
    
    zs=Ns*zn;
    
    zu=Nu*zn;
    
    zt=Nt*zn;
    
    
    J=[xs, ys, zs
        xt, yt, zt
        xu, yu, zu];
    
    % Sarrus
    detJ=xs*yt*zu+ys*zt*xu+zs*xt*yu-zs*yt*xu-ys*xt*zu-xs*zt*yu;
    
    
    % derivate funzioni di forma rispetto a x,y,z
    Nstu=[Ns;Nt;Nu];
    Nxyz=J\Nstu;
     
    % MASS MATRIX
    Mc=Mc+n'*n*rho*detJ;
    % Kiso=Kiso+B'*D*B*detJ;
    V=V+detJ;
    
    if detJ<mindetJ
        mindetJ=detJ;
    end
    
    if GyroFlag >= 1
    % ROTATIONAL EFFECTS
    %--------------------------------------------------------------------------
    %CORIOLIS MATRIX
    %reorganize shape functions matrix

    NG=[n(1,1) 0 0 n(1,2) 0 0 n(1,3) 0 0 n(1,4) 0 0 n(1,5) 0 0 n(1,6) 0 0 n(1,7) 0 0 n(1,8) 0 0;
        0 n(1,1) 0 0 n(1,2) 0 0 n(1,3) 0 0 n(1,4) 0 0 n(1,5) 0 0 n(1,6) 0 0 n(1,7) 0 0 n(1,8) 0;
        0 0 n(1,1) 0 0 n(1,2) 0 0 n(1,3) 0 0 n(1,4) 0 0 n(1,5) 0 0 n(1,6) 0 0 n(1,7) 0 0 n(1,8)];
    
    %Omega matrix around Z axis only for Corioli's components
    OMC=2*[0 1 0;
        -1 0 0;
        0 0 0];
    %Corioli's matrix integration
    CDc=CDc+rho*(NG'*OMC*NG)*detJ;  %Coriolis matrix

    %Gyroscopic matrix integration
     NZ=Nxyz(3,:);
    [Ge]=GyroscopicMatrix(n(1,1:8),NZ,XX,YY);
    GDe=GDe+rho*Ge*detJ;
    
    elseif GyroFlag < 1
    GDe=0;
    CDc=0;     
    else 
    GDe=0;
    CDc=0;

    end
end

% element distortion
DP=mindetJ*8/V;

end