function [ General_stress,General_strain,Delta_stress, Delta_HardeningParameter_q_SmallStrain,G_tangent]= SmallStrain(General_strain, Delta_General_strain ,G0,G_elastic,Gamma_r,General_stress, Elastic_matrix)

%   This function add small strain behavior to a constitutive model with constant shear modulus

%--------------------------------inputs------------------------------------------%
%General_stress                              :stress tensor
%General_strain                              :strain tensor 
%Delta_General_strain                        :change of strain

%Elastic_matrix                              :elastic matrix of the base model
%G_elastic                                   :constant elastic shear modulus
%G0                                          :parameter of small strain behavior
%Gamma_r                                     :parameter of small strain behavior
%--------------------------------------------------------------------------------%



%--------------------------------outputs------------------------------------------%
%General_stress                              :Updated stress tensor
%General_strain                              :Updated strain tensor 

%Delta_General_stress                        :Calculated change of stress
%Delta_HardeningParameter_q_SmallStrain      :Calculated change of hardening parameter
%G_tangent                                   :Calculated tangent modulus
%--------------------------------------------------------------------------------%






%-----------------------------------------------------------Updating deviatoric strain----------------------------------------------------------%
      Epsilon_q_SmallStrain=(2^0.5)*(((General_strain(1)-General_strain(2))^2+(General_strain(2)-General_strain(3))^2+...
          (General_strain(3)-General_strain(1))^2+3* (General_strain(4)^2+ General_strain(5)^2+ General_strain(6)^2))^0.5)/3;  %deviatoric strain
      
      Delta_Epsilon_q_SmallStrain=(2^0.5)*(((Delta_General_strain(1)-Delta_General_strain(2))^2+(Delta_General_strain(2)-...
          Delta_General_strain(3))^2+(Delta_General_strain(3)-Delta_General_strain(1))^2+3* (Delta_General_strain(4)^2+...
          Delta_General_strain(5)^2+ Delta_General_strain(6)^2))^0.5)/3;                                                        %change of deviatoric strain
      
      Epsilon_q_SmallStrain= Epsilon_q_SmallStrain+Delta_Epsilon_q_SmallStrain;                                                 %updated deviatoric strain
%-------------------------------------------------------------------------------------------------------------------------------------------------%


%----------------------------------------------alpha and derivative of yield surface based on alpha ------------------------------------------------%            
     alpha= Epsilon_q_SmallStrain*(1-G0/G_elastic/((1+0.385/Gamma_r*Epsilon_q_SmallStrain)));                        %calculation of internal variable alpha 
     d_alpha=Delta_Epsilon_q_SmallStrain*(1-G0/G_elastic/((1+0.385/Gamma_r*Epsilon_q_SmallStrain)^2));               %change of internal variable alpha
         
             varaiable_1=1+G0/G_elastic+0.385/Gamma_r*(alpha);
             varaiable_2=1-G0/G_elastic-0.385/Gamma_r*(alpha);
             varaiable_3=((varaiable_2)^(2)+4*0.385/Gamma_r*(alpha))^(0.5);
             varaiable_4= varaiable_1/varaiable_3;  
             varaiable_5_SmallStrain=1-varaiable_4;                               
             varaiable_6=(1+0.385/Gamma_r*Epsilon_q_SmallStrain)^2;
             varaiable_7=(1.5*G0/varaiable_6);          
      
     d_q0_d_alpha=varaiable_5_SmallStrain*varaiable_7;                                                           %derivative of hardening parameter based on alpha 
%----------------------------------------------------------------------------------------------------------------------------------------------------%
        
        
%-------------------------------------------- derivative of yield surface based on stress -----------------------------------------------------------%
     q_SmallStrain=(0.5*(General_stress(1)-General_stress(2))^2+0.5*(General_stress(2)-General_stress(3))^2+0.5*...
         (General_stress(3)-General_stress(1))^2+3*(General_stress(4)^2+ General_stress(5)^2+ General_stress(6)^2))^0.5;       %deviatoric stress

     if q_SmallStrain==0                                                
                 dFss_dsigma=[-1; .5; 0.5; 0; 0; 0];                  %for the first step. This is for undrained triaxial test and should be adjusted in other cases
     else
                 dFss_dsigma=-[2*General_stress(1)-General_stress(2)-General_stress(3);
                     2*General_stress(2)-General_stress(1)-General_stress(3); 
                     2*General_stress(3)-General_stress(2)-General_stress(1);
                     6*(General_stress(4)); 
                     6*(General_stress(5));
                     6*(General_stress(6));]/(2*(q_SmallStrain));
     end
%----------------------------------------------------------------------------------------------------------------------------------------------------%
        


%----------------------------------------------------------- Plastic_Multiplier ---------------------------------------------------------------------%
     Plastic_Multiplier_SmallStrainNume=(dFss_dsigma'*Elastic_matrix*Delta_General_strain+d_q0_d_alpha*d_alpha);
     Plastic_Multiplier_SmallStrainDenum=(dFss_dsigma'*Elastic_matrix*dFss_dsigma);
     Plastic_Multiplier_SmallStrain=Plastic_Multiplier_SmallStrainNume/Plastic_Multiplier_SmallStrainDenum;
%----------------------------------------------------------------------------------------------------------------------------------------------------%
        

%----------------------------------------------------------- Stress calculation----------------------------------------------------------------------%
    Delta_beta= Plastic_Multiplier_SmallStrain*dFss_dsigma;                                     %calculation of internal variable beta
    Delta_stress=Elastic_matrix*(Delta_General_strain-Delta_beta);
%----------------------------------------------------------------------------------------------------------------------------------------------------%


%----------------------------------------------------------- Tangent modulus-------------------------------------------------------------------------%
    Delta_q_SmallStrain=(0.5*(Delta_stress(1)-Delta_stress(2))^2+0.5*(Delta_stress(2)-Delta_stress(3))^2+0.5*(Delta_stress(3)-Delta_stress(1))^2+3*(Delta_stress(4)^2+ Delta_stress(5)^2+ Delta_stress(6)^2))^0.5;        
    G_tangent=Delta_q_SmallStrain/3/Delta_Epsilon_q_SmallStrain;
%----------------------------------------------------------------------------------------------------------------------------------------------------%
     
                   
 %----------------------------------------------------------------- outputs---------------------------------------------------------------------------%
    Delta_HardeningParameter_q_SmallStrain=Delta_q_SmallStrain;   
    General_stress=General_stress+ Delta_stress;    
    General_strain=General_strain+Delta_General_strain;
 %----------------------------------------------------------------------------------------------------------------------------------------------------%
   
    
     
    
    
     
   


end

