function [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,stress_trial]=...
    NeoHookean_Compressible_response(NeoHookean_Compressible,material_type_e,F,...
    stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp)


LAM=NeoHookean_Compressible(material_type_e(1,2),1);
MU=NeoHookean_Compressible(material_type_e(1,2),2);

cf=F'*F; % right cauchy tensor

cc=cf^(-1); % Inverse of the right cauchy tensor
dl=eye(3); %Kroneker Delta

X12=1/2; X13=1/3; X23=2/3; X43=4/3; X53=5/3; X89=8/9;

%Vector interpretation of 2nd oreder tensors:
C1=cf(1,1); C2=cf(2,2); C3=cf(3,3); C4=cf(1,2); C5=cf(2,3); C6=cf(1,3);
CC=[cc(1,1) cc(2,2) cc(3,3) cc(1,2) cc(2,3) cc(1,3)]';

%Invarients
I1=C1+C2+C3;
I2=C1*C2+C1*C3+C2*C3-C4^2-C5^2-C6^2;
I3=det(cf);
J3=sqrt(I3);

%derivatives of invarients

I1E=2*eye(3);
I2E=2*[C2+C3, -C4, -C6;
    -C4, C3+C1, -C5;
    -C6, -C5,   C1+C2];
I3E=2*I3*cc;

W1 = I3^(-X13); W2 = X13*I1*I3^(-X43); W3 = I3^(-X23);
W4 = X23*I2*I3^(-X53); W5 = X12*I3^(-X12);

% Tensor Products
I3EdyadI3E = 4*I3^2*(CC*CC'); %2nd order
[CCdyadCC,EYE] = TensorProduct(cc,cc,cc); %[AdyadB,EYE]=TensorProduct(A,B,cc)
[I3EdyadI3E_4th,~]= TensorProduct(I3E,I3E); % fourth order dyadic product

I3EE = [ 0     4*C3  4*C2  0    -4*C5  0      ;
    4*C3  0     4*C1  0     0    -4*C6   ;
    4*C2  4*C1  0    -4*C4  0     0      ;
    0     0    -4*C4 -2*C3  2*C6  2*C5   ;
    -4*C5  0     0     2*C6 -2*C1  2*C4   ;
    0    -4*C6  0     2*C5  2*C4 -2*C2]  ; % 2nd order
I3EE_4th = 4*I3*(CCdyadCC-EYE); %4th order

ef=(cf-eye(3))/2; %Green Lagrangian strain
%Vector interpretation of the 2nd oreder tensors
EF=[ef(1,1) ef(2,2) ef(3,3) ef(1,2) ef(2,3) ef(1,3)]';
epsilon=EF;

%Second order Material stiffness
Dtan_S = -X12*I3^(-2)*(-MU+LAM*log(J3))*I3EdyadI3E+X12*I3^(-2)*(LAM*X12)*I3EdyadI3E+...
    X12*I3^(-1)*(-MU+LAM*log(J3))*I3EE;
Dtan_S_4th = -X12*I3^(-2)*(-MU+LAM*log(J3))*I3EdyadI3E_4th+X12*I3^(-2)*(LAM*X12)*I3EdyadI3E_4th+...
    X12*I3^(-1)*(-MU+LAM*log(J3))*I3EE_4th;

%Second Piola-Kirchhoff stress
secondPiola = X12*I3^(-1)*(-MU+LAM*log(J3))*I3E+X12*MU*I1E;

firstPiola=F*secondPiola;
stress_trial=[firstPiola(1,1);firstPiola(2,2);firstPiola(3,3);firstPiola(1,2);firstPiola(2,1);firstPiola(2,3);firstPiola(3,2);firstPiola(3,1);firstPiola(1,3)];

%2nd piola transformation to the first piola transformation
Dtan_P_4th = Second2First(F,Dtan_S_4th,secondPiola,dl); %4th order

Dtan=zeros(9);
idx1=[1 2 3 1 2 2 3 3 1];
idx2=[1 2 3 2 1 3 2 1 3];

for i=1:9
    for j=1:9
        Dtan(i,j)= Dtan_P_4th(idx1(i),idx2(i),idx1(j),idx2(j));
    end
end

end