function Q = generalized_force(Dir_inT,Rho_inT,F_final,Dir_contact,Rho_contact,F_contact,u)
%{
Function for defining the generalized forces at time t of the integration.
--------------------------------------------------------------------------
Syntax :
Q = generalized_force(Dir_inT,Rho_inT,F_final,Dir_contact,Rho_contact,F_contact,u)
--------------------------------------------------------------------------
Function Description :
This function defines the generalized forces acting on the system based on
the given values for the contact and muscle forces. If at some point hi on
the humerus the contact is lost the relevant contact force will be zero.
Yet, we have a zero force in the formula.
--------------------------------------------------------------------------
%}

F_resultant = zeros(3,1);       % initialize the resultant of forces
M_resultant = zeros(3,1);       % initialize the resultant of moments

% define the resultant force and momemnts at/around HH by the muscles
for j = 1:size(F_final,2) 
       
    F_resultant = F_final(j)*Dir_inT(:,j) + F_resultant;
    M_resultant = cross(Rho_inT(:,j),F_final(j)*Dir_inT(:,j)) + M_resultant;
    
end


% define the resultant force and momemnts at/around HH due to the contact
for i = 1:size(F_contact,1)-1
    for j = 2:size(F_contact,2)-1
        
        F_resultant = F_contact(i,j)*Dir_contact{i,j} + F_resultant;
        M_resultant = cross(Rho_contact{i,j},F_contact(i,j)*Dir_contact{i,j}) + M_resultant;
    end
end

% the contact force associated with the top point
F_resultant = F_contact(1,1)*Dir_contact{1,1} + F_resultant;
M_resultant = cross(Rho_contact{1,1},F_contact(1,1)*Dir_contact{1,1}) + M_resultant;

% the contact force associated with the bottom point
F_resultant = F_contact(size(F_contact,1),size(F_contact,2))*Dir_contact{size(F_contact,1),size(F_contact,2)} + F_resultant;
M_resultant = cross(Rho_contact{size(F_contact,1),size(F_contact,2)},F_contact(size(F_contact,1),size(F_contact,2))*Dir_contact{size(F_contact,1),size(F_contact,2)}) + M_resultant;


% calculate the generalized coordinates at time t of the integration
e = versatile_function('define partial velocities matrix',u); % define the partial velocities matrix
Q = [zeros(3); eye(3)]*F_resultant + [e; zeros(3)]*M_resultant;
return;