syms x y x3 y3 q L1 r1 r2 Dy Dx x21 y21
x3 = r2.*cos(q);
y3 = r2.*sin(q);

F1 = (x-x3).^2+(y-y3).^2-L1.^2;
F2 = (x-Dx).^2+(y-Dy).^2-r1.^2;

[x2 y2] = solve(F1, F2,'x','y');

q = 0; %:pi./20:2.*pi;
r2 = 5.0;
r1 = 10.0;
L1 = 24.0;
Dx = 4.0;
Dy = 26.0;


q=pi/2:pi./10:1.8*pi;


%x21=1./2.*(-1./(4.*r2.^2-8.*Dy.*r2.*sin(q)-8.*r2.*cos(q).*Dx+4.*Dy.^2+4.*Dx.^2).*(4.*Dy.^3-4.*r2.*sin(q).*L1.^2-8.*Dy.*r2.^2.*sin(q).^2-8.*r2.^2.*cos(q).*Dx.*sin(q)-4.*Dy.*r1.^2+4.*Dy.*L1.^2-8.*r2.*cos(q).*Dx.*Dy+4.*r2.^3.*sin(q)+4.*Dx.^2.*Dy+4.*Dx.^2.*r2.*sin(q)+4.*Dy.*r2.^2-4.*Dy.^2.*r2.*sin(q)+4.*r1.^2.*r2.*sin(q)+4.*(8.*r2.^2.*cos(q).*Dx.*sin(q).*Dy.*r1.^2-16.*r2.^2.*cos(q).*Dx.^3.*sin(q).*Dy-16.*r2.^4.*cos(q).*Dx.*sin(q).*Dy+4.*r2.^3.*cos(q).*Dx.*sin(q).^2.*r1.^2-4.*Dy.^2.*r1.^2.*r2.*cos(q).*Dx-4.*Dy.*r1.^2.*Dx.^2.*r2.*sin(q)-4.*r2.*cos(q).*Dx.*r1.^2.*L1.^2-8.*Dy.^3.*r2.^2.*cos(q).*Dx.*sin(q)+4.*r2.^3.*sin(q).^2.*L1.^2.*cos(q).*Dx+8.*r2.^2.*sin(q).*L1.^2.*cos(q).*Dx.*Dy-4.*r2.*sin(q).*L1.^2.*Dx.^2.*Dy+8.*Dy.*r2.^4.*sin(q).^3.*cos(q).*Dx+4.*Dy.^2.*r2.^3.*sin(q).^2.*cos(q).*Dx+2.*r2.*cos(q).*Dx.*L1.^4-8.*r2.*cos(q).*Dx.^3.*r1.^2-8.*r2.*cos(q).*Dx.^3.*L1.^2-8.*r2.^3.*cos(q).*Dx.*r1.^2-8.*r2.^3.*cos(q).*Dx.*L1.^2-4.*r2.^2.*cos(q).^2.*Dx.^2.*Dy.^2+8.*r2.*cos(q).*Dx.^3.*Dy.^2+8.*r2.^3.*cos(q).*Dx.*Dy.^2+16.*r2.^3.*sin(q).*Dx.^2.*Dy+4.*Dx.^4.*Dy.*r2.*sin(q)-2.*Dx.^2.*r2.^2.*sin(q).^2.*r1.^2-4.*Dy.^2.*L1.^2.*r2.*cos(q).*Dx+8.*r2.^3.*cos(q).^2.*Dx.^2.*sin(q).*Dy+4.*Dy.^3.*Dx.^2.*r2.*sin(q)- ...
%   Dx.^2.*r1.^4-Dx.^2.*L1.^4+2.*Dx.^4.*r1.^2+2.*Dx.^4.*L1.^2+r2.^6.*sin(q).^2-Dx.^6+2.*Dx.^2.*r1.^2.*L1.^2-8.*r2.^4.*cos(q).^2.*Dx.^2+r1.^4.*r2.^2.*sin(q).^2+6.*r2.^5.*cos(q).*Dx+20.*r2.^3.*cos(q).*Dx.^3-8.*r2.^2.*cos(q).^2.*Dx.^4+6.*r2.*cos(q).*Dx.^5+6.*r2.^4.*sin(q).^2.*Dx.^2+4.*r2.^5.*sin(q).*Dy-2.*r2.^4.*sin(q).^2.*r1.^2+5.*Dx.^4.*r2.^2.*sin(q).^2+r2.^2.*sin(q).^2.*L1.^4-2.*r2.^4.*sin(q).^2.*L1.^2+4.*Dy.^2.*r2.^4.*sin(q).^4-4.*Dy.*r2.^5.*sin(q).^3+Dy.^4.*r2.^2.*sin(q).^2+4.*Dy.^3.*r2.^3.*sin(q)-2.*Dy.^2.*r2.^4.*sin(q).^2-4.*Dy.^3.*r2.^3.*sin(q).^3+4.*r2.^2.*Dx.^2.*r1.^2+4.*r2.^2.*Dx.^2.*L1.^2+2.*r2.^2.*r1.^2.*L1.^2-8.*Dx.^2.*Dy.^2.*r2.^2+2.*Dy.^4.*r2.*cos(q).*Dx+4.*r2.^3.*sin(q).^3.*L1.^2.*Dy-2.*r2.^2.*sin(q).^2.*L1.^2.*Dx.^2-4.*r2.^3.*sin(q).*L1.^2.*Dy-2.*r2.^2.*sin(q).^2.*L1.^2.*Dy.^2-2.*r2.^2.*sin(q).^2.*L1.^2.*r1.^2-2.*Dy.^2.*r2.^2.*sin(q).^2.*r1.^2+2.*Dy.^2.*r2.^2.*sin(q).^2.*Dx.^2-12.*Dy.*r2.^3.*sin(q).^3.*Dx.^2+4.*Dy.*r2.^3.*sin(q).^3.*r1.^2+4.*r2.^4.*cos(q).^2.*Dx.^2.*sin(q).^2-4.*r2.^5.*cos(q).*Dx.*sin(q).^2-12.*r2.^3.*cos(q).*Dx.^3.*sin(q).^2-4.*Dy.*r1.^2.*r2.^3.*sin(q)+8.*r2.^2.*cos(q).^2.*Dx.^2.*r1.^2+8.*r2.^2.*cos(q).^2.*Dx.^2.*L1.^2+2.*r2.*cos(q).*Dx.*r1.^4-7.*r2.^2.*Dx.^4+2.*r2.^4.*r1.^2-Dy.^4.*Dx.^2-Dy.^4.*r2.^2-2.*Dx.^4.*Dy.^2-2.*Dy.^2.*r2.^4-r2.^2.*r1.^4-r2.^2.*L1.^4-7.*r2.^4.*Dx.^2+2.*r2.^4.*L1.^2-r2.^6+2.*Dy.^2.*r1.^2.*Dx.^2+2.*Dy.^2.*r1.^2.*r2.^2+2.*Dy.^2.*L1.^2.*Dx.^2+2.*Dy.^2.*L1.^2.*r2.^2).^(1./2)).*r2.*sin(q)...
%   -Dy.^2-Dx.^2-L1.^2+r1.^2+r2.^2+1./(4.*r2.^2-8.*Dy.*r2.*sin(q)-8.*r2.*cos(q).*Dx+4.*Dy.^2+4.*Dx.^2).*(4.*Dy.^3-4.*r2.*sin(q).*L1.^2-8.*Dy.*r2.^2.*sin(q).^2-8.*r2.^2.*cos(q).*Dx.*sin(q)-4.*Dy.*r1.^2+4.*Dy.*L1.^2-8.*r2.*cos(q).*Dx.*Dy+4.*r2.^3.*sin(q)+4.*Dx.^2.*Dy+4.*Dx.^2.*r2.*sin(q)+4.*Dy.*r2.^2-4.*Dy.^2.*r2.*sin(q)+4.*r1.^2.*r2.*sin(q)+4.*(8.*r2.^2.*cos(q).*Dx.*sin(q).*Dy.*r1.^2-16.*r2.^2.*cos(q).*Dx.^3.*sin(q).*Dy-16.*r2.^4.*cos(q).*Dx.*sin(q).*Dy+4.*r2.^3.*cos(q).*Dx.*sin(q).^2.*r1.^2-4.*Dy.^2.*r1.^2.*r2.*cos(q).*Dx-4.*Dy.*r1.^2.*Dx.^2.*r2.*sin(q)-4.*r2.*cos(q).*Dx.*r1.^2.*L1.^2-8.*Dy.^3.*r2.^2.*cos(q).*Dx.*sin(q)+4.*r2.^3.*sin(q).^2.*L1.^2.*cos(q).*Dx+8.*r2.^2.*sin(q).*L1.^2.*cos(q).*Dx.*Dy-4.*r2.*sin(q).*L1.^2.*Dx.^2.*Dy+8.*Dy.*r2.^4.*sin(q).^3.*cos(q).*Dx+4.*Dy.^2.*r2.^3.*sin(q).^2.*cos(q).*Dx+2.*r2.*cos(q).*Dx.*L1.^4-8.*r2.*cos(q).*Dx.^3.*r1.^2-8.*r2.*cos(q).*Dx.^3.*L1.^2-8.*r2.^3.*cos(q).*Dx.*r1.^2-8.*r2.^3.*cos(q).*Dx.*L1.^2-4.*r2.^2.*cos(q).^2.*Dx.^2.*Dy.^2+8.*r2.*cos(q).*Dx.^3.*Dy.^2+8.*r2.^3.*cos(q).*Dx.*Dy.^2+16.*r2.^3.*sin(q).*Dx.^2.*Dy+4.*Dx.^4.*Dy.*r2.*sin(q)-2.*Dx.^2.*r2.^2.*sin(q).^2.*r1.^2-4.*Dy.^2.*L1.^2.*r2.*cos(q).*Dx+8.*r2.^3.*cos(q).^2.*Dx.^2.*sin(q).*Dy+4.*Dy.^3 ...
%   .*Dx.^2.*r2.*sin(q)-Dx.^2.*r1.^4-Dx.^2.*L1.^4+2.*Dx.^4.*r1.^2+2.*Dx.^4.*L1.^2+r2.^6.*sin(q).^2-Dx.^6+2.*Dx.^2.*r1.^2.*L1.^2-8.*r2.^4.*cos(q).^2.*Dx.^2+r1.^4.*r2.^2.*sin(q).^2+6.*r2.^5.*cos(q).*Dx+20.*r2.^3.*cos(q).*Dx.^3-8.*r2.^2.*cos(q).^2.*Dx.^4+6.*r2.*cos(q).*Dx.^5+6.*r2.^4.*sin(q).^2.*Dx.^2+4.*r2.^5.*sin(q).*Dy-2.*r2.^4.*sin(q).^2.*r1.^2+5.*Dx.^4.*r2.^2.*sin(q).^2+r2.^2.*sin(q).^2.*L1.^4-2.*r2.^4.*sin(q).^2.*L1.^2+4.*Dy.^2.*r2.^4.*sin(q).^4-4.*Dy.*r2.^5.*sin(q).^3+Dy.^4.*r2.^2.*sin(q).^2+4.*Dy.^3.*r2.^3.*sin(q)-2.*Dy.^2.*r2.^4.*sin(q).^2-4.*Dy.^3.*r2.^3.*sin(q).^3+4.*r2.^2.*Dx.^2.*r1.^2+4.*r2.^2.*Dx.^2.*L1.^2+2.*r2.^2.*r1.^2.*L1.^2-8.*Dx.^2.*Dy.^2.*r2.^2+2.*Dy.^4.*r2.*cos(q).*Dx+4.*r2.^3.*sin(q).^3.*L1.^2.*Dy-2.*r2.^2.*sin(q).^2.*L1.^2.*Dx.^2-4.*r2.^3.*sin(q).*L1.^2.*Dy-2.*r2.^2.*sin(q).^2.*L1.^2.*Dy.^2-2.*r2.^2.*sin(q).^2.*L1.^2.*r1.^2-2.*Dy.^2.*r2.^2.*sin(q).^2.*r1.^2+2.*Dy.^2.*r2.^2.*sin(q).^2.*Dx.^2-12.*Dy.*r2.^3.*sin(q).^3.*Dx.^2+4.*Dy.*r2.^3.*sin(q).^3.*r1.^2+4.*r2.^4.*cos(q).^2.*Dx.^2.*sin(q).^2-4.*r2.^5.*cos(q).*Dx.*sin(q).^2-12.*r2.^3.*cos(q).*Dx.^3.*sin(q).^2-4.*Dy.*r1.^2.*r2.^3.*sin(q)+8.*r2.^2.*cos(q).^2.*Dx.^2.*r1.^2+8.*r2.^2.*cos(q).^2.*Dx.^2.*L1.^2+2.*r2.*cos(q).*Dx.*r1.^4-7.*r2.^2.*Dx.^4+2.*r2.^4.*r1.^2 ...
%   -Dy.^4.*Dx.^2-Dy.^4.*r2.^2-2.*Dx.^4.*Dy.^2-2.*Dy.^2.*r2.^4-r2.^2.*r1.^4-r2.^2.*L1.^4-7.*r2.^4.*Dx.^2+2.*r2.^4.*L1.^2-r2.^6+2.*Dy.^2.*r1.^2.*Dx.^2+2.*Dy.^2.*r1.^2.*r2.^2+2.*Dy.^2.*L1.^2.*Dx.^2+2.*Dy.^2.*L1.^2.*r2.^2).^(1./2)).*Dy)./(-Dx+r2.*cos(q));
   
%y21=1./2./(4.*r2.^2-8.*Dy.*r2.*sin(q)-8.*r2.*cos(q).*Dx+4.*Dy.^2+4.*Dx.^2).*(4.*Dy.^3-4.*r2.*sin(q).*L1.^2-8.*Dy.*r2.^2.*sin(q).^2-8.*r2.^2.*cos(q).*Dx.*sin(q)-4.*Dy.*r1.^2+4.*Dy.*L1.^2-8.*r2.*cos(q).*Dx.*Dy+4.*r2.^3.*sin(q)+4.*Dx.^2.*Dy+4.*Dx.^2.*r2.*sin(q)+4.*Dy.*r2.^2-4.*Dy.^2.*r2.*sin(q)+4.*r1.^2.*r2.*sin(q)+4.*(8.*r2.^2.*cos(q).*Dx.*sin(q).*Dy.*r1.^2-16.*r2.^2.*cos(q).*Dx.^3.*sin(q).*Dy-16.*r2.^4.*cos(q).*Dx.*sin(q).*Dy+4.*r2.^3.*cos(q).*Dx.*sin(q).^2.*r1.^2-4.*Dy.^2.*r1.^2.*r2.*cos(q).*Dx-4.*Dy.*r1.^2.*Dx.^2.*r2.*sin(q)-4.*r2.*cos(q).*Dx.*r1.^2.*L1.^2-8.*Dy.^3.*r2.^2.*cos(q).*Dx.*sin(q)+4.*r2.^3.*sin(q).^2.*L1.^2.*cos(q).*Dx+8.*r2.^2.*sin(q).*L1.^2.*cos(q).*Dx.*Dy-4.*r2.*sin(q).*L1.^2.*Dx.^2.*Dy+8.*Dy.*r2.^4.*sin(q).^3.*cos(q).*Dx+4.*Dy.^2.*r2.^3.*sin(q).^2.*cos(q).*Dx+2.*r2.*cos(q).*Dx.*L1.^4-8.*r2.*cos(q).*Dx.^3.*r1.^2-8.*r2.*cos(q).*Dx.^3.*L1.^2-8.*r2.^3.*cos(q).*Dx.*r1.^2-8.*r2.^3.*cos(q).*Dx.*L1.^2-4.*r2.^2.*cos(q).^2.*Dx.^2.*Dy.^2+8.*r2.*cos(q).*Dx.^3.*Dy.^2+8.*r2.^3.*cos(q).*Dx.*Dy.^2+16.*r2.^3.*sin(q).*Dx.^2.*Dy+4.*Dx.^4.*Dy.*r2.*sin(q)-2.*Dx.^2.*r2.^2.*sin(q).^2.*r1.^2-4.*Dy.^2.*L1.^2.*r2.*cos(q).*Dx+8.*r2.^3.*cos(q).^2.*Dx.^2.*sin(q).*Dy+4.*Dy.^3.*Dx.^2.*r2.*sin(q)-Dx.^2.*r1.^4-Dx.^2.*L1.^4+2.*Dx.^4.*r1.^2+2.*Dx.^4.*L1.^2+r2.^6.*sin(q).^2-Dx.^6+2.*Dx.^2.*r1.^2.*L1.^2-8.*r2.^4.*cos(q).^2.*Dx.^2+r1.^4.*r2.^2.*sin(q).^2+6.*r2.^5.*cos(q).*Dx+20.*r2.^3.*cos(q).*Dx.^3-8.*r2.^2.*cos(q).^2.*Dx.^4+6.*r2.*cos(q).*Dx.^5+6.*r2.^4.*sin(q).^2.*Dx.^2+4.*r2.^5.*sin(q).*Dy-2.*r2.^4.*sin(q).^2.*r1.^2+5.*Dx.^4.*r2.^2.*sin(q).^2+r2.^2.*sin(q).^2.*L1.^4-2.*r2.^4.*sin(q).^2.*L1.^2+4.*Dy.^2.*r2.^4.*sin(q).^4-4.*Dy.*r2.^5.*sin(q).^3+Dy.^4.*r2.^2.*sin(q).^2+4.*Dy.^3.*r2.^3.*sin(q)-2.*Dy.^2.*r2.^4.*sin(q).^2-4.*Dy.^3.*r2.^3.*sin(q).^3+4.*r2.^2.*Dx.^2.*r1.^2+4.*r2.^2.*Dx.^2.*L1.^2+2.*r2.^2.*r1.^2.*L1.^2-8.*Dx.^2.*Dy.^2.*r2.^2+2.*Dy.^4.*r2.*cos(q).*Dx+4.*r2.^3.*sin(q).^3.*L1.^2.*Dy-2.*r2.^2.*sin(q).^2.*L1.^2.*Dx.^2-4.*r2.^3.*sin(q).*L1.^2.*Dy-2.*r2.^2.*sin(q).^2.*L1.^2.*Dy.^2-2.*r2.^2.*sin(q).^2.*L1.^2.*r1.^2-2.*Dy.^2.*r2.^2.*sin(q).^2.*r1.^2+2.*Dy.^2.*r2.^2.*sin(q).^2.*Dx.^2-12.*Dy.*r2.^3.*sin(q).^3.*Dx.^2+4.*Dy.*r2.^3.*sin(q).^3.*r1.^2+4.*r2.^4.*cos(q).^2.*Dx.^2.*sin(q).^2-4.*r2.^5.*cos(q).*Dx.*sin(q).^2-12.*r2.^3.*cos(q).*Dx.^3.*sin(q).^2-4.*Dy.*r1.^2.*r2.^3.*sin(q)+8.*r2.^2.*cos(q).^2.*Dx.^2.*r1.^2+8.*r2.^2.*cos(q).^2.*Dx.^2.*L1.^2+2.*r2.*cos(q).*Dx.*r1.^4-7.*r2.^2.*Dx.^4+2.*r2.^4.*r1.^2-Dy.^4.*Dx.^2-Dy.^4.*r2.^2-2.*Dx.^4.*Dy.^2-2.*Dy.^2.*r2.^4-r2.^2.*r1.^4-r2.^2.*L1.^4-7.*r2.^4.*Dx.^2+2.*r2.^4.*L1.^2-r2.^6+2.*Dy.^2.*r1.^2.*Dx.^2+2.*Dy.^2.*r1.^2.*r2.^2+2.*Dy.^2.*L1.^2.*Dx.^2+2.*Dy.^2.*L1.^2.*r2.^2).^(1./2));


x22 = [1/2*(-1/(4*r2^2-8*Dy*r2*sin(q)-8*Dx*r2*cos(q)+4*Dy^2+4*Dx^2)*(-4*Dy^2*r2*sin(q)+4*Dx^2*Dy+4*Dy*L1^2+4*Dx^2*r2*sin(q)+4*Dy*r2^2+4*Dy^3+4*r1^2*r2*sin(q)+4*r2^3*sin(q)-4*r2*sin(q)*L1^2-8*Dy*r2^2*sin(q)^2-4*Dy*r1^2-8*r2^2*cos(q)*Dx*sin(q)-8*r2*cos(q)*Dx*Dy-4*(-2*Dx^4*Dy^2-Dx^2*Dy^4-2*Dy^2*r2^4-Dy^4*r2^2+2*r2^4*r1^2+2*r2^4*L1^2-r2^2*L1^4-7*r2^2*Dx^4-7*r2^4*Dx^2-r2^2*r1^4+2*Dx^4*r1^2+2*Dx^4*L1^2-Dx^2*L1^4-Dx^2*r1^4+r2^6*sin(q)^2+2*Dx^2*Dy^2*L1^2-8*Dx^2*Dy^2*r2^2+2*Dx^2*Dy^2*r1^2+2*Dy^2*L1^2*r2^2+2*Dy^2*r2^2*r1^2+4*r2^2*Dx^2*r1^2+4*r2^2*Dx^2*L1^2+2*r2^2*r1^2*L1^2-12*Dx^3*r2^3*sin(q)^2*cos(q)-4*Dy*r2^3*r1^2*sin(q)+8*Dy^2*r2^3*cos(q)*Dx+2*Dy^4*r2*cos(q)*Dx-2*r1^2*r2^2*sin(q)^2*L1^2+4*r1^2*r2^3*sin(q)^3*Dy+2*Dx^2*r1^2*L1^2+Dy^4*r2^2*sin(q)^2+4*Dy^3*r2^3*sin(q)-2*Dy^2*r2^4*sin(q)^2-4*Dy^3*r2^3*sin(q)^3+5*Dx^4*r2^2*sin(q)^2+6*Dx^2*r2^4*sin(q)^2+4*Dy*r2^5*sin(q)+r1^4*r2^2*sin(q)^2-2*r1^2*r2^4*sin(q)^2-2*r2^4*sin(q)^2*L1^2-4*r2^5*sin(q)^3*Dy+r2^2*sin(q)^2*L1^4+4*Dy^2*r2^4*sin(q)^4+6*r2^5*cos(q)*Dx+20*r2^3*cos(q)*Dx^3-8*Dx^2*r2^4*cos(q)^2-8*Dx^4*r2^2*cos(q)^2+6*Dx^5*r2*cos(q)-8*r2^3*cos(q)*Dx*L1^2+8*Dy*L1^2*r2^2*cos(q)*Dx*sin(q)-4*Dy^2*L1^2*r2*cos(q)*Dx-4*Dy^2*r1^2*r2*cos(q)*Dx+4*r2^3*sin(q)^2*L1^2*cos(q)*Dx+8*Dy*r2^4*sin(q)^3*cos(q)*Dx+8*r2^3*cos(q)^2*Dx^2*sin(q)*Dy-16*Dy*r2^4*cos(q)*Dx*sin(q)+4*r1^2*r2^3*sin(q)^2*cos(q)*Dx+8*r1^2*r2^2*sin(q)*cos(q)*Dx*Dy+4*Dy^3*r2*sin(q)*Dx^2+2*Dy^2*r2^2*sin(q)^2*Dx^2-2*Dy^2*r2^2*sin(q)^2*r1^2-2*Dy^2*r2^2*sin(q)^2*L1^2+4*Dx^4*Dy*r2*sin(q)+16*Dx^2*Dy*r2^3*sin(q)+8*Dx^3*Dy^2*r2*cos(q)-4*Dy*L1^2*r2^3*sin(q)-2*Dx^2*r2^2*sin(q)^2*r1^2-2*Dx^2*r2^2*sin(q)^2*L1^2-12*Dx^2*r2^3*sin(q)^3*Dy-8*Dx^3*r2*cos(q)*r1^2-8*Dx^3*r2*cos(q)*L1^2+8*Dx^2*r2^2*cos(q)^2*L1^2+8*Dx^2*r2^2*cos(q)^2*r1^2+2*Dx*r2*cos(q)*L1^4+2*Dx*r2*cos(q)*r1^4+4*r2^4*cos(q)^2*Dx^2*sin(q)^2-4*r2^2*cos(q)^2*Dx^2*Dy^2-8*r2^3*cos(q)*Dx*r1^2-4*r2^5*sin(q)^2*cos(q)*Dx+...
      4*r2^3*sin(q)^3*L1^2*Dy-r2^6-4*Dx*r2*cos(q)*r1^2*L1^2+4*Dy^2*r2^3*sin(q)^2*cos(q)*Dx-8*Dy^3*r2^2*sin(q)*cos(q)*Dx-4*Dx^2*Dy*r1^2*r2*sin(q)-4*Dx^2*Dy*r2*sin(q)*L1^2-16*Dx^3*Dy*r2^2*cos(q)*sin(q)-Dx^6)^(1/2))*r2*sin(q)-Dy^2-Dx^2-L1^2+r1^2+r2^2+1/(4*r2^2-8*Dy*r2*sin(q)-8*Dx*r2*cos(q)+4*Dy^2+4*Dx^2)*(-4*Dy^2*r2*sin(q)+4*Dx^2*Dy+4*Dy*L1^2+4*Dx^2*r2*sin(q)+4*Dy*r2^2+4*Dy^3+4*r1^2*r2*sin(q)+4*r2^3*sin(q)-4*r2*sin(q)*L1^2-8*Dy*r2^2*sin(q)^2-4*Dy*r1^2-8*r2^2*cos(q)*Dx*sin(q)-8*r2*cos(q)*Dx*Dy-4*(-2*Dx^4*Dy^2-Dx^2*Dy^4-2*Dy^2*r2^4-Dy^4*r2^2+2*r2^4*r1^2+2*r2^4*L1^2-r2^2*L1^4-7*r2^2*Dx^4-7*r2^4*Dx^2-r2^2*r1^4+2*Dx^4*r1^2+2*Dx^4*L1^2-Dx^2*L1^4-Dx^2*r1^4+r2^6*sin(q)^2+2*Dx^2*Dy^2*L1^2-8*Dx^2*Dy^2*r2^2+2*Dx^2*Dy^2*r1^2+2*Dy^2*L1^2*r2^2+2*Dy^2*r2^2*r1^2+4*r2^2*Dx^2*r1^2+4*r2^2*Dx^2*L1^2+2*r2^2*r1^2*L1^2-12*Dx^3*r2^3*sin(q)^2*cos(q)-4*Dy*r2^3*r1^2*sin(q)+8*Dy^2*r2^3*cos(q)*Dx+2*Dy^4*r2*cos(q)*Dx-2*r1^2*r2^2*sin(q)^2*L1^2+4*r1^2*r2^3*sin(q)^3*Dy+2*Dx^2*r1^2*L1^2+Dy^4*r2^2*sin(q)^2+4*Dy^3*r2^3*sin(q)-2*Dy^2*r2^4*sin(q)^2-4*Dy^3*r2^3*sin(q)^3+5*Dx^4*r2^2*sin(q)^2+6*Dx^2*r2^4*sin(q)^2+4*Dy*r2^5*sin(q)+r1^4*r2^2*sin(q)^2-2*r1^2*r2^4*sin(q)^2-2*r2^4*sin(q)^2*L1^2-4*r2^5*sin(q)^3*Dy+r2^2*sin(q)^2*L1^4+4*Dy^2*r2^4*sin(q)^4+6*r2^5*cos(q)*Dx+20*r2^3*cos(q)*Dx^3-8*Dx^2*r2^4*cos(q)^2-8*Dx^4*r2^2*cos(q)^2+6*Dx^5*r2*cos(q)-8*r2^3*cos(q)*Dx*L1^2+8*Dy*L1^2*r2^2*cos(q)*Dx*sin(q)-4*Dy^2*L1^2*r2*cos(q)*Dx-4*Dy^2*r1^2*r2*cos(q)*Dx+4*r2^3*sin(q)^2*L1^2*cos(q)*Dx+8*Dy*r2^4*sin(q)^3*cos(q)*Dx+8*r2^3*cos(q)^2*Dx^2*sin(q)*Dy-16*Dy*r2^4*cos(q)*Dx*sin(q)+4*r1^2*r2^3*sin(q)^2*cos(q)*Dx+8*r1^2*r2^2*sin(q)*cos(q)*Dx*Dy+4*Dy^3*r2*sin(q)*Dx^2+2*Dy^2*r2^2*sin(q)^2*Dx^2-2*Dy^2*r2^2*sin(q)^2*r1^2-2*Dy^2*r2^2*sin(q)^2*L1^2+4*Dx^4*Dy*r2*sin(q)+16*Dx^2*Dy*r2^3*sin(q)+8*Dx^3*Dy^2*r2*cos(q)-4*Dy*L1^2*r2^3*sin(q)-2*Dx^2*r2^2*sin(q)^2*r1^2-2*Dx^2*r2^2*sin(q)^2*L1^2-12*Dx^2*r2^3*sin(q)^3*Dy-8*Dx^3 ...
      *r2*cos(q)*r1^2-8*Dx^3*r2*cos(q)*L1^2+8*Dx^2*r2^2*cos(q)^2*L1^2+8*Dx^2*r2^2*cos(q)^2*r1^2+2*Dx*r2*cos(q)*L1^4+2*Dx*r2*cos(q)*r1^4+4*r2^4*cos(q)^2*Dx^2*sin(q)^2-4*r2^2*cos(q)^2*Dx^2*Dy^2-8*r2^3*cos(q)*Dx*r1^2-4*r2^5*sin(q)^2*cos(q)*Dx+4*r2^3*sin(q)^3*L1^2*Dy-r2^6-4*Dx*r2*cos(q)*r1^2*L1^2+4*Dy^2*r2^3*sin(q)^2*cos(q)*Dx-8*Dy^3*r2^2*sin(q)*cos(q)*Dx-4*Dx^2*Dy*r1^2*r2*sin(q)-4*Dx^2*Dy*r2*sin(q)*L1^2-16*Dx^3*Dy*r2^2*cos(q)*sin(q)-Dx^6)^(1/2))*Dy)/(-Dx+r2*cos(q))];
y22 = [ 1/2/(4*r2^2-8*Dy*r2*sin(q)-8*Dx*r2*cos(q)+4*Dy^2+4*Dx^2)*(-4*Dy^2*r2*sin(q)+4*Dx^2*Dy+4*Dy*L1^2+4*Dx^2*r2*sin(q)+4*Dy*r2^2+4*Dy^3+4*r1^2*r2*sin(q)+4*r2^3*sin(q)-4*r2*sin(q)*L1^2-8*Dy*r2^2*sin(q)^2-4*Dy*r1^2-8*r2^2*cos(q)*Dx*sin(q)-8*r2*cos(q)*Dx*Dy-4*(-2*Dx^4*Dy^2-Dx^2*Dy^4-2*Dy^2*r2^4-Dy^4*r2^2+2*r2^4*r1^2+2*r2^4*L1^2-r2^2*L1^4-7*r2^2*Dx^4-7*r2^4*Dx^2-r2^2*r1^4+2*Dx^4*r1^2+2*Dx^4*L1^2-Dx^2*L1^4-Dx^2*r1^4+r2^6*sin(q)^2+2*Dx^2*Dy^2*L1^2-8*Dx^2*Dy^2*r2^2+2*Dx^2*Dy^2*r1^2+2*Dy^2*L1^2*r2^2+2*Dy^2*r2^2*r1^2+4*r2^2*Dx^2*r1^2+4*r2^2*Dx^2*L1^2+2*r2^2*r1^2*L1^2-12*Dx^3*r2^3*sin(q)^2*cos(q)-4*Dy*r2^3*r1^2*sin(q)+8*Dy^2*r2^3*cos(q)*Dx+2*Dy^4*r2*cos(q)*Dx-2*r1^2*r2^2*sin(q)^2*L1^2+4*r1^2*r2^3*sin(q)^3*Dy+2*Dx^2*r1^2*L1^2+Dy^4*r2^2*sin(q)^2+4*Dy^3*r2^3*sin(q)-2*Dy^2*r2^4*sin(q)^2-4*Dy^3*r2^3*sin(q)^3+5*Dx^4*r2^2*sin(q)^2+6*Dx^2*r2^4*sin(q)^2+4*Dy*r2^5*sin(q)+r1^4*r2^2*sin(q)^2-2*r1^2*r2^4*sin(q)^2-2*r2^4*sin(q)^2*L1^2-4*r2^5*sin(q)^3*Dy+r2^2*sin(q)^2*L1^4+4*Dy^2*r2^4*sin(q)^4+6*r2^5*cos(q)*Dx+20*r2^3*cos(q)*Dx^3-8*Dx^2*r2^4*cos(q)^2-8*Dx^4*r2^2*cos(q)^2+6*Dx^5*r2*cos(q)-8*r2^3*cos(q)*Dx*L1^2+8*Dy*L1^2*r2^2*cos(q)*Dx*sin(q)-4*Dy^2*L1^2*r2*cos(q)*Dx-4*Dy^2*r1^2*r2*cos(q)*Dx+4*r2^3*sin(q)^2*L1^2*cos(q)*Dx+8*Dy*r2^4*sin(q)^3*cos(q)*Dx+8*r2^3*cos(q)^2*Dx^2*sin(q)*Dy-16*Dy*r2^4*cos(q)*Dx*sin(q)+4*r1^2*r2^3*sin(q)^2*cos(q)*Dx+8*r1^2*r2^2*sin(q)*cos(q)*Dx*Dy+4*Dy^3*r2*sin(q)*Dx^2+2*Dy^2*r2^2*sin(q)^2*Dx^2-2*Dy^2*r2^2*sin(q)^2*r1^2-2*Dy^2*r2^2*sin(q)^2*L1^2+4*Dx^4*Dy*r2*sin(q)+16*Dx^2*Dy*r2^3*sin(q)+8*Dx^3*Dy^2*r2*cos(q)-4*Dy*L1^2*r2^3*sin(q)-2*Dx^2*r2^2*sin(q)^2*r1^2-2*Dx^2*r2^2*sin(q)^2*L1^2-12*Dx^2*r2^3*sin(q)^3*Dy-8*Dx^3*r2*cos(q)*r1^2-8*Dx^3*r2*cos(q)*L1^2+8*Dx^2*r2^2*cos(q)^2*L1^2+8*Dx^2*r2^2*cos(q)^2*r1^2+2*Dx*r2*cos(q)*L1^4+2*Dx*r2*cos(q)*r1^4+4*r2^4*cos(q)^2*Dx^2*sin(q)^2-4*r2^2*cos(q)^2*Dx^2*Dy^2-8*r2^3*cos(q)*Dx*r1^2-4*r2^5*sin(q)^2*cos(q)*Dx+4*r2^3*sin(q)^3*L1^2*Dy-r2^6-4*Dx*r2*cos(q)*r1^2*L1^2+4*Dy^2*r2^3*sin(q)^2*cos(q)*Dx-8*Dy^3*r2^2*sin(q)*cos(q)*Dx-4*Dx^2*Dy*r1^2*r2*sin(q)-4*Dx^2*Dy*r2*sin(q)*L1^2-16*Dx^3*Dy*r2^2*cos(q)*sin(q)-Dx^6)^(1/2))]; 

xp= subs(x21);
yp= subs(y21);


clf
axis equal
   
Dxx = Dx*ones(size(xp));
Dyy = Dy*ones(size(yp));
plot([Dxx' xp']', [Dyy' yp']');
   