% dforceRRR.m
% dynamic forces analysis 
% RRR dyad: Ai->R  P->R  Aj->R 
% link (i) with CM at Ci
% link (j) with CM at Cj
function out=dforceRRR...
(rAi_,rAj_,rP_,rCi_,rCj_,Fi_,Mi_,Fj_,Mj_)
% rAi_ is the position vector of Ai
% rAj_ is the position vector of Aj
% rP_  is the position vector of P
% rCi_ is the position vector of CM of (i)
% rCj_ is the position vector of CM of (j)
% Fi_ is the total force on (i) at Ci including
%     the inertia force at Ci and 
%     the gravitational force at Ci
% Mi_ is the total moment on (i) including
%     the moment of inertia on (i)
% Fj_ is the total force on (j) at Cj including
%     the inertia force at Cj and 
%     the gravitational force at Cj
% Mj_ is the total moment on (j) including
%     the moment of inertia on (j)

syms FAix FAiy FAjx FAjy 
% FAix is the reaction force on (i) at Ai on x-axis
% FAiy is the reaction force on (i) at Ai on y-axis
FAi_=[FAix, FAiy, 0]; % unknown joint force  at Ai
% FAjx is the reaction force on (j) at Aj on x-axis
% FAjy is the rection force on (j) at Aj on y-axis
FAj_=[FAjx, FAjy, 0]; % unknown joint force  at Aj

% sum forces of (i) and (j)
% FAi_ + Fi_ + FAj_ + Fj_ = 0_
eqF_ = FAi_+Fi_+FAj_+Fj_;
eqFx = eqF_(1); % (1)
eqFy = eqF_(2); % (2)

rPCi_ = rCi_ - rP_; % position vector PCi_ 
rPCj_ = rCj_ - rP_; % position vector PCi_ 
rPAi_ = rAi_ - rP_; % position vector PAi_ 
rPAj_ = rAj_ - rP_; % position vector PAj_ 

% sum moments of (i) about P
% PCi_ x Fi_ + PAi_ x FAi_ + Mi_ = 0_
eqMi_ = cross(rPCi_,Fi_)+cross(rPAi_,FAi_)+Mi_;
eqMiz = eqMi_(3); % (3)

% sum moments of (j) about P
% PCj_ x Fj_ + PAj_ x FAj_ + Mj_ = 0_
eqMj_ = cross(rPCj_,Fj_)+cross(rPAj_,FAj_)+Mj_;
eqMjz = eqMj_(3); % (4)

% Eqs.(1)-(4) =>
sol=solve(eqFx, eqFy, eqMiz, eqMjz);
FAixs=eval(sol.FAix);
FAiys=eval(sol.FAiy);
FAjxs=eval(sol.FAjx);
FAjys=eval(sol.FAjy);

% numerical values 
FAis_ = [FAixs, FAiys, 0];
FAjs_ = [FAjxs, FAjys, 0];

% sum forces of (i)
% FAi_ + Fi_ + Fji_ = 0_ => 
% Fji_ is the reaction force of (j) on (i) at P
Fji_ = - Fi_ - FAis_;
Fjix = Fji_(1); % x-component of Fji_
Fjiy = Fji_(2); % y-component of Fji_

out = [FAixs,FAiys,FAjxs,FAjys,Fjix,Fjiy];
end