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