function [ip1,ip2,status] = NE_int_2circles(circle1,circle2)
% NE_int_2circles - find intersection points of 2 circles
% On input:
%     circle1 (1x3 vector): center and radius of circle 1
%     circle2 (1x3 vector): center and radius of circle 2
% On output:
%     ip1 (1x2 vector): intersection point 1
%     ip2 (1x2 vector): intersection point 2
%     status (int): 0 if no intersection
%                   1 if 1 point
%                   2 if 2 points;
%                   3 if same circle
% Call:
%     [p1,p2,st] = NE_int_2circles([0,0,2],[1,0,1]);
% Author:
%     T. Henderson
%     UU
%     Summer 2023
%

ZERO_THRESH = 0.001;

ip1 = [];
ip2 = [];
status = 0;

C1 = circle1(1:2);
r1 = circle1(3);
C2 = circle2(1:2);
r2 = circle2(3);
d = norm(C1-C2);

% both great circles through origin
if isinf(r1)&isinf(r2)
    ip1 = [0,0];
    ip2 = [0,0];
    status = 1;
    return
end

% No intersection
if d>r1+r2|d+min(r1,r2)<max(r1,r2)|...  % non-overlapping
        d+r1<r2|d+r2<r1                 % one inside other
    return
end

% Same circle
if norm(C2-C1)<ZERO_THRESH&norm(r2-r1)<ZERO_THRESH
    ip1 = [C1(1)+r1,C1(2)];
    ip2 = [-(C1(1)+r1),C1(2)];
    status = 3;
    return
end

% 1 or 2 intersection points
[ip1,ip2,status] = NE_int_E2_2circles(circle1,circle2);
d1 = norm(ip1);
d2 = norm(ip2);
status = 2;
if d1<1
    ip2 = ip1;
else
    ip1 = ip2;
end
