1 files changed, 77 insertions, 0 deletions
diff --git a/R/intersectionCercles.R b/R/intersectionCercles.R
new file mode 100755
index 0000000..8dd2234
--- /dev/null
+++ b/R/intersectionCercles.R
@@ -0,0 +1,77 @@
+
+##############################
+# 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);
+}
+\ No newline at end of file