function [Dtan,stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,J3,stress_trial] =...
    FerraraPandolfi_Incomp_response(Ferrara_Pandolfi_Macro,material_type_e,F,...
    stress_gp,stress_vmisses_gp,strain_gp,strain_vmisses_gp,Mbar,Nbar1)


KAPPA=Ferrara_Pandolfi_Macro(material_type_e(1,2),1);
MU=Ferrara_Pandolfi_Macro(material_type_e(1,2),2);
K4=Ferrara_Pandolfi_Macro(material_type_e(1,2),3);
KA=Ferrara_Pandolfi_Macro(material_type_e(1,2),4);
GAMA=Ferrara_Pandolfi_Macro(material_type_e(1,2),5);

a4=cosd(GAMA)*Mbar'+sind(GAMA)*Nbar1'; %The orientation of the first family of fibers
a6=-cosd(GAMA)*Mbar'+sind(GAMA)*Nbar1'; %The orientation of the first family of fibers

ma4=a4*a4'; %Structural tensor a4_dyad_a4
ma6=a6*a6'; %Structural tensor a6_dyad_a6

cf=F'*F; % right cauchy tensor

cc=cf^(-1); % Inverse of the right cauchy tensor
ef=(cf-eye(3))/2; %Green Lagrangian strain
dl=eye(3); %Kroneker Delta

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

%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);
EF=[ef(1,1) ef(2,2) ef(3,3) ef(1,2) ef(2,3) ef(1,3)]';
CC=[cc(1,1) cc(2,2) cc(3,3) cc(1,2) cc(2,3) cc(1,3)]';
MA4=[ma4(1,1) ma4(2,2) ma4(3,3) ma4(1,2) ma4(2,3) ma4(1,3)]';
MA6=[ma6(1,1) ma6(2,2) ma6(3,3) ma6(1,2) ma6(2,3) ma6(1,3)]';

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


W1 = I3^(-X13); W2 = X13*I1*I3^(-X43); W3 = I3^(-X23);
W4 = X23*I2*I3^(-X53); W5 = X12*I3^(-X12);
W6 = X13*I4*I3^(-X43); W7 = X13*I6*I3^(-X43);
J4= I3^(-X13)*I4;
J6=I3^(-X13)*I6;
EXP4=exp(KA*((J4-1)^2));
EXP6=exp(KA*((J6-1)^2));

%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;


I4E=2*ma4;
I6E=2*ma6;

%Deviatoric Invarients
J1E = W1*I1E - W2*I3E;
J2E = W3*I2E - W4*I3E;
J3E = W5*I3E;
J4E = W1*I4E + ((-1/3)*(I3^(-4/3))*I3E*I4);
J6E = W1*I6E + ((-1/3)*(I3^(-4/3))*I3E*I6);

%Second Piola-Kirchhoff stress
stress_vol=0.5*KAPPA*log(J3)*(J3^(-2))*I3E;
% % % % %         stress_vol=X12*KAPPA*log(J3)*(I3^(-1))*I3E;
%         stress_vol=KAPPA*log10(J3)*CC;
% % % % %         stress_iso=X12*MU*J1E;

%          stress_iso=X12*MU*((-1/3)*(I3^(-4/3))*I3E*(I1*I3^-X13)+I3^(-1/3)*J1E);
stress_iso=X12*MU*J1E;
if (J4>1 && J6>1)
    GEAR4=1;
    GEAR6=1;
elseif (J4>1 && J6<=1)
    GEAR4=1;
    GEAR6=0;
elseif (J4<=1 && J6>1)
    GEAR4=0;
    GEAR6=1;
elseif (J4<=1 && J6<=1)
    GEAR4=0;
    GEAR6=0;
end

stress_ich=GEAR4*(K4*EXP4*(J4-1)*J4E)+...
    GEAR6*(K4*EXP6*(J6-1)*J6E);



secondPiola = stress_vol+stress_iso+stress_ich;
stress_trial_fiber = stress_ich ;
cauchy=(1/J3)*F*secondPiola*F';
cauchy_fiber=(1/J3)*F*[stress_trial_fiber(1) stress_trial_fiber(4) stress_trial_fiber(6);stress_trial_fiber(4) stress_trial_fiber(2) stress_trial_fiber(5);stress_trial_fiber(6) stress_trial_fiber(5) stress_trial_fiber(3)]*F';
secondPiola_voigt=[secondPiola(1,1) secondPiola(2,2) secondPiola(3,3) secondPiola(1,2) secondPiola(2,3) secondPiola(1,3)]';



% 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
[J1EdyadJ3E_4th,~]= TensorProduct(J1E,J3E); % fourth order dyadic product
[J3EdyadJ1E_4th,~]= TensorProduct(J3E,J1E); % fourth order dyadic product
[J3EdyadJ3E_4th,~]= TensorProduct(J3E,J3E); % fourth order dyadic product
[J3EdyadJ4E_4th,~]= TensorProduct(J3E,J4E); % fourth order dyadic product
[J4EdyadJ3E_4th,~]= TensorProduct(J4E,J3E); % fourth order dyadic product
[J3EdyadJ6E_4th,~]= TensorProduct(J3E,J6E); % fourth order dyadic product
[J6EdyadJ3E_4th,~]= TensorProduct(J6E,J3E); % fourth order dyadic product
[J4EdyadJ4E_4th,~]= TensorProduct(J4E,J4E); % fourth order dyadic product
[J6EdyadJ6E_4th,~]= TensorProduct(J6E,J6E); % fourth order dyadic product

%Second order Material stiffness
I2EE = [0  4  4  0  0  0; 4  0  4  0  0  0; 4  4  0  0  0  0;
    0  0  0 -2  0  0; 0  0  0  0 -2  0; 0  0  0  0  0 -2];
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]  ;

I3EE_4th = 4*I3*(CCdyadCC-EYE); %4th order

W1 = X23*I3^(-X12);    W2 = X89*I1*I3^(-X43); W3 = X13*I1*I3^(-X43);
W4 = X43*I3^(-X12);    W5 = X89*I2*I3^(-X53); W6 = I3^(-X23);
W7 = X23*I2*I3^(-X53); W8 = I3^(-X12);        W9 = X12*I3^(-X12);
W10= X89*I4*I3^(-X43); W11= X89*I6*I3^(-X43); W12= X13*I4*I3^(-X43);
W13= X13*I6*I3^(-X43);

J1EE_4th = -W1*(J1EdyadJ3E_4th + J3EdyadJ1E_4th) + W2*(J3EdyadJ3E_4th) - W3*I3EE_4th;
J4EE_4th = -W1*(J4EdyadJ3E_4th + J3EdyadJ4E_4th) + W10*(J3EdyadJ3E_4th)- W12*I3EE_4th;
J6EE_4th = -W1*(J3EdyadJ6E_4th + J6EdyadJ3E_4th) + W11*(J3EdyadJ3E_4th)- W13*I3EE_4th;

%Second order Material stiffness
Dtan_vol = X12*KAPPA*(-(I3^(-2))*log(J3)*I3EdyadI3E_4th + X12*(I3^(-2))*I3EdyadI3E_4th + (I3^(-1))*log(J3)*I3EE_4th);
Dtan_iso = X12*MU*J1EE_4th;

Dtan_ich = GEAR4*(K4*EXP4*((J4EdyadJ4E_4th)+2*KA*((J4-1)^2)*(J4EdyadJ4E_4th)+(J4-1)*J4EE_4th))...
    + GEAR6*(K4*EXP6*((J6EdyadJ6E_4th)+2*KA*((J6-1)^2)*(J6EdyadJ6E_4th)+(J6-1)*J6EE_4th));

Dtan_S_4th = Dtan_vol + Dtan_iso + Dtan_ich;

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



