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
27 #include <IntAna_QuadQuadGeo.ixx>
29 #include <IntAna_IntConicQuad.hxx>
30 #include <StdFail_NotDone.hxx>
31 #include <Standard_DomainError.hxx>
32 #include <Standard_OutOfRange.hxx>
33 #include <math_DirectPolynomialRoots.hxx>
43 #include <gp_Pnt2d.hxx>
44 #include <gp_Vec2d.hxx>
45 #include <gp_Dir2d.hxx>
49 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D);
51 void RefineDir(gp_Dir& aDir);
53 //=======================================================================
55 //purpose : O p e r a t i o n s D i v e r s e s s u r d e s A x 1
56 //=======================================================================
59 AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2);
61 void Distance(Standard_Real& dist,
62 Standard_Real& Param1,
63 Standard_Real& Param2);
65 gp_Pnt PtIntersect() {
68 Standard_Boolean Coplanar(void) {
71 Standard_Boolean Same(void) {
72 return theparallel && (thedistance<myEPSILON_DISTANCE);
74 Standard_Real Distance(void) {
77 Standard_Boolean Intersect(void) {
78 return (thecoplanar && (!theparallel));
80 Standard_Boolean Parallel(void) {
83 Standard_Boolean Normal(void) {
88 Standard_Real Det33(const Standard_Real a11,
89 const Standard_Real a12,
90 const Standard_Real a13,
91 const Standard_Real a21,
92 const Standard_Real a22,
93 const Standard_Real a23,
94 const Standard_Real a31,
95 const Standard_Real a32,
96 const Standard_Real a33) {
97 Standard_Real theReturn =
98 a11*(a22*a33-a32*a23) - a21*(a12*a33-a32*a13) + a31*(a12*a23-a22*a13) ;
106 Standard_Real thedistance;
107 Standard_Boolean theparallel;
108 Standard_Boolean thecoplanar;
109 Standard_Boolean thenormal;
111 Standard_Real myEPSILON_DISTANCE;
112 Standard_Real myEPSILON_AXES_PARA;
115 //=======================================================================
116 //function : AxeOperator::AxeOperator
118 //=======================================================================
119 AxeOperator::AxeOperator(const gp_Ax1& A1,const gp_Ax1& A2)
121 myEPSILON_DISTANCE=0.00000000000001;
122 myEPSILON_AXES_PARA=0.000000000001;
125 //---------------------------------------------------------------------
126 gp_Dir V1=Axe1.Direction();
127 gp_Dir V2=Axe2.Direction();
128 gp_Pnt P1=Axe1.Location();
129 gp_Pnt P2=Axe2.Location();
133 thecoplanar= Standard_False;
134 thenormal = Standard_False;
136 //--- check if the two axis are parallel
137 theparallel=V1.IsParallel(V2, myEPSILON_AXES_PARA);
138 //--- Distance between the two axis
139 gp_XYZ perp(A1.Direction().XYZ().Crossed(A2.Direction().XYZ()));
142 thedistance = L1.Distance(A2.Location());
145 thedistance = Abs(gp_Vec(perp.Normalized()).Dot(gp_Vec(Axe1.Location(),
148 //--- check if Axis are Coplanar
150 if(thedistance<myEPSILON_DISTANCE) {
151 D33=Det33(V1.X(),V1.Y(),V1.Z()
152 ,V2.X(),V2.Y(),V2.Z()
153 ,P1.X()-P2.X(),P1.Y()-P2.Y(),P1.Z()-P2.Z());
154 if(Abs(D33)<=myEPSILON_DISTANCE) {
155 thecoplanar=Standard_True;
159 thecoplanar=Standard_True;
160 thenormal=(V1.Dot(V2)==0.0)? Standard_True : Standard_False;
162 //--- check if the two axis are concurrent
163 if(thecoplanar && (!theparallel)) {
164 Standard_Real smx=P2.X() - P1.X();
165 Standard_Real smy=P2.Y() - P1.Y();
166 Standard_Real smz=P2.Z() - P1.Z();
167 Standard_Real Det1,Det2,Det3,A;
168 Det1=V1.Y() * V2.X() - V1.X() * V2.Y();
169 Det2=V1.Z() * V2.Y() - V1.Y() * V2.Z();
170 Det3=V1.Z() * V2.X() - V1.X() * V2.Z();
172 if((Det1!=0.0) && ((Abs(Det1) >= Abs(Det2))&&(Abs(Det1) >= Abs(Det3)))) {
173 A=(smy * V2.X() - smx * V2.Y())/Det1;
176 && ((Abs(Det2) >= Abs(Det1))
177 &&(Abs(Det2) >= Abs(Det3)))) {
178 A=(smz * V2.Y() - smy * V2.Z())/Det2;
181 A=(smz * V2.X() - smx * V2.Z())/Det3;
183 ptintersect.SetCoord( P1.X() + A * V1.X()
185 ,P1.Z() + A * V1.Z());
188 ptintersect.SetCoord(0,0,0); //-- Pour eviter des FPE
191 //=======================================================================
192 //function : Distance
194 //=======================================================================
195 void AxeOperator::Distance(Standard_Real& dist,
196 Standard_Real& Param1,
197 Standard_Real& Param2)
199 gp_Vec O1O2(Axe1.Location(),Axe2.Location());
200 gp_Dir U1 = Axe1.Direction(); //-- juste pour voir.
201 gp_Dir U2 = Axe2.Direction();
203 gp_Dir N = U1.Crossed(U2);
204 Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
206 U1.Z(),U2.Z(),N.Z());
208 dist = Det33(U1.X(),U2.X(),O1O2.X(),
209 U1.Y(),U2.Y(),O1O2.Y(),
210 U1.Z(),U2.Z(),O1O2.Z()) / D;
211 Param1 = Det33(O1O2.X(),U2.X(),N.X(),
212 O1O2.Y(),U2.Y(),N.Y(),
213 O1O2.Z(),U2.Z(),N.Z()) / (-D);
214 //------------------------------------------------------------
215 //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
216 //-- soit : Segment perpendiculaire : O1+P1 D1
218 Param2 = Det33(U1.X(),O1O2.X(),N.X(),
219 U1.Y(),O1O2.Y(),N.Y(),
220 U1.Z(),O1O2.Z(),N.Z()) / (D);
223 //=======================================================================
224 //function : DirToAx2
225 //purpose : returns a gp_Ax2 where D is the main direction
226 //=======================================================================
227 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
229 Standard_Real x=D.X(); Standard_Real ax=Abs(x);
230 Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
231 Standard_Real z=D.Z(); Standard_Real az=Abs(z);
232 if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
233 return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
235 else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
236 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
239 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
242 //=======================================================================
243 //function : IntAna_QuadQuadGeo
244 //purpose : Empty constructor
245 //=======================================================================
246 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
247 : done(Standard_False),
249 typeres(IntAna_Empty),
260 myCommonGen(Standard_False),
265 //=======================================================================
266 //function : InitTolerances
268 //=======================================================================
269 void IntAna_QuadQuadGeo::InitTolerances()
271 myEPSILON_DISTANCE = 0.00000000000001;
272 myEPSILON_ANGLE_CONE = 0.000000000001;
273 myEPSILON_MINI_CIRCLE_RADIUS = 0.000000001;
274 myEPSILON_CYLINDER_DELTA_RADIUS = 0.0000000000001;
275 myEPSILON_CYLINDER_DELTA_DISTANCE= 0.0000001;
276 myEPSILON_AXES_PARA = 0.000000000001;
278 //=======================================================================
279 //function : IntAna_QuadQuadGeo
281 //=======================================================================
282 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1,
284 const Standard_Real TolAng,
285 const Standard_Real Tol)
286 : done(Standard_False),
288 typeres(IntAna_Empty),
299 myCommonGen(Standard_False),
303 Perform(P1,P2,TolAng,Tol);
305 //=======================================================================
308 //=======================================================================
309 void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
311 const Standard_Real TolAng,
312 const Standard_Real Tol)
314 Standard_Real A1, B1, C1, D1, A2, B2, C2, D2, dist1, dist2, aMVD;
319 P1.Coefficients(A1,B1,C1,D1);
320 P2.Coefficients(A2,B2,C2,D2);
322 gp_Vec aVN1(A1,B1,C1);
323 gp_Vec aVN2(A2,B2,C2);
324 gp_Vec vd(aVN1.Crossed(aVN2));
326 const gp_Pnt& aLocP1=P1.Location();
327 const gp_Pnt& aLocP2=P2.Location();
329 dist1=A2*aLocP1.X() + B2*aLocP1.Y() + C2*aLocP1.Z() + D2;
330 dist2=A1*aLocP2.X() + B1*aLocP2.Y() + C1*aLocP2.Z() + D1;
334 // normalles are collinear - planes are same or parallel
335 typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same
339 Standard_Real denom, denom2, ddenom, par1, par2;
340 Standard_Real X1, Y1, Z1, X2, Y2, Z2, aEps;
343 denom=A1*A2 + B1*B2 + C1*C2;
344 denom2 = denom*denom;
345 ddenom = 1. - denom2;
347 denom = ( Abs(ddenom) <= aEps ) ? aEps : ddenom;
352 gp_Vec inter1(aVN1.Crossed(vd));
353 gp_Vec inter2(aVN2.Crossed(vd));
355 X1=aLocP1.X() + par1*inter1.X();
356 Y1=aLocP1.Y() + par1*inter1.Y();
357 Z1=aLocP1.Z() + par1*inter1.Z();
358 X2=aLocP2.X() + par2*inter2.X();
359 Y2=aLocP2.Y() + par2*inter2.Y();
360 Z2=aLocP2.Z() + par2*inter2.Z();
362 pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
364 typeres = IntAna_Line;
367 //-------------------------------------------------------
368 // When the value of the angle between the planes is small
369 // the origin of intersection line is computed with error
370 // [ ~0.0001 ] that can not br considered as small one
372 // for {A~=2.e-6, dist1=4.2e-5, dist2==1.e-4} =>
373 // {denom=3.4e-12, par1=12550297.6, par2=32605552.9, etc}
375 // the origin should be refined if it is possible
377 Standard_Real aTreshAng, aTreshDist;
379 aTreshAng=2.e-6; // 1.e-4 deg
382 if (aMVD < aTreshAng) {
383 Standard_Real aDist1, aDist2;
385 aDist1=A1*pt1.X() + B1*pt1.Y() + C1*pt1.Z() + D1;
386 aDist2=A2*pt1.X() + B2*pt1.Y() + C2*pt1.Z() + D2;
388 if (fabs(aDist1)>aTreshDist || fabs(aDist2)>aTreshDist) {
389 Standard_Boolean bIsDone, bIsParallel;
390 IntAna_IntConicQuad aICQ;
394 gp_Lin aL1(pt1, aDN1);
396 aICQ.Perform(aL1, P1, TolAng, Tol);
397 bIsDone=aICQ.IsDone();
402 const gp_Pnt& aPnt1=aICQ.Point(1);
403 //----------------------------------
405 gp_Dir aDL2(dir1.Crossed(aDN1));
406 gp_Lin aL2(aPnt1, aDL2);
408 aICQ.Perform(aL2, P2, TolAng, Tol);
409 bIsDone=aICQ.IsDone();
414 bIsParallel=aICQ.IsParallel();
419 const gp_Pnt& aPnt2=aICQ.Point(1);
427 //=======================================================================
428 //function : IntAna_QuadQuadGeo
429 //purpose : Pln Cylinder
430 //=======================================================================
431 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
432 ,const gp_Cylinder& Cl
433 ,const Standard_Real Tolang
434 ,const Standard_Real Tol
435 ,const Standard_Real H)
436 : done(Standard_False),
438 typeres(IntAna_Empty),
449 myCommonGen(Standard_False),
453 Perform(P,Cl,Tolang,Tol,H);
455 //=======================================================================
458 //=======================================================================
459 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
460 ,const gp_Cylinder& Cl
461 ,const Standard_Real Tolang
462 ,const Standard_Real Tol
463 ,const Standard_Real H)
465 done = Standard_False;
466 Standard_Real dist,radius;
467 Standard_Real A,B,C,D;
469 Standard_Real sint,cost,h;
470 gp_XYZ axex,axey,omega;
474 radius = Cl.Radius();
476 gp_Lin axec(Cl.Axis());
477 gp_XYZ normp(P.Axis().Direction().XYZ());
479 P.Coefficients(A,B,C,D);
480 axec.Location().Coord(X,Y,Z);
481 // la distance axe/plan est evaluee a l origine de l axe.
482 dist = A*X + B*Y + C*Z + D;
484 Standard_Real tolang = Tolang;
485 Standard_Boolean newparams = Standard_False;
487 gp_Vec ldv( axec.Direction() );
489 Standard_Real dA = Abs( ldv.Angle( npv ) );
492 Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
493 Standard_Real dangle = Abs( dang );
494 if( dangle > Tolang )
496 Standard_Real sinda = Abs( Sin( dangle ) );
497 Standard_Real dif = Abs( sinda - Tol );
501 newparams = Standard_True;
507 IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
509 if (inter.IsParallel()) {
510 // Le resultat de l intersection Plan-Cylindre est de type droite.
511 // il y a 1 ou 2 droites
513 typeres = IntAna_Line;
514 omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
516 if (Abs(Abs(dist)-radius) < Tol)
523 gp_XYZ omegaXYZ(X,Y,Z);
524 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
525 Standard_Real Xt,Yt,Zt,distt;
526 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
527 distt = A*Xt + B*Yt + C*Zt + D;
528 gp_XYZ omega1(omegaXYZtrnsl.X()-distt*A,
529 omegaXYZtrnsl.Y()-distt*B,
530 omegaXYZtrnsl.Z()-distt*C );
532 ppt1.SetXYZ( omega1 );
533 gp_Vec vv1(pt1,ppt1);
538 dir1 = axec.Direction();
540 else if (Abs(dist) < radius)
543 h = Sqrt(radius*radius - dist*dist);
544 axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
546 pt1.SetXYZ(omega - h*axey);
547 pt2.SetXYZ(omega + h*axey);
551 gp_XYZ omegaXYZ(X,Y,Z);
552 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
553 Standard_Real Xt,Yt,Zt,distt,ht;
554 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
555 distt = A*Xt + B*Yt + C*Zt + D;
556 // ht = Sqrt(radius*radius - distt*distt);
557 Standard_Real anSqrtArg = radius*radius - distt*distt;
558 ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
560 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A,
561 omegaXYZtrnsl.Y()-distt*B,
562 omegaXYZtrnsl.Z()-distt*C );
564 ppt1.SetXYZ( omega1 - ht*axey);
565 ppt2.SetXYZ( omega1 + ht*axey);
566 gp_Vec vv1(pt1,ppt1);
567 gp_Vec vv2(pt2,ppt2);
575 dir1 = axec.Direction();
576 dir2 = axec.Direction();
581 // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
582 // et ne pas etre seulement supprime...
585 typeres = IntAna_Empty;
588 else { // Il y a un point d intersection. C est le centre du cercle
589 // ou de l ellipse solution.
592 axey = normp.Crossed(axec.Direction().XYZ());
593 sint = axey.Modulus();
595 pt1 = inter.Point(1);
597 if (sint < Tol/radius) {
599 // on construit un cercle avec comme axes X et Y ceux du cylindre
600 typeres = IntAna_Circle;
602 dir1 = axec.Direction(); // axe Z
603 dir2 = Cl.Position().XDirection();
608 // on construit un ellipse
609 typeres = IntAna_Ellipse;
610 cost = Abs(axec.Direction().XYZ().Dot(normp));
611 axex = axey.Crossed(normp);
613 dir1.SetXYZ(normp); //Modif ds ce bloc
616 param1 = radius/cost;
621 done = Standard_True;
623 //=======================================================================
624 //function : IntAna_QuadQuadGeo
626 //=======================================================================
627 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
629 const Standard_Real Tolang,
630 const Standard_Real Tol)
631 : done(Standard_False),
633 typeres(IntAna_Empty),
644 myCommonGen(Standard_False),
648 Perform(P,Co,Tolang,Tol);
650 //=======================================================================
653 //=======================================================================
654 void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
656 const Standard_Real Tolang,
657 const Standard_Real Tol)
660 done = Standard_False;
663 Standard_Real A,B,C,D;
665 Standard_Real dist,sint,cost,sina,cosa,angl,costa;
669 gp_Lin axec(Co.Axis());
670 P.Coefficients(A,B,C,D);
671 gp_Pnt apex(Co.Apex());
674 dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
676 gp_XYZ normp = P.Axis().Direction().XYZ();
677 if(P.Direct()==Standard_False) { //-- lbr le 14 jan 97
681 axey = normp.Crossed(Co.Axis().Direction().XYZ());
682 axex = axey.Crossed(normp);
685 angl = Co.SemiAngle();
688 sina = Abs(Sin(angl));
691 // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
693 sint = axey.Modulus();
694 cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
696 // Le calcul de costa permet de determiner si le plan contient
697 // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
699 costa = cost*cosa - sint*sina; // sin((PI/2 -t)-angl)=cos(t+angl)
700 // avec t ramene entre 0 et pi/2.
702 if (Abs(dist) < Tol) {
703 // on considere que le plan contient le sommet du cone.
704 // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
705 // selon l inclinaison du plan.
707 if (Abs(costa) < Tolang) { // plan parallele a la generatrice
708 typeres = IntAna_Line;
710 gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
711 // point sur l axe du cone cote z positif
713 dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
714 ptonaxe = ptonaxe - dist*normp;
716 dir1.SetXYZ(ptonaxe - pt1.XYZ());
718 else if (cost < sina) { // plan "interieur" au cone
719 typeres = IntAna_Line;
723 dh = Sqrt(sina*sina-cost*cost)/cosa;
724 dir1.SetXYZ(axex + dh*axey);
725 dir2.SetXYZ(axex - dh*axey);
727 else { // plan "exterieur" au cone
728 typeres = IntAna_Point;
734 // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
735 // l inclinaison du plan.
736 Standard_Real deltacenter, distance;
739 // Le plan contient la direction de l axe du cone. La solution est
741 typeres = IntAna_Hyperbola;
743 pt1.SetXYZ(apex.XYZ()-dist*normp);
747 param1 = param2 = Abs(dist/Tan(angl));
748 param1bis = param2bis = Abs(dist);
752 IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
754 gp_Pnt center(inter.Point(1));
756 // En fonction de la position de l intersection par rapport au sommet
757 // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
758 // sur axec est negatif (voir definition du cone)
760 distance = apex.Distance(center);
762 if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
767 if (Abs(costa) < Tolang) { // plan parallele a une generatrice
768 typeres = IntAna_Parabola;
770 deltacenter = distance/2./cosa;
772 pt1.SetXYZ(center.XYZ()-deltacenter*axex);
775 param1 = deltacenter*sina*sina;
777 else if (sint < Tolang) { // plan perpendiculaire a l axe
778 typeres = IntAna_Circle;
781 dir1 = Co.Position().Direction();
782 dir2 = Co.Position().XDirection();
783 param1 = apex.Distance(center)*Abs(Tan(angl));
785 else if (cost < sina ) {
786 typeres = IntAna_Hyperbola;
790 deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
791 pt1.SetXYZ(center.XYZ() - deltacenter*axex);
795 param1 = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
796 param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
799 else { // on a alors cost > sina
800 typeres = IntAna_Ellipse;
802 Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
803 deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
805 pt1.SetXYZ(center.XYZ() + deltacenter*axex);
809 param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
814 //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
815 //-- des hyperboles trop bizarres
816 //-- On retourne False -> Traitement par biparametree
817 static Standard_Real EllipseLimit = 1.0E+9; //OCC513(apo) 1000000
818 static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
819 if(typeres==IntAna_Ellipse && nbint>=1) {
820 if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit) {
825 if(typeres==IntAna_Hyperbola && nbint>=2) {
826 if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit) {
827 done = Standard_False;
831 if(typeres==IntAna_Hyperbola && nbint>=1) {
832 if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit) {
838 done = Standard_True;
841 //=======================================================================
842 //function : IntAna_QuadQuadGeo
843 //purpose : Pln Sphere
844 //=======================================================================
845 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
847 : done(Standard_False),
849 typeres(IntAna_Empty),
860 myCommonGen(Standard_False),
866 //=======================================================================
869 //=======================================================================
870 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
874 done = Standard_False;
875 Standard_Real A,B,C,D,dist, radius;
879 // debug JAG : on met typeres = IntAna_Empty par defaut...
880 typeres = IntAna_Empty;
882 P.Coefficients(A,B,C,D);
883 S.Location().Coord(X,Y,Z);
886 dist = A * X + B * Y + C * Z + D;
888 if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
889 // on a une seule solution : le point projection du centre de la sphere
892 typeres = IntAna_Point;
893 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
895 else if (Abs(dist) < radius) {
896 // on a un cercle solution
898 typeres = IntAna_Circle;
899 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
900 dir1 = P.Axis().Direction();
901 if(P.Direct()==Standard_False) dir1.Reverse();
902 dir2 = P.Position().XDirection();
903 param1 = Sqrt(radius*radius - dist*dist);
905 param2bis=0.0; //-- pour eviter param2bis not used ....
906 done = Standard_True;
909 //=======================================================================
910 //function : IntAna_QuadQuadGeo
911 //purpose : Cylinder - Cylinder
912 //=======================================================================
913 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
914 const gp_Cylinder& Cyl2,
915 const Standard_Real Tol)
916 : done(Standard_False),
918 typeres(IntAna_Empty),
929 myCommonGen(Standard_False),
933 Perform(Cyl1,Cyl2,Tol);
935 //=======================================================================
938 //=======================================================================
939 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
940 const gp_Cylinder& Cyl2,
941 const Standard_Real Tol)
944 //---------------------------- Parallel axes -------------------------
945 AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
946 Standard_Real R1=Cyl1.Radius();
947 Standard_Real R2=Cyl2.Radius();
948 Standard_Real RmR, RmR_Relative;
949 RmR=(R1>R2)? (R1-R2) : (R2-R1);
952 Rmax=(R1>R2)? R1 : R2;
953 RmR_Relative=RmR/Rmax;
956 Standard_Real DistA1A2=A1A2.Distance();
958 if(A1A2.Parallel()) {
964 typeres=IntAna_Empty;
967 else { //-- DistA1A2 > Tol
968 gp_Pnt P1=Cyl1.Location();
969 gp_Pnt P2t=Cyl2.Location();
971 //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
972 gp_Dir DirCyl = Cyl1.Position().Direction();
973 Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
975 P2.SetCoord( P2t.X() - ProjP2OnDirCyl1*DirCyl.X()
976 ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
977 ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
979 Standard_Real R1pR2=R1+R2;
980 if(DistA1A2>(R1pR2+Tol)) {
981 typeres=IntAna_Empty;
984 else if(DistA1A2>(R1pR2)) {
985 //-- 1 Tangent line -------------------------------------OK
990 Standard_Real R1_R1pR2=R1/R1pR2;
991 pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
992 ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
993 ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
996 else if(DistA1A2>RmR) {
997 //-- 2 lines ---------------------------------------------OK
1002 gp_Dir DirA1A2=gp_Dir(P1P2);
1003 gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
1005 Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);
1007 // Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
1008 Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
1009 Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
1011 if((Beta+Beta)<Tol) {
1013 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
1014 ,P1.Y() + Alpha*DirA1A2.Y()
1015 ,P1.Z() + Alpha*DirA1A2.Z());
1018 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
1019 ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
1020 ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
1022 pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
1023 ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
1024 ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
1027 else if(DistA1A2>(RmR-Tol)) {
1028 //-- 1 Tangent ------------------------------------------OK
1029 typeres=IntAna_Line;
1032 Standard_Real R1_RmR=R1/RmR;
1034 if(R1 < R2) R1_RmR = -R1_RmR;
1036 pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
1037 ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
1038 ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
1042 typeres=IntAna_Empty;
1046 else { //-- No Parallel Axis ---------------------------------OK
1047 if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
1048 && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE)) {
1049 //-- PI/2 between the two axis and Intersection
1050 //-- and identical radius
1051 typeres=IntAna_Ellipse;
1053 gp_Dir DirCyl1=Cyl1.Position().Direction();
1054 gp_Dir DirCyl2=Cyl2.Position().Direction();
1055 pt1=pt2=A1A2.PtIntersect();
1057 Standard_Real A=DirCyl1.Angle(DirCyl2);
1059 B=Abs(Sin(0.5*(M_PI-A)));
1062 if(A==0.0 || B==0.0) {
1063 typeres=IntAna_Same;
1068 gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
1069 dir1 = gp_Dir(dircyl1.Added(dircyl2));
1070 dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
1072 param2 = Cyl1.Radius() / A;
1073 param1 = Cyl1.Radius() / B;
1074 param2bis= param1bis = Cyl1.Radius();
1075 if(param1 < param1bis) {
1076 A=param1; param1=param1bis; param1bis=A;
1078 if(param2 < param2bis) {
1079 A=param2; param2=param2bis; param2bis=A;
1083 if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol) {
1084 typeres = IntAna_Point;
1085 Standard_Real d,p1,p2;
1087 gp_Dir D1 = Cyl1.Axis().Direction();
1088 gp_Dir D2 = Cyl2.Axis().Direction();
1089 A1A2.Distance(d,p1,p2);
1090 gp_Pnt P = Cyl1.Axis().Location();
1091 gp_Pnt P1(P.X() - p1*D1.X(),
1094 P = Cyl2.Axis().Location();
1095 gp_Pnt P2(P.X() - p2*D2.X(),
1101 pt1.SetCoord(P1.X() + p1*D1.X(),
1103 P1.Z() + p1*D1.Z());
1107 typeres=IntAna_NoGeometricSolution;
1112 //=======================================================================
1113 //function : IntAna_QuadQuadGeo
1114 //purpose : Cylinder - Cone
1115 //=======================================================================
1116 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1118 const Standard_Real Tol)
1119 : done(Standard_False),
1121 typeres(IntAna_Empty),
1132 myCommonGen(Standard_False),
1136 Perform(Cyl,Con,Tol);
1138 //=======================================================================
1139 //function : Perform
1141 //=======================================================================
1142 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1144 const Standard_Real )
1147 AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1149 gp_Pnt Pt=Con.Apex();
1150 Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1151 gp_Dir dir=Cyl.Position().Direction();
1152 pt1.SetCoord( Pt.X() + dist*dir.X()
1153 ,Pt.Y() + dist*dir.Y()
1154 ,Pt.Z() + dist*dir.Z());
1155 pt2.SetCoord( Pt.X() - dist*dir.X()
1156 ,Pt.Y() - dist*dir.Y()
1157 ,Pt.Z() - dist*dir.Z());
1159 param1=param2=Cyl.Radius();
1161 typeres=IntAna_Circle;
1165 typeres=IntAna_NoGeometricSolution;
1168 //=======================================================================
1170 //purpose : Cylinder - Sphere
1171 //=======================================================================
1172 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1173 const gp_Sphere& Sph,
1174 const Standard_Real Tol)
1175 : done(Standard_False),
1177 typeres(IntAna_Empty),
1188 myCommonGen(Standard_False),
1192 Perform(Cyl,Sph,Tol);
1194 //=======================================================================
1195 //function : Perform
1197 //=======================================================================
1198 void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1199 ,const gp_Sphere& Sph
1200 ,const Standard_Real)
1203 gp_Pnt Pt=Sph.Location();
1204 AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1205 if((A1A2.Intersect() && Pt.Distance(A1A2.PtIntersect())==0.0 )
1207 if(Sph.Radius() < Cyl.Radius()) {
1208 typeres = IntAna_Empty;
1211 Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1212 gp_Dir dir=Cyl.Position().Direction();
1214 typeres=IntAna_Circle;
1215 pt1.SetCoord( Pt.X() + dist*dir.X()
1216 ,Pt.Y() + dist*dir.Y()
1217 ,Pt.Z() + dist*dir.Z());
1219 param1 = Cyl.Radius();
1220 if(dist>RealEpsilon()) {
1221 pt2.SetCoord( Pt.X() - dist*dir.X()
1222 ,Pt.Y() - dist*dir.Y()
1223 ,Pt.Z() - dist*dir.Z());
1224 param2=Cyl.Radius();
1230 typeres=IntAna_NoGeometricSolution;
1234 //=======================================================================
1235 //function : IntAna_QuadQuadGeo
1236 //purpose : Cone - Cone
1237 //=======================================================================
1238 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1239 const gp_Cone& Con2,
1240 const Standard_Real Tol)
1241 : done(Standard_False),
1243 typeres(IntAna_Empty),
1254 myCommonGen(Standard_False),
1258 Perform(Con1,Con2,Tol);
1261 //=======================================================================
1262 //function : Perform
1264 //=======================================================================
1265 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1266 const gp_Cone& Con2,
1267 const Standard_Real Tol)
1271 Standard_Real tg1, tg2, aDA1A2, aTol2;
1272 gp_Pnt aPApex1, aPApex2;
1274 Standard_Real TOL_APEX_CONF = 1.e-10;
1277 tg1=Tan(Con1.SemiAngle());
1278 tg2=Tan(Con2.SemiAngle());
1280 if((tg1 * tg2) < 0.) {
1285 aPApex1=Con1.Apex();
1286 aPApex2=Con2.Apex();
1287 aDA1A2=aPApex1.SquareDistance(aPApex2);
1289 AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1295 gp_Pnt P=Con1.Apex();
1296 gp_Dir D=Con1.Position().Direction();
1297 Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1299 if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
1300 if (fabs(d) < TOL_APEX_CONF) {
1301 typeres = IntAna_Point;
1306 x=(d*tg2)/(tg1+tg2);
1307 pt1.SetCoord( P.X() + x*D.X()
1312 x=(d*tg2)/(tg2-tg1);
1313 pt2.SetCoord( P.X() + x*D.X()
1319 typeres=IntAna_Circle;
1322 if (fabs(d) < TOL_APEX_CONF) {
1323 typeres=IntAna_Same;
1326 typeres=IntAna_Circle;
1329 pt1.SetCoord( P.X() + x*D.X()
1332 param1 = Abs(x * tg1);
1336 } //-- fin A1A2.Same
1338 else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1339 //-- voir AnVer12mai98
1340 Standard_Real DistA1A2=A1A2.Distance();
1341 gp_Dir DA1=Con1.Position().Direction();
1342 gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1343 gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1344 Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1346 gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
1347 O1O2n.Y()-O1O2_DA1*DA1.Y(),
1348 O1O2n.Z()-O1O2_DA1*DA1.Z());
1349 gp_Dir DB1=gp_Dir(O1_Proj_A2);
1351 Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1352 Standard_Real ABSTG1 = Abs(tg1);
1353 Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1354 Standard_Real X1 = X2+yO1O2;
1356 gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1357 Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1358 Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1360 gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1361 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1362 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1363 gp_Vec P1MO1O2(P1,MO1O2);
1365 gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1366 gp_Dir OrthoPln = DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1368 IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1369 if(INTER_QUAD_PLN.IsDone()) {
1370 switch(INTER_QUAD_PLN.TypeInter()) {
1371 case IntAna_Ellipse: {
1372 typeres=IntAna_Ellipse;
1373 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1375 dir1 = E.Position().Direction();
1376 dir2 = E.Position().XDirection();
1377 param1 = E.MajorRadius();
1378 param1bis = E.MinorRadius();
1382 case IntAna_Circle: {
1383 typeres=IntAna_Circle;
1384 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1386 dir1 = C.Position().XDirection();
1387 dir2 = C.Position().YDirection();
1388 param1 = C.Radius();
1392 case IntAna_Hyperbola: {
1393 typeres=IntAna_Hyperbola;
1394 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1395 pt1 = pt2 = H.Location();
1396 dir1 = H.Position().Direction();
1397 dir2 = H.Position().XDirection();
1398 param1 = param2 = H.MajorRadius();
1399 param1bis = param2bis = H.MinorRadius();
1404 typeres=IntAna_Line;
1405 gp_Lin H=INTER_QUAD_PLN.Line(1);
1406 pt1 = pt2 = H.Location();
1407 dir1 = dir2 = H.Position().Direction();
1408 param1 = param2 = 0.0;
1409 param1bis = param2bis = 0.0;
1414 typeres=IntAna_NoGeometricSolution;
1417 }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1419 else if (aDA1A2<aTol2) {
1421 // When apices are coinsided there can be 3 possible cases
1422 // 3.1 - empty solution (iRet=0)
1423 // 3.2 - one line when cone1 touches cone2 (iRet=1)
1424 // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1426 Standard_Integer iRet;
1427 Standard_Real aGamma, aBeta1, aBeta2;
1428 Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1429 Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1430 gp_Pnt2d aP0, aPA1, aP1, aPA2;
1434 // Preliminary analysis. Determination of iRet
1439 aPA1.SetCoord(aD1, 0.);
1440 aP0.SetCoord(0., 0.);
1444 aGamma=aAx1.Angle(aAx2);
1445 if (aGamma>aHalfPI){
1448 aCosGamma=Cos(aGamma);
1449 aSinGamma=Sin(aGamma);
1451 aBeta1=Con1.SemiAngle();
1452 aTgBeta1=Tan(aBeta1);
1453 aTgBeta1=Abs(aTgBeta1);
1455 aBeta2=Con2.SemiAngle();
1456 aTgBeta2=Tan(aBeta2);
1457 aTgBeta2=Abs(aTgBeta2);
1460 aP1.SetCoord(aD1, aR1);
1463 aVAx2.SetCoord(aCosGamma, aSinGamma);
1464 gp_Dir2d aDAx2(aVAx2);
1465 gp_Lin2d aLAx2(aP0, aDAx2);
1467 gp_Vec2d aV(aP0, aP1);
1469 aPA2=aP0.Translated(aDx*aDAx2);
1472 aDx=aPA2.Distance(aP0);
1476 aRD2=aPA2.Distance(aP1);
1478 if (aRD2>(aR2+Tol)) {
1480 typeres=IntAna_Empty; //nothing
1484 iRet=1; //touch case => 1 line
1485 if (aRD2<(aR2-Tol)) {
1486 iRet=2;//intersection => couple of lines
1489 // Finding the solution in 3D
1492 gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1493 gp_Dir aD3Ax1, aD3Ax2;
1495 IntAna_QuadQuadGeo aIntr;
1497 aQApex1=Con1.Apex();
1498 aD3Ax1=aAx1.Direction();
1499 aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1500 aQApex1.Y()+aD1*aD3Ax1.Y(),
1501 aQApex1.Z()+aD1*aD3Ax1.Z());
1503 aDx=aD3Ax1.Dot(aAx2.Direction());
1507 aD3Ax2=aAx2.Direction();
1509 aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1511 aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1512 aQApex1.Y()+aD2*aD3Ax2.Y(),
1513 aQApex1.Z()+aD2*aD3Ax2.Z());
1515 gp_Pln aPln1(aQA1, aD3Ax1);
1516 gp_Pln aPln2(aQA2, aD3Ax2);
1518 aIntr.Perform(aPln1, aPln2, Tol, Tol);
1519 if (!aIntr.IsDone()) {
1520 iRet=-1; // just in case. it must not be so
1521 typeres=IntAna_NoGeometricSolution;
1526 const gp_Dir& aDLin=aLin.Direction();
1527 gp_Vec aVLin(aDLin);
1528 gp_Pnt aOrig=aLin.Location();
1529 gp_Vec aVr(aQA1, aOrig);
1531 aQX=aOrig.Translated(aDx*aVLin);
1535 typeres=IntAna_Line;
1546 gp_Vec aVX(aQApex1, aQX);
1553 aDa=aQA1.Distance(aQX);
1554 aDx=sqrt(aR1*aR1-aDa*aDa);
1555 aQX1=aQX.Translated(aDx*aVLin);
1556 aQX2=aQX.Translated(-aDx*aVLin);
1560 gp_Vec aVX1(aQApex1, aQX1);
1562 gp_Vec aVX2(aQApex1, aQX2);
1565 } //else if (aDA1A2<aTol2) {
1566 //Case when cones have common generatrix
1567 else if(A1A2.Intersect()) {
1568 //Check if apex of one cone belongs another one
1569 Standard_Real u, v, tol2 = Tol*Tol;
1570 ElSLib::Parameters(Con2, aPApex1, u, v);
1571 gp_Pnt p = ElSLib::Value(u, v, Con2);
1572 if(aPApex1.SquareDistance(p) > tol2) {
1573 typeres=IntAna_NoGeometricSolution;
1577 ElSLib::Parameters(Con1, aPApex2, u, v);
1578 p = ElSLib::Value(u, v, Con1);
1579 if(aPApex2.SquareDistance(p) > tol2) {
1580 typeres=IntAna_NoGeometricSolution;
1584 //Cones have a common generatrix passing through apexes
1585 myCommonGen = Standard_True;
1587 //common generatrix of cones
1588 gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1590 //Intersection point of axes
1591 gp_Pnt aPAxeInt = A1A2.PtIntersect();
1593 //Characteristic point of intersection curve
1594 u = ElCLib::Parameter(aGen, aPAxeInt);
1595 myPChar = ElCLib::Value(u, aGen);
1598 //Other generatrixes of cones laying in maximal plane
1599 gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1600 gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
1602 //Intersection point of generatrixes
1603 gp_Dir aN; //solution plane normal
1604 gp_Dir aD1 = aGen1.Direction();
1606 gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1608 if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1609 aN = aD1.Crossed(aD2);
1611 else if(aGen1.SquareDistance(aGen2) > tol2) {
1612 //Something wrong ???
1613 typeres=IntAna_NoGeometricSolution;
1617 gp_Dir D1 = aGen1.Position().Direction();
1618 gp_Dir D2 = aGen2.Position().Direction();
1619 gp_Pnt O1 = aGen1.Location();
1620 gp_Pnt O2 = aGen2.Location();
1621 Standard_Real D1DotD2 = D1.Dot(D2);
1622 Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1623 gp_Vec O1O2 (O1,O2);
1624 Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1626 gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1628 aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1629 aN = aD1.Crossed(aD2);
1631 //Plane that must contain intersection curves
1632 gp_Pln anIntPln(myPChar, aN);
1634 IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1636 if(INTER_QUAD_PLN.IsDone()) {
1637 switch(INTER_QUAD_PLN.TypeInter()) {
1638 case IntAna_Ellipse: {
1639 typeres=IntAna_Ellipse;
1640 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1642 dir1 = E.Position().Direction();
1643 dir2 = E.Position().XDirection();
1644 param1 = E.MajorRadius();
1645 param1bis = E.MinorRadius();
1649 case IntAna_Circle: {
1650 typeres=IntAna_Circle;
1651 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1653 dir1 = C.Position().XDirection();
1654 dir2 = C.Position().YDirection();
1655 param1 = C.Radius();
1659 case IntAna_Parabola: {
1660 typeres=IntAna_Parabola;
1661 gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1662 pt1 = Prb.Location();
1663 dir1 = Prb.Position().Direction();
1664 dir2 = Prb.Position().XDirection();
1665 param1 = Prb.Focal();
1669 case IntAna_Hyperbola: {
1670 typeres=IntAna_Hyperbola;
1671 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1672 pt1 = pt2 = H.Location();
1673 dir1 = H.Position().Direction();
1674 dir2 = H.Position().XDirection();
1675 param1 = param2 = H.MajorRadius();
1676 param1bis = param2bis = H.MinorRadius();
1681 typeres=IntAna_NoGeometricSolution;
1687 typeres=IntAna_NoGeometricSolution;
1690 //=======================================================================
1691 //function : IntAna_QuadQuadGeo
1692 //purpose : Sphere - Cone
1693 //=======================================================================
1694 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1696 const Standard_Real Tol)
1697 : done(Standard_False),
1699 typeres(IntAna_Empty),
1710 myCommonGen(Standard_False),
1714 Perform(Sph,Con,Tol);
1716 //=======================================================================
1717 //function : Perform
1719 //=======================================================================
1720 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1722 const Standard_Real)
1728 AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1729 gp_Pnt Pt=Sph.Location();
1731 if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1733 gp_Pnt ConApex= Con.Apex();
1734 Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1736 if(dApexSphCenter>RealEpsilon()) {
1737 ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1740 ConDir = Con.Position().Direction();
1743 Standard_Real Rad=Sph.Radius();
1744 Standard_Real tga=Tan(Con.SemiAngle());
1748 //-- x: Roots of (x**2 + y**2 = Rad**2)
1749 //-- tga = y / (x+dApexSphCenter)
1750 Standard_Real tgatga = tga * tga;
1751 math_DirectPolynomialRoots Eq( 1.0+tgatga
1752 ,2.0*tgatga*dApexSphCenter
1753 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1755 Standard_Integer nbsol=Eq.NbSolutions();
1757 typeres=IntAna_Empty;
1760 typeres=IntAna_Circle;
1762 Standard_Real x = Eq.Value(1);
1763 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1765 pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1766 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1767 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1768 param1 = tga * dApexSphCenterpx;
1769 param1 = Abs(param1);
1771 if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1772 typeres=IntAna_PointAndCircle;
1777 Standard_Real x=Eq.Value(2);
1778 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1780 pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1781 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1782 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1783 param2 = tga * dApexSphCenterpx;
1784 param2 = Abs(param2);
1786 if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1787 typeres=IntAna_PointAndCircle;
1794 done=Standard_False;
1798 typeres=IntAna_NoGeometricSolution;
1802 //=======================================================================
1803 //function : IntAna_QuadQuadGeo
1804 //purpose : Sphere - Sphere
1805 //=======================================================================
1806 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1807 ,const gp_Sphere& Sph2
1808 ,const Standard_Real Tol)
1809 : done(Standard_False),
1811 typeres(IntAna_Empty),
1822 myCommonGen(Standard_False),
1826 Perform(Sph1,Sph2,Tol);
1828 //=======================================================================
1829 //function : Perform
1831 //=======================================================================
1832 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1833 const gp_Sphere& Sph2,
1834 const Standard_Real Tol)
1837 gp_Pnt O1=Sph1.Location();
1838 gp_Pnt O2=Sph2.Location();
1839 Standard_Real dO1O2=O1.Distance(O2);
1840 Standard_Real R1=Sph1.Radius();
1841 Standard_Real R2=Sph2.Radius();
1842 Standard_Real Rmin,Rmax;
1843 typeres=IntAna_Empty;
1844 param2bis=0.0; //-- pour eviter param2bis not used ....
1846 if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1848 if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1849 typeres = IntAna_Same;
1852 if(dO1O2<=Tol) { return; }
1853 gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1854 Standard_Real t = Rmax - dO1O2 - Rmin;
1856 //----------------------------------------------------------------------
1857 //-- |----------------- R1 --------------------|
1858 //-- |----dO1O2-----|-----------R2----------|
1861 //-- |------ R1 ------|---------dO1O2----------|
1862 //-- |-------------------R2-----------------------|
1864 //----------------------------------------------------------------------
1865 if(t >= 0.0 && t <=Tol) {
1866 typeres = IntAna_Point;
1869 if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1870 else t2=(-R1+(dO1O2-R2))*0.5;
1872 pt1.SetCoord( O1.X() + t2*Dir.X()
1873 ,O1.Y() + t2*Dir.Y()
1874 ,O1.Z() + t2*Dir.Z());
1877 //-----------------------------------------------------------------
1878 //-- |----------------- dO1O2 --------------------|
1879 //-- |----R1-----|-----------R2----------|-Tol-|
1881 //-- |----------------- Rmax --------------------|
1882 //-- |----Rmin----|-------dO1O2-------|-Tol-|
1884 //-----------------------------------------------------------------
1885 if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1886 typeres=IntAna_Empty;
1889 //---------------------------------------------------------------
1892 //---------------------------------------------------------------
1893 Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1894 Standard_Real Beta = R1*R1-Alpha*Alpha;
1895 Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1897 if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1898 typeres = IntAna_Point;
1899 Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1902 typeres = IntAna_Circle;
1906 pt1.SetCoord( O1.X() + Alpha*Dir.X()
1907 ,O1.Y() + Alpha*Dir.Y()
1908 ,O1.Z() + Alpha*Dir.Z());
1916 //=======================================================================
1917 //function : IntAna_QuadQuadGeo
1918 //purpose : Plane - Torus
1919 //=======================================================================
1920 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
1921 const gp_Torus& Tor,
1922 const Standard_Real Tol)
1923 : done(Standard_False),
1925 typeres(IntAna_Empty),
1936 myCommonGen(Standard_False),
1940 Perform(Pln,Tor,Tol);
1942 //=======================================================================
1943 //function : Perform
1945 //=======================================================================
1946 void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
1947 const gp_Torus& Tor,
1948 const Standard_Real Tol)
1950 done = Standard_True;
1952 Standard_Real aRMin, aRMaj;
1954 aRMin = Tor.MinorRadius();
1955 aRMaj = Tor.MajorRadius();
1956 if (aRMin >= aRMaj) {
1957 typeres = IntAna_NoGeometricSolution;
1961 const gp_Ax1 aPlnAx = Pln.Axis();
1962 const gp_Ax1 aTorAx = Tor.Axis();
1964 Standard_Boolean bParallel, bNormal;
1966 bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
1967 bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
1968 if (!bNormal && !bParallel) {
1969 typeres = IntAna_NoGeometricSolution;
1973 Standard_Real aDist;
1975 gp_Pnt aTorLoc = aTorAx.Location();
1977 Standard_Real aDt, X, Y, Z, A, B, C, D;
1979 Pln.Coefficients(A,B,C,D);
1980 aTorLoc.Coord(X,Y,Z);
1981 aDist = A*X + B*Y + C*Z + D;
1983 if ((Abs(aDist) - aRMin) > Tol) {
1984 typeres=IntAna_Empty;
1988 typeres = IntAna_Circle;
1990 pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
1991 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
1992 param1 = aRMaj + aDt;
1993 dir1 = aTorAx.Direction();
1995 if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
1997 param2 = aRMaj - aDt;
2004 aDist = Pln.Distance(aTorLoc);
2005 if (aDist > myEPSILON_DISTANCE) {
2006 typeres = IntAna_NoGeometricSolution;
2010 typeres = IntAna_Circle;
2011 param2 = param1 = aRMin;
2012 dir2 = dir1 = aPlnAx.Direction();
2015 gp_Dir aDir = aTorAx.Direction()^dir1;
2016 pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
2017 pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
2021 //=======================================================================
2022 //function : IntAna_QuadQuadGeo
2023 //purpose : Cylinder - Torus
2024 //=======================================================================
2025 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
2026 const gp_Torus& Tor,
2027 const Standard_Real Tol)
2028 : done(Standard_False),
2030 typeres(IntAna_Empty),
2041 myCommonGen(Standard_False),
2045 Perform(Cyl,Tor,Tol);
2047 //=======================================================================
2048 //function : Perform
2050 //=======================================================================
2051 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
2052 const gp_Torus& Tor,
2053 const Standard_Real Tol)
2055 done = Standard_True;
2057 Standard_Real aRMin, aRMaj;
2059 aRMin = Tor.MinorRadius();
2060 aRMaj = Tor.MajorRadius();
2061 if (aRMin >= aRMaj) {
2062 typeres = IntAna_NoGeometricSolution;
2066 const gp_Ax1 aCylAx = Cyl.Axis();
2067 const gp_Ax1 aTorAx = Tor.Axis();
2069 const gp_Lin aLin(aTorAx);
2070 const gp_Pnt aLocCyl = Cyl.Location();
2072 if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
2073 (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
2074 typeres = IntAna_NoGeometricSolution;
2078 Standard_Real aRCyl;
2080 aRCyl = Cyl.Radius();
2081 if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2082 typeres = IntAna_Empty;
2086 typeres = IntAna_Circle;
2088 Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2089 gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2091 dir1 = aTorAx.Direction();
2092 pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2095 if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2096 (aRCyl < (aRMaj + aRMin))) {
2098 pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2104 //=======================================================================
2105 //function : IntAna_QuadQuadGeo
2106 //purpose : Cone - Torus
2107 //=======================================================================
2108 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2109 const gp_Torus& Tor,
2110 const Standard_Real Tol)
2111 : done(Standard_False),
2113 typeres(IntAna_Empty),
2124 myCommonGen(Standard_False),
2128 Perform(Con,Tor,Tol);
2130 //=======================================================================
2131 //function : Perform
2133 //=======================================================================
2134 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2135 const gp_Torus& Tor,
2136 const Standard_Real Tol)
2138 done = Standard_True;
2140 Standard_Real aRMin, aRMaj;
2142 aRMin = Tor.MinorRadius();
2143 aRMaj = Tor.MajorRadius();
2144 if (aRMin >= aRMaj) {
2145 typeres = IntAna_NoGeometricSolution;
2149 const gp_Ax1 aConAx = Con.Axis();
2150 const gp_Ax1 aTorAx = Tor.Axis();
2152 const gp_Lin aLin(aTorAx);
2153 const gp_Pnt aConApex = Con.Apex();
2155 if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2156 (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2157 typeres = IntAna_NoGeometricSolution;
2161 Standard_Real anAngle, aDist, aParam[4], aDt;
2163 gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2166 anAngle = Con.SemiAngle();
2167 aTorLoc = aTorAx.Location();
2169 aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2170 gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2171 gp_Ax1 anAxCLRot(aConApex, aDN);
2172 gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2173 gp_Dir aDL = aConL.Position().Direction();
2174 gp_Dir aXDir = Tor.XAxis().Direction();
2176 typeres = IntAna_Empty;
2178 for (i = 0; i < 2; ++i) {
2182 aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2184 aDist = aConL.Distance(aPCT);
2185 if (aDist > aRMin+Tol) {
2189 typeres = IntAna_Circle;
2191 gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
2192 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2195 gp_XYZ aDVal = aDt*aDL.XYZ();
2196 aP.SetXYZ(aPh + aDVal);
2197 aParam[nbint] = aLin.Distance(aP);
2198 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2199 aDir[nbint] = aTorAx.Direction();
2201 if ((aDist < aRMin) && (aDt > Tol)) {
2202 aP.SetXYZ(aPh - aDVal);
2203 aParam[nbint] = aLin.Distance(aP);
2204 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2205 aDir[nbint] = aDir[nbint-1];
2210 for (i = 0; i < nbint; ++i) {
2242 //=======================================================================
2243 //function : IntAna_QuadQuadGeo
2244 //purpose : Sphere - Torus
2245 //=======================================================================
2246 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2247 const gp_Torus& Tor,
2248 const Standard_Real Tol)
2249 : done(Standard_False),
2251 typeres(IntAna_Empty),
2262 myCommonGen(Standard_False),
2266 Perform(Sph,Tor,Tol);
2268 //=======================================================================
2269 //function : Perform
2271 //=======================================================================
2272 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2273 const gp_Torus& Tor,
2274 const Standard_Real Tol)
2276 done = Standard_True;
2278 Standard_Real aRMin, aRMaj;
2280 aRMin = Tor.MinorRadius();
2281 aRMaj = Tor.MajorRadius();
2282 if (aRMin >= aRMaj) {
2283 typeres = IntAna_NoGeometricSolution;
2287 const gp_Ax1 aTorAx = Tor.Axis();
2288 const gp_Lin aLin(aTorAx);
2289 const gp_Pnt aSphLoc = Sph.Location();
2291 if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2292 typeres = IntAna_NoGeometricSolution;
2296 Standard_Real aRSph, aDist;
2299 gp_Dir aXDir = Tor.XAxis().Direction();
2300 aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2301 aRSph = Sph.Radius();
2303 gp_Vec aVec12(aTorLoc, aSphLoc);
2304 aDist = aVec12.Magnitude();
2305 if (((aDist - Tol) > (aRMin + aRSph)) ||
2306 ((aDist + Tol) < Abs(aRMin - aRSph))) {
2307 typeres = IntAna_Empty;
2311 typeres = IntAna_Circle;
2313 Standard_Real anAlpha, aBeta;
2315 anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2316 aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2318 gp_Dir aDir12(aVec12);
2319 gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2320 gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2323 gp_XYZ aDVal = aBeta*aDC.XYZ();
2324 aP.SetXYZ(aPh + aDVal);
2325 param1 = aLin.Distance(aP);
2326 pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2327 dir1 = aTorAx.Direction();
2329 if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) &&
2330 (aDVal.Modulus() > Tol)) {
2331 aP.SetXYZ(aPh - aDVal);
2332 param2 = aLin.Distance(aP);
2333 pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2339 //=======================================================================
2340 //function : IntAna_QuadQuadGeo
2341 //purpose : Torus - Torus
2342 //=======================================================================
2343 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2344 const gp_Torus& Tor2,
2345 const Standard_Real Tol)
2346 : done(Standard_False),
2348 typeres(IntAna_Empty),
2359 myCommonGen(Standard_False),
2363 Perform(Tor1,Tor2,Tol);
2365 //=======================================================================
2366 //function : Perform
2368 //=======================================================================
2369 void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2370 const gp_Torus& Tor2,
2371 const Standard_Real Tol)
2373 done = Standard_True;
2375 Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2377 aRMin1 = Tor1.MinorRadius();
2378 aRMaj1 = Tor1.MajorRadius();
2379 aRMin2 = Tor2.MinorRadius();
2380 aRMaj2 = Tor2.MajorRadius();
2381 if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2382 typeres = IntAna_NoGeometricSolution;
2386 const gp_Ax1 anAx1 = Tor1.Axis();
2387 const gp_Ax1 anAx2 = Tor2.Axis();
2390 if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2391 (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2392 typeres = IntAna_NoGeometricSolution;
2396 gp_Pnt aLoc1, aLoc2;
2398 aLoc1 = anAx1.Location();
2399 aLoc2 = anAx2.Location();
2401 if (aLoc1.IsEqual(aLoc2, Tol) &&
2402 (Abs(aRMin1 - aRMin2) <= Tol) &&
2403 (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2404 typeres = IntAna_Same;
2408 Standard_Real aDist;
2411 gp_Dir aXDir1 = Tor1.XAxis().Direction();
2412 aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2413 aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2415 gp_Vec aV12(aP1, aP2);
2416 aDist = aV12.Magnitude();
2417 if (((aDist - Tol) > (aRMin1 + aRMin2)) ||
2418 ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2419 typeres = IntAna_Empty;
2423 typeres = IntAna_Circle;
2425 Standard_Real anAlpha, aBeta;
2427 anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2428 aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2430 gp_Dir aDir12(aV12);
2431 gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2432 gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2435 gp_XYZ aDVal = aBeta*aDC.XYZ();
2436 aP.SetXYZ(aPh + aDVal);
2437 param1 = aL1.Distance(aP);
2438 pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2439 dir1 = anAx1.Direction();
2441 if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) &&
2442 aDVal.Modulus() > Tol) {
2443 aP.SetXYZ(aPh - aDVal);
2444 param2 = aL1.Distance(aP);
2445 pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2451 //=======================================================================
2453 //purpose : Returns a Point
2454 //=======================================================================
2455 gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
2457 if(!done) { StdFail_NotDone::Raise(); }
2458 if(n>nbint || n<1) { Standard_DomainError::Raise(); }
2459 if(typeres==IntAna_PointAndCircle) {
2460 if(n!=1) { Standard_DomainError::Raise(); }
2461 if(param1==0.0) return(pt1);
2464 else if(typeres==IntAna_Point) {
2465 if(n==1) return(pt1);
2469 // WNT (what can you expect from MicroSoft ?)
2470 return gp_Pnt(0,0,0);
2472 //=======================================================================
2474 //purpose : Returns a Line
2475 //=======================================================================
2476 gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
2478 if(!done) { StdFail_NotDone::Raise(); }
2479 if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2480 Standard_DomainError::Raise();
2482 if(n==1) { return(gp_Lin(pt1,dir1)); }
2483 else { return(gp_Lin(pt2,dir2)); }
2485 //=======================================================================
2487 //purpose : Returns a Circle
2488 //=======================================================================
2489 gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
2491 if(!done) { StdFail_NotDone::Raise(); }
2492 if(typeres==IntAna_PointAndCircle) {
2493 if(n!=1) { Standard_DomainError::Raise(); }
2494 if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2495 return(gp_Circ(DirToAx2(pt2,dir2),param2));
2497 else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2498 Standard_DomainError::Raise();
2500 if (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2501 else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2502 else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2503 else { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
2506 //=======================================================================
2507 //function : Ellipse
2508 //purpose : Returns a Elips
2509 //=======================================================================
2510 gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2512 if(!done) { StdFail_NotDone::Raise(); }
2513 if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2514 Standard_DomainError::Raise();
2518 Standard_Real R1=param1, R2=param1bis, aTmp;
2520 aTmp=R1; R1=R2; R2=aTmp;
2522 gp_Ax2 anAx2(pt1, dir1 ,dir2);
2523 gp_Elips anElips (anAx2, R1, R2);
2527 Standard_Real R1=param2, R2=param2bis, aTmp;
2529 aTmp=R1; R1=R2; R2=aTmp;
2531 gp_Ax2 anAx2(pt2, dir2 ,dir1);
2532 gp_Elips anElips (anAx2, R1, R2);
2536 //=======================================================================
2537 //function : Parabola
2538 //purpose : Returns a Parabola
2539 //=======================================================================
2540 gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
2543 StdFail_NotDone::Raise();
2545 if (typeres!=IntAna_Parabola) {
2546 Standard_DomainError::Raise();
2548 if((n>nbint) || (n!=1)) {
2549 Standard_OutOfRange::Raise();
2551 return(gp_Parab(gp_Ax2( pt1
2556 //=======================================================================
2557 //function : Hyperbola
2558 //purpose : Returns a Hyperbola
2559 //=======================================================================
2560 gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
2563 StdFail_NotDone::Raise();
2565 if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2566 Standard_DomainError::Raise();
2569 return(gp_Hypr(gp_Ax2( pt1
2572 ,param1,param1bis));
2575 return(gp_Hypr(gp_Ax2( pt2
2578 ,param2,param2bis));
2581 //=======================================================================
2582 //function : HasCommonGen
2584 //=======================================================================
2585 Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2589 //=======================================================================
2592 //=======================================================================
2593 const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2597 //=======================================================================
2598 //function : RefineDir
2600 //=======================================================================
2601 void RefineDir(gp_Dir& aDir)
2603 Standard_Integer k, m, n;
2604 Standard_Real aC[3];
2606 aDir.Coord(aC[0], aC[1], aC[2]);
2610 for (k=0; k<3; ++k) {
2611 if (aC[k]==1. || aC[k]==-1.) {
2614 else if (aC[k]!=0.) {
2620 Standard_Real aEps, aR1, aR2, aNum;
2626 for (k=0; k<3; ++k) {
2628 aNum=(m)? aC[k] : -aC[k];
2629 if (aNum>aR1 && aNum<aR2) {
2642 aDir.SetCoord(aC[0], aC[1], aC[2]);