src/IntAna/IntAna_QuadQuadGeo.cxx

   1 // Created on: 1992-08-06
   2 // Created by: Laurent BUCHARD
   3 // Copyright (c) 1992-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 //----------------------------------------------------------------------
  18 //-- Purposse: Geometric Intersection between two Natural Quadric
  19 //--          If the intersection is not a conic,
  20 //--          analytical methods must be called.
  21 //----------------------------------------------------------------------
  22 #ifndef OCCT_DEBUG
  23 #define No_Standard_RangeError
  24 #define No_Standard_OutOfRange
  25 #endif
  26
  27 #include <IntAna_QuadQuadGeo.ixx>
  28
  29 #include <IntAna_IntConicQuad.hxx>
  30 #include <StdFail_NotDone.hxx>
  31 #include <Standard_DomainError.hxx>
  32 #include <Standard_OutOfRange.hxx>
  33 #include <math_DirectPolynomialRoots.hxx>
  34
  35 #include <gp.hxx>
  36 #include <gp_Pln.hxx>
  37 #include <gp_Vec.hxx>
  38 #include <ElSLib.hxx>
  39 #include <ElCLib.hxx>
  40
  41 #include <gp_Dir.hxx>
  42 #include <gp_XYZ.hxx>
  43 #include <gp_Pnt2d.hxx>
  44 #include <gp_Vec2d.hxx>
  45 #include <gp_Dir2d.hxx>
  46
  47
  48 static
  49   gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
  50 static
  51   void RefineDir(gp_Dir& aDir);
  52
  53 //=======================================================================
  54 //class :  AxeOperator
  55 //purpose  : O p e r a t i o n s   D i v e r s e s  s u r   d e s   A x 1
  56 //=======================================================================
  57 class AxeOperator {
  58  public:
  59   AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
  60
  61   void Distance(Standard_Real& dist,
  62                 Standard_Real& Param1,
  63                 Standard_Real& Param2);
  64
  65   gp_Pnt PtIntersect()              {
  66     return ptintersect;
  67   }
  68   Standard_Boolean Coplanar(void)   {
  69     return thecoplanar;
  70   }
  71   Standard_Boolean Same(void)       {
  72     return theparallel && (thedistance<myEPSILON_DISTANCE);
  73   }
  74   Standard_Real Distance(void)      {
  75     return thedistance ;
  76   }
  77   Standard_Boolean Intersect(void)  {
  78     return (thecoplanar && (!theparallel));
  79   }
  80   Standard_Boolean Parallel(void)   {
  81     return theparallel;
  82   }
  83   Standard_Boolean Normal(void)     {
  84     return thenormal;
  85   }
  86
  87  protected:
  88   Standard_Real Det33(const Standard_Real a11,
  89                       const Standard_Real a12,
  90                       const Standard_Real a13,
  91                       const Standard_Real a21,
  92                       const Standard_Real a22,
  93                       const Standard_Real a23,
  94                       const Standard_Real a31,
  95                       const Standard_Real a32,
  96                       const Standard_Real a33) {
  97     Standard_Real theReturn =
  98       a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;
  99     return theReturn ;
 100   }
 101
 102  private:
 103   gp_Pnt ptintersect;
 104   gp_Ax1 Axe1;
 105   gp_Ax1 Axe2;
 106   Standard_Real thedistance;
 107   Standard_Boolean theparallel;
 108   Standard_Boolean thecoplanar;
 109   Standard_Boolean thenormal;
 110   //
 111   Standard_Real myEPSILON_DISTANCE;
 112   Standard_Real myEPSILON_AXES_PARA;
 113 };
 114
 115 //=======================================================================
 116 //function : AxeOperator::AxeOperator
 117 //purpose  :
 118 //=======================================================================
 119 AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2)
 120 {
 121   myEPSILON_DISTANCE=0.00000000000001;
 122   myEPSILON_AXES_PARA=0.000000000001;
 123   Axe1=A1;
 124   Axe2=A2;
 125   //---------------------------------------------------------------------
 126   gp_Dir V1=Axe1.Direction();
 127   gp_Dir V2=Axe2.Direction();
 128   gp_Pnt P1=Axe1.Location();
 129   gp_Pnt P2=Axe2.Location();
 130   //
 131   RefineDir(V1);
 132   RefineDir(V2);
 133   thecoplanar= Standard_False;
 134   thenormal  = Standard_False;
 135
 136   //--- check if the two axis are parallel
 137   theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);
 138   //--- Distance between the two axis
 139   gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
 140   if (theparallel) {
 141     gp_Lin L1(A1);
 142     thedistance = L1.Distance(A2.Location());
 143   }
 144   else {
 145     thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
 146                                                            Axe2.Location())));
 147   }
 148   //--- check if Axis are Coplanar
 149   Standard_Real D33;
 150   if(thedistance<myEPSILON_DISTANCE) {
 151     D33=Det33(V1.X(),V1.Y(),V1.Z()
 152               ,V2.X(),V2.Y(),V2.Z()
 153               ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
 154     if(Abs(D33)<=myEPSILON_DISTANCE) {
 155       thecoplanar=Standard_True;
 156     }
 157   }
 158   else {
 159     thecoplanar=Standard_True;
 160     thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
 161   }
 162   //--- check if the two axis are concurrent
 163   if(thecoplanar && (!theparallel)) {
 164     Standard_Real smx=P2.X() - P1.X();
 165     Standard_Real smy=P2.Y() - P1.Y();
 166     Standard_Real smz=P2.Z() - P1.Z();
 167     Standard_Real Det1,Det2,Det3,A;
 168     Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
 169     Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
 170     Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
 171
 172     if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
 173       A=(smy * V2.X() - smx * V2.Y())/Det1;
 174     }
 175     else if((Det2!=0.0)
 176             && ((Abs(Det2) >= Abs(Det1))
 177                 &&(Abs(Det2) >= Abs(Det3)))) {
 178       A=(smz * V2.Y() - smy * V2.Z())/Det2;
 179     }
 180     else {
 181       A=(smz * V2.X() - smx * V2.Z())/Det3;
 182     }
 183     ptintersect.SetCoord( P1.X() + A * V1.X()
 184                          ,P1.Y() + A * V1.Y()
 185                          ,P1.Z() + A * V1.Z());
 186   }
 187   else {
 188     ptintersect.SetCoord(0,0,0);  //-- Pour eviter des FPE
 189   }
 190 }
 191 //=======================================================================
 192 //function : Distance
 193 //purpose  :
 194 //=======================================================================
 195 void AxeOperator::Distance(Standard_Real& dist,
 196                            Standard_Real& Param1,
 197                            Standard_Real& Param2)
 198  {
 199   gp_Vec O1O2(Axe1.Location(),Axe2.Location());
 200   gp_Dir U1 = Axe1.Direction();   //-- juste pour voir.
 201   gp_Dir U2 = Axe2.Direction();
 202
 203   gp_Dir N  = U1.Crossed(U2);
 204   Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
 205                           U1.Y(),U2.Y(),N.Y(),
 206                           U1.Z(),U2.Z(),N.Z());
 207   if(D) {
 208     dist = Det33(U1.X(),U2.X(),O1O2.X(),
 209                  U1.Y(),U2.Y(),O1O2.Y(),
 210                  U1.Z(),U2.Z(),O1O2.Z()) / D;
 211     Param1 = Det33(O1O2.X(),U2.X(),N.X(),
 212                    O1O2.Y(),U2.Y(),N.Y(),
 213                    O1O2.Z(),U2.Z(),N.Z()) / (-D);
 214     //------------------------------------------------------------
 215     //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
 216     //-- soit : Segment perpendiculaire : O1+P1 D1
 217     //--                                  O2-P2 D2
 218     Param2 = Det33(U1.X(),O1O2.X(),N.X(),
 219                    U1.Y(),O1O2.Y(),N.Y(),
 220                    U1.Z(),O1O2.Z(),N.Z()) / (D);
 221   }
 222 }
 223 //=======================================================================
 224 //function : DirToAx2
 225 //purpose  : returns a gp_Ax2 where D is the main direction
 226 //=======================================================================
 227 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
 228 {
 229   Standard_Real x=D.X(); Standard_Real ax=Abs(x);
 230   Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
 231   Standard_Real z=D.Z(); Standard_Real az=Abs(z);
 232   if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
 233     return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
 234   }
 235   else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
 236     return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
 237   }
 238   else {
 239     return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
 240   }
 241 }
 242 //=======================================================================
 243 //function : IntAna_QuadQuadGeo
 244 //purpose  : Empty constructor
 245 //=======================================================================
 246 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
 247     : done(Standard_False),
 248       nbint(0),
 249       typeres(IntAna_Empty),
 250       pt1(0,0,0),
 251       pt2(0,0,0),
 252       pt3(0,0,0),
 253       pt4(0,0,0),
 254       param1(0),
 255       param2(0),
 256       param3(0),
 257       param4(0),
 258       param1bis(0),
 259       param2bis(0),
 260       myCommonGen(Standard_False),
 261       myPChar(0,0,0)
 262 {
 263   InitTolerances();
 264 }
 265 //=======================================================================
 266 //function : InitTolerances
 267 //purpose  :
 268 //=======================================================================
 269 void IntAna_QuadQuadGeo::InitTolerances()
 270 {
 271   myEPSILON_DISTANCE               = 0.00000000000001;
 272   myEPSILON_ANGLE_CONE             = 0.000000000001;
 273   myEPSILON_MINI_CIRCLE_RADIUS     = 0.000000001;
 274   myEPSILON_CYLINDER_DELTA_RADIUS  = 0.0000000000001;
 275   myEPSILON_CYLINDER_DELTA_DISTANCE= 0.0000001;
 276   myEPSILON_AXES_PARA              = 0.000000000001;
 277 }
 278 //=======================================================================
 279 //function : IntAna_QuadQuadGeo
 280 //purpose  : Pln  Pln
 281 //=======================================================================
 282 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1,
 283                                        const gp_Pln& P2,
 284                                        const Standard_Real TolAng,
 285                                        const Standard_Real Tol)
 286 : done(Standard_False),
 287   nbint(0),
 288   typeres(IntAna_Empty),
 289   pt1(0,0,0),
 290   pt2(0,0,0),
 291   pt3(0,0,0),
 292   pt4(0,0,0),
 293   param1(0),
 294   param2(0),
 295   param3(0),
 296   param4(0),
 297   param1bis(0),
 298   param2bis(0),
 299   myCommonGen(Standard_False),
 300   myPChar(0,0,0)
 301 {
 302   InitTolerances();
 303   Perform(P1,P2,TolAng,Tol);
 304 }
 305 //=======================================================================
 306 //function : Perform
 307 //purpose  :
 308 //=======================================================================
 309 void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
 310                                   const gp_Pln& P2,
 311                                   const Standard_Real TolAng,
 312                                   const Standard_Real Tol)
 313 {
 314   Standard_Real A1, B1, C1, D1, A2, B2, C2, D2, dist1, dist2, aMVD;
 315   //
 316   done=Standard_False;
 317   param2bis=0.;
 318   //
 319   P1.Coefficients(A1,B1,C1,D1);
 320   P2.Coefficients(A2,B2,C2,D2);
 321   //
 322   gp_Vec aVN1(A1,B1,C1);
 323   gp_Vec aVN2(A2,B2,C2);
 324   gp_Vec vd(aVN1.Crossed(aVN2));
 325   //
 326   const gp_Pnt& aLocP1=P1.Location();
 327   const gp_Pnt& aLocP2=P2.Location();
 328   //
 329   dist1=A2*aLocP1.X() + B2*aLocP1.Y() + C2*aLocP1.Z() + D2;
 330   dist2=A1*aLocP2.X() + B1*aLocP2.Y() + C1*aLocP2.Z() + D1;
 331   //
 332   aMVD=vd.Magnitude();
 333   if(aMVD <=TolAng) {
 334     // normalles are collinear - planes are same or parallel
 335     typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same
 336       : IntAna_Empty;
 337   }
 338   else {
 339     Standard_Real denom, denom2, ddenom, par1, par2;
 340     Standard_Real X1, Y1, Z1, X2, Y2, Z2, aEps;
 341     //
 342     aEps=1.e-16;
 343     denom=A1*A2 + B1*B2 + C1*C2;
 344     denom2 = denom*denom;
 345     ddenom = 1. - denom2;
 346
 347     denom = ( Abs(ddenom) <= aEps ) ? aEps : ddenom;
 348
 349     par1 = dist1/denom;
 350     par2 = -dist2/denom;
 351
 352     gp_Vec inter1(aVN1.Crossed(vd));
 353     gp_Vec inter2(aVN2.Crossed(vd));
 354
 355     X1=aLocP1.X() + par1*inter1.X();
 356     Y1=aLocP1.Y() + par1*inter1.Y();
 357     Z1=aLocP1.Z() + par1*inter1.Z();
 358     X2=aLocP2.X() + par2*inter2.X();
 359     Y2=aLocP2.Y() + par2*inter2.Y();
 360     Z2=aLocP2.Z() + par2*inter2.Z();
 361
 362     pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
 363     dir1 = gp_Dir(vd);
 364     typeres = IntAna_Line;
 365     nbint = 1;
 366     //
 367     //-------------------------------------------------------
 368     // When the value of the angle between the planes is small
 369     // the origin of intersection line is computed with error
 370     // [ ~0.0001 ] that can not br considered as small one
 371     // e.g.
 372     // for {A~=2.e-6, dist1=4.2e-5, dist2==1.e-4} =>
 373     // {denom=3.4e-12, par1=12550297.6, par2=32605552.9, etc}
 374     // So,
 375     // the origin should be refined if it is possible
 376     //
 377     Standard_Real aTreshAng, aTreshDist;
 378     //
 379     aTreshAng=2.e-6; // 1.e-4 deg
 380     aTreshDist=1.e-12;
 381     //
 382     if (aMVD < aTreshAng) {
 383       Standard_Real aDist1, aDist2;
 384       //
 385       aDist1=A1*pt1.X() + B1*pt1.Y() + C1*pt1.Z() + D1;
 386       aDist2=A2*pt1.X() + B2*pt1.Y() + C2*pt1.Z() + D2;
 387       //
 388       if (fabs(aDist1)>aTreshDist || fabs(aDist2)>aTreshDist) {
 389         Standard_Boolean bIsDone, bIsParallel;
 390         IntAna_IntConicQuad aICQ;
 391         //
 392         // 1.
 393         gp_Dir aDN1(aVN1);
 394         gp_Lin aL1(pt1, aDN1);
 395         //
 396         aICQ.Perform(aL1, P1, TolAng, Tol);
 397         bIsDone=aICQ.IsDone();
 398         if (!bIsDone) {
 399           return;
 400         }
 401         //
 402         const gp_Pnt& aPnt1=aICQ.Point(1);
 403         //----------------------------------
 404         // 2.
 405         gp_Dir aDL2(dir1.Crossed(aDN1));
 406         gp_Lin aL2(aPnt1, aDL2);
 407         //
 408         aICQ.Perform(aL2, P2, TolAng, Tol);
 409         bIsDone=aICQ.IsDone();
 410         if (!bIsDone) {
 411           return;
 412         }
 413         //
 414         bIsParallel=aICQ.IsParallel();
 415         if (bIsParallel) {
 416           return;
 417         }
 418         //
 419         const gp_Pnt& aPnt2=aICQ.Point(1);
 420         //
 421         pt1=aPnt2;
 422       }
 423     }
 424   }
 425   done=Standard_True;
 426 }
 427 //=======================================================================
 428 //function : IntAna_QuadQuadGeo
 429 //purpose  : Pln Cylinder
 430 //=======================================================================
 431 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
 432                                        ,const gp_Cylinder& Cl
 433                                        ,const Standard_Real Tolang
 434                                        ,const Standard_Real Tol
 435                                        ,const Standard_Real H)
 436      : done(Standard_False),
 437        nbint(0),
 438        typeres(IntAna_Empty),
 439        pt1(0,0,0),
 440        pt2(0,0,0),
 441        pt3(0,0,0),
 442        pt4(0,0,0),
 443        param1(0),
 444        param2(0),
 445        param3(0),
 446        param4(0),
 447        param1bis(0),
 448        param2bis(0),
 449        myCommonGen(Standard_False),
 450        myPChar(0,0,0)
 451 {
 452   InitTolerances();
 453   Perform(P,Cl,Tolang,Tol,H);
 454 }
 455 //=======================================================================
 456 //function : Perform
 457 //purpose  :
 458 //=======================================================================
 459   void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
 460                                    ,const gp_Cylinder& Cl
 461                                    ,const Standard_Real Tolang
 462                                    ,const Standard_Real Tol
 463                                    ,const Standard_Real H)
 464 {
 465   done = Standard_False;
 466   Standard_Real dist,radius;
 467   Standard_Real A,B,C,D;
 468   Standard_Real X,Y,Z;
 469   Standard_Real sint,cost,h;
 470   gp_XYZ axex,axey,omega;
 471
 472
 473   param2bis=0.0;
 474   radius = Cl.Radius();
 475
 476   gp_Lin axec(Cl.Axis());
 477   gp_XYZ normp(P.Axis().Direction().XYZ());
 478
 479   P.Coefficients(A,B,C,D);
 480   axec.Location().Coord(X,Y,Z);
 481   // la distance axe/plan est evaluee a l origine de l axe.
 482   dist = A*X + B*Y + C*Z + D;
 483
 484   Standard_Real tolang = Tolang;
 485   Standard_Boolean newparams = Standard_False;
 486
 487   gp_Vec ldv( axec.Direction() );
 488   gp_Vec npv( normp );
 489   Standard_Real dA = Abs( ldv.Angle( npv ) );
 490   if( dA > (M_PI/4.) )
 491     {
 492       Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
 493       Standard_Real dangle = Abs( dang );
 494       if( dangle > Tolang )
 495         {
 496           Standard_Real sinda = Abs( Sin( dangle ) );
 497           Standard_Real dif = Abs( sinda - Tol );
 498           if( dif < Tol )
 499             {
 500               tolang = sinda * 2.;
 501               newparams = Standard_True;
 502             }
 503         }
 504     }
 505
 506   nbint = 0;
 507   IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
 508
 509   if (inter.IsParallel()) {
 510     // Le resultat de l intersection Plan-Cylindre est de type droite.
 511     // il y a 1 ou 2 droites
 512
 513     typeres = IntAna_Line;
 514     omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
 515
 516     if (Abs(Abs(dist)-radius) < Tol)
 517       {
 518         nbint = 1;
 519         pt1.SetXYZ(omega);
 520
 521         if( newparams )
 522           {
 523             gp_XYZ omegaXYZ(X,Y,Z);
 524             gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
 525             Standard_Real Xt,Yt,Zt,distt;
 526             omegaXYZtrnsl.Coord(Xt,Yt,Zt);
 527             distt = A*Xt + B*Yt + C*Zt + D;
 528             gp_XYZ omega1(omegaXYZtrnsl.X()-distt*A,
 529                           omegaXYZtrnsl.Y()-distt*B,
 530                           omegaXYZtrnsl.Z()-distt*C );
 531             gp_Pnt ppt1;
 532             ppt1.SetXYZ( omega1 );
 533             gp_Vec vv1(pt1,ppt1);
 534             gp_Dir dd1( vv1 );
 535             dir1 = dd1;
 536           }
 537         else
 538           dir1 = axec.Direction();
 539     }
 540     else if (Abs(dist) < radius)
 541       {
 542         nbint = 2;
 543         h = Sqrt(radius*radius - dist*dist);
 544         axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
 545
 546         pt1.SetXYZ(omega - h*axey);
 547         pt2.SetXYZ(omega + h*axey);
 548
 549         if( newparams )
 550           {
 551             gp_XYZ omegaXYZ(X,Y,Z);
 552             gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
 553             Standard_Real Xt,Yt,Zt,distt,ht;
 554             omegaXYZtrnsl.Coord(Xt,Yt,Zt);
 555             distt = A*Xt + B*Yt + C*Zt + D;
 556             //             ht = Sqrt(radius*radius - distt*distt);
 557             Standard_Real anSqrtArg = radius*radius - distt*distt;
 558             ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
 559
 560             gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A,
 561                           omegaXYZtrnsl.Y()-distt*B,
 562                           omegaXYZtrnsl.Z()-distt*C );
 563             gp_Pnt ppt1,ppt2;
 564             ppt1.SetXYZ( omega1 - ht*axey);
 565             ppt2.SetXYZ( omega1 + ht*axey);
 566             gp_Vec vv1(pt1,ppt1);
 567             gp_Vec vv2(pt2,ppt2);
 568             gp_Dir dd1( vv1 );
 569             gp_Dir dd2( vv2 );
 570             dir1 = dd1;
 571             dir2 = dd2;
 572           }
 573         else
 574           {
 575             dir1 = axec.Direction();
 576             dir2 = axec.Direction();
 577           }
 578     }
 579     //  else nbint = 0
 580
 581     // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
 582     // et ne pas etre seulement supprime...
 583
 584     else {
 585       typeres = IntAna_Empty;
 586     }
 587   }
 588   else {     // Il y a un point d intersection. C est le centre du cercle
 589              // ou de l ellipse solution.
 590
 591     nbint = 1;
 592     axey = normp.Crossed(axec.Direction().XYZ());
 593     sint = axey.Modulus();
 594
 595     pt1 = inter.Point(1);
 596
 597     if (sint < Tol/radius) {
 598
 599       // on construit un cercle avec comme axes X et Y ceux du cylindre
 600       typeres = IntAna_Circle;
 601
 602       dir1 = axec.Direction(); // axe Z
 603       dir2 = Cl.Position().XDirection();
 604       param1 = radius;
 605     }
 606     else {
 607
 608       // on construit un ellipse
 609       typeres = IntAna_Ellipse;
 610       cost = Abs(axec.Direction().XYZ().Dot(normp));
 611       axex = axey.Crossed(normp);
 612
 613       dir1.SetXYZ(normp);         //Modif ds ce bloc
 614       dir2.SetXYZ(axex);
 615
 616       param1    = radius/cost;
 617       param1bis = radius;
 618     }
 619   }
 620
 621   done = Standard_True;
 622 }
 623 //=======================================================================
 624 //function : IntAna_QuadQuadGeo
 625 //purpose  : Pln Cone
 626 //=======================================================================
 627   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
 628                                          const gp_Cone& Co,
 629                                          const Standard_Real Tolang,
 630                                          const Standard_Real Tol)
 631 : done(Standard_False),
 632   nbint(0),
 633   typeres(IntAna_Empty),
 634   pt1(0,0,0),
 635   pt2(0,0,0),
 636   pt3(0,0,0),
 637   pt4(0,0,0),
 638   param1(0),
 639   param2(0),
 640   param3(0),
 641   param4(0),
 642   param1bis(0),
 643   param2bis(0),
 644   myCommonGen(Standard_False),
 645   myPChar(0,0,0)
 646 {
 647   InitTolerances();
 648   Perform(P,Co,Tolang,Tol);
 649 }
 650 //=======================================================================
 651 //function : Perform
 652 //purpose  :
 653 //=======================================================================
 654   void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
 655                                    const gp_Cone& Co,
 656                                    const Standard_Real Tolang,
 657                                    const Standard_Real Tol)
 658 {
 659
 660   done = Standard_False;
 661   nbint = 0;
 662
 663   Standard_Real A,B,C,D;
 664   Standard_Real X,Y,Z;
 665   Standard_Real dist,sint,cost,sina,cosa,angl,costa;
 666   Standard_Real dh;
 667   gp_XYZ axex,axey;
 668
 669   gp_Lin axec(Co.Axis());
 670   P.Coefficients(A,B,C,D);
 671   gp_Pnt apex(Co.Apex());
 672
 673   apex.Coord(X,Y,Z);
 674   dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
 675
 676   gp_XYZ normp = P.Axis().Direction().XYZ();
 677   if(P.Direct()==Standard_False) {  //-- lbr le 14 jan 97
 678     normp.Reverse();
 679   }
 680
 681   axey = normp.Crossed(Co.Axis().Direction().XYZ());
 682   axex = axey.Crossed(normp);
 683
 684
 685   angl = Co.SemiAngle();
 686
 687   cosa = Cos(angl);
 688   sina = Abs(Sin(angl));
 689
 690
 691   // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
 692
 693   sint = axey.Modulus();
 694   cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
 695
 696   // Le calcul de costa permet de determiner si le plan contient
 697   // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
 698
 699   costa = cost*cosa - sint*sina;  // sin((PI/2 -t)-angl)=cos(t+angl)
 700                                   // avec  t ramene entre 0 et pi/2.
 701
 702   if (Abs(dist) < Tol) {
 703     // on considere que le plan contient le sommet du cone.
 704     // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
 705     // selon l inclinaison du plan.
 706
 707     if (Abs(costa) < Tolang) {          // plan parallele a la generatrice
 708       typeres = IntAna_Line;
 709       nbint = 1;
 710       gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
 711       // point sur l axe du cone cote z positif
 712
 713       dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
 714       ptonaxe = ptonaxe - dist*normp;
 715       pt1 = apex;
 716       dir1.SetXYZ(ptonaxe - pt1.XYZ());
 717     }
 718     else if (cost < sina) {   // plan "interieur" au cone
 719       typeres = IntAna_Line;
 720       nbint = 2;
 721       pt1 = apex;
 722       pt2 = apex;
 723       dh = Sqrt(sina*sina-cost*cost)/cosa;
 724       dir1.SetXYZ(axex + dh*axey);
 725       dir2.SetXYZ(axex - dh*axey);
 726     }
 727     else {                              // plan "exterieur" au cone
 728       typeres = IntAna_Point;
 729       nbint = 1;
 730       pt1 = apex;
 731     }
 732   }
 733   else {
 734     // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
 735     // l inclinaison du plan.
 736     Standard_Real deltacenter, distance;
 737
 738     if (cost < Tolang) {
 739       // Le plan contient la direction de l axe du cone. La solution est
 740       // l hyperbole
 741       typeres = IntAna_Hyperbola;
 742       nbint = 2;
 743       pt1.SetXYZ(apex.XYZ()-dist*normp);
 744       pt2 = pt1;
 745       dir1=normp;
 746       dir2.SetXYZ(axex);
 747       param1     = param2 = Abs(dist/Tan(angl));
 748       param1bis  = param2bis = Abs(dist);
 749     }
 750     else {
 751
 752       IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
 753
 754       gp_Pnt center(inter.Point(1));
 755
 756       // En fonction de la position de l intersection par rapport au sommet
 757       // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
 758       // sur axec est negatif (voir definition du cone)
 759
 760       distance = apex.Distance(center);
 761
 762       if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
 763         axex.Reverse();
 764         axey.Reverse();
 765       }
 766
 767       if (Abs(costa) < Tolang) {  // plan parallele a une generatrice
 768         typeres = IntAna_Parabola;
 769         nbint = 1;
 770         deltacenter = distance/2./cosa;
 771         axex.Normalize();
 772         pt1.SetXYZ(center.XYZ()-deltacenter*axex);
 773         dir1 = normp;
 774         dir2.SetXYZ(axex);
 775         param1 = deltacenter*sina*sina;
 776       }
 777       else if (sint  < Tolang) {            // plan perpendiculaire a l axe
 778         typeres = IntAna_Circle;
 779         nbint = 1;
 780         pt1 = center;
 781         dir1 = Co.Position().Direction();
 782         dir2 = Co.Position().XDirection();
 783         param1 = apex.Distance(center)*Abs(Tan(angl));
 784       }
 785       else if (cost < sina ) {
 786         typeres = IntAna_Hyperbola;
 787         nbint = 2;
 788         axex.Normalize();
 789
 790         deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
 791         pt1.SetXYZ(center.XYZ() - deltacenter*axex);
 792         pt2 = pt1;
 793         dir1 = normp;
 794         dir2.SetXYZ(axex);
 795         param1    = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
 796         param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
 797
 798       }
 799       else {                    // on a alors cost > sina
 800         typeres = IntAna_Ellipse;
 801         nbint = 1;
 802         Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
 803         deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
 804         axex.Normalize();
 805         pt1.SetXYZ(center.XYZ() + deltacenter*axex);
 806         dir1 = normp;
 807         dir2.SetXYZ(axex);
 808         param1    = radius;
 809         param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
 810       }
 811     }
 812   }
 813
 814   //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
 815   //-- des hyperboles trop bizarres
 816   //-- On retourne False -> Traitement par biparametree
 817   static Standard_Real EllipseLimit   = 1.0E+9; //OCC513(apo) 1000000
 818   static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
 819   if(typeres==IntAna_Ellipse && nbint>=1) {
 820     if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit)  {
 821       done=Standard_False;
 822       return;
 823     }
 824   }
 825   if(typeres==IntAna_Hyperbola && nbint>=2) {
 826     if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit)  {
 827       done = Standard_False;
 828       return;
 829     }
 830   }
 831   if(typeres==IntAna_Hyperbola && nbint>=1) {
 832     if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit)  {
 833       done=Standard_False;
 834       return;
 835     }
 836   }
 837
 838   done = Standard_True;
 839 }
 840
 841 //=======================================================================
 842 //function : IntAna_QuadQuadGeo
 843 //purpose  : Pln Sphere
 844 //=======================================================================
 845 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
 846                                        const gp_Sphere& S)
 847 : done(Standard_False),
 848   nbint(0),
 849   typeres(IntAna_Empty),
 850   pt1(0,0,0),
 851   pt2(0,0,0),
 852   pt3(0,0,0),
 853   pt4(0,0,0),
 854   param1(0),
 855   param2(0),
 856   param3(0),
 857   param4(0),
 858   param1bis(0),
 859   param2bis(0),
 860   myCommonGen(Standard_False),
 861   myPChar(0,0,0)
 862 {
 863   InitTolerances();
 864   Perform(P,S);
 865 }
 866 //=======================================================================
 867 //function : Perform
 868 //purpose  :
 869 //=======================================================================
 870 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
 871                                  ,const gp_Sphere& S)
 872 {
 873
 874   done = Standard_False;
 875   Standard_Real A,B,C,D,dist, radius;
 876   Standard_Real X,Y,Z;
 877
 878   nbint = 0;
 879 // debug JAG : on met typeres = IntAna_Empty par defaut...
 880   typeres = IntAna_Empty;
 881
 882   P.Coefficients(A,B,C,D);
 883   S.Location().Coord(X,Y,Z);
 884   radius = S.Radius();
 885
 886   dist = A * X + B * Y + C * Z + D;
 887
 888   if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
 889     // on a une seule solution : le point projection du centre de la sphere
 890     // sur le plan
 891     nbint = 1;
 892     typeres = IntAna_Point;
 893     pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
 894   }
 895   else if (Abs(dist) < radius) {
 896     // on a un cercle solution
 897     nbint = 1;
 898     typeres = IntAna_Circle;
 899     pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
 900     dir1 = P.Axis().Direction();
 901     if(P.Direct()==Standard_False) dir1.Reverse();
 902     dir2 = P.Position().XDirection();
 903     param1 = Sqrt(radius*radius - dist*dist);
 904   }
 905   param2bis=0.0; //-- pour eviter param2bis not used ....
 906   done = Standard_True;
 907 }
 908
 909 //=======================================================================
 910 //function : IntAna_QuadQuadGeo
 911 //purpose  : Cylinder - Cylinder
 912 //=======================================================================
 913 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
 914                                        const gp_Cylinder& Cyl2,
 915                                        const Standard_Real Tol)
 916 : done(Standard_False),
 917   nbint(0),
 918   typeres(IntAna_Empty),
 919   pt1(0,0,0),
 920   pt2(0,0,0),
 921   pt3(0,0,0),
 922   pt4(0,0,0),
 923   param1(0),
 924   param2(0),
 925   param3(0),
 926   param4(0),
 927   param1bis(0),
 928   param2bis(0),
 929   myCommonGen(Standard_False),
 930   myPChar(0,0,0)
 931 {
 932   InitTolerances();
 933   Perform(Cyl1,Cyl2,Tol);
 934 }
 935 //=======================================================================
 936 //function : Perform
 937 //purpose  :
 938 //=======================================================================
 939 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
 940                                  const gp_Cylinder& Cyl2,
 941                                  const Standard_Real Tol)
 942 {
 943   done=Standard_True;
 944   //---------------------------- Parallel axes -------------------------
 945   AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
 946   Standard_Real R1=Cyl1.Radius();
 947   Standard_Real R2=Cyl2.Radius();
 948   Standard_Real RmR, RmR_Relative;
 949   RmR=(R1>R2)? (R1-R2) : (R2-R1);
 950   {
 951     Standard_Real Rmax;
 952     Rmax=(R1>R2)? R1 : R2;
 953     RmR_Relative=RmR/Rmax;
 954   }
 955
 956   Standard_Real DistA1A2=A1A2.Distance();
 957
 958   if(A1A2.Parallel())
 959   {
 960     if(DistA1A2<=Tol)
 961     {
 962       if(RmR<=Tol)
 963       {
 964         typeres=IntAna_Same;
 965       }
 966       else
 967       {
 968         typeres=IntAna_Empty;
 969       }
 970     }
 971     else
 972     {  //-- DistA1A2 > Tol
 973       gp_Pnt P1=Cyl1.Location();
 974       gp_Pnt P2t=Cyl2.Location();
 975       gp_Pnt P2;
 976       //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
 977       gp_Dir DirCyl = Cyl1.Position().Direction();
 978       Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
 979
 980       //P2 is a projection the location of the 2nd cylinder on the base
 981       //of the 1st cylinder
 982       P2.SetCoord(P2t.X() - ProjP2OnDirCyl1*DirCyl.X(),
 983                   P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y(),
 984                   P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
 985       //--
 986       Standard_Real R1pR2=R1+R2;
 987       if(DistA1A2>(R1pR2+Tol))
 988       {
 989         typeres=IntAna_Empty;
 990         nbint=0;
 991       }
 992       else if((R1pR2 - DistA1A2) <= RealSmall())
 993       {
 994         //-- 1 Tangent line -------------------------------------OK
 995         typeres=IntAna_Line;
 996
 997         nbint=1;
 998         dir1=DirCyl;
 999         Standard_Real R1_R1pR2=R1/R1pR2;
1000         pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X()),
1001                       P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y()),
1002                       P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
1003       }
1004       else if(DistA1A2>RmR)
1005       {
1006         //-- 2 lines ---------------------------------------------OK
1007         typeres=IntAna_Line;
1008         nbint=2;
1009         dir1=DirCyl;
1010         dir2=dir1;
1011
1012         const Standard_Real aR1R1 = R1*R1;
1013
1014         /*
1015                       P1
1016                       o
1017                     * | *
1018                   * O1|   *
1019               A o-----o-----o B
1020                   *   |   *
1021                     * | *
1022                       o
1023                       P2
1024
1025           Two cylinders have axes collinear. Therefore, problem can be reformulated as
1026           to find intersection point of two circles (the bases of the cylinders) on
1027           the plane: 1st circle has center P1 and radius R1 (the radius of the
1028           1st cylinder) and 2nd circle has center P2 and radius R2 (the radius of the
1029           2nd cylinder). The plane is the base of the 1st cylinder. Points A and B
1030           are intersection point of these circles. Distance P1P2 is equal to DistA1A2.
1031           O1 is the intersection point of P1P2 and AB segments.
1032
1033           At that, if distance AB < Tol we consider that the circles are tangent and
1034           has only one intersection point.
1035
1036             AB = 2*R1*sin(angle AP1P2).
1037           Accordingly,
1038             AB^2 < Tol^2 => 4*R1*R1*sin(angle AP1P2)^2 < Tol*Tol.
1039         */
1040
1041
1042         //Cosine and Square of Sine of the A-P1-P2 angle
1043         const Standard_Real aCos = 0.5*(aR1R1-R2*R2+DistA1A2*DistA1A2)/(R1*DistA1A2);
1044         const Standard_Real aSin2 = 1-aCos*aCos;
1045
1046         const Standard_Boolean isTangent =((4.0*aR1R1*aSin2) < Tol*Tol);
1047
1048         //Normalized vector P1P2
1049         const gp_Vec DirA1A2((P2.XYZ() - P1.XYZ())/DistA1A2);
1050
1051         if(isTangent)
1052         {
1053           //Intercept the segment from P1 point along P1P2 direction
1054           //and having |P1O1| length
1055           nbint=1;
1056           pt1.SetXYZ(P1.XYZ() + DirA1A2.XYZ()*R1*aCos);
1057         }
1058         else
1059         {
1060           //Sine of the A-P1-P2 angle (if aSin2 < 0 then isTangent == TRUE =>
1061           //go to another branch)
1062           const Standard_Real aSin = sqrt(aSin2);
1063
1064           //1. Rotate P1P2 to the angle A-P1-P2 relative to P1
1065           //(clockwise and anticlockwise for getting
1066           //two intersection points).
1067           //2. Intercept the segment from P1 along direction,
1068           //determined in the preview paragraph and having R1 length
1069           const gp_Dir  &aXDir = Cyl1.Position().XDirection(),
1070                         &aYDir = Cyl1.Position().YDirection();
1071           const gp_Vec  aR1Xdir = R1*aXDir.XYZ(),
1072                         aR1Ydir = R1*aYDir.XYZ();
1073
1074           //Source 2D-coordinates of the P1P2 vector normalized
1075           //in coordinate system, based on the X- and Y-directions
1076           //of the 1st cylinder in the plane of the 1st cylinder base
1077           //(P1 is the origin of the coordinate system).
1078           const Standard_Real aDx = DirA1A2.Dot(aXDir),
1079                               aDy = DirA1A2.Dot(aYDir);
1080
1081           //New coordinate (after rotation) of the P1P2 vector normalized.
1082           Standard_Real aNewDx = aDx*aCos - aDy*aSin,
1083                         aNewDy = aDy*aCos + aDx*aSin;
1084           pt1.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
1085
1086           aNewDx = aDx*aCos + aDy*aSin;
1087           aNewDy = aDy*aCos - aDx*aSin;
1088           pt2.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
1089         }
1090       }
1091       else if(DistA1A2>(RmR-Tol))
1092       {
1093         //-- 1 Tangent ------------------------------------------OK
1094         typeres=IntAna_Line;
1095         nbint=1;
1096         dir1=DirCyl;
1097         Standard_Real R1_RmR=R1/RmR;
1098
1099         if(R1 < R2)
1100           R1_RmR = -R1_RmR;
1101
1102         pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X()),
1103                       P1.Y() + R1_RmR * (P2.Y()-P1.Y()),
1104                       P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
1105       }
1106       else {
1107         nbint=0;
1108         typeres=IntAna_Empty;
1109       }
1110     }
1111   }
1112   else { //-- No Parallel Axis ---------------------------------OK
1113     if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
1114        && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE))
1115     {
1116       //-- PI/2 between the two axis   and   Intersection
1117       //-- and identical radius
1118       typeres=IntAna_Ellipse;
1119       nbint=2;
1120       gp_Dir DirCyl1=Cyl1.Position().Direction();
1121       gp_Dir DirCyl2=Cyl2.Position().Direction();
1122       pt1=pt2=A1A2.PtIntersect();
1123
1124       Standard_Real A=DirCyl1.Angle(DirCyl2);
1125       Standard_Real B;
1126       B=Abs(Sin(0.5*(M_PI-A)));
1127       A=Abs(Sin(0.5*A));
1128
1129       if(A==0.0 || B==0.0)
1130       {
1131         typeres=IntAna_Same;
1132         return;
1133       }
1134
1135       gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
1136       dir1 = gp_Dir(dircyl1.Added(dircyl2));
1137       dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
1138
1139       param2   = Cyl1.Radius() / A;
1140       param1   = Cyl1.Radius() / B;
1141       param2bis= param1bis = Cyl1.Radius();
1142       if(param1 < param1bis)
1143       {
1144         A=param1;
1145         param1=param1bis;
1146         param1bis=A;
1147       }
1148
1149       if(param2 < param2bis)
1150       {
1151         A=param2;
1152         param2=param2bis;
1153         param2bis=A;
1154       }
1155     }
1156     else
1157     {
1158       if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol)
1159       {
1160         typeres = IntAna_Point;
1161         Standard_Real d,p1,p2;
1162
1163         gp_Dir D1 = Cyl1.Axis().Direction();
1164         gp_Dir D2 = Cyl2.Axis().Direction();
1165         A1A2.Distance(d,p1,p2);
1166         gp_Pnt P = Cyl1.Axis().Location();
1167         gp_Pnt P1(P.X() - p1*D1.X(),
1168                   P.Y() - p1*D1.Y(),
1169                   P.Z() - p1*D1.Z());
1170
1171         P = Cyl2.Axis().Location();
1172
1173         gp_Pnt P2(P.X() - p2*D2.X(),
1174                   P.Y() - p2*D2.Y(),
1175                   P.Z() - p2*D2.Z());
1176
1177         gp_Vec P1P2(P1,P2);
1178         D1=gp_Dir(P1P2);
1179         p1=Cyl1.Radius();
1180
1181         pt1.SetCoord(P1.X() + p1*D1.X(),
1182                      P1.Y() + p1*D1.Y(),
1183                      P1.Z() + p1*D1.Z());
1184         nbint = 1;
1185       }
1186       else
1187       {
1188         typeres=IntAna_NoGeometricSolution;
1189       }
1190     }
1191   }
1192 }
1193 //=======================================================================
1194 //function : IntAna_QuadQuadGeo
1195 //purpose  : Cylinder - Cone
1196 //=======================================================================
1197 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1198                                        const gp_Cone& Con,
1199                                        const Standard_Real Tol)
1200 : done(Standard_False),
1201   nbint(0),
1202   typeres(IntAna_Empty),
1203   pt1(0,0,0),
1204   pt2(0,0,0),
1205   pt3(0,0,0),
1206   pt4(0,0,0),
1207   param1(0),
1208   param2(0),
1209   param3(0),
1210   param4(0),
1211   param1bis(0),
1212   param2bis(0),
1213   myCommonGen(Standard_False),
1214   myPChar(0,0,0)
1215 {
1216   InitTolerances();
1217   Perform(Cyl,Con,Tol);
1218 }
1219 //=======================================================================
1220 //function : Perform
1221 //purpose  :
1222 //=======================================================================
1223   void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1224                                    const gp_Cone& Con,
1225                                    const Standard_Real )
1226 {
1227   done=Standard_True;
1228   AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1229   if(A1A2.Same()) {
1230     gp_Pnt Pt=Con.Apex();
1231     Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1232     gp_Dir dir=Cyl.Position().Direction();
1233     pt1.SetCoord( Pt.X() + dist*dir.X()
1234                  ,Pt.Y() + dist*dir.Y()
1235                  ,Pt.Z() + dist*dir.Z());
1236     pt2.SetCoord( Pt.X() - dist*dir.X()
1237                  ,Pt.Y() - dist*dir.Y()
1238                  ,Pt.Z() - dist*dir.Z());
1239     dir1=dir2=dir;
1240     param1=param2=Cyl.Radius();
1241     nbint=2;
1242     typeres=IntAna_Circle;
1243
1244   }
1245   else {
1246     typeres=IntAna_NoGeometricSolution;
1247   }
1248 }
1249 //=======================================================================
1250 //function :
1251 //purpose  : Cylinder - Sphere
1252 //=======================================================================
1253   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1254                                          const gp_Sphere& Sph,
1255                                          const Standard_Real Tol)
1256 : done(Standard_False),
1257   nbint(0),
1258   typeres(IntAna_Empty),
1259   pt1(0,0,0),
1260   pt2(0,0,0),
1261   pt3(0,0,0),
1262   pt4(0,0,0),
1263   param1(0),
1264   param2(0),
1265   param3(0),
1266   param4(0),
1267   param1bis(0),
1268   param2bis(0),
1269   myCommonGen(Standard_False),
1270   myPChar(0,0,0)
1271 {
1272   InitTolerances();
1273   Perform(Cyl,Sph,Tol);
1274 }
1275 //=======================================================================
1276 //function : Perform
1277 //purpose  :
1278 //=======================================================================
1279   void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1280                                    ,const gp_Sphere& Sph
1281                                    ,const Standard_Real)
1282 {
1283   done=Standard_True;
1284   gp_Pnt Pt=Sph.Location();
1285   AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1286   if((A1A2.Intersect()  && Pt.Distance(A1A2.PtIntersect())==0.0 )
1287      || (A1A2.Same()))      {
1288     if(Sph.Radius() < Cyl.Radius()) {
1289       typeres = IntAna_Empty;
1290     }
1291     else {
1292       Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1293       gp_Dir dir=Cyl.Position().Direction();
1294       dir1 = dir2 = dir;
1295       typeres=IntAna_Circle;
1296       pt1.SetCoord( Pt.X() + dist*dir.X()
1297                    ,Pt.Y() + dist*dir.Y()
1298                    ,Pt.Z() + dist*dir.Z());
1299       nbint=1;
1300       param1 = Cyl.Radius();
1301       if(dist>RealEpsilon()) {
1302         pt2.SetCoord( Pt.X() - dist*dir.X()
1303                      ,Pt.Y() - dist*dir.Y()
1304                      ,Pt.Z() - dist*dir.Z());
1305         param2=Cyl.Radius();
1306         nbint=2;
1307       }
1308     }
1309   }
1310   else {
1311     typeres=IntAna_NoGeometricSolution;
1312   }
1313 }
1314
1315 //=======================================================================
1316 //function : IntAna_QuadQuadGeo
1317 //purpose  : Cone - Cone
1318 //=======================================================================
1319   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1320                                          const gp_Cone& Con2,
1321                                          const Standard_Real Tol)
1322 : done(Standard_False),
1323   nbint(0),
1324   typeres(IntAna_Empty),
1325   pt1(0,0,0),
1326   pt2(0,0,0),
1327   pt3(0,0,0),
1328   pt4(0,0,0),
1329   param1(0),
1330   param2(0),
1331   param3(0),
1332   param4(0),
1333   param1bis(0),
1334   param2bis(0),
1335   myCommonGen(Standard_False),
1336   myPChar(0,0,0)
1337 {
1338   InitTolerances();
1339   Perform(Con1,Con2,Tol);
1340 }
1341 //
1342 //=======================================================================
1343 //function : Perform
1344 //purpose  :
1345 //=======================================================================
1346   void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1347                                    const gp_Cone& Con2,
1348                                    const Standard_Real Tol)
1349 {
1350   done=Standard_True;
1351   //
1352   Standard_Real tg1, tg2, aDA1A2, aTol2;
1353   gp_Pnt aPApex1, aPApex2;
1354
1355   Standard_Real TOL_APEX_CONF = 1.e-10;
1356
1357   //
1358   tg1=Tan(Con1.SemiAngle());
1359   tg2=Tan(Con2.SemiAngle());
1360
1361   if((tg1 * tg2) < 0.) {
1362     tg2 = -tg2;
1363   }
1364   //
1365   aTol2=Tol*Tol;
1366   aPApex1=Con1.Apex();
1367   aPApex2=Con2.Apex();
1368   aDA1A2=aPApex1.SquareDistance(aPApex2);
1369   //
1370   AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1371   //
1372   // 1
1373   if(A1A2.Same()) {
1374     //-- two circles
1375     Standard_Real x;
1376     gp_Pnt P=Con1.Apex();
1377     gp_Dir D=Con1.Position().Direction();
1378     Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1379
1380     if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
1381       if (fabs(d) < TOL_APEX_CONF) {
1382         typeres = IntAna_Point;
1383         nbint = 1;
1384         pt1 = P;
1385         return;
1386       }
1387       x=(d*tg2)/(tg1+tg2);
1388       pt1.SetCoord( P.X() + x*D.X()
1389                    ,P.Y() + x*D.Y()
1390                    ,P.Z() + x*D.Z());
1391       param1=Abs(x*tg1);
1392
1393       x=(d*tg2)/(tg2-tg1);
1394       pt2.SetCoord( P.X() + x*D.X()
1395                    ,P.Y() + x*D.Y()
1396                    ,P.Z() + x*D.Z());
1397       param2=Abs(x*tg1);
1398       dir1 = dir2 = D;
1399       nbint=2;
1400       typeres=IntAna_Circle;
1401     }
1402     else {
1403       if (fabs(d) < TOL_APEX_CONF) {
1404         typeres=IntAna_Same;
1405       }
1406       else {
1407         typeres=IntAna_Circle;
1408         nbint=1;
1409         x=d*0.5;
1410         pt1.SetCoord( P.X() + x*D.X()
1411                      ,P.Y() + x*D.Y()
1412                      ,P.Z() + x*D.Z());
1413         param1 = Abs(x * tg1);
1414         dir1 = D;
1415       }
1416     }
1417   } //-- fin A1A2.Same
1418   // 2
1419   else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1420     //-- voir AnVer12mai98
1421     Standard_Real DistA1A2=A1A2.Distance();
1422     gp_Dir DA1=Con1.Position().Direction();
1423     gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1424     gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1425     Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1426
1427     gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
1428                       O1O2n.Y()-O1O2_DA1*DA1.Y(),
1429                       O1O2n.Z()-O1O2_DA1*DA1.Z());
1430     gp_Dir DB1=gp_Dir(O1_Proj_A2);
1431
1432     Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1433     Standard_Real ABSTG1 = Abs(tg1);
1434     Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1435     Standard_Real X1 = X2+yO1O2;
1436
1437     gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1438               Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1439               Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1440
1441     gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1442                  0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1443                  0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1444     gp_Vec P1MO1O2(P1,MO1O2);
1445
1446     gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1447     gp_Dir OrthoPln =  DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1448
1449     IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1450       if(INTER_QUAD_PLN.IsDone()) {
1451       switch(INTER_QUAD_PLN.TypeInter()) {
1452       case IntAna_Ellipse:         {
1453         typeres=IntAna_Ellipse;
1454         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1455         pt1 = E.Location();
1456         dir1 = E.Position().Direction();
1457         dir2 = E.Position().XDirection();
1458         param1 = E.MajorRadius();
1459         param1bis = E.MinorRadius();
1460         nbint = 1;
1461         break;
1462       }
1463       case IntAna_Circle: {
1464         typeres=IntAna_Circle;
1465         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1466         pt1 = C.Location();
1467         dir1 = C.Position().XDirection();
1468         dir2 = C.Position().YDirection();
1469         param1 = C.Radius();
1470         nbint = 1;
1471         break;
1472       }
1473       case IntAna_Hyperbola: {
1474         typeres=IntAna_Hyperbola;
1475         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1476         pt1 = pt2 = H.Location();
1477         dir1 = H.Position().Direction();
1478         dir2 = H.Position().XDirection();
1479         param1 = param2 = H.MajorRadius();
1480         param1bis = param2bis = H.MinorRadius();
1481         nbint = 2;
1482         break;
1483       }
1484       case IntAna_Line: {
1485         typeres=IntAna_Line;
1486         gp_Lin H=INTER_QUAD_PLN.Line(1);
1487         pt1 = pt2 = H.Location();
1488         dir1 = dir2 = H.Position().Direction();
1489         param1 = param2 = 0.0;
1490         param1bis = param2bis = 0.0;
1491         nbint = 2;
1492         break;
1493       }
1494       default:
1495         typeres=IntAna_NoGeometricSolution;
1496       }
1497     }
1498   }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1499   // 3
1500   else if (aDA1A2<aTol2) {
1501     //
1502     // When apices are coinsided there can be 3 possible cases
1503     // 3.1 - empty solution (iRet=0)
1504     // 3.2 - one line  when cone1 touches cone2 (iRet=1)
1505     // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1506     //
1507     Standard_Integer iRet;
1508     Standard_Real aGamma, aBeta1, aBeta2;
1509     Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1510     Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1511     gp_Pnt2d aP0, aPA1, aP1, aPA2;
1512     gp_Vec2d aVAx2;
1513     gp_Ax1 aAx1, aAx2;
1514     //
1515     // Preliminary analysis. Determination of iRet
1516     //
1517     iRet=0;
1518     aHalfPI=0.5*M_PI;
1519     aD1=1.;
1520     aPA1.SetCoord(aD1, 0.);
1521     aP0.SetCoord(0., 0.);
1522     //
1523     aAx1=Con1.Axis();
1524     aAx2=Con2.Axis();
1525     aGamma=aAx1.Angle(aAx2);
1526     if (aGamma>aHalfPI){
1527       aGamma=M_PI-aGamma;
1528     }
1529     aCosGamma=Cos(aGamma);
1530     aSinGamma=Sin(aGamma);
1531     //
1532     aBeta1=Con1.SemiAngle();
1533     aTgBeta1=Tan(aBeta1);
1534     aTgBeta1=Abs(aTgBeta1);
1535     //
1536     aBeta2=Con2.SemiAngle();
1537     aTgBeta2=Tan(aBeta2);
1538     aTgBeta2=Abs(aTgBeta2);
1539     //
1540     aR1=aD1*aTgBeta1;
1541     aP1.SetCoord(aD1, aR1);
1542     //
1543     // PA2
1544     aVAx2.SetCoord(aCosGamma, aSinGamma);
1545     gp_Dir2d aDAx2(aVAx2);
1546     gp_Lin2d aLAx2(aP0, aDAx2);
1547     //
1548     gp_Vec2d aV(aP0, aP1);
1549     aDx=aVAx2.Dot(aV);
1550     aPA2=aP0.Translated(aDx*aDAx2);
1551     //
1552     // aR2
1553     aDx=aPA2.Distance(aP0);
1554     aR2=aDx*aTgBeta2;
1555     //
1556     // aRD2
1557     aRD2=aPA2.Distance(aP1);
1558     //
1559     if (aRD2>(aR2+Tol)) {
1560       iRet=0;
1561       typeres=IntAna_Empty; //nothing
1562       return;
1563     }
1564     //
1565     iRet=1; //touch case => 1 line
1566     if (aRD2<(aR2-Tol)) {
1567       iRet=2;//intersection => couple of lines
1568     }
1569     //
1570     // Finding the solution in 3D
1571     //
1572     Standard_Real aDa;
1573     gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1574     gp_Dir aD3Ax1, aD3Ax2;
1575     gp_Lin aLin;
1576     IntAna_QuadQuadGeo aIntr;
1577     //
1578     aQApex1=Con1.Apex();
1579     aD3Ax1=aAx1.Direction();
1580     aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1581                   aQApex1.Y()+aD1*aD3Ax1.Y(),
1582                   aQApex1.Z()+aD1*aD3Ax1.Z());
1583     //
1584     aDx=aD3Ax1.Dot(aAx2.Direction());
1585     if (aDx<0.) {
1586       aAx2.Reverse();
1587     }
1588     aD3Ax2=aAx2.Direction();
1589     //
1590     aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1591     //
1592     aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1593                   aQApex1.Y()+aD2*aD3Ax2.Y(),
1594                   aQApex1.Z()+aD2*aD3Ax2.Z());
1595     //
1596     gp_Pln aPln1(aQA1, aD3Ax1);
1597     gp_Pln aPln2(aQA2, aD3Ax2);
1598     //
1599     aIntr.Perform(aPln1, aPln2, Tol, Tol);
1600     if (!aIntr.IsDone() || 0 == aIntr.NbSolutions()) {
1601       iRet=-1; // just in case. it must not be so
1602       typeres=IntAna_NoGeometricSolution;
1603       return;
1604     }
1605     //
1606     aLin=aIntr.Line(1);
1607     const gp_Dir& aDLin=aLin.Direction();
1608     gp_Vec aVLin(aDLin);
1609     gp_Pnt aOrig=aLin.Location();
1610     gp_Vec aVr(aQA1, aOrig);
1611     aDx=aVLin.Dot(aVr);
1612     aQX=aOrig.Translated(aDx*aVLin);
1613     //
1614     // Final part
1615     //
1616     typeres=IntAna_Line;
1617     //
1618     param1=0.;
1619     param2 =0.;
1620     param1bis=0.;
1621     param2bis=0.;
1622     //
1623     if (iRet==1) {
1624       // one line
1625       nbint=1;
1626       pt1=aQApex1;
1627       gp_Vec aVX(aQApex1, aQX);
1628       dir1=gp_Dir(aVX);
1629     }
1630
1631     else {//iRet=2
1632       // two lines
1633       nbint=2;
1634       aDa=aQA1.Distance(aQX);
1635       aDx=sqrt(aR1*aR1-aDa*aDa);
1636       aQX1=aQX.Translated(aDx*aVLin);
1637       aQX2=aQX.Translated(-aDx*aVLin);
1638       //
1639       pt1=aQApex1;
1640       pt2=aQApex1;
1641       gp_Vec aVX1(aQApex1, aQX1);
1642       dir1=gp_Dir(aVX1);
1643       gp_Vec aVX2(aQApex1, aQX2);
1644       dir2=gp_Dir(aVX2);
1645     }
1646   } //else if (aDA1A2<aTol2) {
1647   //Case when cones have common generatrix
1648   else if(A1A2.Intersect()) {
1649     //Check if apex of one cone belongs another one
1650     Standard_Real u, v, tol2 = Tol*Tol;
1651     ElSLib::Parameters(Con2, aPApex1, u, v);
1652     gp_Pnt p = ElSLib::Value(u, v, Con2);
1653     if(aPApex1.SquareDistance(p) > tol2) {
1654       typeres=IntAna_NoGeometricSolution;
1655       return;
1656     }
1657     //
1658     ElSLib::Parameters(Con1, aPApex2, u, v);
1659     p = ElSLib::Value(u, v, Con1);
1660     if(aPApex2.SquareDistance(p) > tol2) {
1661       typeres=IntAna_NoGeometricSolution;
1662       return;
1663     }
1664
1665     //Cones have a common generatrix passing through apexes
1666     myCommonGen = Standard_True;
1667
1668     //common generatrix of cones
1669     gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1670
1671     //Intersection point of axes
1672     gp_Pnt aPAxeInt = A1A2.PtIntersect();
1673
1674     //Characteristic point of intersection curve
1675     u = ElCLib::Parameter(aGen, aPAxeInt);
1676     myPChar = ElCLib::Value(u, aGen);
1677
1678
1679     //Other generatrixes of cones laying in maximal plane
1680     gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1681     gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
1682     //
1683     //Intersection point of generatrixes
1684     gp_Dir aN; //solution plane normal
1685     gp_Dir aD1 = aGen1.Direction();
1686
1687     gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1688
1689     if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1690       aN = aD1.Crossed(aD2);
1691     }
1692     else if(aGen1.SquareDistance(aGen2) > tol2) {
1693       //Something wrong ???
1694       typeres=IntAna_NoGeometricSolution;
1695       return;
1696     }
1697     else {
1698       gp_Dir D1 = aGen1.Position().Direction();
1699       gp_Dir D2 = aGen2.Position().Direction();
1700       gp_Pnt O1 = aGen1.Location();
1701       gp_Pnt O2 = aGen2.Location();
1702       Standard_Real D1DotD2 = D1.Dot(D2);
1703       Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1704       gp_Vec O1O2 (O1,O2);
1705       Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1706       U2 /= aSin;
1707       gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1708
1709       aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1710       aN = aD1.Crossed(aD2);
1711     }
1712     //Plane that must contain intersection curves
1713     gp_Pln anIntPln(myPChar, aN);
1714
1715     IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1716
1717       if(INTER_QUAD_PLN.IsDone()) {
1718       switch(INTER_QUAD_PLN.TypeInter()) {
1719       case IntAna_Ellipse:         {
1720         typeres=IntAna_Ellipse;
1721         gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1722         pt1 = E.Location();
1723         dir1 = E.Position().Direction();
1724         dir2 = E.Position().XDirection();
1725         param1 = E.MajorRadius();
1726         param1bis = E.MinorRadius();
1727         nbint = 1;
1728         break;
1729       }
1730       case IntAna_Circle: {
1731         typeres=IntAna_Circle;
1732         gp_Circ C=INTER_QUAD_PLN.Circle(1);
1733         pt1 = C.Location();
1734         dir1 = C.Position().XDirection();
1735         dir2 = C.Position().YDirection();
1736         param1 = C.Radius();
1737         nbint = 1;
1738         break;
1739       }
1740       case IntAna_Parabola: {
1741         typeres=IntAna_Parabola;
1742         gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1743         pt1 = Prb.Location();
1744         dir1 = Prb.Position().Direction();
1745         dir2 = Prb.Position().XDirection();
1746         param1 = Prb.Focal();
1747         nbint = 1;
1748         break;
1749       }
1750       case IntAna_Hyperbola: {
1751         typeres=IntAna_Hyperbola;
1752         gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1753         pt1 = pt2 = H.Location();
1754         dir1 = H.Position().Direction();
1755         dir2 = H.Position().XDirection();
1756         param1 = param2 = H.MajorRadius();
1757         param1bis = param2bis = H.MinorRadius();
1758         nbint = 2;
1759         break;
1760       }
1761       default:
1762         typeres=IntAna_NoGeometricSolution;
1763       }
1764     }
1765   }
1766
1767   else {
1768     typeres=IntAna_NoGeometricSolution;
1769   }
1770 }
1771 //=======================================================================
1772 //function : IntAna_QuadQuadGeo
1773 //purpose  : Sphere - Cone
1774 //=======================================================================
1775   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1776                                          const gp_Cone& Con,
1777                                          const Standard_Real Tol)
1778 : done(Standard_False),
1779   nbint(0),
1780   typeres(IntAna_Empty),
1781   pt1(0,0,0),
1782   pt2(0,0,0),
1783   pt3(0,0,0),
1784   pt4(0,0,0),
1785   param1(0),
1786   param2(0),
1787   param3(0),
1788   param4(0),
1789   param1bis(0),
1790   param2bis(0),
1791   myCommonGen(Standard_False),
1792   myPChar(0,0,0)
1793 {
1794   InitTolerances();
1795   Perform(Sph,Con,Tol);
1796 }
1797 //=======================================================================
1798 //function : Perform
1799 //purpose  :
1800 //=======================================================================
1801   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1802                                    const gp_Cone& Con,
1803                                    const Standard_Real)
1804 {
1805
1806   //
1807   done=Standard_True;
1808   //
1809   AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1810   gp_Pnt Pt=Sph.Location();
1811   //
1812   if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1813      || A1A2.Same()) {
1814     gp_Pnt ConApex= Con.Apex();
1815     Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1816     gp_Dir ConDir;
1817     if(dApexSphCenter>RealEpsilon()) {
1818       ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1819     }
1820     else {
1821       ConDir = Con.Position().Direction();
1822     }
1823
1824     Standard_Real Rad=Sph.Radius();
1825     Standard_Real tga=Tan(Con.SemiAngle());
1826
1827
1828     //-- 2 circles
1829     //-- x: Roots of    (x**2 + y**2 = Rad**2)
1830     //--                tga = y / (x+dApexSphCenter)
1831     Standard_Real tgatga = tga * tga;
1832     math_DirectPolynomialRoots Eq( 1.0+tgatga
1833                                   ,2.0*tgatga*dApexSphCenter
1834                                   ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1835     if(Eq.IsDone()) {
1836       Standard_Integer nbsol=Eq.NbSolutions();
1837       if(nbsol==0) {
1838         typeres=IntAna_Empty;
1839       }
1840       else {
1841         typeres=IntAna_Circle;
1842         if(nbsol>=1) {
1843           Standard_Real x = Eq.Value(1);
1844           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1845           nbint=1;
1846           pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1847                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1848                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1849           param1 = tga * dApexSphCenterpx;
1850           param1 = Abs(param1);
1851           dir1 = ConDir;
1852           if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1853             typeres=IntAna_PointAndCircle;
1854             param1=0.0;
1855           }
1856         }
1857         if(nbsol>=2) {
1858           Standard_Real x=Eq.Value(2);
1859           Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1860           nbint=2;
1861           pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1862                        ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1863                        ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1864           param2 = tga * dApexSphCenterpx;
1865           param2 = Abs(param2);
1866           dir2=ConDir;
1867           if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1868             typeres=IntAna_PointAndCircle;
1869             param2=0.0;
1870           }
1871         }
1872       }
1873     }
1874     else {
1875       done=Standard_False;
1876     }
1877   }
1878   else {
1879     typeres=IntAna_NoGeometricSolution;
1880   }
1881 }
1882
1883 //=======================================================================
1884 //function : IntAna_QuadQuadGeo
1885 //purpose  : Sphere - Sphere
1886 //=======================================================================
1887   IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1888                                          ,const gp_Sphere& Sph2
1889                                          ,const Standard_Real Tol)
1890 : done(Standard_False),
1891   nbint(0),
1892   typeres(IntAna_Empty),
1893   pt1(0,0,0),
1894   pt2(0,0,0),
1895   pt3(0,0,0),
1896   pt4(0,0,0),
1897   param1(0),
1898   param2(0),
1899   param3(0),
1900   param4(0),
1901   param1bis(0),
1902   param2bis(0),
1903   myCommonGen(Standard_False),
1904   myPChar(0,0,0)
1905 {
1906   InitTolerances();
1907   Perform(Sph1,Sph2,Tol);
1908 }
1909 //=======================================================================
1910 //function : Perform
1911 //purpose  :
1912 //=======================================================================
1913   void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1914                                    const gp_Sphere& Sph2,
1915                                    const Standard_Real Tol)
1916 {
1917   done=Standard_True;
1918   gp_Pnt O1=Sph1.Location();
1919   gp_Pnt O2=Sph2.Location();
1920   Standard_Real dO1O2=O1.Distance(O2);
1921   Standard_Real R1=Sph1.Radius();
1922   Standard_Real R2=Sph2.Radius();
1923   Standard_Real Rmin,Rmax;
1924   typeres=IntAna_Empty;
1925   param2bis=0.0; //-- pour eviter param2bis not used ....
1926
1927   if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1928
1929   if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1930     typeres = IntAna_Same;
1931   }
1932   else {
1933     if(dO1O2<=Tol) { return; }
1934     gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1935     Standard_Real t = Rmax - dO1O2 - Rmin;
1936
1937     //----------------------------------------------------------------------
1938     //--        |----------------- R1 --------------------|
1939     //--        |----dO1O2-----|-----------R2----------|
1940     //--                                            --->--<-- t
1941     //--
1942     //--        |------ R1 ------|---------dO1O2----------|
1943     //--     |-------------------R2-----------------------|
1944     //--  --->--<-- t
1945     //----------------------------------------------------------------------
1946     if(t >= 0.0  && t <=Tol) {
1947       typeres = IntAna_Point;
1948       nbint = 1;
1949       Standard_Real t2;
1950       if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1951       else         t2=(-R1+(dO1O2-R2))*0.5;
1952
1953       pt1.SetCoord( O1.X() + t2*Dir.X()
1954                    ,O1.Y() + t2*Dir.Y()
1955                    ,O1.Z() + t2*Dir.Z());
1956     }
1957     else  {
1958       //-----------------------------------------------------------------
1959       //--        |----------------- dO1O2 --------------------|
1960       //--        |----R1-----|-----------R2----------|-Tol-|
1961       //--
1962       //--        |----------------- Rmax --------------------|
1963       //--        |----Rmin----|-------dO1O2-------|-Tol-|
1964       //--
1965       //-----------------------------------------------------------------
1966       if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1967         typeres=IntAna_Empty;
1968       }
1969       else {
1970         //---------------------------------------------------------------
1971         //--
1972         //--
1973         //---------------------------------------------------------------
1974         Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1975         Standard_Real Beta = R1*R1-Alpha*Alpha;
1976         Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1977
1978         if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1979           typeres = IntAna_Point;
1980           Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1981         }
1982         else {
1983           typeres = IntAna_Circle;
1984           dir1 = Dir;
1985           param1 = Beta;
1986         }
1987         pt1.SetCoord( O1.X() + Alpha*Dir.X()
1988                      ,O1.Y() + Alpha*Dir.Y()
1989                      ,O1.Z() + Alpha*Dir.Z());
1990
1991         nbint=1;
1992       }
1993     }
1994   }
1995 }
1996
1997 //=======================================================================
1998 //function : IntAna_QuadQuadGeo
1999 //purpose  : Plane - Torus
2000 //=======================================================================
2001 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
2002                                        const gp_Torus& Tor,
2003                                        const Standard_Real Tol)
2004 : done(Standard_False),
2005   nbint(0),
2006   typeres(IntAna_Empty),
2007   pt1(0,0,0),
2008   pt2(0,0,0),
2009   pt3(0,0,0),
2010   pt4(0,0,0),
2011   param1(0),
2012   param2(0),
2013   param3(0),
2014   param4(0),
2015   param1bis(0),
2016   param2bis(0),
2017   myCommonGen(Standard_False),
2018   myPChar(0,0,0)
2019 {
2020   InitTolerances();
2021   Perform(Pln,Tor,Tol);
2022 }
2023 //=======================================================================
2024 //function : Perform
2025 //purpose  :
2026 //=======================================================================
2027 void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
2028                                  const gp_Torus& Tor,
2029                                  const Standard_Real Tol)
2030 {
2031   done = Standard_True;
2032   //
2033   Standard_Real aRMin, aRMaj;
2034   //
2035   aRMin = Tor.MinorRadius();
2036   aRMaj = Tor.MajorRadius();
2037   if (aRMin >= aRMaj) {
2038     typeres = IntAna_NoGeometricSolution;
2039     return;
2040   }
2041   //
2042   const gp_Ax1 aPlnAx = Pln.Axis();
2043   const gp_Ax1 aTorAx = Tor.Axis();
2044   //
2045   Standard_Boolean bParallel, bNormal;
2046   //
2047   bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
2048   bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
2049   if (!bNormal && !bParallel) {
2050     typeres = IntAna_NoGeometricSolution;
2051     return;
2052   }
2053   //
2054   Standard_Real aDist;
2055   //
2056   gp_Pnt aTorLoc = aTorAx.Location();
2057   if (bParallel) {
2058     Standard_Real aDt, X, Y, Z, A, B, C, D;
2059     //
2060     Pln.Coefficients(A,B,C,D);
2061     aTorLoc.Coord(X,Y,Z);
2062     aDist = A*X + B*Y + C*Z + D;
2063     //
2064     if ((Abs(aDist) - aRMin) > Tol) {
2065       typeres=IntAna_Empty;
2066       return;
2067     }
2068     //
2069     typeres = IntAna_Circle;
2070     //
2071     pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
2072     aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2073     param1 = aRMaj + aDt;
2074     dir1 = aTorAx.Direction();
2075     nbint = 1;
2076     if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
2077       pt2 = pt1;
2078       param2 = aRMaj - aDt;
2079       dir2 = dir1;
2080       nbint = 2;
2081     }
2082   }
2083   //
2084   else {
2085     aDist = Pln.Distance(aTorLoc);
2086     if (aDist > myEPSILON_DISTANCE) {
2087       typeres = IntAna_NoGeometricSolution;
2088       return;
2089     }
2090     //
2091     typeres = IntAna_Circle;
2092     param2 = param1 = aRMin;
2093     dir2 = dir1 = aPlnAx.Direction();
2094     nbint = 2;
2095     //
2096     gp_Dir aDir = aTorAx.Direction()^dir1;
2097     pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
2098     pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
2099   }
2100 }
2101
2102 //=======================================================================
2103 //function : IntAna_QuadQuadGeo
2104 //purpose  : Cylinder - Torus
2105 //=======================================================================
2106 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
2107                                        const gp_Torus& Tor,
2108                                        const Standard_Real Tol)
2109 : done(Standard_False),
2110   nbint(0),
2111   typeres(IntAna_Empty),
2112   pt1(0,0,0),
2113   pt2(0,0,0),
2114   pt3(0,0,0),
2115   pt4(0,0,0),
2116   param1(0),
2117   param2(0),
2118   param3(0),
2119   param4(0),
2120   param1bis(0),
2121   param2bis(0),
2122   myCommonGen(Standard_False),
2123   myPChar(0,0,0)
2124 {
2125   InitTolerances();
2126   Perform(Cyl,Tor,Tol);
2127 }
2128 //=======================================================================
2129 //function : Perform
2130 //purpose  :
2131 //=======================================================================
2132 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
2133                                  const gp_Torus& Tor,
2134                                  const Standard_Real Tol)
2135 {
2136   done = Standard_True;
2137   //
2138   Standard_Real aRMin, aRMaj;
2139   //
2140   aRMin = Tor.MinorRadius();
2141   aRMaj = Tor.MajorRadius();
2142   if (aRMin >= aRMaj) {
2143     typeres = IntAna_NoGeometricSolution;
2144     return;
2145   }
2146   //
2147   const gp_Ax1 aCylAx = Cyl.Axis();
2148   const gp_Ax1 aTorAx = Tor.Axis();
2149   //
2150   const gp_Lin aLin(aTorAx);
2151   const gp_Pnt aLocCyl = Cyl.Location();
2152   //
2153   if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
2154       (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2155     typeres = IntAna_NoGeometricSolution;
2156     return;
2157   }
2158   //
2159   Standard_Real aRCyl;
2160   //
2161   aRCyl = Cyl.Radius();
2162   if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2163     typeres = IntAna_Empty;
2164     return;
2165   }
2166   //
2167   typeres = IntAna_Circle;
2168   //
2169   Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2170   gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2171   //
2172   dir1 = aTorAx.Direction();
2173   pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2174   param1 = aRCyl;
2175   nbint = 1;
2176   if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2177                        (aRCyl < (aRMaj + aRMin))) {
2178     dir2 = dir1;
2179     pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2180     param2 = param1;
2181     nbint = 2;
2182   }
2183 }
2184
2185 //=======================================================================
2186 //function : IntAna_QuadQuadGeo
2187 //purpose  : Cone - Torus
2188 //=======================================================================
2189 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2190                                        const gp_Torus& Tor,
2191                                        const Standard_Real Tol)
2192 : done(Standard_False),
2193   nbint(0),
2194   typeres(IntAna_Empty),
2195   pt1(0,0,0),
2196   pt2(0,0,0),
2197   pt3(0,0,0),
2198   pt4(0,0,0),
2199   param1(0),
2200   param2(0),
2201   param3(0),
2202   param4(0),
2203   param1bis(0),
2204   param2bis(0),
2205   myCommonGen(Standard_False),
2206   myPChar(0,0,0)
2207 {
2208   InitTolerances();
2209   Perform(Con,Tor,Tol);
2210 }
2211 //=======================================================================
2212 //function : Perform
2213 //purpose  :
2214 //=======================================================================
2215 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2216                                  const gp_Torus& Tor,
2217                                  const Standard_Real Tol)
2218 {
2219   done = Standard_True;
2220   //
2221   Standard_Real aRMin, aRMaj;
2222   //
2223   aRMin = Tor.MinorRadius();
2224   aRMaj = Tor.MajorRadius();
2225   if (aRMin >= aRMaj) {
2226     typeres = IntAna_NoGeometricSolution;
2227     return;
2228   }
2229   //
2230   const gp_Ax1 aConAx = Con.Axis();
2231   const gp_Ax1 aTorAx = Tor.Axis();
2232   //
2233   const gp_Lin aLin(aTorAx);
2234   const gp_Pnt aConApex = Con.Apex();
2235   //
2236   if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2237       (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2238     typeres = IntAna_NoGeometricSolution;
2239     return;
2240   }
2241   //
2242   Standard_Real anAngle, aDist, aParam[4], aDt;
2243   Standard_Integer i;
2244   gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2245   gp_Dir aDir[4];
2246   //
2247   anAngle = Con.SemiAngle();
2248   aTorLoc = aTorAx.Location();
2249   //
2250   aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2251   gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2252   gp_Ax1 anAxCLRot(aConApex, aDN);
2253   gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2254   gp_Dir aDL = aConL.Position().Direction();
2255   gp_Dir aXDir = Tor.XAxis().Direction();
2256   //
2257   typeres = IntAna_Empty;
2258   //
2259   for (i = 0; i < 2; ++i) {
2260     if (i) {
2261       aXDir.Reverse();
2262     }
2263     aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2264     //
2265     aDist = aConL.Distance(aPCT);
2266     if (aDist > aRMin+Tol) {
2267       continue;
2268     }
2269     //
2270     typeres = IntAna_Circle;
2271     //
2272     gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
2273     aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2274     //
2275     gp_Pnt aP;
2276     gp_XYZ aDVal = aDt*aDL.XYZ();
2277     aP.SetXYZ(aPh + aDVal);
2278     aParam[nbint] = aLin.Distance(aP);
2279     aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2280     aDir[nbint] = aTorAx.Direction();
2281     ++nbint;
2282     if ((aDist < aRMin) && (aDt > Tol)) {
2283       aP.SetXYZ(aPh - aDVal);
2284       aParam[nbint] = aLin.Distance(aP);
2285       aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2286       aDir[nbint] = aDir[nbint-1];
2287       ++nbint;
2288     }
2289   }
2290   //
2291   for (i = 0; i < nbint; ++i) {
2292     switch (i) {
2293     case 0:{
2294       pt1 = aPt[i];
2295       param1 = aParam[i];
2296       dir1 = aDir[i];
2297       break;
2298     }
2299     case 1:{
2300       pt2 = aPt[i];
2301       param2 = aParam[i];
2302       dir2 = aDir[i];
2303       break;
2304     }
2305     case 2:{
2306       pt3 = aPt[i];
2307       param3 = aParam[i];
2308       dir3 = aDir[i];
2309       break;
2310     }
2311     case 3:{
2312       pt4 = aPt[i];
2313       param4 = aParam[i];
2314       dir4 = aDir[i];
2315       break;
2316     }
2317     default:
2318       break;
2319     }
2320   }
2321 }
2322
2323 //=======================================================================
2324 //function : IntAna_QuadQuadGeo
2325 //purpose  : Sphere - Torus
2326 //=======================================================================
2327 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2328                                        const gp_Torus& Tor,
2329                                        const Standard_Real Tol)
2330 : done(Standard_False),
2331   nbint(0),
2332   typeres(IntAna_Empty),
2333   pt1(0,0,0),
2334   pt2(0,0,0),
2335   pt3(0,0,0),
2336   pt4(0,0,0),
2337   param1(0),
2338   param2(0),
2339   param3(0),
2340   param4(0),
2341   param1bis(0),
2342   param2bis(0),
2343   myCommonGen(Standard_False),
2344   myPChar(0,0,0)
2345 {
2346   InitTolerances();
2347   Perform(Sph,Tor,Tol);
2348 }
2349 //=======================================================================
2350 //function : Perform
2351 //purpose  :
2352 //=======================================================================
2353 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2354                                  const gp_Torus& Tor,
2355                                  const Standard_Real Tol)
2356 {
2357   done = Standard_True;
2358   //
2359   Standard_Real aRMin, aRMaj;
2360   //
2361   aRMin = Tor.MinorRadius();
2362   aRMaj = Tor.MajorRadius();
2363   if (aRMin >= aRMaj) {
2364     typeres = IntAna_NoGeometricSolution;
2365     return;
2366   }
2367   //
2368   const gp_Ax1 aTorAx = Tor.Axis();
2369   const gp_Lin aLin(aTorAx);
2370   const gp_Pnt aSphLoc = Sph.Location();
2371   //
2372   if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2373     typeres = IntAna_NoGeometricSolution;
2374     return;
2375   }
2376   //
2377   Standard_Real aRSph, aDist;
2378   gp_Pnt aTorLoc;
2379   //
2380   gp_Dir aXDir = Tor.XAxis().Direction();
2381   aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2382   aRSph = Sph.Radius();
2383   //
2384   gp_Vec aVec12(aTorLoc, aSphLoc);
2385   aDist = aVec12.Magnitude();
2386   if (((aDist - Tol) > (aRMin + aRSph)) ||
2387       ((aDist + Tol) < Abs(aRMin - aRSph))) {
2388     typeres = IntAna_Empty;
2389     return;
2390   }
2391   //
2392   typeres = IntAna_Circle;
2393   //
2394   Standard_Real anAlpha, aBeta;
2395   //
2396   anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2397   aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2398   //
2399   gp_Dir aDir12(aVec12);
2400   gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2401   gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2402   //
2403   gp_Pnt aP;
2404   gp_XYZ aDVal = aBeta*aDC.XYZ();
2405   aP.SetXYZ(aPh + aDVal);
2406   param1 = aLin.Distance(aP);
2407   pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2408   dir1 = aTorAx.Direction();
2409   nbint = 1;
2410   if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) &&
2411       (aDVal.Modulus() > Tol)) {
2412     aP.SetXYZ(aPh - aDVal);
2413     param2 = aLin.Distance(aP);
2414     pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2415     dir2 = dir1;
2416     nbint = 2;
2417   }
2418 }
2419
2420 //=======================================================================
2421 //function : IntAna_QuadQuadGeo
2422 //purpose  : Torus - Torus
2423 //=======================================================================
2424 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2425                                        const gp_Torus& Tor2,
2426                                        const Standard_Real Tol)
2427 : done(Standard_False),
2428   nbint(0),
2429   typeres(IntAna_Empty),
2430   pt1(0,0,0),
2431   pt2(0,0,0),
2432   pt3(0,0,0),
2433   pt4(0,0,0),
2434   param1(0),
2435   param2(0),
2436   param3(0),
2437   param4(0),
2438   param1bis(0),
2439   param2bis(0),
2440   myCommonGen(Standard_False),
2441   myPChar(0,0,0)
2442 {
2443   InitTolerances();
2444   Perform(Tor1,Tor2,Tol);
2445 }
2446 //=======================================================================
2447 //function : Perform
2448 //purpose  :
2449 //=======================================================================
2450 void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2451                                  const gp_Torus& Tor2,
2452                                  const Standard_Real Tol)
2453 {
2454   done = Standard_True;
2455   //
2456   Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2457   //
2458   aRMin1 = Tor1.MinorRadius();
2459   aRMaj1 = Tor1.MajorRadius();
2460   aRMin2 = Tor2.MinorRadius();
2461   aRMaj2 = Tor2.MajorRadius();
2462   if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2463     typeres = IntAna_NoGeometricSolution;
2464     return;
2465   }
2466   //
2467   const gp_Ax1 anAx1 = Tor1.Axis();
2468   const gp_Ax1 anAx2 = Tor2.Axis();
2469   //
2470   gp_Lin aL1(anAx1);
2471   if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2472       (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2473     typeres = IntAna_NoGeometricSolution;
2474     return;
2475   }
2476   //
2477   gp_Pnt aLoc1, aLoc2;
2478   //
2479   aLoc1 = anAx1.Location();
2480   aLoc2 = anAx2.Location();
2481   //
2482   if (aLoc1.IsEqual(aLoc2, Tol) &&
2483       (Abs(aRMin1 - aRMin2) <= Tol) &&
2484       (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2485     typeres = IntAna_Same;
2486     return;
2487   }
2488   //
2489   Standard_Real aDist;
2490   gp_Pnt aP1, aP2;
2491   //
2492   gp_Dir aXDir1 = Tor1.XAxis().Direction();
2493   aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2494   aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2495   //
2496   gp_Vec aV12(aP1, aP2);
2497   aDist = aV12.Magnitude();
2498   if (((aDist - Tol) > (aRMin1 + aRMin2)) ||
2499       ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2500     typeres = IntAna_Empty;
2501     return;
2502   }
2503   //
2504   typeres = IntAna_Circle;
2505   //
2506   Standard_Real anAlpha, aBeta;
2507   //
2508   anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2509   aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2510   //
2511   gp_Dir aDir12(aV12);
2512   gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2513   gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2514   //
2515   gp_Pnt aP;
2516   gp_XYZ aDVal = aBeta*aDC.XYZ();
2517   aP.SetXYZ(aPh + aDVal);
2518   param1 = aL1.Distance(aP);
2519   pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2520   dir1 = anAx1.Direction();
2521   nbint = 1;
2522   if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) &&
2523       aDVal.Modulus() > Tol) {
2524     aP.SetXYZ(aPh - aDVal);
2525     param2 = aL1.Distance(aP);
2526     pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2527     dir2 = dir1;
2528     nbint = 2;
2529   }
2530 }
2531
2532 //=======================================================================
2533 //function : Point
2534 //purpose  : Returns a Point
2535 //=======================================================================
2536   gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
2537 {
2538   if(!done)          {    StdFail_NotDone::Raise();        }
2539   if(n>nbint || n<1) {    Standard_DomainError::Raise();   }
2540   if(typeres==IntAna_PointAndCircle) {
2541     if(n!=1) { Standard_DomainError::Raise();  }
2542     if(param1==0.0) return(pt1);
2543     return(pt2);
2544   }
2545   else if(typeres==IntAna_Point) {
2546     if(n==1) return(pt1);
2547     return(pt2);
2548   }
2549
2550   // WNT (what can you expect from MicroSoft ?)
2551   return gp_Pnt(0,0,0);
2552 }
2553 //=======================================================================
2554 //function : Line
2555 //purpose  : Returns a Line
2556 //=======================================================================
2557   gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
2558 {
2559   if(!done)        {   StdFail_NotDone::Raise();   }
2560   if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2561     Standard_DomainError::Raise();
2562     }
2563   if(n==1) {  return(gp_Lin(pt1,dir1));   }
2564   else {      return(gp_Lin(pt2,dir2));   }
2565 }
2566 //=======================================================================
2567 //function : Circle
2568 //purpose  : Returns a Circle
2569 //=======================================================================
2570   gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
2571 {
2572   if(!done) {    StdFail_NotDone::Raise();     }
2573   if(typeres==IntAna_PointAndCircle) {
2574     if(n!=1) { Standard_DomainError::Raise();  }
2575     if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2576     return(gp_Circ(DirToAx2(pt2,dir2),param2));
2577   }
2578   else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2579     Standard_DomainError::Raise();
2580     }
2581   if      (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2582   else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2583   else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2584   else           { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
2585 }
2586
2587 //=======================================================================
2588 //function : Ellipse
2589 //purpose  : Returns a Elips
2590 //=======================================================================
2591   gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2592 {
2593   if(!done) {     StdFail_NotDone::Raise();     }
2594   if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2595     Standard_DomainError::Raise();
2596   }
2597
2598   if(n==1) {
2599     Standard_Real R1=param1, R2=param1bis, aTmp;
2600     if (R1<R2) {
2601       aTmp=R1; R1=R2; R2=aTmp;
2602     }
2603     gp_Ax2 anAx2(pt1, dir1 ,dir2);
2604     gp_Elips anElips (anAx2, R1, R2);
2605     return anElips;
2606   }
2607   else {
2608     Standard_Real R1=param2, R2=param2bis, aTmp;
2609     if (R1<R2) {
2610       aTmp=R1; R1=R2; R2=aTmp;
2611     }
2612     gp_Ax2 anAx2(pt2, dir2 ,dir1);
2613     gp_Elips anElips (anAx2, R1, R2);
2614     return anElips;
2615   }
2616 }
2617 //=======================================================================
2618 //function : Parabola
2619 //purpose  : Returns a Parabola
2620 //=======================================================================
2621   gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
2622 {
2623   if(!done) {
2624     StdFail_NotDone::Raise();
2625     }
2626   if (typeres!=IntAna_Parabola) {
2627     Standard_DomainError::Raise();
2628   }
2629   if((n>nbint) || (n!=1)) {
2630     Standard_OutOfRange::Raise();
2631   }
2632   return(gp_Parab(gp_Ax2( pt1
2633                          ,dir1
2634                          ,dir2)
2635                   ,param1));
2636 }
2637 //=======================================================================
2638 //function : Hyperbola
2639 //purpose  : Returns a Hyperbola
2640 //=======================================================================
2641   gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
2642 {
2643   if(!done) {
2644     StdFail_NotDone::Raise();
2645     }
2646   if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2647     Standard_DomainError::Raise();
2648     }
2649   if(n==1) {
2650     return(gp_Hypr(gp_Ax2( pt1
2651                           ,dir1
2652                           ,dir2)
2653                    ,param1,param1bis));
2654   }
2655   else {
2656     return(gp_Hypr(gp_Ax2( pt2
2657                           ,dir1
2658                           ,dir2.Reversed())
2659                    ,param2,param2bis));
2660   }
2661 }
2662 //=======================================================================
2663 //function : HasCommonGen
2664 //purpose  :
2665 //=======================================================================
2666 Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2667 {
2668   return myCommonGen;
2669 }
2670 //=======================================================================
2671 //function : PChar
2672 //purpose  :
2673 //=======================================================================
2674 const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2675 {
2676   return myPChar;
2677 }
2678 //=======================================================================
2679 //function : RefineDir
2680 //purpose  :
2681 //=======================================================================
2682 void RefineDir(gp_Dir& aDir)
2683 {
2684   Standard_Integer k, m, n;
2685   Standard_Real aC[3];
2686   //
2687   aDir.Coord(aC[0], aC[1], aC[2]);
2688   //
2689   m=0;
2690   n=0;
2691   for (k=0; k<3; ++k) {
2692     if (aC[k]==1. || aC[k]==-1.) {
2693       ++m;
2694     }
2695     else if (aC[k]!=0.) {
2696       ++n;
2697     }
2698   }
2699   //
2700   if (m && n) {
2701     Standard_Real aEps, aR1, aR2, aNum;
2702     //
2703     aEps=RealEpsilon();
2704     aR1=1.-aEps;
2705     aR2=1.+aEps;
2706     //
2707     for (k=0; k<3; ++k) {
2708       m=(aC[k]>0.);
2709       aNum=(m)? aC[k] : -aC[k];
2710       if (aNum>aR1 && aNum<aR2) {
2711         if (m) {
2712           aC[k]=1.;
2713         }
2714         else {
2715           aC[k]=-1.;
2716         }
2717         //
2718         aC[(k+1)%3]=0.;
2719         aC[(k+2)%3]=0.;
2720         break;
2721       }
2722     }
2723     aDir.SetCoord(aC[0], aC[1], aC[2]);
2724   }
2725 }
2726
2727
2728