r1 = normal(S(1,1),M1(1,1),S(1,2),M1(1,2));
r2 = normal(S(1,1),M2(1,1),S(1,2),M2(1,2));
r3 = normal(S(1,1),M3(1,1),S(1,2),M3(1,2));
Z1 =intersect2(x,y,M1(1),M1(2),r1);
Z2=intersect2(x,y,M2(1),M2(2),r2);
Z3=intersect2(x,y,M3(1),M3(2),r3);
n1(i) = norm(Z2(i)-Z1(i));
n2(i) = norm(Z3(i)-Z2(i));
n3(i) = norm(Z1(i)-Z3(i));
[xm1,ym1] = circle(M1(1,1),M1(1,2),r1);
[xm2,ym2] = circle(M2(1,1),M2(1,2),r2);
[xm3,ym3] = circle(M3(1,1),M3(1,2),r3);
X = reshape(ZZ(:,1),[v,v]);
Y = reshape(ZZ(:,2),[v,v]);
contour(X,Y,Z,30,'ShowText','off')
function [xunit, yunit] = circle(x,y,r)
function d = normal(x1,x2,y1,y2)
d = sqrt((x2-x1).^2+(y2-y1).^2);
function [Z] = intersect2(x,y,cx,cy,r)
theta = atan2((y-cy) ,(x-cx));