function [Mc,V,DP,GDe,CDc]=Masses20(rho,coord,n27,Ns27,Nt27,Nu27,GyroFlag)

Mc=zeros(60);

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

V=0;
mindetJ=+inf;


weights=reshape([0.171467764060357d0,0.274348422496571d0,0.171467764060357d0,
    0.274348422496571d0,0.438957475994513d0,0.274348422496571d0,
    0.171467764060357d0,0.274348422496571d0,0.171467764060357d0,
    0.274348422496571d0,0.438957475994513d0,0.274348422496571d0,
    0.438957475994513d0,0.702331961591221d0,0.438957475994513d0,
    0.274348422496571d0,0.438957475994513d0,0.274348422496571d0,
    0.171467764060357d0,0.274348422496571d0,0.171467764060357d0,
    0.274348422496571d0,0.438957475994513d0,0.274348422496571d0,
    0.171467764060357d0,0.274348422496571d0,0.171467764060357d0],27,1);

if GyroFlag >=1
    
    GDe=zeros(60);
    CDc=zeros(60);
    % 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
        0,-pos,-pos;
        pos,0,-pos;
        0,pos,-pos;
        -pos,0,-pos;
        0,-pos,pos;
        pos,0,pos;
        0,pos,pos;
        -pos,0,pos;
        -pos,-pos,0;
        pos,-pos,0;
        pos,pos,0;
        -pos,pos,0];
    
    XX=zeros(20);
    YY=zeros(20);
    
    for tt=1:20
        
        
        xi=Bpos(tt,1);
        et=Bpos(tt,2);
        ze=Bpos(tt,3);
        omg=1-xi;
        omh=1-et;
        omr=1-ze;
        opg=1+xi;
        oph=1+et;
        opr=1+ze;
        tpgphpr=opg+oph+ze;
        tmgphpr=omg+oph+ze;
        tmgmhpr=omg+omh+ze;
        tpgmhpr=opg+omh+ze;
        tpgphmr=opg+oph-ze;
        tmgphmr=omg+oph-ze;
        tmgmhmr=omg+omh-ze;
        tpgmhmr=opg+omh-ze;
        omgopg=omg*opg/4;
        omhoph=omh*oph/4;
        omropr=omr*opr/4;
        omgmopg=(omg-opg)/4;
        omhmoph=(omh-oph)/4;
        omrmopr=(omr-opr)/4;
        
        %     shape functions
        
        Nt(1,1)=-omg*omh*omr*tpgphpr/8;
        Nt(1,2)=-opg*omh*omr*tmgphpr/8;
        Nt(1,3)=-opg*oph*omr*tmgmhpr/8;
        Nt(1,4)=-omg*oph*omr*tpgmhpr/8;
        Nt(1,5)=-omg*omh*opr*tpgphmr/8;
        Nt(1,6)=-opg*omh*opr*tmgphmr/8;
        Nt(1,7)=-opg*oph*opr*tmgmhmr/8;
        Nt(1,8)=-omg*oph*opr*tpgmhmr/8;
        Nt(1,9)=omgopg*omh*omr;
        Nt(1,10)=omhoph*opg*omr;
        Nt(1,11)=omgopg*oph*omr;
        Nt(1,12)=omhoph*omg*omr;
        Nt(1,13)=omgopg*omh*opr;
        Nt(1,14)=omhoph*opg*opr;
        Nt(1,15)=omgopg*oph*opr;
        Nt(1,16)=omhoph*omg*opr;
        Nt(1,17)=omropr*omg*omh;
        Nt(1,18)=omropr*opg*omh;
        Nt(1,19)=omropr*opg*oph;
        Nt(1,20)=omropr*omg*oph;
        Bt=Nt*coord;
        
        for rr=1:20
            
            
            xi=Bpos(rr,1);
            et=Bpos(rr,2);
            ze=Bpos(rr,3);
            omg=1-xi;
            omh=1-et;
            omr=1-ze;
            opg=1+xi;
            oph=1+et;
            opr=1+ze;
            tpgphpr=opg+oph+ze;
            tmgphpr=omg+oph+ze;
            tmgmhpr=omg+omh+ze;
            tpgmhpr=opg+omh+ze;
            tpgphmr=opg+oph-ze;
            tmgphmr=omg+oph-ze;
            tmgmhmr=omg+omh-ze;
            tpgmhmr=opg+omh-ze;
            omgopg=omg*opg/4;
            omhoph=omh*oph/4;
            omropr=omr*opr/4;
            omgmopg=(omg-opg)/4;
            omhmoph=(omh-oph)/4;
            omrmopr=(omr-opr)/4;
            
            %     shape functions
            
            Nr(1,1)=-omg*omh*omr*tpgphpr/8;
            Nr(1,2)=-opg*omh*omr*tmgphpr/8;
            Nr(1,3)=-opg*oph*omr*tmgmhpr/8;
            Nr(1,4)=-omg*oph*omr*tpgmhpr/8;
            Nr(1,5)=-omg*omh*opr*tpgphmr/8;
            Nr(1,6)=-opg*omh*opr*tmgphmr/8;
            Nr(1,7)=-opg*oph*opr*tmgmhmr/8;
            Nr(1,8)=-omg*oph*opr*tpgmhmr/8;
            Nr(1,9)=omgopg*omh*omr;
            Nr(1,10)=omhoph*opg*omr;
            Nr(1,11)=omgopg*oph*omr;
            Nr(1,12)=omhoph*omg*omr;
            Nr(1,13)=omgopg*omh*opr;
            Nr(1,14)=omhoph*opg*opr;
            Nr(1,15)=omgopg*oph*opr;
            Nr(1,16)=omhoph*omg*opr;
            Nr(1,17)=omropr*omg*omh;
            Nr(1,18)=omropr*opg*omh;
            Nr(1,19)=omropr*opg*oph;
            Nr(1,20)=omropr*omg*oph;
            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


for ii=1:27
  
    n=n27(:,:,ii);
    
    Nxi=Ns27(ii,:);
    Net=Nt27(ii,:);
    Nze=Nu27(ii,:);
    
    xs=Nxi*xn;
    
    xu=Nze*xn;
    
    xt=Net*xn;
    
    ys=Nxi*yn;
    
    yu=Nze*yn;
    
    yt=Net*yn;
    
    zs=Nxi*zn;
    
    zu=Nze*zn;
    
    zt=Net*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=[Nxi;Net;Nze];
    Nxyz=J\Nstu;
    
    
    % MASS MATRIX
    Mc=Mc+n'*n*rho*detJ*weights(ii,1);
    % Kiso=Kiso+B'*D*B*detJ;
    V=V+detJ*weights(ii,1);
    
    if detJ<mindetJ
        mindetJ=detJ;
    end
    
    if GyroFlag >= 1
        % ROTATIONAL EFFECTS
        %--------------------------------------------------------------------------
        %CORIOLIS MATRIX
        %reorganize shape functions matrix
        
        ord=1:3:60;
        NG=zeros(3,60);
        NG(1,ord)=n(1,1:20);
        NG(2,ord+1)=n(1,1:20);
        NG(3,ord+2)=n(1,1:20);
        
        %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*weights(ii,1); %Coriolis matrix
        
        %Gyroscopic matrix integration
        NZ=Nxyz(3,:);
        [Ge]=GyroscopicMatrix2ord(n(1,1:20),NZ,XX,YY);
        GDe=GDe+rho*Ge*detJ*weights(ii,1); % Gyroscopic matrix
        
    elseif GyroFlag < 1
        GDe=0;
        CDc=0;
    else
        GDe=0;
        CDc=0;
        
    end
end

DP=mindetJ*8/V;

end