xp1:=px-cos(x);
xp2:=py-sin(x);
lp:=sqrt((xp1)^2+(xp2)^2);
ep1:=xp1/lp; ep2:=xp2/lp;                  # unit vector XP

xq1:=a-cos(x);
xq2:=-sin(x);
lq:=sqrt((xq1)^2+(xq2)^2);
eq1:=xq1/lq; eq2:=xq2/lq;                  # unit vector XQ

f:=(ep1+eq1)*sin(x)-(ep2+eq2)*cos(x);