function [I1, I2, status] = circlesIntersectionPoints(P1, r1, P2, r2)
% circlesIntersectionPoints  Compute the intersection points of two circle
%
%   [I1, I2, status] = circlesIntersectionPoints(P1, r1, P2, r2)
%
%Input parameters
%   
%   P1 is the center of circle 1
%   r1 is the radius of circle 1
%   P2 is the center of circle 2
%   r2 is the radius of circle 1
%
%   
%Output parameters
%   
%   I1 and I2 are the coordinates of intersection points
%   status is the status of the intersection:
%       * -1 : Circles are too far apart and do not intersect;
%       * -2 : One circle is inside the other and do not intersect;
%       * -3 : Circles are merged , there are an infinite number of points
%       * 1  : Single intersection point
%       * 2  : two intersection points
%
%Usage example:
%
%   [i1, i2, s] = circlesIntersectionPoints([3,4], 3.4, [-1,2], 2.6)
%
%
%   Detailed explaination on https://lucidar.me/en/mathematics/how-to-calculate-the-intersection-points-of-two-circles/
%   See more on https://lucidar.me/en/matlab/how-to-calculate-intersection-points-between-two-circles/
%

I1=0;
I2=0;

%% Distance between centers
d = sqrt( (P2(1)-P1(1))*(P2(1)-P1(1)) + (P2(2)-P1(2))*(P2(2)-P1(2)) );

%% Check cases
if (d>r1+r2)
    % Circles are too far apart and do not intersect;
    status = -1;
    return
end

if (d<abs(r1-r2))
    % One circle is inside the other and do not intersect;
    status = -2;
    return
end

if (d==0 && r1==r2)
    % Circles are merged 
    status = -3;
    return
end

if (d==r1+r2)
    % Single intersection point;
    status = 1;
end

if (d<r1+r2)
    % Two intersection points.
    status = 1;
end

%% a and b
a = (r1*r1 - r2*r2 + d*d) / (2*d);
b = (r2*r2 - r1*r1 + d*d) / (2*d);

%% h
h=sqrt(r1*r1 - a*a);

%% P5
P5 = [ P1(1) + (a/d)*(P2(1)-P1(1)) , P1(2) + (a/d)*(P2(2)-P1(2)) ];

%% Intersection points
I1 = [ P5(1) - h*(P2(2)-P1(2))/d ; P5(2) + h*(P2(1)-P1(1))/d];
I2 = [ P5(1) + h*(P2(2)-P1(2))/d ; P5(2) - h*(P2(1)-P1(1))/d];
end