/* circle_circle.c find intersections of two circles      */
/* known x0, y0, r0 for circle 0, x1, y1, r1 for circle 1 */
/* compute d, length of line between centers              */
/* d = sqrt((x1-x0)^2 + (y1-y0)^2)                        */
/* quit if d > r0+r1,  if d < |r0-r1| and other cases     */
/* assume a point 2, x2, y2, on line between centers      */
/* distance a from (x0,y0) to (x2,y2)                     */
/* a = (r0^2 - r1^2 + d^2)/(2*d)                          */
/* x2 = x0 + a*(x1-x0)/d                                  */
/* y2 = y0 + a*(y1-y0)/d                                  */
/* distance h from intersection to a on (x2,y2)           */
/* h = sqrt(r0^2 - a^2)                                   */
/* solution, intersection (x3,y3) is                      */
/* x3 = x2 +/- h*(y1-y0)/d    h dx/dy                     */
/* y3 = y3 +/- h*(x1-x0)/d    h dy/dx                     */
/* check 4 possible cases against circle 0 and circle 1   */
#include <stdio.h>
#include <math.h>
#undef abs
#define abs(x) ((x)<0.0?(-(x)):(x))

int main(int argc, char *argv[])
{
  double x0=-0.1, y0=-0.2, r0=2.0;
  double x1=2.0, y1=1.0, r1=2.1;
  double x2, y2, x30, y30, x31, y31, a, d, h;
  double x3p, x3m, y3p, y3m, err;
  double b;

  printf("find circle to circle intersection \n");
  printf("circle 0 at (%f,%f) with radius %f \n", x0, y0, r0);
  printf("circle 1 at (%f,%f) with radius %f \n", x1, y1, r1);

  d = sqrt((x1-x0)*(x1-x0) + (y1-y0)*(y1-y0));
  printf("distance d between centers %f \n", d);

  if(d > r0+r1 || d <= abs(r0-r1))
    printf("no intersection, thus expect no solution \n");
  /* continue, assume solution exists */

  a = (r0*r0 - r1*r1 + d*d)/(2.0*d);
  printf("distance a (x0,y0) to (x2,y2) %f \n", a);

  x2 = x0 + a*(x1-x0)/d;
  y2 = y0 + a*(y1-y0)/d;
  printf("point (x2,y2) on line between centers (%f,%f) \n", x2, y2);

  h = sqrt(r0*r0 - a*a);
  printf("distance h from intersection to a on (x2,y2) %f \n", h);

  /* just for checking 
  err = a*a + h*h - r0*r0;
  printf("check a^2+h^2=r0^2 err=%e \n", err);
  b = d - a;
  err = b*b + h*h - r1*r1;
  printf("b=%f, check b^2+h^2-r1^2 err=%e \n", b, err);
  */

  x3p = x2 + h*(y1-y0)/d;
  x3m = x2 - h*(y1-y0)/d;
  y3p = y2 + h*(x1-x0)/d;
  y3m = y2 - h*(x1-x0)/d;

  /* check 4 cases */
  x30 = x3p-x0; /* relative to center */
  y30 = y3p-y0;
  x31 = x3p-x1;
  y31 = y3p-y1;
  err = abs((x30*x30 + y30*y30)-r0*r0)+abs((x31*x31 + y31*y31)-r1*r1);
  if(err>1.0e-5)
    printf("try (x3p,y3p) (%f,%f) err=%e \n", x3p, y3p, err);
  else
    printf(" a solution, intersection (x3,y3) is (%f,%f) \n", x3p, y3p);

  x30 = x3p-x0;
  y30 = y3m-y0;
  x31 = x3p-x1;
  y31 = y3m-y1;
  err = abs((x30*x30 + y30*y30)-r0*r0)+abs((x31*x31 + y31*y31)-r1*r1);
  if(err>1.0e-5)
    printf("try (x3p,y3m) (%f,%f) err=%e \n", x3p, y3m, err);
  else
    printf(" a solution, intersection (x3,y3) is (%f,%f) \n", x3p, y3m);

  x30 = x3m-x0;
  y30 = y3p-y0;
  x31 = x3m-x1;
  y31 = y3p-y1;
  err = abs((x30*x30 + y30*y30)-r0*r0)+abs((x31*x31 + y31*y31)-r1*r1);
  if(err>1.0e-5)
    printf("try (x3m,y3p) (%f,%f) err=%e \n", x3m, y3p, err);
  else
    printf(" a solution, intersection (x3,y3) is (%f,%f) \n", x3m, y3p);

  x30 = x3m-x0;
  y30 = y3m-y0;
  x31 = x3m-x1;
  y31 = y3m-y1;
  err = abs((x30*x30 + y30*y30)-r0*r0)+abs((x31*x31 + y31*y31)-r1*r1);
  if(err>1.0e-5)
    printf("try (x3m,y3m) (%f,%f) err=%e \n", x3m, y3m, err);
  else
    printf(" a solution, intersection (x3,y3) is (%f,%f) \n", x3m, y3m);

return 0;
} /* end circle_to_circle.c */