function [STRESS, BU, ALPHA, EP, D]=MULPLAST(MP,D,L,B,ALPHA,EP,LTAN)
%******************************************************************************
% Multiplicative plasticity with linear combined hardening                    %
%******************************************************************************
%%                                                                            %
 IDEN=[1 1 1]'; TWO3=2/3; STWO3=sqrt(TWO3); EPS=1E-12*[0 1 2]';     % Constants
 MU=MP(2); BETA=MP(3); H=MP(4); Y0=MP(5);                 % Material properties
 R=inv(eye(3)-L);                                   % Inc. deformation gradient
 BM=[B(1) B(4) B(6);B(4) B(2) B(5);B(6) B(5) B(3)];  % Previous left C-G tensor
 BM=R*BM*R';                                    % Trial elastic left C-G tensor
 B=[BM(1,1) BM(2,2) BM(3,3) BM(1,2) BM(2,3) BM(1,3)]';         % Voigt notation
 [V,P]=eig(BM);                                    % Eigenvalues & eigenvectors
 % EIGEN=sort(real([P(1,1) P(2,2) P(3,3)]))';               % Principal stretches
 % EIGEN=sort(EIGEN+EPS*(rand(1,3)-0.5));        % Perturb duplicated eigenvalues
 % M=zeros(6,3);                                  % Eigenvector matrices M = N*N'
 % for K=1:3                       % Store 3 eigenvect matrices in Voigt notation
 %  KB=1+mod(K,3);
 %  KC=1+mod(KB,3);
 %  DA=1/((EIGEN(KB)-EIGEN(K))*(EIGEN(KC)-EIGEN(K)));
 %  M(1,K)=((B(1)-EIGEN(KB))*(B(1)-EIGEN(KC))+B(4)*B(4)+B(6)*B(6))*DA;
 %  M(2,K)=((B(2)-EIGEN(KB))*(B(2)-EIGEN(KC))+B(4)*B(4)+B(5)*B(5))*DA;
 %  M(3,K)=((B(3)-EIGEN(KB))*(B(3)-EIGEN(KC))+B(5)*B(5)+B(6)*B(6))*DA;
 %  M(4,K)=(B(4)*(B(1)-EIGEN(KB)+B(2)-EIGEN(KC))+B(5)*B(6))*DA;
 %  M(5,K)=(B(5)*(B(2)-EIGEN(KB)+B(3)-EIGEN(KC))+B(4)*B(6))*DA;
 %  M(6,K)=(B(6)*(B(3)-EIGEN(KB)+B(1)-EIGEN(KC))+B(4)*B(5))*DA;
 % end
 EIGEN=[P(1,1); P(2,2); P(3,3)]+EPS;            % Square of principal stretches
 M=[V(1,:).^2; V(2,:).^2; V(3,:).^2;            % Eigenvector matrices M = N*N'
	V(1,:).*V(2,:); V(2,:).*V(3,:); V(1,:).*V(3,:)];                          %
 DEPS=0.5*log(EIGEN);                           % Logarithmic principal strains
 SIGTR=D*DEPS;                                         % Trial principal stress
 ETA=SIGTR - ALPHA - sum(SIGTR)*IDEN/3;                        % Shifted stress
 ETAT=norm(ETA);                                       % Norm of shifted stress
 FYLD=ETAT - STWO3*(Y0+(1-BETA)*H*EP);                    %trial yield function
 if FYLD<Y0*1E-6                                          % Material is elastic
  SIG=SIGTR;                                         % Updated principal stress 
  STRESS=M*SIG;                                          % Updated Caucy stress
  BU=B;                                       % Updated elastic left C-G tensor
 else                                                     % Material is plastic
  GAMMA=FYLD/(2*MU + TWO3*H);                       % Plastic consistency param
  EP=EP + GAMMA*STWO3;                            % Updated eff. plastic strain
  N=ETA/ETAT;                                         % Unit vector normal to f
  DEPS=DEPS - GAMMA*N;                                 % Updated elastic strain
  SIG=SIGTR - 2*MU*GAMMA*N;                          % Updated principal stress
  ALPHA=ALPHA + TWO3*BETA*H*GAMMA*N;                      % Updated back stress
  STRESS=M*SIG;                                         % Updated Cauchy stress
  BU=M*exp(2*DEPS);                           % Updated elastic left C-G tensor
  if LTAN                          % Tangent stiffness in pricipal stress space
   V1=4*MU^2/(2*MU+TWO3*H); V2=4*MU^2*GAMMA/ETAT;                % Coefficients
   D=D-(V1-V2)*(N*N')+V2*((IDEN*IDEN')/3-eye(3,3));         % Tangent stiffness
  end                                                      % End of LTAN branch
 end                                                    % End of plastic branch
 if LTAN                                          % Calculate tangent stiffness
  J1 = sum(EIGEN); J3 = prod(EIGEN);                               % Invariants
  I2=[1 1 1 0 0 0]';                  % Rank-2 identity tensor (Voigt notation)
  I4=eye(6);I4(4,4)=.5;I4(5,5)=.5;I4(6,6)=.5;          % Rank-4 identity tensor
  IBB=[0,B(4)^2-B(1)*B(2),B(6)^2-B(1)*B(3),0,B(4)*B(6)-B(1)*B(5),0;           %
       B(4)*B(4)-B(1)*B(2),0,B(5)^2-B(2)*B(3),0,0,B(4)*B(5)-B(2)*B(6);        %
       B(6)^2-B(1)*B(3),B(5)^2-B(2)*B(3),0,B(5)*B(6)-B(3)*B(4),0,0;           %
       0,0,B(5)*B(6)-B(3)*B(4),(B(1)*B(2)-B(4)^2)/2,...                       %
	      (B(2)*B(6)-B(4)*B(5))/2,(B(1)*B(5)-B(4)*B(6))/2;                    %
       B(4)*B(6)-B(1)*B(5),0,0,(B(2)*B(6)-B(4)*B(5))/2,...                    %
	      (B(2)*B(3)-B(5)^2)/2,(B(3)*B(4)-B(5)*B(6))/2;                       %
       0,B(4)*B(5)-B(2)*B(6),0,(B(1)*B(5)-B(4)*B(6))/2,...                    %
	      (B(3)*B(4)-B(5)*B(6))/2,(B(1)*B(3)-B(6)^2)/2];                      %
  D1=1/((EIGEN(2)-EIGEN(1))*(EIGEN(3)-EIGEN(1)));                % Coefficients
  D2=1/((EIGEN(3)-EIGEN(2))*(EIGEN(1)-EIGEN(2)));                             %
  D3=1/((EIGEN(1)-EIGEN(3))*(EIGEN(2)-EIGEN(3)));                             %
  T11=-J3*D1/EIGEN(1); T12=-J3*D2/EIGEN(2); T13=-J3*D3/EIGEN(3); % Coefficients
  T21=D1*EIGEN(1);     T22=D2*EIGEN(2);     T23=D3*EIGEN(3);                  %
  T31=T21*(J1-4*EIGEN(1));T32=T22*(J1-4*EIGEN(2));T33=T23*(J1-4*EIGEN(3));    %
  C1=D1*IBB+T11*(I4-(I2-M(:,1))*(I2-M(:,1))')+T21*(B*M(:,1)'+M(:,1)*B')...    %
	+T31*M(:,1)*M(:,1)';                                                      %
  C2=D2*IBB+T12*(I4-(I2-M(:,2))*(I2-M(:,2))')+T22*(B*M(:,2)'+M(:,2)*B')...    %
	+T32*M(:,2)*M(:,2)';                                                      %
  C3=D3*IBB+T13*(I4-(I2-M(:,3))*(I2-M(:,3))')+T23*(B*M(:,3)'+M(:,3)*B')...    %
	+T33*M(:,3)*M(:,3)';                                                      %
  D = M*D*M' + 2*(C1*SIG(1)+C2*SIG(2)+C3*SIG(3));           % Tangent stiffness
 end                                                       % End of LTAN branch
end                                                                           %