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
6 // This file is part of Open CASCADE Technology software library.
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
17 //----------------------------------------------------------------------
18 //-- Purposse: Geometric Intersection between two Natural Quadric
19 //-- If the intersection is not a conic,
20 //-- analytical methods must be called.
21 //----------------------------------------------------------------------
23 #define No_Standard_RangeError
24 #define No_Standard_OutOfRange
31 #include <gp_Circ.hxx>
32 #include <gp_Cone.hxx>
33 #include <gp_Cylinder.hxx>
35 #include <gp_Dir2d.hxx>
36 #include <gp_Elips.hxx>
37 #include <gp_Hypr.hxx>
39 #include <gp_Parab.hxx>
42 #include <gp_Pnt2d.hxx>
43 #include <gp_Sphere.hxx>
44 #include <gp_Torus.hxx>
46 #include <gp_Vec2d.hxx>
48 #include <IntAna_IntConicQuad.hxx>
49 #include <IntAna_QuadQuadGeo.hxx>
50 #include <math_DirectPolynomialRoots.hxx>
51 #include <Standard_DomainError.hxx>
52 #include <Standard_OutOfRange.hxx>
53 #include <StdFail_NotDone.hxx>
56 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
58 void RefineDir(gp_Dir& aDir);
60 //=======================================================================
62 //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
63 //=======================================================================
66 AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
68 void Distance(Standard_Real& dist,
69 Standard_Real& Param1,
70 Standard_Real& Param2);
72 gp_Pnt PtIntersect() {
75 Standard_Boolean Coplanar(void) {
78 Standard_Boolean Same(void) {
79 return theparallel && (thedistance<myEPSILON_DISTANCE);
81 Standard_Real Distance(void) {
84 Standard_Boolean Intersect(void) {
85 return (thecoplanar && (!theparallel));
87 Standard_Boolean Parallel(void) {
90 Standard_Boolean Normal(void) {
95 Standard_Real Det33(const Standard_Real a11,
96 const Standard_Real a12,
97 const Standard_Real a13,
98 const Standard_Real a21,
99 const Standard_Real a22,
100 const Standard_Real a23,
101 const Standard_Real a31,
102 const Standard_Real a32,
103 const Standard_Real a33) {
104 Standard_Real theReturn =
105 a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;
113 Standard_Real thedistance;
114 Standard_Boolean theparallel;
115 Standard_Boolean thecoplanar;
116 Standard_Boolean thenormal;
118 Standard_Real myEPSILON_DISTANCE;
119 Standard_Real myEPSILON_AXES_PARA;
122 //=======================================================================
123 //function : AxeOperator::AxeOperator
125 //=======================================================================
126 AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2)
128 myEPSILON_DISTANCE=1.0e-14;
129 myEPSILON_AXES_PARA=Precision::Angular();
132 //---------------------------------------------------------------------
133 gp_Dir V1=Axe1.Direction();
134 gp_Dir V2=Axe2.Direction();
135 gp_Pnt P1=Axe1.Location();
136 gp_Pnt P2=Axe2.Location();
140 thecoplanar= Standard_False;
141 thenormal = Standard_False;
143 //--- check if the two axis are parallel
144 theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);
145 //--- Distance between the two axis
146 gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
149 thedistance = L1.Distance(A2.Location());
152 thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
155 //--- check if Axis are Coplanar
157 if(thedistance<myEPSILON_DISTANCE) {
158 D33=Det33(V1.X(),V1.Y(),V1.Z()
159 ,V2.X(),V2.Y(),V2.Z()
160 ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
161 if(Abs(D33)<=myEPSILON_DISTANCE) {
162 thecoplanar=Standard_True;
166 thecoplanar=Standard_True;
167 thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
169 //--- check if the two axis are concurrent
170 if(thecoplanar && (!theparallel)) {
171 Standard_Real smx=P2.X() - P1.X();
172 Standard_Real smy=P2.Y() - P1.Y();
173 Standard_Real smz=P2.Z() - P1.Z();
174 Standard_Real Det1,Det2,Det3,A;
175 Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
176 Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
177 Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
179 if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
180 A=(smy * V2.X() - smx * V2.Y())/Det1;
183 && ((Abs(Det2) >= Abs(Det1))
184 &&(Abs(Det2) >= Abs(Det3)))) {
185 A=(smz * V2.Y() - smy * V2.Z())/Det2;
188 A=(smz * V2.X() - smx * V2.Z())/Det3;
190 ptintersect.SetCoord( P1.X() + A * V1.X()
192 ,P1.Z() + A * V1.Z());
195 ptintersect.SetCoord(0,0,0); //-- Pour eviter des FPE
198 //=======================================================================
199 //function : Distance
201 //=======================================================================
202 void AxeOperator::Distance(Standard_Real& dist,
203 Standard_Real& Param1,
204 Standard_Real& Param2)
206 gp_Vec O1O2(Axe1.Location(),Axe2.Location());
207 gp_Dir U1 = Axe1.Direction(); //-- juste pour voir.
208 gp_Dir U2 = Axe2.Direction();
210 gp_Dir N = U1.Crossed(U2);
211 Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
213 U1.Z(),U2.Z(),N.Z());
215 dist = Det33(U1.X(),U2.X(),O1O2.X(),
216 U1.Y(),U2.Y(),O1O2.Y(),
217 U1.Z(),U2.Z(),O1O2.Z()) / D;
218 Param1 = Det33(O1O2.X(),U2.X(),N.X(),
219 O1O2.Y(),U2.Y(),N.Y(),
220 O1O2.Z(),U2.Z(),N.Z()) / (-D);
221 //------------------------------------------------------------
222 //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
223 //-- soit : Segment perpendiculaire : O1+P1 D1
225 Param2 = Det33(U1.X(),O1O2.X(),N.X(),
226 U1.Y(),O1O2.Y(),N.Y(),
227 U1.Z(),O1O2.Z(),N.Z()) / (D);
230 //=======================================================================
231 //function : DirToAx2
232 //purpose : returns a gp_Ax2 where D is the main direction
233 //=======================================================================
234 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
236 Standard_Real x=D.X(); Standard_Real ax=Abs(x);
237 Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
238 Standard_Real z=D.Z(); Standard_Real az=Abs(z);
239 if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
240 return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
242 else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
243 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
246 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
249 //=======================================================================
250 //function : IntAna_QuadQuadGeo
251 //purpose : Empty constructor
252 //=======================================================================
253 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
254 : done(Standard_False),
256 typeres(IntAna_Empty),
267 myCommonGen(Standard_False),
272 //=======================================================================
273 //function : InitTolerances
275 //=======================================================================
276 void IntAna_QuadQuadGeo::InitTolerances()
278 myEPSILON_DISTANCE = 1.0e-14;
279 myEPSILON_ANGLE_CONE = Precision::Angular();
280 myEPSILON_MINI_CIRCLE_RADIUS = 0.01*Precision::Confusion();
281 myEPSILON_CYLINDER_DELTA_RADIUS = 1.0e-13;
282 myEPSILON_CYLINDER_DELTA_DISTANCE= Precision::Confusion();
283 myEPSILON_AXES_PARA = Precision::Angular();
285 //=======================================================================
286 //function : IntAna_QuadQuadGeo
288 //=======================================================================
289 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1,
291 const Standard_Real TolAng,
292 const Standard_Real Tol)
293 : done(Standard_False),
295 typeres(IntAna_Empty),
306 myCommonGen(Standard_False),
310 Perform(P1,P2,TolAng,Tol);
312 //=======================================================================
315 //=======================================================================
316 void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
318 const Standard_Real TolAng,
319 const Standard_Real Tol)
321 Standard_Real A1, B1, C1, D1, A2, B2, C2, D2, dist1, dist2, aMVD;
326 P1.Coefficients(A1,B1,C1,D1);
327 P2.Coefficients(A2,B2,C2,D2);
329 gp_Vec aVN1(A1,B1,C1);
330 gp_Vec aVN2(A2,B2,C2);
331 gp_Vec vd(aVN1.Crossed(aVN2));
333 const gp_Pnt& aLocP1=P1.Location();
334 const gp_Pnt& aLocP2=P2.Location();
336 dist1=A2*aLocP1.X() + B2*aLocP1.Y() + C2*aLocP1.Z() + D2;
337 dist2=A1*aLocP2.X() + B1*aLocP2.Y() + C1*aLocP2.Z() + D1;
341 // normalles are collinear - planes are same or parallel
342 typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same
346 Standard_Real denom, denom2, ddenom, par1, par2;
347 Standard_Real X1, Y1, Z1, X2, Y2, Z2, aEps;
350 denom=A1*A2 + B1*B2 + C1*C2;
351 denom2 = denom*denom;
352 ddenom = 1. - denom2;
354 denom = ( Abs(ddenom) <= aEps ) ? aEps : ddenom;
359 gp_Vec inter1(aVN1.Crossed(vd));
360 gp_Vec inter2(aVN2.Crossed(vd));
362 X1=aLocP1.X() + par1*inter1.X();
363 Y1=aLocP1.Y() + par1*inter1.Y();
364 Z1=aLocP1.Z() + par1*inter1.Z();
365 X2=aLocP2.X() + par2*inter2.X();
366 Y2=aLocP2.Y() + par2*inter2.Y();
367 Z2=aLocP2.Z() + par2*inter2.Z();
369 pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
371 typeres = IntAna_Line;
374 //-------------------------------------------------------
375 // When the value of the angle between the planes is small
376 // the origin of intersection line is computed with error
377 // [ ~0.0001 ] that can not br considered as small one
379 // for {A~=2.e-6, dist1=4.2e-5, dist2==1.e-4} =>
380 // {denom=3.4e-12, par1=12550297.6, par2=32605552.9, etc}
382 // the origin should be refined if it is possible
384 Standard_Real aTreshAng, aTreshDist;
386 aTreshAng=2.e-6; // 1.e-4 deg
389 if (aMVD < aTreshAng) {
390 Standard_Real aDist1, aDist2;
392 aDist1=A1*pt1.X() + B1*pt1.Y() + C1*pt1.Z() + D1;
393 aDist2=A2*pt1.X() + B2*pt1.Y() + C2*pt1.Z() + D2;
395 if (fabs(aDist1)>aTreshDist || fabs(aDist2)>aTreshDist) {
396 Standard_Boolean bIsDone, bIsParallel;
397 IntAna_IntConicQuad aICQ;
401 gp_Lin aL1(pt1, aDN1);
403 aICQ.Perform(aL1, P1, TolAng, Tol);
404 bIsDone=aICQ.IsDone();
409 const gp_Pnt& aPnt1=aICQ.Point(1);
410 //----------------------------------
412 gp_Dir aDL2(dir1.Crossed(aDN1));
413 gp_Lin aL2(aPnt1, aDL2);
415 aICQ.Perform(aL2, P2, TolAng, Tol);
416 bIsDone=aICQ.IsDone();
421 bIsParallel=aICQ.IsParallel();
426 const gp_Pnt& aPnt2=aICQ.Point(1);
434 //=======================================================================
435 //function : IntAna_QuadQuadGeo
436 //purpose : Pln Cylinder
437 //=======================================================================
438 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
439 ,const gp_Cylinder& Cl
440 ,const Standard_Real Tolang
441 ,const Standard_Real Tol
442 ,const Standard_Real H)
443 : done(Standard_False),
445 typeres(IntAna_Empty),
456 myCommonGen(Standard_False),
460 Perform(P,Cl,Tolang,Tol,H);
462 //=======================================================================
465 //=======================================================================
466 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
467 ,const gp_Cylinder& Cl
468 ,const Standard_Real Tolang
469 ,const Standard_Real Tol
470 ,const Standard_Real H)
472 done = Standard_False;
473 Standard_Real dist,radius;
474 Standard_Real A,B,C,D;
476 Standard_Real sint,cost,h;
477 gp_XYZ axex,axey,omega;
481 radius = Cl.Radius();
483 gp_Lin axec(Cl.Axis());
484 gp_XYZ normp(P.Axis().Direction().XYZ());
486 P.Coefficients(A,B,C,D);
487 axec.Location().Coord(X,Y,Z);
488 // la distance axe/plan est evaluee a l origine de l axe.
489 dist = A*X + B*Y + C*Z + D;
491 Standard_Real tolang = Tolang;
492 Standard_Boolean newparams = Standard_False;
494 gp_Vec ldv( axec.Direction() );
496 Standard_Real dA = Abs( ldv.Angle( npv ) );
499 Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
500 Standard_Real dangle = Abs( dang );
501 if( dangle > Tolang )
503 Standard_Real sinda = Abs( Sin( dangle ) );
504 Standard_Real dif = Abs( sinda - Tol );
508 newparams = Standard_True;
514 IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
516 if (inter.IsParallel()) {
517 // Le resultat de l intersection Plan-Cylindre est de type droite.
518 // il y a 1 ou 2 droites
520 typeres = IntAna_Line;
521 omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
523 if (Abs(Abs(dist)-radius) < Tol)
530 gp_XYZ omegaXYZ(X,Y,Z);
531 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
532 Standard_Real Xt,Yt,Zt,distt;
533 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
534 distt = A*Xt + B*Yt + C*Zt + D;
535 gp_XYZ omega1(omegaXYZtrnsl.X()-distt*A,
536 omegaXYZtrnsl.Y()-distt*B,
537 omegaXYZtrnsl.Z()-distt*C );
539 ppt1.SetXYZ( omega1 );
540 gp_Vec vv1(pt1,ppt1);
545 dir1 = axec.Direction();
547 else if (Abs(dist) < radius)
550 h = Sqrt(radius*radius - dist*dist);
551 axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
553 pt1.SetXYZ(omega - h*axey);
554 pt2.SetXYZ(omega + h*axey);
558 gp_XYZ omegaXYZ(X,Y,Z);
559 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
560 Standard_Real Xt,Yt,Zt,distt,ht;
561 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
562 distt = A*Xt + B*Yt + C*Zt + D;
563 // ht = Sqrt(radius*radius - distt*distt);
564 Standard_Real anSqrtArg = radius*radius - distt*distt;
565 ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
567 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A,
568 omegaXYZtrnsl.Y()-distt*B,
569 omegaXYZtrnsl.Z()-distt*C );
571 ppt1.SetXYZ( omega1 - ht*axey);
572 ppt2.SetXYZ( omega1 + ht*axey);
573 gp_Vec vv1(pt1,ppt1);
574 gp_Vec vv2(pt2,ppt2);
582 dir1 = axec.Direction();
583 dir2 = axec.Direction();
588 // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
589 // et ne pas etre seulement supprime...
592 typeres = IntAna_Empty;
595 else { // Il y a un point d intersection. C est le centre du cercle
596 // ou de l ellipse solution.
599 axey = normp.Crossed(axec.Direction().XYZ());
600 sint = axey.Modulus();
602 pt1 = inter.Point(1);
604 if (sint < Tol/radius) {
606 // on construit un cercle avec comme axes X et Y ceux du cylindre
607 typeres = IntAna_Circle;
609 dir1 = axec.Direction(); // axe Z
610 dir2 = Cl.Position().XDirection();
615 // on construit un ellipse
616 typeres = IntAna_Ellipse;
617 cost = Abs(axec.Direction().XYZ().Dot(normp));
618 axex = axey.Crossed(normp);
620 dir1.SetXYZ(normp); //Modif ds ce bloc
623 param1 = radius/cost;
628 done = Standard_True;
630 //=======================================================================
631 //function : IntAna_QuadQuadGeo
633 //=======================================================================
634 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
636 const Standard_Real Tolang,
637 const Standard_Real Tol)
638 : done(Standard_False),
640 typeres(IntAna_Empty),
651 myCommonGen(Standard_False),
655 Perform(P,Co,Tolang,Tol);
657 //=======================================================================
660 //=======================================================================
661 void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
663 const Standard_Real Tolang,
664 const Standard_Real Tol)
667 done = Standard_False;
670 Standard_Real A,B,C,D;
672 Standard_Real dist,sint,cost,sina,cosa,angl,costa;
676 gp_Lin axec(Co.Axis());
677 P.Coefficients(A,B,C,D);
678 gp_Pnt apex(Co.Apex());
681 dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
683 gp_XYZ normp = P.Axis().Direction().XYZ();
684 if(P.Direct()==Standard_False) { //-- lbr le 14 jan 97
688 axey = normp.Crossed(Co.Axis().Direction().XYZ());
689 axex = axey.Crossed(normp);
692 angl = Co.SemiAngle();
695 sina = Abs(Sin(angl));
698 // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
700 sint = axey.Modulus();
701 cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
703 // Le calcul de costa permet de determiner si le plan contient
704 // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
706 costa = cost*cosa - sint*sina; // sin((PI/2 -t)-angl)=cos(t+angl)
707 // avec t ramene entre 0 et pi/2.
709 if (Abs(dist) < Tol) {
710 // on considere que le plan contient le sommet du cone.
711 // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
712 // selon l inclinaison du plan.
714 if (Abs(costa) < Tolang) { // plan parallele a la generatrice
715 typeres = IntAna_Line;
717 gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
718 // point sur l axe du cone cote z positif
720 dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
721 ptonaxe = ptonaxe - dist*normp;
723 dir1.SetXYZ(ptonaxe - pt1.XYZ());
725 else if (cost < sina) { // plan "interieur" au cone
726 typeres = IntAna_Line;
730 dh = Sqrt(sina*sina-cost*cost)/cosa;
731 dir1.SetXYZ(axex + dh*axey);
732 dir2.SetXYZ(axex - dh*axey);
734 else { // plan "exterieur" au cone
735 typeres = IntAna_Point;
741 // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
742 // l inclinaison du plan.
743 Standard_Real deltacenter, distance;
746 // Le plan contient la direction de l axe du cone. La solution est
748 typeres = IntAna_Hyperbola;
750 pt1.SetXYZ(apex.XYZ()-dist*normp);
754 param1 = param2 = Abs(dist/Tan(angl));
755 param1bis = param2bis = Abs(dist);
759 IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
761 gp_Pnt center(inter.Point(1));
763 // En fonction de la position de l intersection par rapport au sommet
764 // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
765 // sur axec est negatif (voir definition du cone)
767 distance = apex.Distance(center);
769 if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
774 if (Abs(costa) < Tolang) { // plan parallele a une generatrice
775 typeres = IntAna_Parabola;
777 deltacenter = distance/2./cosa;
779 pt1.SetXYZ(center.XYZ()-deltacenter*axex);
782 param1 = deltacenter*sina*sina;
784 else if (sint < Tolang) { // plan perpendiculaire a l axe
785 typeres = IntAna_Circle;
788 dir1 = Co.Position().Direction();
789 dir2 = Co.Position().XDirection();
790 param1 = apex.Distance(center)*Abs(Tan(angl));
792 else if (cost < sina ) {
793 typeres = IntAna_Hyperbola;
797 deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
798 pt1.SetXYZ(center.XYZ() - deltacenter*axex);
802 param1 = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
803 param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
806 else { // on a alors cost > sina
807 typeres = IntAna_Ellipse;
809 Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
810 deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
812 pt1.SetXYZ(center.XYZ() + deltacenter*axex);
816 param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
821 //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
822 //-- des hyperboles trop bizarres
823 //-- On retourne False -> Traitement par biparametree
824 static Standard_Real EllipseLimit = 1.0E+9; //OCC513(apo) 1000000
825 static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
826 if(typeres==IntAna_Ellipse && nbint>=1) {
827 if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit) {
832 if(typeres==IntAna_Hyperbola && nbint>=2) {
833 if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit) {
834 done = Standard_False;
838 if(typeres==IntAna_Hyperbola && nbint>=1) {
839 if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit) {
845 done = Standard_True;
848 //=======================================================================
849 //function : IntAna_QuadQuadGeo
850 //purpose : Pln Sphere
851 //=======================================================================
852 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
854 : done(Standard_False),
856 typeres(IntAna_Empty),
867 myCommonGen(Standard_False),
873 //=======================================================================
876 //=======================================================================
877 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
881 done = Standard_False;
882 Standard_Real A,B,C,D,dist, radius;
886 // debug JAG : on met typeres = IntAna_Empty par defaut...
887 typeres = IntAna_Empty;
889 P.Coefficients(A,B,C,D);
890 S.Location().Coord(X,Y,Z);
893 dist = A * X + B * Y + C * Z + D;
895 if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
896 // on a une seule solution : le point projection du centre de la sphere
899 typeres = IntAna_Point;
900 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
902 else if (Abs(dist) < radius) {
903 // on a un cercle solution
905 typeres = IntAna_Circle;
906 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
907 dir1 = P.Axis().Direction();
908 if(P.Direct()==Standard_False) dir1.Reverse();
909 dir2 = P.Position().XDirection();
910 param1 = Sqrt(radius*radius - dist*dist);
912 param2bis=0.0; //-- pour eviter param2bis not used ....
913 done = Standard_True;
916 //=======================================================================
917 //function : IntAna_QuadQuadGeo
918 //purpose : Cylinder - Cylinder
919 //=======================================================================
920 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
921 const gp_Cylinder& Cyl2,
922 const Standard_Real Tol)
923 : done(Standard_False),
925 typeres(IntAna_Empty),
936 myCommonGen(Standard_False),
940 Perform(Cyl1,Cyl2,Tol);
942 //=======================================================================
945 //=======================================================================
946 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
947 const gp_Cylinder& Cyl2,
948 const Standard_Real Tol)
951 //---------------------------- Parallel axes -------------------------
952 AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
953 Standard_Real R1=Cyl1.Radius();
954 Standard_Real R2=Cyl2.Radius();
955 Standard_Real RmR, RmR_Relative;
956 RmR=(R1>R2)? (R1-R2) : (R2-R1);
959 Rmax=(R1>R2)? R1 : R2;
960 RmR_Relative=RmR/Rmax;
963 Standard_Real DistA1A2=A1A2.Distance();
975 typeres=IntAna_Empty;
979 { //-- DistA1A2 > Tol
980 gp_Pnt P1=Cyl1.Location();
981 gp_Pnt P2t=Cyl2.Location();
983 //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
984 gp_Dir DirCyl = Cyl1.Position().Direction();
985 Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
987 //P2 is a projection the location of the 2nd cylinder on the base
988 //of the 1st cylinder
989 P2.SetCoord(P2t.X() - ProjP2OnDirCyl1*DirCyl.X(),
990 P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y(),
991 P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
993 Standard_Real R1pR2=R1+R2;
994 if(DistA1A2>(R1pR2+Tol))
996 typeres=IntAna_Empty;
999 else if((R1pR2 - DistA1A2) <= RealSmall())
1001 //-- 1 Tangent line -------------------------------------OK
1002 typeres=IntAna_Line;
1006 Standard_Real R1_R1pR2=R1/R1pR2;
1007 pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X()),
1008 P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y()),
1009 P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
1011 else if(DistA1A2>RmR)
1013 //-- 2 lines ---------------------------------------------OK
1014 typeres=IntAna_Line;
1019 const Standard_Real aR1R1 = R1*R1;
1032 Two cylinders have axes collinear. Therefore, problem can be reformulated as
1033 to find intersection point of two circles (the bases of the cylinders) on
1034 the plane: 1st circle has center P1 and radius R1 (the radius of the
1035 1st cylinder) and 2nd circle has center P2 and radius R2 (the radius of the
1036 2nd cylinder). The plane is the base of the 1st cylinder. Points A and B
1037 are intersection point of these circles. Distance P1P2 is equal to DistA1A2.
1038 O1 is the intersection point of P1P2 and AB segments.
1040 At that, if distance AB < Tol we consider that the circles are tangent and
1041 has only one intersection point.
1043 AB = 2*R1*sin(angle AP1P2).
1045 AB^2 < Tol^2 => 4*R1*R1*sin(angle AP1P2)^2 < Tol*Tol.
1049 //Cosine and Square of Sine of the A-P1-P2 angle
1050 const Standard_Real aCos = 0.5*(aR1R1-R2*R2+DistA1A2*DistA1A2)/(R1*DistA1A2);
1051 const Standard_Real aSin2 = 1-aCos*aCos;
1053 const Standard_Boolean isTangent =((4.0*aR1R1*aSin2) < Tol*Tol);
1055 //Normalized vector P1P2
1056 const gp_Vec DirA1A2((P2.XYZ() - P1.XYZ())/DistA1A2);
1060 //Intercept the segment from P1 point along P1P2 direction
1061 //and having |P1O1| length
1063 pt1.SetXYZ(P1.XYZ() + DirA1A2.XYZ()*R1*aCos);
1067 //Sine of the A-P1-P2 angle (if aSin2 < 0 then isTangent == TRUE =>
1068 //go to another branch)
1069 const Standard_Real aSin = sqrt(aSin2);
1071 //1. Rotate P1P2 to the angle A-P1-P2 relative to P1
1072 //(clockwise and anticlockwise for getting
1073 //two intersection points).
1074 //2. Intercept the segment from P1 along direction,
1075 //determined in the preview paragraph and having R1 length
1076 const gp_Dir &aXDir = Cyl1.Position().XDirection(),
1077 &aYDir = Cyl1.Position().YDirection();
1078 const gp_Vec aR1Xdir = R1*aXDir.XYZ(),
1079 aR1Ydir = R1*aYDir.XYZ();
1081 //Source 2D-coordinates of the P1P2 vector normalized
1082 //in coordinate system, based on the X- and Y-directions
1083 //of the 1st cylinder in the plane of the 1st cylinder base
1084 //(P1 is the origin of the coordinate system).
1085 const Standard_Real aDx = DirA1A2.Dot(aXDir),
1086 aDy = DirA1A2.Dot(aYDir);
1088 //New coordinate (after rotation) of the P1P2 vector normalized.
1089 Standard_Real aNewDx = aDx*aCos - aDy*aSin,
1090 aNewDy = aDy*aCos + aDx*aSin;
1091 pt1.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
1093 aNewDx = aDx*aCos + aDy*aSin;
1094 aNewDy = aDy*aCos - aDx*aSin;
1095 pt2.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
1098 else if(DistA1A2>(RmR-Tol))
1100 //-- 1 Tangent ------------------------------------------OK
1101 typeres=IntAna_Line;
1104 Standard_Real R1_RmR=R1/RmR;
1109 pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X()),
1110 P1.Y() + R1_RmR * (P2.Y()-P1.Y()),
1111 P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
1115 typeres=IntAna_Empty;
1119 else { //-- No Parallel Axis ---------------------------------OK
1120 if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
1121 && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE))
1123 //-- PI/2 between the two axis and Intersection
1124 //-- and identical radius
1125 typeres=IntAna_Ellipse;
1127 gp_Dir DirCyl1=Cyl1.Position().Direction();
1128 gp_Dir DirCyl2=Cyl2.Position().Direction();
1129 pt1=pt2=A1A2.PtIntersect();
1131 Standard_Real A=DirCyl1.Angle(DirCyl2);
1133 B=Abs(Sin(0.5*(M_PI-A)));
1136 if(A==0.0 || B==0.0)
1138 typeres=IntAna_Same;
1142 gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
1143 dir1 = gp_Dir(dircyl1.Added(dircyl2));
1144 dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
1146 param2 = Cyl1.Radius() / A;
1147 param1 = Cyl1.Radius() / B;
1148 param2bis= param1bis = Cyl1.Radius();
1149 if(param1 < param1bis)
1156 if(param2 < param2bis)
1165 if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol)
1167 typeres = IntAna_Point;
1168 Standard_Real d,p1,p2;
1170 gp_Dir D1 = Cyl1.Axis().Direction();
1171 gp_Dir D2 = Cyl2.Axis().Direction();
1172 A1A2.Distance(d,p1,p2);
1173 gp_Pnt P = Cyl1.Axis().Location();
1174 gp_Pnt P1(P.X() - p1*D1.X(),
1178 P = Cyl2.Axis().Location();
1180 gp_Pnt P2(P.X() - p2*D2.X(),
1188 pt1.SetCoord(P1.X() + p1*D1.X(),
1190 P1.Z() + p1*D1.Z());
1195 typeres=IntAna_NoGeometricSolution;
1200 //=======================================================================
1201 //function : IntAna_QuadQuadGeo
1202 //purpose : Cylinder - Cone
1203 //=======================================================================
1204 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1206 const Standard_Real Tol)
1207 : done(Standard_False),
1209 typeres(IntAna_Empty),
1220 myCommonGen(Standard_False),
1224 Perform(Cyl,Con,Tol);
1226 //=======================================================================
1227 //function : Perform
1229 //=======================================================================
1230 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1232 const Standard_Real )
1235 AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1237 gp_Pnt Pt=Con.Apex();
1238 Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1239 gp_Dir dir=Cyl.Position().Direction();
1240 pt1.SetCoord( Pt.X() + dist*dir.X()
1241 ,Pt.Y() + dist*dir.Y()
1242 ,Pt.Z() + dist*dir.Z());
1243 pt2.SetCoord( Pt.X() - dist*dir.X()
1244 ,Pt.Y() - dist*dir.Y()
1245 ,Pt.Z() - dist*dir.Z());
1247 param1=param2=Cyl.Radius();
1249 typeres=IntAna_Circle;
1253 typeres=IntAna_NoGeometricSolution;
1256 //=======================================================================
1258 //purpose : Cylinder - Sphere
1259 //=======================================================================
1260 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1261 const gp_Sphere& Sph,
1262 const Standard_Real Tol)
1263 : done(Standard_False),
1265 typeres(IntAna_Empty),
1276 myCommonGen(Standard_False),
1280 Perform(Cyl,Sph,Tol);
1282 //=======================================================================
1283 //function : Perform
1285 //=======================================================================
1286 void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1287 ,const gp_Sphere& Sph
1288 ,const Standard_Real)
1291 gp_Pnt Pt=Sph.Location();
1292 AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1293 if((A1A2.Intersect() && Pt.Distance(A1A2.PtIntersect())==0.0 )
1295 if(Sph.Radius() < Cyl.Radius()) {
1296 typeres = IntAna_Empty;
1299 Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1300 gp_Dir dir=Cyl.Position().Direction();
1302 typeres=IntAna_Circle;
1303 pt1.SetCoord( Pt.X() + dist*dir.X()
1304 ,Pt.Y() + dist*dir.Y()
1305 ,Pt.Z() + dist*dir.Z());
1307 param1 = Cyl.Radius();
1308 if(dist>RealEpsilon()) {
1309 pt2.SetCoord( Pt.X() - dist*dir.X()
1310 ,Pt.Y() - dist*dir.Y()
1311 ,Pt.Z() - dist*dir.Z());
1312 param2=Cyl.Radius();
1318 typeres=IntAna_NoGeometricSolution;
1322 //=======================================================================
1323 //function : IntAna_QuadQuadGeo
1324 //purpose : Cone - Cone
1325 //=======================================================================
1326 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1327 const gp_Cone& Con2,
1328 const Standard_Real Tol)
1329 : done(Standard_False),
1331 typeres(IntAna_Empty),
1342 myCommonGen(Standard_False),
1346 Perform(Con1,Con2,Tol);
1349 //=======================================================================
1350 //function : Perform
1352 //=======================================================================
1353 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1354 const gp_Cone& Con2,
1355 const Standard_Real Tol)
1359 Standard_Real tg1, tg2, aDA1A2, aTol2;
1360 gp_Pnt aPApex1, aPApex2;
1362 Standard_Real TOL_APEX_CONF = 1.e-10;
1365 tg1=Tan(Con1.SemiAngle());
1366 tg2=Tan(Con2.SemiAngle());
1368 if((tg1 * tg2) < 0.) {
1373 aPApex1=Con1.Apex();
1374 aPApex2=Con2.Apex();
1375 aDA1A2=aPApex1.SquareDistance(aPApex2);
1377 AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1383 gp_Pnt P=Con1.Apex();
1384 gp_Dir D=Con1.Position().Direction();
1385 Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1387 if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
1388 if (fabs(d) < TOL_APEX_CONF) {
1389 typeres = IntAna_Point;
1394 x=(d*tg2)/(tg1+tg2);
1395 pt1.SetCoord( P.X() + x*D.X()
1400 x=(d*tg2)/(tg2-tg1);
1401 pt2.SetCoord( P.X() + x*D.X()
1407 typeres=IntAna_Circle;
1410 if (fabs(d) < TOL_APEX_CONF) {
1411 typeres=IntAna_Same;
1414 typeres=IntAna_Circle;
1417 pt1.SetCoord( P.X() + x*D.X()
1420 param1 = Abs(x * tg1);
1424 } //-- fin A1A2.Same
1426 else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1427 //-- voir AnVer12mai98
1428 Standard_Real DistA1A2=A1A2.Distance();
1429 gp_Dir DA1=Con1.Position().Direction();
1430 gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1431 gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1432 Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1434 gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
1435 O1O2n.Y()-O1O2_DA1*DA1.Y(),
1436 O1O2n.Z()-O1O2_DA1*DA1.Z());
1437 gp_Dir DB1=gp_Dir(O1_Proj_A2);
1439 Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1440 Standard_Real ABSTG1 = Abs(tg1);
1441 Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1442 Standard_Real X1 = X2+yO1O2;
1444 gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1445 Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1446 Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1448 gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1449 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1450 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1451 gp_Vec P1MO1O2(P1,MO1O2);
1453 gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1454 gp_Dir OrthoPln = DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1456 IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1457 if(INTER_QUAD_PLN.IsDone()) {
1458 switch(INTER_QUAD_PLN.TypeInter()) {
1459 case IntAna_Ellipse: {
1460 typeres=IntAna_Ellipse;
1461 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1463 dir1 = E.Position().Direction();
1464 dir2 = E.Position().XDirection();
1465 param1 = E.MajorRadius();
1466 param1bis = E.MinorRadius();
1470 case IntAna_Circle: {
1471 typeres=IntAna_Circle;
1472 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1474 dir1 = C.Position().XDirection();
1475 dir2 = C.Position().YDirection();
1476 param1 = C.Radius();
1480 case IntAna_Hyperbola: {
1481 typeres=IntAna_Hyperbola;
1482 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1483 pt1 = pt2 = H.Location();
1484 dir1 = H.Position().Direction();
1485 dir2 = H.Position().XDirection();
1486 param1 = param2 = H.MajorRadius();
1487 param1bis = param2bis = H.MinorRadius();
1492 typeres=IntAna_Line;
1493 gp_Lin H=INTER_QUAD_PLN.Line(1);
1494 pt1 = pt2 = H.Location();
1495 dir1 = dir2 = H.Position().Direction();
1496 param1 = param2 = 0.0;
1497 param1bis = param2bis = 0.0;
1502 typeres=IntAna_NoGeometricSolution;
1505 }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1507 else if (aDA1A2<aTol2) {
1509 // When apices are coinsided there can be 3 possible cases
1510 // 3.1 - empty solution (iRet=0)
1511 // 3.2 - one line when cone1 touches cone2 (iRet=1)
1512 // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1514 Standard_Integer iRet;
1515 Standard_Real aGamma, aBeta1, aBeta2;
1516 Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1517 Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1518 gp_Pnt2d aP0, aPA1, aP1, aPA2;
1522 // Preliminary analysis. Determination of iRet
1527 aPA1.SetCoord(aD1, 0.);
1528 aP0.SetCoord(0., 0.);
1532 aGamma=aAx1.Angle(aAx2);
1533 if (aGamma>aHalfPI){
1536 aCosGamma=Cos(aGamma);
1537 aSinGamma=Sin(aGamma);
1539 aBeta1=Con1.SemiAngle();
1540 aTgBeta1=Tan(aBeta1);
1541 aTgBeta1=Abs(aTgBeta1);
1543 aBeta2=Con2.SemiAngle();
1544 aTgBeta2=Tan(aBeta2);
1545 aTgBeta2=Abs(aTgBeta2);
1548 aP1.SetCoord(aD1, aR1);
1551 aVAx2.SetCoord(aCosGamma, aSinGamma);
1552 gp_Dir2d aDAx2(aVAx2);
1553 gp_Lin2d aLAx2(aP0, aDAx2);
1555 gp_Vec2d aV(aP0, aP1);
1557 aPA2=aP0.Translated(aDx*aDAx2);
1560 aDx=aPA2.Distance(aP0);
1564 aRD2=aPA2.Distance(aP1);
1566 if (aRD2>(aR2+Tol)) {
1568 typeres=IntAna_Empty; //nothing
1572 iRet=1; //touch case => 1 line
1573 if (aRD2<(aR2-Tol)) {
1574 iRet=2;//intersection => couple of lines
1577 // Finding the solution in 3D
1580 gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1581 gp_Dir aD3Ax1, aD3Ax2;
1583 IntAna_QuadQuadGeo aIntr;
1585 aQApex1=Con1.Apex();
1586 aD3Ax1=aAx1.Direction();
1587 aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1588 aQApex1.Y()+aD1*aD3Ax1.Y(),
1589 aQApex1.Z()+aD1*aD3Ax1.Z());
1591 aDx=aD3Ax1.Dot(aAx2.Direction());
1595 aD3Ax2=aAx2.Direction();
1597 aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1599 aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1600 aQApex1.Y()+aD2*aD3Ax2.Y(),
1601 aQApex1.Z()+aD2*aD3Ax2.Z());
1603 gp_Pln aPln1(aQA1, aD3Ax1);
1604 gp_Pln aPln2(aQA2, aD3Ax2);
1606 aIntr.Perform(aPln1, aPln2, Tol, Tol);
1607 if (!aIntr.IsDone() || 0 == aIntr.NbSolutions()) {
1608 iRet=-1; // just in case. it must not be so
1609 typeres=IntAna_NoGeometricSolution;
1614 const gp_Dir& aDLin=aLin.Direction();
1615 gp_Vec aVLin(aDLin);
1616 gp_Pnt aOrig=aLin.Location();
1617 gp_Vec aVr(aQA1, aOrig);
1619 aQX=aOrig.Translated(aDx*aVLin);
1623 typeres=IntAna_Line;
1634 gp_Vec aVX(aQApex1, aQX);
1641 aDa=aQA1.Distance(aQX);
1642 aDx=sqrt(aR1*aR1-aDa*aDa);
1643 aQX1=aQX.Translated(aDx*aVLin);
1644 aQX2=aQX.Translated(-aDx*aVLin);
1648 gp_Vec aVX1(aQApex1, aQX1);
1650 gp_Vec aVX2(aQApex1, aQX2);
1653 } //else if (aDA1A2<aTol2) {
1654 //Case when cones have common generatrix
1655 else if(A1A2.Intersect()) {
1656 //Check if apex of one cone belongs another one
1657 Standard_Real u, v, tol2 = Tol*Tol;
1658 ElSLib::Parameters(Con2, aPApex1, u, v);
1659 gp_Pnt p = ElSLib::Value(u, v, Con2);
1660 if(aPApex1.SquareDistance(p) > tol2) {
1661 typeres=IntAna_NoGeometricSolution;
1665 ElSLib::Parameters(Con1, aPApex2, u, v);
1666 p = ElSLib::Value(u, v, Con1);
1667 if(aPApex2.SquareDistance(p) > tol2) {
1668 typeres=IntAna_NoGeometricSolution;
1672 //Cones have a common generatrix passing through apexes
1673 myCommonGen = Standard_True;
1675 //common generatrix of cones
1676 gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1678 //Intersection point of axes
1679 gp_Pnt aPAxeInt = A1A2.PtIntersect();
1681 //Characteristic point of intersection curve
1682 u = ElCLib::Parameter(aGen, aPAxeInt);
1683 myPChar = ElCLib::Value(u, aGen);
1686 //Other generatrixes of cones laying in maximal plane
1687 gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1688 gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
1690 //Intersection point of generatrixes
1691 gp_Dir aN; //solution plane normal
1692 gp_Dir aD1 = aGen1.Direction();
1694 gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1696 if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1697 aN = aD1.Crossed(aD2);
1699 else if(aGen1.SquareDistance(aGen2) > tol2) {
1700 //Something wrong ???
1701 typeres=IntAna_NoGeometricSolution;
1705 gp_Dir D1 = aGen1.Position().Direction();
1706 gp_Dir D2 = aGen2.Position().Direction();
1707 gp_Pnt O1 = aGen1.Location();
1708 gp_Pnt O2 = aGen2.Location();
1709 Standard_Real D1DotD2 = D1.Dot(D2);
1710 Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1711 gp_Vec O1O2 (O1,O2);
1712 Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1714 gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1716 aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1717 aN = aD1.Crossed(aD2);
1719 //Plane that must contain intersection curves
1720 gp_Pln anIntPln(myPChar, aN);
1722 IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1724 if(INTER_QUAD_PLN.IsDone()) {
1725 switch(INTER_QUAD_PLN.TypeInter()) {
1726 case IntAna_Ellipse: {
1727 typeres=IntAna_Ellipse;
1728 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1730 dir1 = E.Position().Direction();
1731 dir2 = E.Position().XDirection();
1732 param1 = E.MajorRadius();
1733 param1bis = E.MinorRadius();
1737 case IntAna_Circle: {
1738 typeres=IntAna_Circle;
1739 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1741 dir1 = C.Position().XDirection();
1742 dir2 = C.Position().YDirection();
1743 param1 = C.Radius();
1747 case IntAna_Parabola: {
1748 typeres=IntAna_Parabola;
1749 gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1750 pt1 = Prb.Location();
1751 dir1 = Prb.Position().Direction();
1752 dir2 = Prb.Position().XDirection();
1753 param1 = Prb.Focal();
1757 case IntAna_Hyperbola: {
1758 typeres=IntAna_Hyperbola;
1759 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1760 pt1 = pt2 = H.Location();
1761 dir1 = H.Position().Direction();
1762 dir2 = H.Position().XDirection();
1763 param1 = param2 = H.MajorRadius();
1764 param1bis = param2bis = H.MinorRadius();
1769 typeres=IntAna_NoGeometricSolution;
1775 typeres=IntAna_NoGeometricSolution;
1778 //=======================================================================
1779 //function : IntAna_QuadQuadGeo
1780 //purpose : Sphere - Cone
1781 //=======================================================================
1782 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1784 const Standard_Real Tol)
1785 : done(Standard_False),
1787 typeres(IntAna_Empty),
1798 myCommonGen(Standard_False),
1802 Perform(Sph,Con,Tol);
1804 //=======================================================================
1805 //function : Perform
1807 //=======================================================================
1808 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1810 const Standard_Real)
1816 AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1817 gp_Pnt Pt=Sph.Location();
1819 if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1821 gp_Pnt ConApex= Con.Apex();
1822 Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1824 if(dApexSphCenter>RealEpsilon()) {
1825 ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1828 ConDir = Con.Position().Direction();
1831 Standard_Real Rad=Sph.Radius();
1832 Standard_Real tga=Tan(Con.SemiAngle());
1836 //-- x: Roots of (x**2 + y**2 = Rad**2)
1837 //-- tga = y / (x+dApexSphCenter)
1838 Standard_Real tgatga = tga * tga;
1839 math_DirectPolynomialRoots Eq( 1.0+tgatga
1840 ,2.0*tgatga*dApexSphCenter
1841 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1843 Standard_Integer nbsol=Eq.NbSolutions();
1845 typeres=IntAna_Empty;
1848 typeres=IntAna_Circle;
1850 Standard_Real x = Eq.Value(1);
1851 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1853 pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1854 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1855 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1856 param1 = tga * dApexSphCenterpx;
1857 param1 = Abs(param1);
1859 if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1860 typeres=IntAna_PointAndCircle;
1865 Standard_Real x=Eq.Value(2);
1866 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1868 pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1869 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1870 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1871 param2 = tga * dApexSphCenterpx;
1872 param2 = Abs(param2);
1874 if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1875 typeres=IntAna_PointAndCircle;
1882 done=Standard_False;
1886 typeres=IntAna_NoGeometricSolution;
1890 //=======================================================================
1891 //function : IntAna_QuadQuadGeo
1892 //purpose : Sphere - Sphere
1893 //=======================================================================
1894 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1895 ,const gp_Sphere& Sph2
1896 ,const Standard_Real Tol)
1897 : done(Standard_False),
1899 typeres(IntAna_Empty),
1910 myCommonGen(Standard_False),
1914 Perform(Sph1,Sph2,Tol);
1916 //=======================================================================
1917 //function : Perform
1919 //=======================================================================
1920 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1921 const gp_Sphere& Sph2,
1922 const Standard_Real Tol)
1925 gp_Pnt O1=Sph1.Location();
1926 gp_Pnt O2=Sph2.Location();
1927 Standard_Real dO1O2=O1.Distance(O2);
1928 Standard_Real R1=Sph1.Radius();
1929 Standard_Real R2=Sph2.Radius();
1930 Standard_Real Rmin,Rmax;
1931 typeres=IntAna_Empty;
1932 param2bis=0.0; //-- pour eviter param2bis not used ....
1934 if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1936 if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1937 typeres = IntAna_Same;
1940 if(dO1O2<=Tol) { return; }
1941 gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1942 Standard_Real t = Rmax - dO1O2 - Rmin;
1944 //----------------------------------------------------------------------
1945 //-- |----------------- R1 --------------------|
1946 //-- |----dO1O2-----|-----------R2----------|
1949 //-- |------ R1 ------|---------dO1O2----------|
1950 //-- |-------------------R2-----------------------|
1952 //----------------------------------------------------------------------
1953 if(t >= 0.0 && t <=Tol) {
1954 typeres = IntAna_Point;
1957 if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1958 else t2=(-R1+(dO1O2-R2))*0.5;
1960 pt1.SetCoord( O1.X() + t2*Dir.X()
1961 ,O1.Y() + t2*Dir.Y()
1962 ,O1.Z() + t2*Dir.Z());
1965 //-----------------------------------------------------------------
1966 //-- |----------------- dO1O2 --------------------|
1967 //-- |----R1-----|-----------R2----------|-Tol-|
1969 //-- |----------------- Rmax --------------------|
1970 //-- |----Rmin----|-------dO1O2-------|-Tol-|
1972 //-----------------------------------------------------------------
1973 if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1974 typeres=IntAna_Empty;
1977 //---------------------------------------------------------------
1980 //---------------------------------------------------------------
1981 Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1982 Standard_Real Beta = R1*R1-Alpha*Alpha;
1983 Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1985 if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1986 typeres = IntAna_Point;
1987 Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1990 typeres = IntAna_Circle;
1994 pt1.SetCoord( O1.X() + Alpha*Dir.X()
1995 ,O1.Y() + Alpha*Dir.Y()
1996 ,O1.Z() + Alpha*Dir.Z());
2004 //=======================================================================
2005 //function : IntAna_QuadQuadGeo
2006 //purpose : Plane - Torus
2007 //=======================================================================
2008 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
2009 const gp_Torus& Tor,
2010 const Standard_Real Tol)
2011 : done(Standard_False),
2013 typeres(IntAna_Empty),
2024 myCommonGen(Standard_False),
2028 Perform(Pln,Tor,Tol);
2030 //=======================================================================
2031 //function : Perform
2033 //=======================================================================
2034 void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
2035 const gp_Torus& Tor,
2036 const Standard_Real Tol)
2038 done = Standard_True;
2040 Standard_Real aRMin, aRMaj;
2042 aRMin = Tor.MinorRadius();
2043 aRMaj = Tor.MajorRadius();
2044 if (aRMin >= aRMaj) {
2045 typeres = IntAna_NoGeometricSolution;
2049 const gp_Ax1 aPlnAx = Pln.Axis();
2050 const gp_Ax1 aTorAx = Tor.Axis();
2052 Standard_Boolean bParallel, bNormal;
2054 bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
2055 bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
2056 if (!bNormal && !bParallel) {
2057 typeres = IntAna_NoGeometricSolution;
2061 Standard_Real aDist;
2063 gp_Pnt aTorLoc = aTorAx.Location();
2065 Standard_Real aDt, X, Y, Z, A, B, C, D, aDR, aTolNum;
2067 aTolNum=myEPSILON_CYLINDER_DELTA_RADIUS;
2069 Pln.Coefficients(A,B,C,D);
2070 aTorLoc.Coord(X,Y,Z);
2071 aDist = A*X + B*Y + C*Z + D;
2073 aDR=Abs(aDist) - aRMin;
2074 if (aDR > aTolNum) {
2075 typeres=IntAna_Empty;
2079 if (Abs(aDR) < aTolNum) {
2083 typeres = IntAna_Circle;
2085 pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
2086 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2087 param1 = aRMaj + aDt;
2088 dir1 = aTorAx.Direction();
2090 if ((aDR < -aTolNum) && (aDt > Tol)) {
2092 param2 = aRMaj - aDt;
2099 aDist = Pln.Distance(aTorLoc);
2100 if (aDist > myEPSILON_DISTANCE) {
2101 typeres = IntAna_NoGeometricSolution;
2105 typeres = IntAna_Circle;
2106 param2 = param1 = aRMin;
2107 dir2 = dir1 = aPlnAx.Direction();
2110 gp_Dir aDir = aTorAx.Direction()^dir1;
2111 pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
2112 pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
2116 //=======================================================================
2117 //function : IntAna_QuadQuadGeo
2118 //purpose : Cylinder - Torus
2119 //=======================================================================
2120 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
2121 const gp_Torus& Tor,
2122 const Standard_Real Tol)
2123 : done(Standard_False),
2125 typeres(IntAna_Empty),
2136 myCommonGen(Standard_False),
2140 Perform(Cyl,Tor,Tol);
2142 //=======================================================================
2143 //function : Perform
2145 //=======================================================================
2146 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
2147 const gp_Torus& Tor,
2148 const Standard_Real Tol)
2150 done = Standard_True;
2152 Standard_Real aRMin, aRMaj;
2154 aRMin = Tor.MinorRadius();
2155 aRMaj = Tor.MajorRadius();
2156 if (aRMin >= aRMaj) {
2157 typeres = IntAna_NoGeometricSolution;
2161 const gp_Ax1 aCylAx = Cyl.Axis();
2162 const gp_Ax1 aTorAx = Tor.Axis();
2164 const gp_Lin aLin(aTorAx);
2165 const gp_Pnt aLocCyl = Cyl.Location();
2167 if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
2168 (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2169 typeres = IntAna_NoGeometricSolution;
2173 Standard_Real aRCyl;
2175 aRCyl = Cyl.Radius();
2176 if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2177 typeres = IntAna_Empty;
2181 typeres = IntAna_Circle;
2183 Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2184 gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2186 dir1 = aTorAx.Direction();
2187 pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2190 if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2191 (aRCyl < (aRMaj + aRMin))) {
2193 pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2199 //=======================================================================
2200 //function : IntAna_QuadQuadGeo
2201 //purpose : Cone - Torus
2202 //=======================================================================
2203 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2204 const gp_Torus& Tor,
2205 const Standard_Real Tol)
2206 : done(Standard_False),
2208 typeres(IntAna_Empty),
2219 myCommonGen(Standard_False),
2223 Perform(Con,Tor,Tol);
2225 //=======================================================================
2226 //function : Perform
2228 //=======================================================================
2229 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2230 const gp_Torus& Tor,
2231 const Standard_Real Tol)
2233 done = Standard_True;
2235 Standard_Real aRMin, aRMaj;
2237 aRMin = Tor.MinorRadius();
2238 aRMaj = Tor.MajorRadius();
2239 if (aRMin >= aRMaj) {
2240 typeres = IntAna_NoGeometricSolution;
2244 const gp_Ax1 aConAx = Con.Axis();
2245 const gp_Ax1 aTorAx = Tor.Axis();
2247 const gp_Lin aLin(aTorAx);
2248 const gp_Pnt aConApex = Con.Apex();
2250 if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2251 (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2252 typeres = IntAna_NoGeometricSolution;
2256 Standard_Real anAngle, aDist, aParam[4], aDt;
2258 gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2261 anAngle = Con.SemiAngle();
2262 aTorLoc = aTorAx.Location();
2264 aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2265 gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2266 gp_Ax1 anAxCLRot(aConApex, aDN);
2267 gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2268 gp_Dir aDL = aConL.Position().Direction();
2269 gp_Dir aXDir = Tor.XAxis().Direction();
2271 typeres = IntAna_Empty;
2273 for (i = 0; i < 2; ++i) {
2277 aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2279 aDist = aConL.Distance(aPCT);
2280 if (aDist > aRMin+Tol) {
2284 typeres = IntAna_Circle;
2286 gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
2287 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2290 gp_XYZ aDVal = aDt*aDL.XYZ();
2291 aP.SetXYZ(aPh + aDVal);
2292 aParam[nbint] = aLin.Distance(aP);
2293 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2294 aDir[nbint] = aTorAx.Direction();
2296 if ((aDist < aRMin) && (aDt > Tol)) {
2297 aP.SetXYZ(aPh - aDVal);
2298 aParam[nbint] = aLin.Distance(aP);
2299 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2300 aDir[nbint] = aDir[nbint-1];
2305 for (i = 0; i < nbint; ++i) {
2337 //=======================================================================
2338 //function : IntAna_QuadQuadGeo
2339 //purpose : Sphere - Torus
2340 //=======================================================================
2341 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2342 const gp_Torus& Tor,
2343 const Standard_Real Tol)
2344 : done(Standard_False),
2346 typeres(IntAna_Empty),
2357 myCommonGen(Standard_False),
2361 Perform(Sph,Tor,Tol);
2363 //=======================================================================
2364 //function : Perform
2366 //=======================================================================
2367 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2368 const gp_Torus& Tor,
2369 const Standard_Real Tol)
2371 done = Standard_True;
2373 Standard_Real aRMin, aRMaj;
2375 aRMin = Tor.MinorRadius();
2376 aRMaj = Tor.MajorRadius();
2377 if (aRMin >= aRMaj) {
2378 typeres = IntAna_NoGeometricSolution;
2382 const gp_Ax1 aTorAx = Tor.Axis();
2383 const gp_Lin aLin(aTorAx);
2384 const gp_Pnt aSphLoc = Sph.Location();
2386 if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2387 typeres = IntAna_NoGeometricSolution;
2391 Standard_Real aRSph, aDist;
2394 gp_Dir aXDir = Tor.XAxis().Direction();
2395 aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2396 aRSph = Sph.Radius();
2398 gp_Vec aVec12(aTorLoc, aSphLoc);
2399 aDist = aVec12.Magnitude();
2400 if (((aDist - Tol) > (aRMin + aRSph)) ||
2401 ((aDist + Tol) < Abs(aRMin - aRSph))) {
2402 typeres = IntAna_Empty;
2406 typeres = IntAna_Circle;
2408 Standard_Real anAlpha, aBeta;
2410 anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2411 aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2413 gp_Dir aDir12(aVec12);
2414 gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2415 gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2418 gp_XYZ aDVal = aBeta*aDC.XYZ();
2419 aP.SetXYZ(aPh + aDVal);
2420 param1 = aLin.Distance(aP);
2421 pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2422 dir1 = aTorAx.Direction();
2424 if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) &&
2425 (aDVal.Modulus() > Tol)) {
2426 aP.SetXYZ(aPh - aDVal);
2427 param2 = aLin.Distance(aP);
2428 pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2434 //=======================================================================
2435 //function : IntAna_QuadQuadGeo
2436 //purpose : Torus - Torus
2437 //=======================================================================
2438 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2439 const gp_Torus& Tor2,
2440 const Standard_Real Tol)
2441 : done(Standard_False),
2443 typeres(IntAna_Empty),
2454 myCommonGen(Standard_False),
2458 Perform(Tor1,Tor2,Tol);
2460 //=======================================================================
2461 //function : Perform
2463 //=======================================================================
2464 void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2465 const gp_Torus& Tor2,
2466 const Standard_Real Tol)
2468 done = Standard_True;
2470 Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2472 aRMin1 = Tor1.MinorRadius();
2473 aRMaj1 = Tor1.MajorRadius();
2474 aRMin2 = Tor2.MinorRadius();
2475 aRMaj2 = Tor2.MajorRadius();
2477 const gp_Ax1& anAx1 = Tor1.Axis();
2478 const gp_Ax1& anAx2 = Tor2.Axis();
2480 const gp_Pnt& aLoc1 = anAx1.Location();
2481 const gp_Pnt& aLoc2 = anAx2.Location();
2484 if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2485 (aL1.Distance(aLoc2) > myEPSILON_DISTANCE)) {
2486 typeres = IntAna_NoGeometricSolution;
2490 if (aLoc1.IsEqual(aLoc2, Tol) &&
2491 (Abs(aRMin1 - aRMin2) <= Tol) &&
2492 (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2493 typeres = IntAna_Same;
2497 if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2498 typeres = IntAna_NoGeometricSolution;
2502 Standard_Real aDist;
2505 gp_Dir aXDir1 = Tor1.XAxis().Direction();
2506 aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2507 aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2509 gp_Vec aV12(aP1, aP2);
2510 aDist = aV12.Magnitude();
2511 if (((aDist - Tol) > (aRMin1 + aRMin2)) ||
2512 ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2513 typeres = IntAna_Empty;
2517 typeres = IntAna_Circle;
2519 Standard_Real anAlpha, aBeta;
2521 anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2522 aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2524 gp_Dir aDir12(aV12);
2525 gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2526 gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2529 gp_XYZ aDVal = aBeta*aDC.XYZ();
2530 aP.SetXYZ(aPh + aDVal);
2531 param1 = aL1.Distance(aP);
2532 pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2533 dir1 = anAx1.Direction();
2535 if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) &&
2536 aDVal.Modulus() > Tol) {
2537 aP.SetXYZ(aPh - aDVal);
2538 param2 = aL1.Distance(aP);
2539 pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2545 //=======================================================================
2547 //purpose : Returns a Point
2548 //=======================================================================
2549 gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
2551 if(!done) { throw StdFail_NotDone(); }
2552 if(n>nbint || n<1) { throw Standard_DomainError(); }
2553 if(typeres==IntAna_PointAndCircle) {
2554 if(n!=1) { throw Standard_DomainError(); }
2555 if(param1==0.0) return(pt1);
2558 else if(typeres==IntAna_Point) {
2559 if(n==1) return(pt1);
2563 // WNT (what can you expect from MicroSoft ?)
2564 return gp_Pnt(0,0,0);
2566 //=======================================================================
2568 //purpose : Returns a Line
2569 //=======================================================================
2570 gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
2572 if(!done) { throw StdFail_NotDone(); }
2573 if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2574 throw Standard_DomainError();
2576 if(n==1) { return(gp_Lin(pt1,dir1)); }
2577 else { return(gp_Lin(pt2,dir2)); }
2579 //=======================================================================
2581 //purpose : Returns a Circle
2582 //=======================================================================
2583 gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
2585 if(!done) { throw StdFail_NotDone(); }
2586 if(typeres==IntAna_PointAndCircle) {
2587 if(n!=1) { throw Standard_DomainError(); }
2588 if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2589 return(gp_Circ(DirToAx2(pt2,dir2),param2));
2591 else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2592 throw Standard_DomainError();
2594 if (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2595 else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2596 else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2597 else { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
2600 //=======================================================================
2601 //function : Ellipse
2602 //purpose : Returns a Elips
2603 //=======================================================================
2604 gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2606 if(!done) { throw StdFail_NotDone(); }
2607 if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2608 throw Standard_DomainError();
2612 Standard_Real R1=param1, R2=param1bis, aTmp;
2614 aTmp=R1; R1=R2; R2=aTmp;
2616 gp_Ax2 anAx2(pt1, dir1 ,dir2);
2617 gp_Elips anElips (anAx2, R1, R2);
2621 Standard_Real R1=param2, R2=param2bis, aTmp;
2623 aTmp=R1; R1=R2; R2=aTmp;
2625 gp_Ax2 anAx2(pt2, dir2 ,dir1);
2626 gp_Elips anElips (anAx2, R1, R2);
2630 //=======================================================================
2631 //function : Parabola
2632 //purpose : Returns a Parabola
2633 //=======================================================================
2634 gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
2637 throw StdFail_NotDone();
2639 if (typeres!=IntAna_Parabola) {
2640 throw Standard_DomainError();
2642 if((n>nbint) || (n!=1)) {
2643 throw Standard_OutOfRange();
2645 return(gp_Parab(gp_Ax2( pt1
2650 //=======================================================================
2651 //function : Hyperbola
2652 //purpose : Returns a Hyperbola
2653 //=======================================================================
2654 gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
2657 throw StdFail_NotDone();
2659 if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2660 throw Standard_DomainError();
2663 return(gp_Hypr(gp_Ax2( pt1
2666 ,param1,param1bis));
2669 return(gp_Hypr(gp_Ax2( pt2
2672 ,param2,param2bis));
2675 //=======================================================================
2676 //function : HasCommonGen
2678 //=======================================================================
2679 Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2683 //=======================================================================
2686 //=======================================================================
2687 const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2691 //=======================================================================
2692 //function : RefineDir
2694 //=======================================================================
2695 void RefineDir(gp_Dir& aDir)
2697 Standard_Integer k, m, n;
2698 Standard_Real aC[3];
2700 aDir.Coord(aC[0], aC[1], aC[2]);
2704 for (k=0; k<3; ++k) {
2705 if (aC[k]==1. || aC[k]==-1.) {
2708 else if (aC[k]!=0.) {
2714 Standard_Real aEps, aR1, aR2, aNum;
2720 for (k=0; k<3; ++k) {
2722 aNum=(m)? aC[k] : -aC[k];
2723 if (aNum>aR1 && aNum<aR2) {
2736 aDir.SetCoord(aC[0], aC[1], aC[2]);