function PPDATA = postprocessing_contact_model(TimeStep,Rot_H2T,DRot_H2T,r_HH,Dr_HH,CMDATA, PPDATA)

%{
Function for defining the contacts data at each joint configuration.
--------------------------------------------------------------------------
Syntax :
PPDATA = postprocessing_contact_model(TimeStep,Rot_H2T,DRot_H2T,r_HH,Dr_HH,CMDATA, PPDATA)
--------------------------------------------------------------------------
Function Description :
This function is responsible for postprocessing the achieved results and
contructs the PPDATA structure.
The function receives the current joint configuration (rotation matrix and
the humeral head translation vector as well as their associated rate of
changes) and defines the contact force at each point on the humeral head
surface based on whether the point is in contact or not. We already defined
the points in the humerus body-fixed coordinate frame through
initilize_contact_model(.) saved in CMDATA.
This function is a replication of the "contact_model" to reproduce the
contact data based on the resulted generalized coordinates and velocities.
Here are the outputs of the function:
            1- deformation: based on the points of the humeral head.

            2- F_contact: contact forces at each points of the humeral
            head.

            3- Dir_contact: directions of the contact forces at each
            points of the humeral head.

            4- F_resultant: the resultant force of all the forces acting on
            each of the points on the humeral head. A single vector for
            each time step.

            5- F_resultant_magnitude: the magnitude of the resultant force
            at each time step. The magnitude of the force achieved form 4.

            6- h_inT: the positions of the points on the humeral head for
            which the penetration is not zero. These positions are defined
            in the thorax coordinate frame. Those points that are not in
            contact were set to zero.

            7- cont_area: the contact area of the points that are in
            contact. The other values are set to zero.

            8- max_deformation: the maximum cartilage deformation at each
            time step.
--------------------------------------------------------------------------
%}
h = CMDATA.h;           % receive h_i s in the frame {H}
a = CMDATA.a;           % receive associtaed surface area of the h_i s

for i = 1:size(h,1)-1          % associated with longtitude (theta)          
    for j = 2:size(h,2)-1      % associated with latitude (phi)   
        
        h_inT{i,j} = r_HH + Rot_H2T*h{i,j}; % hi in the {T} frame
        
        u(i,j) = norm(h_inT{i,j},2) - CMDATA.r_g; % the penetration of the point hi
                                                  % in the glenoid fossa
        
        Dh_inT{i,j} = Dr_HH + DRot_H2T*h{i,j}; % rate of change of hi 
        Du(i,j) = dot(h_inT{i,j},Dh_inT{i,j})/norm(h_inT{i,j},2); % rate of change of the ui
        
            if gt(u(i,j),0) & lt(u(i,j),0.98*CMDATA.b) % contact condition of the points
                F_contact(i,j)   = -a(i,j)*CMDATA.s*log(1-u(i,j)/CMDATA.b) + CMDATA.c*Du(i,j);
                Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
                Rho_contact{i,j} = Rot_H2T*h{i,j};
                h_inT{i,j} = h_inT{i,j};
                cont_area(i,j) = a(i,j);
            elseif ge(u(i,j),0.98*CMDATA.b)
                %u(i,j) = 0.98*b;
                F_contact(i,j)   = -a(i,j)*CMDATA.s*log(1-0.98) + CMDATA.c*Du(i,j);
                Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
                Rho_contact{i,j} = Rot_H2T*h{i,j};
                u(i,j) = 0.98*CMDATA.b;
                h_inT{i,j} = h_inT{i,j};
                cont_area(i,j) = a(i,j);
                %disp(['u' num2str(i) ',' num2str(i) '=' num2str(u(i,j)) '>=0.98b']);
            else % and if there is no contact there will be no contact force
                F_contact(i,j) = 0; 
                Dir_contact{i,j} = -h_inT{i,j}/norm(h_inT{i,j},2);
                Rho_contact{i,j} = Rot_H2T*h{i,j};
                u(i,j) = 0;
                h_inT{i,j} = 0;
                cont_area(i,j) = 0;
            end
    end
         
end

% the same calculations for the top and bottom points
% THE TOP POINTS

        h_inT{1,1} = r_HH + Rot_H2T*h{1,1}; % hi in the {T} frame
        
        u(1,1) = norm(h_inT{1,1},2) - CMDATA.r_g; % the penetration of the point hi
                                                  % in the glenoid fossa

        Dh_inT{1,1} = Dr_HH + DRot_H2T*h{1,1}; % rate of change of h1 
        Du(1,1) = dot(h_inT{1,1},Dh_inT{1,1})/norm(h_inT{1,1},2); % rate of change of the ui
                                          
                                                  
            if gt(u(1,1),0) & lt(u(1,1),0.98*CMDATA.b) % contact condition of the points
                F_contact(1,1)   = -a(1,1)*CMDATA.s*log(1-u(1,1)/CMDATA.b) + CMDATA.c*Du(1,1);
                Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
                Rho_contact{1,1} = Rot_H2T*h{1,1};
                h_inT{1,1} = h_inT{1,1};
                cont_area(1,1) = a(1,1);
            elseif ge(u(1,1),0.98*CMDATA.b)
                %u(i,j) = 0.98*b;
                F_contact(1,1)   = -a(1,1)*CMDATA.s*log(1-0.98) + CMDATA.c*Du(1,1);
                Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
                Rho_contact{1,1} = Rot_H2T*h{1,1};
                u(1,1) = 0.98*CMDATA.b;
                h_inT{1,1} = h_inT{1,1};
                cont_area(1,1) = a(1,1);
                %disp(['u' num2str(1) ',' num2str(1) '=' num2str(u(1,1)) '>=0.98b']);
            else % and if there is no contact there will be no contact force
                F_contact(1,1) = 0; 
                Dir_contact{1,1} = -h_inT{1,1}/norm(h_inT{1,1},2);
                Rho_contact{1,1} = Rot_H2T*h{1,1};
                u(1,1) = 0;
                h_inT{1,1} = 0;
                cont_area(1,1) = 0;
            end

% THE BOTTOM POINT
k = size(h,1); l = size(h,2);

        h_inT{k,l} = r_HH + Rot_H2T*h{k,l}; % hi in the {T} frame
        
        u(k,l) = norm(h_inT{k,l},2) - CMDATA.r_g; % the penetration of the point hi
                                           % in the glenoid fossa
        Dh_inT{k,l} = Dr_HH + DRot_H2T*h{k,l}; % rate of change of h bottom 
        Du(k,l) = dot(h_inT{k,l},Dh_inT{k,l})/norm(h_inT{k,l},2); % rate of change of the ui
                                          
            if gt(u(k,l),0) & lt(u(k,l),0.98*CMDATA.b) % contact condition of the points
                F_contact(k,l)   = -a(k,l)*CMDATA.s*log(1-u(k,l)/CMDATA.b) + CMDATA.c*Du(k,l);
                Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
                Rho_contact{k,l} = Rot_H2T*h{k,l};
                h_inT{k,l} = h_inT{k,l};
                cont_area(k,l) = a(k,l);
            elseif ge(u(k,l),0.98*CMDATA.b)
                %u(k,l) = = 0.98*b;
                F_contact(k,l)   = -a(k,l)*CMDATA.s*log(1-0.98) + CMDATA.c*Du(k,l);
                Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
                Rho_contact{k,l} = Rot_H2T*h{k,l};
                u(k,l) = 0.98*CMDATA.b;
                h_inT{k,l} = h_inT{k,l};
                cont_area(k,l) = a(k,l);
                %disp(['u' num2str(k) ',' num2str(l) '=' num2str(u(k,l)) '>=0.98b']);
            else % and if there is no contact there will be no contact force
                F_contact(k,l) = 0; 
                Dir_contact{k,l} = -h_inT{k,l}/norm(h_inT{k,l},2);
                Rho_contact{k,l} = Rot_H2T*h{k,l};
                u(k,l) = 0;
                h_inT{k,l} = 0;
                cont_area(k,l) = 0;
            end
            
% NOTE: these values are defined for the inputted time step. We make the structure
% containing the information for the whole time horizon in the for loop in
% the main function of the postprocessing.
 
 PPDATA.deformation{TimeStep} = u; % deformations of all the points on the humeral head
 PPDATA.max_deformation(TimeStep) = max(u(:)); % maximum deformations at each time step
 PPDATA.F_contact{TimeStep} = F_contact;
 PPDATA.Dir_contact{TimeStep} = Dir_contact;
 PPDATA.h_inT{TimeStep} = h_inT;
 PPDATA.cont_area(TimeStep) = 1e6*sum(cont_area(:)); % contact area [mm2]
 % % define the resultant contact force
 F_resultant = postprocessing_contact_force(Dir_contact,F_contact); % function to deine the resultant contact force
 PPDATA.F_resultant{TimeStep} = F_resultant;
 PPDATA.F_resultant_magnitude(TimeStep) = norm(F_resultant,2);% magnitude of the resultant muscle force
 
 
% deformation = u;            % deformations of all the points on the humeral head            
% F_resultant = postprocessing_contact_force(Dir_contact,F_contact); % function to deine the resultant contact force
% F_resultant_magnitude = norm(F_resultant,2); % magnitude of the resultant muscle force
 

return;