##############################
# Calcule les ou le point d'intersection de deux cercle :
# Cercle 1 de coordonnés (a1,b1) de rayon d1
# Cercle 2 de coordonnés (a2,b2) de rayon d2
# Precision = nombre de chiffre aprés la virgule
##############################
getIntersection=function(a1,b1,d1,a2,b2,d2,precision=5){
  # Compute E (cf paper)
  computeE=function(a1,b1,a2,b2){
    numerateur=-2*a1+2*a2;
    denominateur=-2*b1+2*b2;
    return(numerateur/denominateur);
  }
  
  # Compute F (cf paper)
  computeF=function(a1,b1,d1,a2,b2,d2){
    numerateur=a1^2+b1^2-a2^2-b2^2-d1^2+d2^2;
    denominateur=-2*b1+2*b2;
    return(numerateur/denominateur);
  }
  
  # Compute intersections if b1 != b2
  computeYNotEqual=function(a1,b1,d1,a2,b2,d2){
    E=computeE(a1,b1,a2,b2);
    F=computeF(a1,b1,d1,a2,b2,d2);
    A=E^2+1;
    B=(-2*a1+2*E*F+2*E*b1);
    C=2*F*b1+b1^2+a1^2+F^2-d1^2;
    
    values=polyroot(c(C,B,A));
    x=c();
    y=c();
    for(i in 1:length(values)){
      if(round(Im(values[i]),digit=precision)==0){
        x=c(x,Re(values[i]));
        y=c(y,-Re(values[i])*E-F);
      }
    }
    if(length(x)==0){
      return(NULL);
    }
    return(matrix(c(x,y),nrow=length(x),ncol=2));
  }
  
  # Compute intersections if b1 == b2 and a1 != a2
  computeYEqualAndXNotEqual=function(a1,b1,d1,a2,b2,d2){
    G=-((a1^2-a2^2-d1^2+d2^2)/(-2*a1+2*a2));
    A=1;
    B=-2*b1;
    C=G^2-2*G*a1+a1^2+b1^2-d1^2;
    values=polyroot(c(C,B,A));
    x=c();
    y=c();
    for(i in 1:length(values)){
      if(round(Im(values[i]),digit=precision)==0){
        x=c(x,G);
        y=c(y,Re(values[i]));
      }
    }
    if(length(y)==0){
      return(NULL);
    }
    return(matrix(c(x,y),nrow=length(y),ncol=2));
    
  }
  
  # Compute intersections
  if(b1!=b2){
    return(computeYNotEqual(a1,b1,d1,a2,b2,d2));
  }
  else if(a1!=a2){
    return(computeYEqualAndXNotEqual(a1,b1,d1,a2,b2,d2));
  }
  # No intersection found
  return(NULL);
}