function [Dir_contact,Rho_contact,F_contact]= contact_model_portugease_1(r_HH,Dr_HH,CMDATA)

%{
Function for defining the contact force applying at the humeral head center.
--------------------------------------------------------------------------
Syntax :
[Dir_contact,Rho_contact,F_contact]= contact_model_portugease_1(r_HH,Dr_HH,CMDATA)
--------------------------------------------------------------------------
Function Description :
Based on the current configuration of the joint this function defines the
force applying at the humeral head center. It's like as if we conected the
humeral head center to the center of the sphere approximate the glenoid
fossa by a spring. The direction of this force is therefore in the
direction of the vector that contects the humeral head center to the
glenoid fossa center (sphere center).
The edge that this contact model has over the primer way of modeling
contact is its satiety to spaning the whole surface of the humeral head to
define the contact points and their associated forces. This can however be
the methods backslide as it does not take into account the lever arm
associated with the contact force on the rotational degrees of freedom. in
other words, this approach considers the humeral head as a point rather a
sphere. Whether or not the contact happens, is defined based on the
position of the humeral head center. Where on the humeral head the contact
is happening? It is not possible with this model to define. Maybe a
combination of this model and my post-processing routine can define this
issue. How correct the results would be in this case is arguable. on
top of this, the model fails short to define the contact area and the
associated cartilage displacement. But, it's a trade off that we make to
reimburse the computational cost. 
--------------------------------------------------------------------------
%}
sympref('HeavisideAtOrigin',0); % to assigne value (0,0) to the heaviside function
        
% the pentration of the humeral head in the fossa (joint for the center of HH)
u = norm(r_HH,2) - CMDATA.clearance;        
Du = dot(r_HH, Dr_HH)/norm(r_HH);       
% force calculation
% F_contact = CMDATA.stiffness*u^(2/3)*heaviside(u);                       % contact without damping
% F_contact = -CMDATA.s*log(1-u/CMDATA.b)*heaviside(u);
 F_contact = (CMDATA.stiffness*u^(2/3) + CMDATA.damping*Du*u^(2/3))*heaviside(u); % contact with damping
Dir_contact = -r_HH/norm(r_HH,2);
Rho_contact = zeros(3,1);


           
return