function MASS2D(NE,NDOF,XYZ,LE,DYN)
%******************************************************************************
%  MAIN PROGRAM COMPUTING MASS MATRIX AND INERTIAL VECTOR FOR                 %
%  2D PLANE-STRAIN QUADRILATERAL ELEMENTS                                     %
%******************************************************************************
%%                                                                            %
 global ACC FORCE GKF                                        % Global variables
 XG=[-1/sqrt(3), 1/sqrt(3)]; WGT=[1.0, 1.0];     % Integration points & weights
 for IE=1:NE                             % Loop over elements. Assemble K and F
  ELXY=XYZ(LE(IE,:),:);                             % Element nodal coordinates
  IDOF=zeros(1,NDOF*size(LE,2));            % Local to global DOF mapping index
  for I=1:size(LE,2)                                   % Loop for element nodes
   IDOF(I*NDOF+[-1 0])=LE(IE,I)*NDOF+[-1 0];                % DOF mapping index
  end                                               % End of element nodes loop
  for LX=1:2, for LY=1:2                          %LOOP OVER INTEGRATION POINTS
   E1=XG(LX); E2=XG(LY);                                    % Integration point
   [SHP,~,DET]=SHAPEL([E1 E2], ELXY);                   % Det & shape functions
   FAC=DYN.RHO*WGT(LX)*WGT(LY)*DET;                        % Integraton weights
   FORCE(IDOF) = FORCE(IDOF)-FAC*SHP*(SHP'*ACC(IDOF));        % Sbtract Ma term
   GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+FAC*(SHP*SHP');                  % Mass Matrix 
  end, end                                     % End of integration points loop
 end                                                      % End of element loop
end                                                                           %