src/Extrema/Extrema_ExtElCS.cxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and/or modify it under
   7 // the terms of the GNU Lesser General Public License version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 //  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134
  16
  17 #include <ElCLib.hxx>
  18 #include <ElSLib.hxx>
  19 #include <Extrema_ExtElC.hxx>
  20 #include <Extrema_ExtElCS.hxx>
  21 #include <Extrema_ExtPElC.hxx>
  22 #include <Extrema_ExtPElS.hxx>
  23 #include <Extrema_POnCurv.hxx>
  24 #include <Extrema_POnSurf.hxx>
  25 #include <gp_Circ.hxx>
  26 #include <gp_Cone.hxx>
  27 #include <gp_Cylinder.hxx>
  28 #include <gp_Hypr.hxx>
  29 #include <gp_Lin.hxx>
  30 #include <gp_Pln.hxx>
  31 #include <gp_Sphere.hxx>
  32 #include <gp_Torus.hxx>
  33 #include <gp_Vec.hxx>
  34 #include <IntAna_IntConicQuad.hxx>
  35 #include <IntAna_Quadric.hxx>
  36 #include <IntAna_QuadQuadGeo.hxx>
  37 #include <Precision.hxx>
  38 #include <Standard_NotImplemented.hxx>
  39 #include <Standard_OutOfRange.hxx>
  40 #include <StdFail_NotDone.hxx>
  41 #include <TColStd_ListOfInteger.hxx>
  42
  43 Extrema_ExtElCS::Extrema_ExtElCS()
  44 {
  45   myDone = Standard_False;
  46   myIsPar = Standard_False;
  47   myNbExt = 0;
  48 }
  49
  50
  51 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
  52                                  const gp_Pln& S)
  53 {
  54   Perform(C, S);
  55 }
  56
  57
  58
  59 void Extrema_ExtElCS::Perform(const gp_Lin& C,
  60                               const gp_Pln& S)
  61 {
  62   myDone = Standard_True;
  63   myIsPar = Standard_False;
  64   myNbExt = 0;
  65
  66   if (C.Direction().IsNormal(S.Axis().Direction(),
  67                             Precision::Angular()))
  68   {
  69     mySqDist = new TColStd_HArray1OfReal(1, 1);
  70     mySqDist->SetValue(1, S.SquareDistance(C));
  71     myIsPar = Standard_True;
  72     myNbExt = 1;
  73   }
  74 }
  75
  76
  77 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
  78                                  const gp_Cylinder& S)
  79 {
  80   Perform(C, S);
  81 }
  82
  83
  84
  85 void Extrema_ExtElCS::Perform(const gp_Lin& C,
  86                               const gp_Cylinder& S)
  87 {
  88   myDone = Standard_False;
  89   myNbExt = 0;
  90   myIsPar = Standard_False;
  91
  92   gp_Ax3 Pos = S.Position();
  93
  94   Standard_Boolean isParallel = Standard_False;
  95
  96   Standard_Real radius = S.Radius();
  97   Extrema_ExtElC Extrem(gp_Lin(Pos.Axis()), C, Precision::Angular());
  98   if (Extrem.IsParallel())
  99   {
 100     isParallel = Standard_True;
 101   }
 102   else
 103   {
 104     Standard_Integer i, aStartIdx = 0;
 105
 106     Extrema_POnCurv myPOnC1, myPOnC2;
 107     Extrem.Points(1, myPOnC1, myPOnC2);
 108     gp_Pnt PonAxis = myPOnC1.Value();
 109     gp_Pnt PC = myPOnC2.Value();
 110
 111     // line intersects the cylinder
 112     if (radius - PonAxis.Distance(PC) > Precision::PConfusion())
 113     {
 114       IntAna_Quadric theQuadric(S);
 115       IntAna_IntConicQuad Inters(C, theQuadric);
 116       if (Inters.IsDone() && Inters.IsInQuadric())
 117       {
 118         isParallel = Standard_True;
 119       }
 120       else if (Inters.IsDone())
 121       {
 122         myNbExt = Inters.NbPoints();
 123         aStartIdx = myNbExt;
 124         if (myNbExt > 0)
 125         {
 126           // Not more than 2 additional points from perpendiculars.
 127           mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
 128           myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
 129           myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
 130           Standard_Real u, v, w;
 131           for (i = 1; i <= myNbExt; i++)
 132           {
 133             mySqDist->SetValue(i, 0.);
 134             gp_Pnt P_int = Inters.Point(i);
 135             w = Inters.ParamOnConic(i);
 136             Extrema_POnCurv PonC(w, P_int);
 137             myPoint1->SetValue(i, PonC);
 138             ElSLib::CylinderParameters(Pos, radius, P_int, u, v);
 139             Extrema_POnSurf PonS(u, v, P_int);
 140             myPoint2->SetValue(i, PonS);
 141           }
 142         }
 143       }
 144     }
 145     else
 146     {
 147       // line is tangent or outside of the cylinder
 148       Extrema_ExtPElS ExPS(PC, S, Precision::Confusion());
 149       if (ExPS.IsDone())
 150       {
 151         if (aStartIdx == 0)
 152         {
 153           myNbExt = ExPS.NbExt();
 154           mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
 155           myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
 156           myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
 157         }
 158         else
 159           myNbExt += ExPS.NbExt();
 160
 161         for (i = aStartIdx + 1; i <= myNbExt; i++)
 162         {
 163           myPoint1->SetValue(i, myPOnC2);
 164           myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
 165           mySqDist->SetValue(i, (myPOnC2.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
 166         }
 167       }
 168     }
 169
 170     myDone = Standard_True;
 171   }
 172
 173   if (isParallel)
 174   {
 175     // Line direction is similar to cylinder axis of rotation.
 176     // The case is possible when either extrema returned parallel status
 177     // or Intersection tool returned infinite number of solutions.
 178     // This is possible due to Intersection algorithm uses more precise
 179     // characteristics to consider given geometries parallel.
 180     // In the latter case there may be several extremas, thus we look for
 181     // the one with the lowest distance and use it as a final solution.
 182     mySqDist = new TColStd_HArray1OfReal(1, 1);
 183     Standard_Real aDist = Extrem.SquareDistance(1);
 184     const Standard_Integer aNbExt = Extrem.NbExt();
 185     for (Standard_Integer i = 2; i <= aNbExt; i++)
 186     {
 187       const Standard_Real aD = Extrem.SquareDistance(i);
 188       if (aD < aDist)
 189       {
 190         aDist = aD;
 191       }
 192     }
 193
 194     aDist = sqrt(aDist) - radius;
 195     mySqDist->SetValue(1, aDist * aDist);
 196     myDone = Standard_True;
 197     myIsPar = Standard_True;
 198     myNbExt = 1;
 199   }
 200 }
 201
 202
 203
 204 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
 205                                  const gp_Cone& S)
 206 {  Perform(C, S);}
 207
 208
 209
 210 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
 211 //                            const gp_Cone& S)
 212 void Extrema_ExtElCS::Perform(const gp_Lin& ,
 213                               const gp_Cone& )
 214 {
 215   throw Standard_NotImplemented();
 216
 217 }
 218
 219
 220
 221 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
 222                                  const gp_Sphere& S)
 223 {
 224   Perform(C, S);
 225 }
 226
 227
 228
 229 void Extrema_ExtElCS::Perform(const gp_Lin& C,
 230                               const gp_Sphere& S)
 231 {
 232   // In case of intersection - return four points:
 233   // 2 intersection points and 2 perpendicular.
 234   // No intersection - only min and max.
 235
 236   myDone = Standard_False;
 237   myNbExt = 0;
 238   myIsPar = Standard_False;
 239   Standard_Integer aStartIdx = 0;
 240
 241   gp_Pnt aCenter = S.Location();
 242
 243   Extrema_ExtPElC Extrem(aCenter, C, Precision::Angular(), RealFirst(), RealLast());
 244
 245   Standard_Integer i;
 246   if (Extrem.IsDone() &&
 247       Extrem.NbExt() > 0)
 248   {
 249     Extrema_POnCurv myPOnC1 =  Extrem.Point(1);
 250     if (myPOnC1.Value().Distance(aCenter) <= S.Radius())
 251     {
 252       IntAna_IntConicQuad aLinSphere(C, S);
 253       if (aLinSphere.IsDone())
 254       {
 255         myNbExt = aLinSphere.NbPoints();
 256         aStartIdx = myNbExt;
 257         // Not more than 2 additional points from perpendiculars.
 258         mySqDist = new TColStd_HArray1OfReal(1, myNbExt + 2);
 259         myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt + 2);
 260         myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt + 2);
 261
 262         for (i = 1; i <= myNbExt; i++)
 263         {
 264           Extrema_POnCurv aCPnt(aLinSphere.ParamOnConic(i), aLinSphere.Point(i));
 265
 266           Standard_Real u,v;
 267           ElSLib::Parameters(S, aLinSphere.Point(i), u, v);
 268           Extrema_POnSurf aSPnt(u, v, aLinSphere.Point(i));
 269
 270           myPoint1->SetValue(i, aCPnt);
 271           myPoint2->SetValue(i, aSPnt);
 272           mySqDist->SetValue(i,(aCPnt.Value()).SquareDistance(aSPnt.Value()));
 273         }
 274       }
 275     }
 276
 277     Extrema_ExtPElS ExPS(myPOnC1.Value(), S, Precision::Confusion());
 278     if (ExPS.IsDone())
 279     {
 280       if (aStartIdx == 0)
 281       {
 282         myNbExt = ExPS.NbExt();
 283         mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
 284         myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
 285         myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
 286       }
 287       else
 288         myNbExt += ExPS.NbExt();
 289
 290       for (i = aStartIdx + 1; i <= myNbExt; i++)
 291       {
 292         myPoint1->SetValue(i, myPOnC1);
 293         myPoint2->SetValue(i, ExPS.Point(i - aStartIdx));
 294         mySqDist->SetValue(i,(myPOnC1.Value()).SquareDistance(ExPS.Point(i - aStartIdx).Value()));
 295       }
 296     }
 297   }
 298   myDone = Standard_True;
 299 }
 300
 301
 302 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Lin& C,
 303                                  const gp_Torus& S)
 304 {  Perform(C, S);}
 305
 306
 307
 308 //void Extrema_ExtElCS::Perform(const gp_Lin& C,
 309 //                            const gp_Torus& S)
 310 void Extrema_ExtElCS::Perform(const gp_Lin& ,
 311                               const gp_Torus& )
 312 {
 313   throw Standard_NotImplemented();
 314
 315 }
 316
 317
 318 //        Circle-?
 319
 320 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
 321                                  const gp_Pln& S)
 322 {
 323   Perform(C, S);
 324 }
 325
 326
 327
 328 void Extrema_ExtElCS::Perform(const gp_Circ& C,
 329                               const gp_Pln& S)
 330 {
 331   myDone = Standard_True;
 332   myIsPar = Standard_False;
 333   myNbExt = 0;
 334
 335   gp_Ax2 Pos = C.Position();
 336   gp_Dir NCirc = Pos.Direction();
 337   gp_Dir NPln = S.Axis().Direction();
 338
 339   Standard_Boolean isParallel = Standard_False;
 340
 341   if (NCirc.IsParallel(NPln, Precision::Angular())) {
 342     isParallel = Standard_True;
 343   }
 344   else {
 345
 346     gp_Dir ExtLine = NCirc ^ NPln;
 347     ExtLine = ExtLine ^ NCirc;
 348     //
 349     gp_Dir XDir = Pos.XDirection();
 350     Standard_Real T[2];
 351     T[0] = XDir.AngleWithRef(ExtLine, NCirc);
 352     if(T[0] < 0.)
 353     {
 354       //Put in period
 355       T[0] += M_PI;
 356     }
 357     T[1] = T[0] + M_PI;
 358     //
 359     myNbExt = 2;
 360     //Check intersection
 361     IntAna_IntConicQuad anInter(C, S,
 362                                 Precision::Angular(),
 363                                 Precision::Confusion());
 364
 365     if (anInter.IsDone() && anInter.IsInQuadric())
 366     {
 367       isParallel = Standard_True;
 368     }
 369     else if (anInter.IsDone())
 370     {
 371       if(anInter.NbPoints() > 1)
 372       {
 373         myNbExt += anInter.NbPoints();
 374       }
 375     }
 376
 377     if (!isParallel)
 378     {
 379       myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
 380       mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
 381       myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
 382
 383       Standard_Integer i;
 384       gp_Pnt PC, PP;
 385       Standard_Real U, V;
 386       Extrema_POnCurv POnC;
 387       Extrema_POnSurf POnS;
 388       for (i = 0; i < 2; ++i)
 389       {
 390         PC = ElCLib::CircleValue(T[i], C.Position(), C.Radius());
 391         POnC.SetValues(T[i], PC);
 392         myPoint1->SetValue(i + 1, POnC);
 393         ElSLib::PlaneParameters(S.Position(), PC, U, V);
 394         PP = ElSLib::PlaneValue(U, V, S.Position());
 395         POnS.SetParameters(U, V, PP);
 396         myPoint2->SetValue(i + 1, POnS);
 397         mySqDist->SetValue(i + 1, PC.SquareDistance(PP));
 398       }
 399       //
 400       if (myNbExt > 2)
 401       {
 402         //Add intersection points
 403         for (i = 1; i <= anInter.NbPoints(); ++i)
 404         {
 405           Standard_Real t = anInter.ParamOnConic(i);
 406           PC = ElCLib::CircleValue(t, C.Position(), C.Radius());
 407           POnC.SetValues(t, PC);
 408           myPoint1->SetValue(i + 2, POnC);
 409           ElSLib::PlaneParameters(S.Position(), PC, U, V);
 410           PP = ElSLib::PlaneValue(U, V, S.Position());
 411           POnS.SetParameters(U, V, PP);
 412           myPoint2->SetValue(i + 2, POnS);
 413           mySqDist->SetValue(i + 2, PC.SquareDistance(PP));
 414         }
 415       }
 416     }
 417   }
 418
 419   if (isParallel)
 420   {
 421     mySqDist = new TColStd_HArray1OfReal(1, 1);
 422     mySqDist->SetValue(1, S.SquareDistance(C.Location()));
 423     myIsPar = Standard_True;
 424     myNbExt = 1;
 425   }
 426 }
 427
 428 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
 429                                  const gp_Cylinder& S)
 430 {
 431   Perform(C, S);
 432 }
 433
 434
 435
 436 //  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134 Begin
 437 // Implementation of the method.
 438 void Extrema_ExtElCS::Perform(const gp_Circ& C,
 439                               const gp_Cylinder& S)
 440 {
 441   myDone  = Standard_False;
 442   myIsPar = Standard_False;
 443   myNbExt = 0;
 444
 445   // Get an axis line of the cylinder.
 446   gp_Lin anAxis(S.Axis());
 447
 448   // Compute extrema between the circle and the line.
 449   Extrema_ExtElC anExtC(anAxis, C, 0.);
 450
 451   if (!anExtC.IsDone())
 452     return;
 453
 454   Standard_Boolean isParallel = Standard_False;
 455
 456   if (anExtC.IsParallel()) {
 457     isParallel = Standard_True;
 458   } else {
 459     Standard_Integer aNbExt   = anExtC.NbExt();
 460     Standard_Integer i;
 461     Standard_Integer aCurI    = 1;
 462     Standard_Real    aTolConf = Precision::Confusion();
 463     Standard_Real    aCylRad  = S.Radius();
 464
 465     // Check whether two objects have intersection points
 466     IntAna_Quadric aCylQuad(S);
 467     IntAna_IntConicQuad aCircCylInter(C, aCylQuad);
 468     Standard_Integer aNbInter = 0;
 469     if (aCircCylInter.IsDone() && aCircCylInter.IsInQuadric())
 470     {
 471       isParallel = Standard_True;
 472     }
 473     else if (aCircCylInter.IsDone())
 474     {
 475       aNbInter = aCircCylInter.NbPoints();
 476     }
 477
 478     if (!isParallel)
 479     {
 480       // Compute the extremas.
 481       myNbExt = 2 * aNbExt + aNbInter;
 482       mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
 483       myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
 484       myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
 485
 486       for (i = 1; i <= aNbExt; i++) {
 487         Extrema_POnCurv aPOnAxis;
 488         Extrema_POnCurv aPOnCirc;
 489         Standard_Real   aSqDist = anExtC.SquareDistance(i);
 490         Standard_Real   aDist = sqrt(aSqDist);
 491
 492         anExtC.Points(i, aPOnAxis, aPOnCirc);
 493
 494         if (aSqDist <= (aTolConf * aTolConf)) {
 495           myNbExt -= 2;
 496           continue;
 497         }
 498
 499         gp_Dir aDir(aPOnAxis.Value().XYZ().Subtracted(aPOnCirc.Value().XYZ()));
 500         Standard_Real aShift[2] = {aDist + aCylRad, aDist - aCylRad};
 501         Standard_Integer j;
 502
 503         for (j = 0; j < 2; j++) {
 504           gp_Vec aVec(aDir);
 505           gp_Pnt aPntOnCyl;
 506
 507           aVec.Multiply(aShift[j]);
 508           aPntOnCyl = aPOnCirc.Value().Translated(aVec);
 509
 510           Standard_Real aU;
 511           Standard_Real aV;
 512
 513           ElSLib::Parameters(S, aPntOnCyl, aU, aV);
 514
 515           Extrema_POnSurf aPOnSurf(aU, aV, aPntOnCyl);
 516
 517           myPoint1->SetValue(aCurI, aPOnCirc);
 518           myPoint2->SetValue(aCurI, aPOnSurf);
 519           mySqDist->SetValue(aCurI++, aShift[j] * aShift[j]);
 520         }
 521       }
 522
 523       // Adding intersection points to the list of extremas
 524       for (i=1; i<=aNbInter; i++)
 525       {
 526         Standard_Real aU;
 527         Standard_Real aV;
 528
 529         gp_Pnt aInterPnt = aCircCylInter.Point(i);
 530
 531         aU = ElCLib::Parameter(C, aInterPnt);
 532         Extrema_POnCurv aPOnCirc(aU, aInterPnt);
 533
 534         ElSLib::Parameters(S, aInterPnt, aU, aV);
 535         Extrema_POnSurf aPOnCyl(aU, aV, aInterPnt);
 536         myPoint1->SetValue(aCurI, aPOnCirc);
 537         myPoint2->SetValue(aCurI, aPOnCyl);
 538         mySqDist->SetValue(aCurI++, 0.0);
 539       }
 540     }
 541   }
 542
 543   myDone = Standard_True;
 544
 545   if (isParallel)
 546   {
 547     // The case is possible when either extrema returned parallel status
 548     // or Intersection tool returned infinite number of solutions.
 549     // This is possible due to Intersection algorithm uses more precise
 550     // characteristics to consider given geometries parallel.
 551     // In the latter case there may be several extremas, thus we look for
 552     // the one with the lowest distance and use it as a final solution.
 553
 554     myIsPar = Standard_True;
 555     myNbExt = 1;
 556     mySqDist = new TColStd_HArray1OfReal(1, 1);
 557     Standard_Real aDist = anExtC.SquareDistance(1);
 558
 559     const Standard_Integer aNbExt = anExtC.NbExt();
 560     for (Standard_Integer i = 2; i <= aNbExt; i++)
 561     {
 562       const Standard_Real aD = anExtC.SquareDistance(i);
 563       if (aD < aDist)
 564       {
 565         aDist = aD;
 566       }
 567     }
 568
 569     aDist = sqrt(aDist) - S.Radius();
 570     mySqDist->SetValue(1, aDist * aDist);
 571   }
 572 }
 573 //  Modified by skv - Thu Jul  7 14:37:05 2005 OCC9134 End
 574
 575
 576
 577 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
 578                                  const gp_Cone& S)
 579 {  Perform(C, S);}
 580
 581
 582
 583 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
 584 //                       const gp_Cone& S)
 585 void Extrema_ExtElCS::Perform(const gp_Circ& ,
 586                          const gp_Cone& )
 587 {
 588   throw Standard_NotImplemented();
 589
 590 }
 591
 592
 593
 594 //=======================================================================
 595 //function : Extrema_ExtElCS
 596 //purpose  : Circle/Sphere
 597 //=======================================================================
 598 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
 599                                  const gp_Sphere& S)
 600 {
 601   Perform(C, S);
 602 }
 603
 604 //=======================================================================
 605 //function : Perform
 606 //purpose  : Circle/Sphere
 607 //=======================================================================
 608 void Extrema_ExtElCS::Perform(const gp_Circ& C,
 609                               const gp_Sphere& S)
 610 {
 611   myDone = Standard_False;
 612   myIsPar = Standard_False;
 613   myNbExt = 0;
 614
 615   if (gp_Lin(C.Axis()).SquareDistance(S.Location()) < Precision::SquareConfusion())
 616   {
 617     // Circle and sphere are parallel
 618     myIsPar = Standard_True;
 619     myDone = Standard_True;
 620     myNbExt = 1;
 621
 622     // Compute distance from circle to the sphere
 623     Standard_Real aSqDistLoc = C.Location().SquareDistance(S.Location());
 624     Standard_Real aSqDist = aSqDistLoc + C.Radius() * C.Radius();
 625     Standard_Real aDist = sqrt(aSqDist) - S.Radius();
 626     mySqDist = new TColStd_HArray1OfReal(1, 1);
 627     mySqDist->SetValue(1, aDist * aDist);
 628     return;
 629   }
 630
 631   // Intersect sphere with circle's plane
 632   gp_Pln CPln(C.Location(), C.Axis().Direction());
 633   IntAna_QuadQuadGeo anInter(CPln, S);
 634   if (!anInter.IsDone())
 635     // not done
 636     return;
 637
 638   if (anInter.TypeInter() != IntAna_Circle)
 639   {
 640     // Intersection is empty or just a point.
 641     // The parallel case has already been considered,
 642     // thus, here we have to find only one minimal solution
 643     myNbExt = 1;
 644     myDone = Standard_True;
 645
 646     mySqDist = new TColStd_HArray1OfReal(1, 1);
 647     myPoint1 = new Extrema_HArray1OfPOnCurv(1, 1);
 648     myPoint2 = new Extrema_HArray1OfPOnSurf(1, 1);
 649
 650     // Compute parameter on circle
 651     const Standard_Real aT = ElCLib::Parameter(C, S.Location());
 652     // Compute point on circle
 653     gp_Pnt aPOnC = ElCLib::Value(aT, C);
 654
 655     // Compute parameters on sphere
 656     Standard_Real aU, aV;
 657     ElSLib::Parameters(S, aPOnC, aU, aV);
 658     // Compute point on sphere
 659     gp_Pnt aPOnS = ElSLib::Value(aU, aV, S);
 660
 661     // Save solution
 662     myPoint1->SetValue(1, Extrema_POnCurv(aT, aPOnC));
 663     myPoint2->SetValue(1, Extrema_POnSurf(aU, aV, aPOnS));
 664     mySqDist->SetValue(1, aPOnC.SquareDistance(aPOnS));
 665     return;
 666   }
 667
 668   // Here, the intersection is a circle
 669
 670   // Intersection circle
 671   gp_Circ aCInt = anInter.Circle(1);
 672
 673   // Perform intersection of the input circle with the intersection circle
 674   Extrema_ExtElC anExtC(C, aCInt);
 675   Standard_Boolean isExtremaCircCircValid =  anExtC.IsDone() // Check if intersection is done
 676                                           && !anExtC.IsParallel() // Parallel case has already been considered
 677                                           && anExtC.NbExt() > 0; // Check that some solutions have been found
 678   if (!isExtremaCircCircValid)
 679     // not done
 680     return;
 681
 682   myDone = Standard_True;
 683
 684   // Few solutions
 685   Standard_Real aNbExt = anExtC.NbExt();
 686   // Find the minimal distance
 687   Standard_Real aMinSqDist = ::RealLast();
 688   for (Standard_Integer i = 1; i <= aNbExt; ++i)
 689   {
 690     Standard_Real aSqDist = anExtC.SquareDistance(i);
 691     if (aSqDist < aMinSqDist)
 692       aMinSqDist = aSqDist;
 693   }
 694
 695   // Collect all solutions close to the minimal one
 696   TColStd_ListOfInteger aSols;
 697   for (Standard_Integer i = 1; i <= aNbExt; ++i)
 698   {
 699     Standard_Real aDiff = anExtC.SquareDistance(i) - aMinSqDist;
 700     if (aDiff < Precision::SquareConfusion())
 701       aSols.Append(i);
 702   }
 703
 704   // Save all minimal solutions
 705   myNbExt = aSols.Extent();
 706
 707   mySqDist = new TColStd_HArray1OfReal(1, myNbExt);
 708   myPoint1 = new Extrema_HArray1OfPOnCurv(1, myNbExt);
 709   myPoint2 = new Extrema_HArray1OfPOnSurf(1, myNbExt);
 710
 711   TColStd_ListIteratorOfListOfInteger it(aSols);
 712   for (Standard_Integer iSol = 1; it.More(); it.Next(), ++iSol)
 713   {
 714     Extrema_POnCurv P1, P2;
 715     anExtC.Points(it.Value(), P1, P2);
 716
 717     // Compute parameters on sphere
 718     Standard_Real aU, aV;
 719     ElSLib::Parameters(S, P1.Value(), aU, aV);
 720     // Compute point on sphere
 721     gp_Pnt aPOnS = ElSLib::Value(aU, aV, S);
 722
 723     // Save solution
 724     myPoint1->SetValue(iSol, P1);
 725     myPoint2->SetValue(iSol, Extrema_POnSurf(aU, aV, aPOnS));
 726     mySqDist->SetValue(iSol, P1.Value().SquareDistance(aPOnS));
 727   }
 728 }
 729
 730 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Circ& C,
 731                                  const gp_Torus& S)
 732 {  Perform(C, S);}
 733
 734
 735
 736 //void Extrema_ExtElCS::Perform(const gp_Circ& C,
 737 //                            const gp_Torus& S)
 738 void Extrema_ExtElCS::Perform(const gp_Circ& ,
 739                               const gp_Torus& )
 740 {
 741   throw Standard_NotImplemented();
 742
 743 }
 744
 745 Extrema_ExtElCS::Extrema_ExtElCS(const gp_Hypr& C,
 746                                  const gp_Pln& S)
 747 {
 748   Perform(C, S);
 749 }
 750
 751
 752
 753 void Extrema_ExtElCS::Perform(const gp_Hypr& C,
 754                               const gp_Pln& S)
 755 {
 756   myDone = Standard_True;
 757   myIsPar = Standard_False;
 758   myNbExt = 0;
 759
 760   gp_Ax2 Pos = C.Position();
 761   gp_Dir NHypr = Pos.Direction();
 762   gp_Dir NPln = S.Axis().Direction();
 763
 764   if (NHypr.IsParallel(NPln, Precision::Angular())) {
 765
 766     mySqDist = new TColStd_HArray1OfReal(1, 1);
 767     mySqDist->SetValue(1, S.SquareDistance(C.Location()));
 768     myIsPar = Standard_True;
 769     myNbExt = 1;
 770   }
 771   else {
 772
 773     gp_Dir XDir = Pos.XDirection();
 774     gp_Dir YDir = Pos.YDirection();
 775
 776     Standard_Real A = C.MinorRadius()*(NPln.Dot(YDir));
 777     Standard_Real B = C.MajorRadius()*(NPln.Dot(XDir));
 778
 779     if(Abs(B) > Abs(A)) {
 780       Standard_Real T = -0.5 * Log((A+B)/(B-A));
 781       gp_Pnt Ph = ElCLib::HyperbolaValue(T, Pos, C.MajorRadius(), C.MinorRadius());
 782       Extrema_POnCurv PC(T, Ph);
 783       myPoint1 = new Extrema_HArray1OfPOnCurv(1,1);
 784       myPoint1->SetValue(1, PC);
 785
 786       mySqDist = new TColStd_HArray1OfReal(1, 1);
 787       mySqDist->SetValue(1, S.SquareDistance(Ph));
 788
 789       Standard_Real U, V;
 790       ElSLib::PlaneParameters(S.Position(), Ph, U, V);
 791       gp_Pnt Pp = ElSLib::PlaneValue(U, V, S.Position());
 792       Extrema_POnSurf PS(U, V, Pp);
 793       myPoint2 = new Extrema_HArray1OfPOnSurf(1,1);
 794       myPoint2->SetValue(1, PS);
 795
 796       myNbExt = 1;
 797     }
 798   }
 799 }
 800
 801
 802 Standard_Boolean Extrema_ExtElCS::IsDone() const
 803 {
 804   return myDone;
 805 }
 806
 807
 808 Standard_Integer Extrema_ExtElCS::NbExt() const
 809 {
 810   if (!IsDone()) throw StdFail_NotDone();
 811   return myNbExt;
 812 }
 813
 814 Standard_Real Extrema_ExtElCS::SquareDistance(const Standard_Integer N) const
 815 {
 816   if (N < 1 || N > NbExt())
 817   {
 818     throw Standard_OutOfRange();
 819   }
 820
 821   return mySqDist->Value(N);
 822 }
 823
 824
 825 void Extrema_ExtElCS::Points(const Standard_Integer N,
 826                              Extrema_POnCurv& P1,
 827                              Extrema_POnSurf& P2) const
 828 {
 829   if (N < 1 || N > NbExt())
 830   {
 831     throw Standard_OutOfRange();
 832   }
 833
 834   P1 = myPoint1->Value(N);
 835   P2 = myPoint2->Value(N);
 836 }
 837
 838
 839 Standard_Boolean Extrema_ExtElCS::IsParallel() const
 840 {
 841   if (!IsDone())
 842   {
 843     throw StdFail_NotDone();
 844   }
 845   return myIsPar;
 846 }