% velaccRRR.m
% velocity and acceleration RRR dyad 
function out=velaccRRR...
(rAi_,rAj_,rP_,vAi_,vAj_,aAi_,aAj_)
% rAi_ is the position vector of Ai
% rAj_ is the position vector of Aj
% rP_  is the position vector of P
% vAi_ is the velocity vector of Ai
% vAj_ is the velocity vector of Aj
% aAi_ is the acceleration vector of Ai
% aAj_ is the acceleration vector of Aj

% components of position, velocity, acceleration of Ai
xAi = rAi_(1);
yAi = rAi_(2);
vAix = vAi_(1);
vAiy = vAi_(2);
aAix = aAi_(1);
aAiy = aAi_(2);

% components of position, velocity, acceleration of Aj
xAj = rAj_(1);
yAj = rAj_(2);
vAjx = vAj_(1);
vAjy = vAj_(2);
aAjx = aAj_(1);
aAjy = aAj_(2);

% components of position vector of P
xP = rP_(1);
yP = rP_(2);

% velocity of P
syms vPxs vPys
eqRRR1v =...
(xAi-xP)*(vAix-vPxs)+(yAi-yP)*(vAiy-vPys);
eqRRR2v =...
(xAj-xP)*(vAjx-vPxs)+(yAj-yP)*(vAjy-vPys);

solRRRv=solve(eqRRR1v, eqRRR2v, vPxs, vPys);
vPx = eval(solRRRv.vPxs);
vPy = eval(solRRRv.vPys);

% acceleration of P
syms aPxs aPys
eqRRR1a =...
(xAi-xP)*(aAix-aPxs)+(vAix-vPx)^2+...
(yAi-yP)*(aAiy-aPys)+(vAiy-vPy)^2;
eqRRR2a =...
(xAj-xP)*(aAjx-aPxs)+(vAjx-vPx)^2+...
(yAj-yP)*(aAjy-aPys)+(vAjy-vPy)^2;

solRRRa=solve(eqRRR1a, eqRRR2a, aPxs, aPys);
aPx = eval(solRRRa.aPxs);
aPy = eval(solRRRa.aPys);

% components of velocity and acceleration of P
out = [vPx vPy aPx aPy];
end