function PLAST2D(MID,PROP,ETAN,NE,NDOF,XYZ,LE,DYNF)
%******************************************************************************
% MAIN PROGRAM COMPUTING GLOBAL STIFFNESS MATRIX RESIDUAL FORCE               %
% FOR 2D PLANE-STRAIN ELASTOPLASTIC MATERIAL                                  %
%******************************************************************************
%%                                                                            %
 global LTAN UPDATE DISPDD DISPTD FORCE GKF SIGMA XQ          % Gloal variables
 XG=[-1/sqrt(3), 1/sqrt(3)]; WGT=[1.0, 1.0];     % Integration points & weights
 INTN=0;                    % Index for history variables (each integration pt)
 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
  DSP=reshape(DISPTD(IDOF),NDOF,4);    % NDOFxNENODE element total displacement
  DSPD=reshape(DISPDD(IDOF),NDOF,4);         % element incremental displacement
  for LX=1:2, for LY=1:2                          %LOOP OVER INTEGRATION POINTS
   E1=XG(LX); E2=XG(LY);                                    % Integration point
   INTN = INTN + 1;                        % Integration point sequential index
   [~,SHPD,DET]=SHAPEL([E1 E2], ELXY);       % Det & shape function derivatives
   FAC=WGT(LX)*WGT(LY)*DET;                                % Integraton weights
   if MID>20, NALPHA=3; STRESSN=SIGMA(7:12,INTN);    % Restore converged stress
   else,      NALPHA=6; STRESSN=SIGMA(1:6,INTN); end       % from previous time
   ALPHAN=XQ(1:NALPHA,INTN);                              % Restore back stress
   EPN=XQ(NALPHA+1,INTN);                    % Restore effective plastic strain
   if MID>10, F=eye(2)+DSP*SHPD'; SHPD=inv(F)'*SHPD; end   % Updated Lagrangian
   DEPS=zeros(3,3); DEPS(1:2,1:2)=DSPD*SHPD';              % Incremental strain
   DDEPS=[DEPS(1,1) DEPS(2,2) 0 DEPS(1,2)+DEPS(2,1) 0 0]';     % Voigt notation
   if MID==1                                         % Infinitesimal plasticity
    [STRESS,ALPHA,EP,DTAN]=COMBHARD(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN,LTAN); %
   elseif MID==11                                  % Finite rotation plasticity
    FAC=FAC*det(F);                            % updated Lagrangian integration
    [STRESSN,ALPHAN]=ROTATESTRESS(DEPS,STRESSN,ALPHAN);         % Rotate stress
    [STRESS,ALPHA,EP,DTAN]=COMBHARD(PROP,ETAN,DDEPS,STRESSN,ALPHAN,EPN,LTAN); %
   elseif MID==21                                   % Multiplicative plasticity
    [STRESS,B,ALPHA,EP,DTAN]=MULPLAST(PROP,ETAN,DEPS,STRESSN,ALPHAN,EPN,LTAN);%
   end                                         % End of plasticity model branch
   if UPDATE                                         % Update plastic variables
    SIGMA(1:6,INTN)=STRESS;                                      % Store stress
    XQ(:,INTN)= [ALPHA; EP];     % Store back stress & effective plastic strain
    if MID>20, SIGMA(7:12,INTN)=B; end        % Store right Cauchy-Green tensor
    continue;                                        % Do not calculate F and K
   end                                                   % End of UPDATE branch
   BM=zeros(4,8); BG=zeros(6,8);                   % displacement-strain matrix
   for I=1:4                                           % Loop for element nodes
    COL=I*NDOF+[-1 0];                                       % DOF column index
    BM(:,COL)=[SHPD(1,I) 0; 0 SHPD(2,I); 0 0; SHPD(2,I) SHPD(1,I)]; % BN matrix
    BG(:,COL)=[SHPD(1,I) 0; SHPD(2,I) 0; 0 0; 0 SHPD(1,I); 0 SHPD(2,I); 0 0]; %
   end                                              % End of element nodes loop
   FORCE(IDOF)=FORCE(IDOF)-FAC*BM'*STRESS(1:4);      % Residual internal forces
   if LTAN                                        % Calculate tangent stiffness
    SIG=[STRESS(1) STRESS(4) STRESS(6);                         % Stress tensor
         STRESS(4) STRESS(2) STRESS(5);                                       %
         STRESS(6) STRESS(5) STRESS(3)];                                      %
    if MID==1                                        % Infinitesimal plasticity
     EKF = BM'*DTAN(1:4,1:4)*BM;                  % Calculate tangent stiffness
    elseif MID==11                                 % Finite rotation plasticity
     CTAN=[-STRESS(1) STRESS(1) STRESS(1) -STRESS(4);  % Finite rotation effect
           STRESS(2) -STRESS(2) STRESS(2) -STRESS(4);                         %
           STRESS(3) STRESS(3) -STRESS(3) 0;                                  %
           -STRESS(4) -STRESS(4) 0 -0.5*(STRESS(1)+STRESS(2))];               %
     SHEAD=kron(eye(2),SIG);                         % Initial stiffness matrix
     EKF=BM'*(DTAN(1:4,1:4)+CTAN)*BM + BG'*SHEAD*BG;  % elem. tangent stiffness
    elseif MID==21                                  % Multiplicative plasticity
     SHEAD=kron(eye(2),SIG);                         % Initial stiffness matrix
     EKF = BM'*DTAN(1:4,1:4)*BM + BG'*SHEAD*BG;   % Calculate tangent stiffness
    end                                              % End of plasticity branch
    GKF(IDOF,IDOF)=GKF(IDOF,IDOF)+DYNF*FAC*EKF;          % Assemble to global K
   end                                                     % End of LTAN branch
  end, end                                     % End of integration points loop
 end                                                      % End of element loop
end                                                                           %