circle_intersection(10)

f(-1,-1,-2,-2,5);

function []=f(x1,y1,x2,y2,a)

if(x1>x2)
    xp=x1;x1=x2;x2=xp;
end

if(y1>y2)
    yp=y1;y1=y2;y2=yp;
end

x=x1-a:0.1:x2+a;
y=y1-a:0.1:y2+a;;
f=@(x,y) 4.*((x-x1).*(x-x1)+(y-y1).*(y-y1)).*((x-x2).*(x-x2)+(y-y2).*(y-y2))-(a.*a-((x-x1).*(x-x1)+(y-y1).*(y-y1))-((x-x2).*(x-x2)+(y-y2).*(y-y2))).^2
fimplicit(f);

%%%Try with rotation +  translation

% a=a/2.;
%
% %%Pretpostavka x1,y1 su blize kooridnatnom pocetku
%
% if(x1^2+y1^2>x2^2+y2^2)
%     xp=x1;x1=x2;x2=xp;
%     yp=y1;y1=y2;y2=yp;
% end
% translacija_do_centra=[0-x1,0-y1];
% disp(translacija_do_centra);
% x1=x1+translacija_do_centra(1);
% y1=y1+translacija_do_centra(2);
% x2=x2+translacija_do_centra(1);
% y2=y2+translacija_do_centra(2);
% %%Now we will rotate y2,x2 to x-axis
%
% % plot(x1,y1,'*');
% % hold on;
% % plot(x2,y2,'*');
% % hold off;
%
%
% ugao=atan(y2/x2);%%Pocetni nagib
% if(y2<0 && x2<0)
%     ugao=-ugao;
% end
% udaljenost=sqrt(x2^2+y2^2);
% x2=sqrt(x2^2+y2^2);%%Sad su na x:osi
%
%
% disp(ugao);
%
% %%Translirajmo da bude simetricno u odnosu na y osu:
%
% x1=x1-udaljenost/2.;
% x2=x2-udaljenost/2.;
%
% % plot(x1,y1,'*');
% % hold on;
% % plot(x2,y2,'*');
% % hold off;
%
%
% e=udaljenost/2.;
% disp(e);
% disp(a);
% b=sqrt(a^2-e^2);
% if(a^2-e^2<=0)fprintf("Ne posoje takve tacke");return;end
%
% fi=-pi:0.01:pi;
%
% x=a*cos(fi);
% y=b*sin(fi);
%
% %%Vracajmo unazad
%
% x1=x1+udaljenost/2.;
% x2=x2+udaljenost/2.;
% x=x+udaljenost/2.;
% x2=sqrt(x2^2+y2^2)/cos(ugao);
% disp(atan(y.\x))
%
% R  = [cos(ugao) -sin(ugao); ...
%       sin(ugao)  cos(ugao)];
% rCoords = R*[x ; y];
% xr = rCoords(1,:)';
% yr = rCoords(2,:)';
%
% xr=xr-translacija_do_centra(1);
% yr=yr-translacija_do_centra(2);
%
% x1=x1-translacija_do_centra(1);
% y1=y1-translacija_do_centra(2);
% x2=x2-translacija_do_centra(1);
% y2=y2-translacija_do_centra(2);
%
%
% plot(x1,y1,'*');
% grid on;
% hold on;
% plot(x2,y2,'*');
% plot(xr,yr);

end