From 3f16d970465aafed6d7332d01e313fc92d5e8db9 Mon Sep 17 00:00:00 2001
From: pdn <pdn@opencascade.com>
Date: Thu, 13 Jun 2013 15:28:55 +0400
Subject: [PATCH] 0023939: Incorrect circle parameter in IntAna

Fix for circle circle intersection in case of one point touching
Test command added
---
 .../GeomliteTest_API2dCommands.cxx            | 48 +++++++++++++++++++
 src/IntAna2d/IntAna2d_AnaIntersection_2.cxx   | 12 ++++-
 src/IntTools/IntTools_FaceFace.cxx            |  1 +
 3 files changed, 59 insertions(+), 2 deletions(-)

diff --git a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx
index 83b40f1d41..dc4821ddab 100755
--- a/src/GeomliteTest/GeomliteTest_API2dCommands.cxx
+++ b/src/GeomliteTest/GeomliteTest_API2dCommands.cxx
@@ -42,6 +42,9 @@
 #include <TColStd_Array1OfReal.hxx>
 #include <GeomAbs_Shape.hxx>
 #include <Precision.hxx>
+#include <Geom2d_Circle.hxx>
+#include <IntAna2d_AnaIntersection.hxx>
+#include <IntAna2d_IntPoint.hxx>
 
 #include <stdio.h>
 #ifdef WNT
@@ -336,6 +339,48 @@ static Standard_Integer intersect(Draw_Interpretor& di, Standard_Integer n, cons
   return 0;
 }
 
+//=======================================================================
+//function : intersect
+//purpose  : 
+//=======================================================================
+
+static Standard_Integer intersect_ana(Draw_Interpretor& di, Standard_Integer n, const char** a)
+{
+  if( n < 2) 
+  {
+    cout<< "2dintana circle circle "<<endl;
+    return 1;
+  }
+  
+  Handle(Geom2d_Curve) C1 = DrawTrSurf::GetCurve2d(a[1]);
+  if ( C1.IsNull() && !C1->IsKind(STANDARD_TYPE(Geom2d_Circle))) 
+    return 1;
+
+  Handle(Geom2d_Curve) C2 = DrawTrSurf::GetCurve2d(a[2]);
+  if ( C2.IsNull() && !C2->IsKind(STANDARD_TYPE(Geom2d_Circle)))
+    return 1;
+
+  Handle(Geom2d_Circle) aCir1 = Handle(Geom2d_Circle)::DownCast(C1);
+  Handle(Geom2d_Circle) aCir2 = Handle(Geom2d_Circle)::DownCast(C2);
+
+  IntAna2d_AnaIntersection Intersector(aCir1->Circ2d(), aCir2->Circ2d());
+
+  Standard_Integer i;
+
+  for ( i = 1; i <= Intersector.NbPoints(); i++) {
+    gp_Pnt2d P = Intersector.Point(i).Value();
+    di<<"Intersection point "<<i<<" : "<<P.X()<<" "<<P.Y()<<"\n";
+	di<<"parameter on the fist: "<<Intersector.Point(i).ParamOnFirst();
+	di<<" parameter on the second: "<<Intersector.Point(i).ParamOnSecond()<<"\n";
+    Handle(Draw_Marker2D) mark = new Draw_Marker2D( P, Draw_X, Draw_vert); 
+    dout << mark;
+  }
+  dout.Flush();
+
+  return 0;
+}
+
+
 
 void GeomliteTest::API2dCommands(Draw_Interpretor& theCommands)
 {
@@ -365,4 +410,7 @@ void GeomliteTest::API2dCommands(Draw_Interpretor& theCommands)
 
   theCommands.Add("2dintersect", "intersect curve curve [Tol]",__FILE__,
 		  intersect,g);
+
+  theCommands.Add("2dintanalytical", "intersect curve curve using IntAna",__FILE__,
+		  intersect_ana,g);
 }
diff --git a/src/IntAna2d/IntAna2d_AnaIntersection_2.cxx b/src/IntAna2d/IntAna2d_AnaIntersection_2.cxx
index a0b97aa2ab..74694c9c1f 100755
--- a/src/IntAna2d/IntAna2d_AnaIntersection_2.cxx
+++ b/src/IntAna2d/IntAna2d_AnaIntersection_2.cxx
@@ -66,7 +66,7 @@ void IntAna2d_AnaIntersection::Perform (const gp_Circ2d& C1,
     if (ang1<0) {ang1=2*M_PI+ang1;}                // On revient entre 0 et 2PI
     lpnt[0].SetValue(XS,YS,ang1,ang2);
   }
-  else if (((sum-d)>Epsilon(d)) && ((d-dif)>Epsilon(d))) {
+  else if (((sum-d)>Epsilon(sum)) && ((d-dif)>Epsilon(sum))) {
     empt=Standard_False;
     para=Standard_False;
     iden=Standard_False;
@@ -78,6 +78,11 @@ void IntAna2d_AnaIntersection::Perform (const gp_Circ2d& C1,
     Standard_Real ref2=Ox2.Angle(ax);                       // Resultat entre -PI et +PI
 
     Standard_Real l1=(d*d + R1*R1 -R2*R2)/(2.0*d);
+    Standard_Real aDet = R1*R1-l1*l1;
+    if(aDet < 0.) {
+      aDet = 0.;
+      l1 = (l1 > 0 ? R1 : - R1);
+    }
     Standard_Real h= Sqrt(R1*R1-l1*l1);
 
     Standard_Real XS1= C1.Location().X() + l1*ax.X()/d - h*ax.Y()/d;
@@ -143,12 +148,15 @@ void IntAna2d_AnaIntersection::Perform (const gp_Circ2d& C1,
     lpnt[0].SetValue(XS1,YS1,ang11,ang21);
     lpnt[1].SetValue(XS2,YS2,ang12,ang22);
   }
-  else if (Abs(d-dif)<=Epsilon(d)) {    // Cercles tangents interieurs
+  else if (Abs(d-dif)<=Epsilon(sum)) {    // Cercles tangents interieurs
     empt=Standard_False;
     para=Standard_False;
     iden=Standard_False;
     nbp=1;
     gp_Vec2d ax(C1.Location(),C2.Location());
+    if(C1.Radius() < C2.Radius())
+      ax.Reverse();
+
     gp_Vec2d Ox1(C1.XAxis().Direction());
     gp_Vec2d Ox2(C2.XAxis().Direction());
     Standard_Real ang1=Ox1.Angle(ax);                       // Resultat entre -PI et +PI
diff --git a/src/IntTools/IntTools_FaceFace.cxx b/src/IntTools/IntTools_FaceFace.cxx
index 76f81f487e..2ef5436edd 100755
--- a/src/IntTools/IntTools_FaceFace.cxx
+++ b/src/IntTools/IntTools_FaceFace.cxx
@@ -330,6 +330,7 @@ static
 //=======================================================================
 IntTools_FaceFace::IntTools_FaceFace()
 {
+  myIsDone=Standard_False;
   myTangentFaces=Standard_False;
   //
   myHS1 = new GeomAdaptor_HSurface ();