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
9 // under the terms of the GNU Lesser General Public 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,Standard_Real& Param1,Standard_Real& Param2)
197 gp_Vec O1O2(Axe1.Location(),Axe2.Location()); //-----------------------------
198 gp_Dir U1 = Axe1.Direction(); //-- juste pour voir.
199 gp_Dir U2 = Axe2.Direction();
201 gp_Dir N = U1.Crossed(U2);
202 Standard_Real D = Det33(U1.X(),U2.X(),N.X(),
204 U1.Z(),U2.Z(),N.Z());
206 dist = Det33(U1.X(),U2.X(),O1O2.X(),
207 U1.Y(),U2.Y(),O1O2.Y(),
208 U1.Z(),U2.Z(),O1O2.Z()) / D;
209 Param1 = Det33(O1O2.X(),U2.X(),N.X(),
210 O1O2.Y(),U2.Y(),N.Y(),
211 O1O2.Z(),U2.Z(),N.Z()) / (-D);
212 //------------------------------------------------------------
213 //-- On resout P1 * Dir1 + P2 * Dir2 + d * N = O1O2
214 //-- soit : Segment perpendiculaire : O1+P1 D1
216 Param2 = Det33(U1.X(),O1O2.X(),N.X(),
217 U1.Y(),O1O2.Y(),N.Y(),
218 U1.Z(),O1O2.Z(),N.Z()) / (D);
221 //=======================================================================
222 //function : DirToAx2
223 //purpose : returns a gp_Ax2 where D is the main direction
224 //=======================================================================
225 gp_Ax2 DirToAx2(const gp_Pnt& P,const gp_Dir& D)
227 Standard_Real x=D.X(); Standard_Real ax=Abs(x);
228 Standard_Real y=D.Y(); Standard_Real ay=Abs(y);
229 Standard_Real z=D.Z(); Standard_Real az=Abs(z);
230 if( (ax==0.0) || ((ax<ay) && (ax<az)) ) {
231 return(gp_Ax2(P,D,gp_Dir(gp_Vec(0.0,-z,y))));
233 else if( (ay==0.0) || ((ay<ax) && (ay<az)) ) {
234 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-z,0.0,x))));
237 return(gp_Ax2(P,D,gp_Dir(gp_Vec(-y,x,0.0))));
240 //=======================================================================
241 //function : IntAna_QuadQuadGeo
242 //purpose : Empty constructor
243 //=======================================================================
244 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(void)
245 : done(Standard_False),
247 typeres(IntAna_Empty),
258 myCommonGen(Standard_False),
263 //=======================================================================
264 //function : InitTolerances
266 //=======================================================================
267 void IntAna_QuadQuadGeo::InitTolerances()
269 myEPSILON_DISTANCE = 0.00000000000001;
270 myEPSILON_ANGLE_CONE = 0.000000000001;
271 myEPSILON_MINI_CIRCLE_RADIUS = 0.000000001;
272 myEPSILON_CYLINDER_DELTA_RADIUS = 0.0000000000001;
273 myEPSILON_CYLINDER_DELTA_DISTANCE= 0.0000001;
274 myEPSILON_AXES_PARA = 0.000000000001;
276 //=======================================================================
277 //function : IntAna_QuadQuadGeo
279 //=======================================================================
280 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P1,
282 const Standard_Real TolAng,
283 const Standard_Real Tol)
284 : done(Standard_False),
286 typeres(IntAna_Empty),
297 myCommonGen(Standard_False),
301 Perform(P1,P2,TolAng,Tol);
303 //=======================================================================
306 //=======================================================================
307 void IntAna_QuadQuadGeo::Perform (const gp_Pln& P1,
309 const Standard_Real TolAng,
310 const Standard_Real Tol)
316 Standard_Real A1 = 0., B1 = 0., C1 = 0., D1 = 0., A2 = 0., B2 = 0., C2 = 0., D2 = 0.;
317 P1.Coefficients(A1,B1,C1,D1);
318 P2.Coefficients(A2,B2,C2,D2);
320 gp_Vec vd(gp_Vec(A1,B1,C1).Crossed(gp_Vec(A2,B2,C2)));
321 Standard_Real dist1= A2*P1.Location().X() + B2*P1.Location().Y() + C2*P1.Location().Z() + D2;
322 Standard_Real dist2= A1*P2.Location().X() + B1*P2.Location().Y() + C1*P2.Location().Z() + D1;
324 if(vd.Magnitude() <=TolAng) {
325 // normalles are collinear - planes are same or parallel
326 typeres = (Abs(dist1) <= Tol && Abs(dist2) <= Tol) ? IntAna_Same : IntAna_Empty;
329 Standard_Real denom=A1*A2 + B1*B2 + C1*C2;
331 Standard_Real denom2 = denom*denom;
332 Standard_Real ddenom = 1. - denom2;
333 //denom = ( Abs(ddenom) <= 1.e-9 ) ? 1.e-9 : ddenom;
334 denom = ( Abs(ddenom) <= 1.e-16 ) ? 1.e-16 : ddenom;
336 Standard_Real par1 = dist1/denom;
337 Standard_Real par2 = -dist2/denom;
339 gp_Vec inter1(gp_Vec(A1,B1,C1).Crossed(vd));
340 gp_Vec inter2(gp_Vec(A2,B2,C2).Crossed(vd));
342 Standard_Real X1=P1.Location().X() + par1*inter1.X();
343 Standard_Real Y1=P1.Location().Y() + par1*inter1.Y();
344 Standard_Real Z1=P1.Location().Z() + par1*inter1.Z();
345 Standard_Real X2=P2.Location().X() + par2*inter2.X();
346 Standard_Real Y2=P2.Location().Y() + par2*inter2.Y();
347 Standard_Real Z2=P2.Location().Z() + par2*inter2.Z();
349 pt1=gp_Pnt((X1+X2)*0.5, (Y1+Y2)*0.5, (Z1+Z2)*0.5);
351 typeres = IntAna_Line;
357 //=======================================================================
358 //function : IntAna_QuadQuadGeo
359 //purpose : Pln Cylinder
360 //=======================================================================
361 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Pln& P
362 ,const gp_Cylinder& Cl
363 ,const Standard_Real Tolang
364 ,const Standard_Real Tol
365 ,const Standard_Real H)
366 : done(Standard_False),
368 typeres(IntAna_Empty),
379 myCommonGen(Standard_False),
383 Perform(P,Cl,Tolang,Tol,H);
385 //=======================================================================
388 //=======================================================================
389 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
390 ,const gp_Cylinder& Cl
391 ,const Standard_Real Tolang
392 ,const Standard_Real Tol
393 ,const Standard_Real H)
395 done = Standard_False;
396 Standard_Real dist,radius;
397 Standard_Real A,B,C,D;
399 Standard_Real sint,cost,h;
400 gp_XYZ axex,axey,omega;
404 radius = Cl.Radius();
406 gp_Lin axec(Cl.Axis());
407 gp_XYZ normp(P.Axis().Direction().XYZ());
409 P.Coefficients(A,B,C,D);
410 axec.Location().Coord(X,Y,Z);
411 dist = A*X + B*Y + C*Z + D; // la distance axe/plan est evaluee a l origine de l axe.
413 Standard_Real tolang = Tolang;
414 Standard_Boolean newparams = Standard_False;
416 gp_Vec ldv( axec.Direction() );
418 Standard_Real dA = Abs( ldv.Angle( npv ) );
421 Standard_Real dang = Abs( ldv.Angle( npv ) ) - M_PI/2.;
422 Standard_Real dangle = Abs( dang );
423 if( dangle > Tolang )
425 Standard_Real sinda = Abs( Sin( dangle ) );
426 Standard_Real dif = Abs( sinda - Tol );
430 newparams = Standard_True;
436 IntAna_IntConicQuad inter(axec,P,tolang,Tol,H);
438 if (inter.IsParallel()) {
439 // Le resultat de l intersection Plan-Cylindre est de type droite.
440 // il y a 1 ou 2 droites
442 typeres = IntAna_Line;
443 omega.SetCoord(X-dist*A,Y-dist*B,Z-dist*C);
445 if (Abs(Abs(dist)-radius) < Tol)
452 gp_XYZ omegaXYZ(X,Y,Z);
453 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
454 Standard_Real Xt,Yt,Zt,distt;
455 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
456 distt = A*Xt + B*Yt + C*Zt + D;
457 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
459 ppt1.SetXYZ( omega1 );
460 gp_Vec vv1(pt1,ppt1);
465 dir1 = axec.Direction();
467 else if (Abs(dist) < radius)
470 h = Sqrt(radius*radius - dist*dist);
471 axey = axec.Direction().XYZ().Crossed(normp); // axey est normalise
473 pt1.SetXYZ(omega - h*axey);
474 pt2.SetXYZ(omega + h*axey);
478 gp_XYZ omegaXYZ(X,Y,Z);
479 gp_XYZ omegaXYZtrnsl( omegaXYZ + 100.*axec.Direction().XYZ() );
480 Standard_Real Xt,Yt,Zt,distt,ht;
481 omegaXYZtrnsl.Coord(Xt,Yt,Zt);
482 distt = A*Xt + B*Yt + C*Zt + D;
483 // ht = Sqrt(radius*radius - distt*distt);
484 Standard_Real anSqrtArg = radius*radius - distt*distt;
485 ht = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
487 gp_XYZ omega1( omegaXYZtrnsl.X()-distt*A, omegaXYZtrnsl.Y()-distt*B, omegaXYZtrnsl.Z()-distt*C );
489 ppt1.SetXYZ( omega1 - ht*axey);
490 ppt2.SetXYZ( omega1 + ht*axey);
491 gp_Vec vv1(pt1,ppt1);
492 gp_Vec vv2(pt2,ppt2);
500 dir1 = axec.Direction();
501 dir2 = axec.Direction();
506 // debug JAG : le nbint = 0 doit etre remplace par typeres = IntAna_Empty
507 // et ne pas etre seulement supprime...
510 typeres = IntAna_Empty;
513 else { // Il y a un point d intersection. C est le centre du cercle
514 // ou de l ellipse solution.
517 axey = normp.Crossed(axec.Direction().XYZ());
518 sint = axey.Modulus();
520 pt1 = inter.Point(1);
522 if (sint < Tol/radius) {
524 // on construit un cercle avec comme axes X et Y ceux du cylindre
525 typeres = IntAna_Circle;
527 dir1 = axec.Direction(); // axe Z
528 dir2 = Cl.Position().XDirection();
533 // on construit un ellipse
534 typeres = IntAna_Ellipse;
535 cost = Abs(axec.Direction().XYZ().Dot(normp));
536 axex = axey.Crossed(normp);
538 dir1.SetXYZ(normp); //Modif ds ce bloc
541 param1 = radius/cost;
546 done = Standard_True;
548 //=======================================================================
549 //function : IntAna_QuadQuadGeo
551 //=======================================================================
552 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
554 const Standard_Real Tolang,
555 const Standard_Real Tol)
556 : done(Standard_False),
558 typeres(IntAna_Empty),
569 myCommonGen(Standard_False),
573 Perform(P,Co,Tolang,Tol);
575 //=======================================================================
578 //=======================================================================
579 void IntAna_QuadQuadGeo::Perform(const gp_Pln& P,
581 const Standard_Real Tolang,
582 const Standard_Real Tol)
585 done = Standard_False;
588 Standard_Real A,B,C,D;
590 Standard_Real dist,sint,cost,sina,cosa,angl,costa;
594 gp_Lin axec(Co.Axis());
595 P.Coefficients(A,B,C,D);
596 gp_Pnt apex(Co.Apex());
599 dist = A*X + B*Y + C*Z + D; // distance signee sommet du cone/ Plan
601 gp_XYZ normp = P.Axis().Direction().XYZ();
602 if(P.Direct()==Standard_False) { //-- lbr le 14 jan 97
606 axey = normp.Crossed(Co.Axis().Direction().XYZ());
607 axex = axey.Crossed(normp);
610 angl = Co.SemiAngle();
613 sina = Abs(Sin(angl));
616 // Angle entre la normale au plan et l axe du cone, ramene entre 0. et PI/2.
618 sint = axey.Modulus();
619 cost = Abs(Co.Axis().Direction().XYZ().Dot(normp));
621 // Le calcul de costa permet de determiner si le plan contient
622 // un generatrice du cone : on calcul Sin((PI/2. - t) - angl)
624 costa = cost*cosa - sint*sina; // sin((PI/2 -t)-angl)=cos(t+angl)
625 // avec t ramene entre 0 et pi/2.
627 if (Abs(dist) < Tol) {
628 // on considere que le plan contient le sommet du cone.
629 // les solutions possibles sont donc : 1 point, 1 droite, 2 droites
630 // selon l inclinaison du plan.
632 if (Abs(costa) < Tolang) { // plan parallele a la generatrice
633 typeres = IntAna_Line;
635 gp_XYZ ptonaxe(apex.XYZ() + 10.*(Co.Axis().Direction().XYZ()));
636 // point sur l axe du cone cote z positif
638 dist = A*ptonaxe.X() + B*ptonaxe.Y() + C*ptonaxe.Z() + D;
639 ptonaxe = ptonaxe - dist*normp;
641 dir1.SetXYZ(ptonaxe - pt1.XYZ());
643 else if (cost < sina) { // plan "interieur" au cone
644 typeres = IntAna_Line;
648 dh = Sqrt(sina*sina-cost*cost)/cosa;
649 dir1.SetXYZ(axex + dh*axey);
650 dir2.SetXYZ(axex - dh*axey);
652 else { // plan "exterieur" au cone
653 typeres = IntAna_Point;
659 // Solutions possibles : cercle, ellipse, parabole, hyperbole selon
660 // l inclinaison du plan.
661 Standard_Real deltacenter, distance;
664 // Le plan contient la direction de l axe du cone. La solution est
666 typeres = IntAna_Hyperbola;
668 pt1.SetXYZ(apex.XYZ()-dist*normp);
672 param1 = param2 = Abs(dist/Tan(angl));
673 param1bis = param2bis = Abs(dist);
677 IntAna_IntConicQuad inter(axec,P,Tolang); // on a necessairement 1 point.
679 gp_Pnt center(inter.Point(1));
681 // En fonction de la position de l intersection par rapport au sommet
682 // du cone, on change l axe x en -x et y en -y. Le parametre du sommet
683 // sur axec est negatif (voir definition du cone)
685 distance = apex.Distance(center);
687 if (inter.ParamOnConic(1) + Co.RefRadius()/Tan(angl) < 0.) {
692 if (Abs(costa) < Tolang) { // plan parallele a une generatrice
693 typeres = IntAna_Parabola;
695 deltacenter = distance/2./cosa;
697 pt1.SetXYZ(center.XYZ()-deltacenter*axex);
700 param1 = deltacenter*sina*sina;
702 else if (sint < Tolang) { // plan perpendiculaire a l axe
703 typeres = IntAna_Circle;
706 dir1 = Co.Position().Direction();
707 dir2 = Co.Position().XDirection();
708 param1 = apex.Distance(center)*Abs(Tan(angl));
710 else if (cost < sina ) {
711 typeres = IntAna_Hyperbola;
715 deltacenter = sint*sina*sina*distance/(sina*sina - cost*cost);
716 pt1.SetXYZ(center.XYZ() - deltacenter*axex);
720 param1 = param2 = cost*sina*cosa*distance /(sina*sina-cost*cost);
721 param1bis = param2bis = cost*sina*distance / Sqrt(sina*sina-cost*cost);
724 else { // on a alors cost > sina
725 typeres = IntAna_Ellipse;
727 Standard_Real radius = cost*sina*cosa*distance/(cost*cost-sina*sina);
728 deltacenter = sint*sina*sina*distance/(cost*cost-sina*sina);
730 pt1.SetXYZ(center.XYZ() + deltacenter*axex);
734 param1bis = cost*sina*distance/ Sqrt(cost*cost - sina*sina);
739 //-- On a du mal a gerer plus loin (Value ProjLib, Params ... )
740 //-- des hyperboles trop bizarres
741 //-- On retourne False -> Traitement par biparametree
742 static Standard_Real EllipseLimit = 1.0E+9; //OCC513(apo) 1000000
743 static Standard_Real HyperbolaLimit = 2.0E+6; //OCC537(apo) 50000
744 if(typeres==IntAna_Ellipse && nbint>=1) {
745 if(Abs(param1) > EllipseLimit || Abs(param1bis) > EllipseLimit) {
750 if(typeres==IntAna_Hyperbola && nbint>=2) {
751 if(Abs(param2) > HyperbolaLimit || Abs(param2bis) > HyperbolaLimit) {
752 done = Standard_False;
756 if(typeres==IntAna_Hyperbola && nbint>=1) {
757 if(Abs(param1) > HyperbolaLimit || Abs(param1bis) > HyperbolaLimit) {
763 done = Standard_True;
766 //=======================================================================
767 //function : IntAna_QuadQuadGeo
768 //purpose : Pln Sphere
769 //=======================================================================
770 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& P,
772 : done(Standard_False),
774 typeres(IntAna_Empty),
785 myCommonGen(Standard_False),
791 //=======================================================================
794 //=======================================================================
795 void IntAna_QuadQuadGeo::Perform( const gp_Pln& P
799 done = Standard_False;
800 Standard_Real A,B,C,D,dist, radius;
804 // debug JAG : on met typeres = IntAna_Empty par defaut...
805 typeres = IntAna_Empty;
807 P.Coefficients(A,B,C,D);
808 S.Location().Coord(X,Y,Z);
811 dist = A * X + B * Y + C * Z + D;
813 if (Abs( Abs(dist) - radius) < Epsilon(radius)) {
814 // on a une seule solution : le point projection du centre de la sphere
817 typeres = IntAna_Point;
818 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
820 else if (Abs(dist) < radius) {
821 // on a un cercle solution
823 typeres = IntAna_Circle;
824 pt1.SetCoord(X - dist*A, Y - dist*B, Z - dist*C);
825 dir1 = P.Axis().Direction();
826 if(P.Direct()==Standard_False) dir1.Reverse();
827 dir2 = P.Position().XDirection();
828 param1 = Sqrt(radius*radius - dist*dist);
830 param2bis=0.0; //-- pour eviter param2bis not used ....
831 done = Standard_True;
834 //=======================================================================
835 //function : IntAna_QuadQuadGeo
836 //purpose : Cylinder - Cylinder
837 //=======================================================================
838 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl1,
839 const gp_Cylinder& Cyl2,
840 const Standard_Real Tol)
841 : done(Standard_False),
843 typeres(IntAna_Empty),
854 myCommonGen(Standard_False),
858 Perform(Cyl1,Cyl2,Tol);
860 //=======================================================================
863 //=======================================================================
864 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
865 const gp_Cylinder& Cyl2,
866 const Standard_Real Tol)
869 //---------------------------- Parallel axes -------------------------
870 AxeOperator A1A2(Cyl1.Axis(),Cyl2.Axis());
871 Standard_Real R1=Cyl1.Radius();
872 Standard_Real R2=Cyl2.Radius();
873 Standard_Real RmR, RmR_Relative;
874 RmR=(R1>R2)? (R1-R2) : (R2-R1);
877 Rmax=(R1>R2)? R1 : R2;
878 RmR_Relative=RmR/Rmax;
881 Standard_Real DistA1A2=A1A2.Distance();
883 if(A1A2.Parallel()) {
889 typeres=IntAna_Empty;
892 else { //-- DistA1A2 > Tol
893 gp_Pnt P1=Cyl1.Location();
894 gp_Pnt P2t=Cyl2.Location();
896 //-- P2t is projected on the plane (P1,DirCylX,DirCylY)
897 gp_Dir DirCyl = Cyl1.Position().Direction();
898 Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
900 P2.SetCoord( P2t.X() - ProjP2OnDirCyl1*DirCyl.X()
901 ,P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y()
902 ,P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
904 Standard_Real R1pR2=R1+R2;
905 if(DistA1A2>(R1pR2+Tol)) {
906 typeres=IntAna_Empty;
909 else if(DistA1A2>(R1pR2)) {
910 //-- 1 Tangent line -------------------------------------OK
915 Standard_Real R1_R1pR2=R1/R1pR2;
916 pt1.SetCoord( P1.X() + R1_R1pR2 * (P2.X()-P1.X())
917 ,P1.Y() + R1_R1pR2 * (P2.Y()-P1.Y())
918 ,P1.Z() + R1_R1pR2 * (P2.Z()-P1.Z()));
921 else if(DistA1A2>RmR) {
922 //-- 2 lines ---------------------------------------------OK
927 gp_Dir DirA1A2=gp_Dir(P1P2);
928 gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
930 Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);
932 // Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
933 Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
934 Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
936 if((Beta+Beta)<Tol) {
938 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
939 ,P1.Y() + Alpha*DirA1A2.Y()
940 ,P1.Z() + Alpha*DirA1A2.Z());
943 pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X()
944 ,P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y()
945 ,P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z() );
947 pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X()
948 ,P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y()
949 ,P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
952 else if(DistA1A2>(RmR-Tol)) {
953 //-- 1 Tangent ------------------------------------------OK
957 Standard_Real R1_RmR=R1/RmR;
959 if(R1 < R2) R1_RmR = -R1_RmR;
961 pt1.SetCoord( P1.X() + R1_RmR * (P2.X()-P1.X())
962 ,P1.Y() + R1_RmR * (P2.Y()-P1.Y())
963 ,P1.Z() + R1_RmR * (P2.Z()-P1.Z()));
967 typeres=IntAna_Empty;
971 else { //-- No Parallel Axis ---------------------------------OK
972 if((RmR_Relative<=myEPSILON_CYLINDER_DELTA_RADIUS)
973 && (DistA1A2 <= myEPSILON_CYLINDER_DELTA_DISTANCE)) {
974 //-- PI/2 between the two axis and Intersection
975 //-- and identical radius
976 typeres=IntAna_Ellipse;
978 gp_Dir DirCyl1=Cyl1.Position().Direction();
979 gp_Dir DirCyl2=Cyl2.Position().Direction();
980 pt1=pt2=A1A2.PtIntersect();
982 Standard_Real A=DirCyl1.Angle(DirCyl2);
984 B=Abs(Sin(0.5*(M_PI-A)));
987 if(A==0.0 || B==0.0) {
993 gp_Vec dircyl1(DirCyl1);gp_Vec dircyl2(DirCyl2);
994 dir1 = gp_Dir(dircyl1.Added(dircyl2));
995 dir2 = gp_Dir(dircyl1.Subtracted(dircyl2));
997 param2 = Cyl1.Radius() / A;
998 param1 = Cyl1.Radius() / B;
999 param2bis= param1bis = Cyl1.Radius();
1000 if(param1 < param1bis) {
1001 A=param1; param1=param1bis; param1bis=A;
1003 if(param2 < param2bis) {
1004 A=param2; param2=param2bis; param2bis=A;
1008 if(Abs(DistA1A2-Cyl1.Radius()-Cyl2.Radius())<Tol) {
1009 typeres = IntAna_Point;
1010 Standard_Real d,p1,p2;
1012 gp_Dir D1 = Cyl1.Axis().Direction();
1013 gp_Dir D2 = Cyl2.Axis().Direction();
1014 A1A2.Distance(d,p1,p2);
1015 gp_Pnt P = Cyl1.Axis().Location();
1016 gp_Pnt P1(P.X() - p1*D1.X(),
1019 P = Cyl2.Axis().Location();
1020 gp_Pnt P2(P.X() - p2*D2.X(),
1026 pt1.SetCoord(P1.X() + p1*D1.X(),
1028 P1.Z() + p1*D1.Z());
1032 typeres=IntAna_NoGeometricSolution;
1037 //=======================================================================
1038 //function : IntAna_QuadQuadGeo
1039 //purpose : Cylinder - Cone
1040 //=======================================================================
1041 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1043 const Standard_Real Tol)
1044 : done(Standard_False),
1046 typeres(IntAna_Empty),
1057 myCommonGen(Standard_False),
1061 Perform(Cyl,Con,Tol);
1063 //=======================================================================
1064 //function : Perform
1066 //=======================================================================
1067 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1069 const Standard_Real )
1072 AxeOperator A1A2(Cyl.Axis(),Con.Axis());
1074 gp_Pnt Pt=Con.Apex();
1075 Standard_Real dist=Cyl.Radius()/(Tan(Con.SemiAngle()));
1076 gp_Dir dir=Cyl.Position().Direction();
1077 pt1.SetCoord( Pt.X() + dist*dir.X()
1078 ,Pt.Y() + dist*dir.Y()
1079 ,Pt.Z() + dist*dir.Z());
1080 pt2.SetCoord( Pt.X() - dist*dir.X()
1081 ,Pt.Y() - dist*dir.Y()
1082 ,Pt.Z() - dist*dir.Z());
1084 param1=param2=Cyl.Radius();
1086 typeres=IntAna_Circle;
1090 typeres=IntAna_NoGeometricSolution;
1093 //=======================================================================
1095 //purpose : Cylinder - Sphere
1096 //=======================================================================
1097 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1098 const gp_Sphere& Sph,
1099 const Standard_Real Tol)
1100 : done(Standard_False),
1102 typeres(IntAna_Empty),
1113 myCommonGen(Standard_False),
1117 Perform(Cyl,Sph,Tol);
1119 //=======================================================================
1120 //function : Perform
1122 //=======================================================================
1123 void IntAna_QuadQuadGeo::Perform( const gp_Cylinder& Cyl
1124 ,const gp_Sphere& Sph
1125 ,const Standard_Real)
1128 gp_Pnt Pt=Sph.Location();
1129 AxeOperator A1A2(Cyl.Axis(),Sph.Position().Axis());
1130 if((A1A2.Intersect() && Pt.Distance(A1A2.PtIntersect())==0.0 )
1132 if(Sph.Radius() < Cyl.Radius()) {
1133 typeres = IntAna_Empty;
1136 Standard_Real dist=Sqrt( Sph.Radius() * Sph.Radius() - Cyl.Radius() * Cyl.Radius() );
1137 gp_Dir dir=Cyl.Position().Direction();
1139 typeres=IntAna_Circle;
1140 pt1.SetCoord( Pt.X() + dist*dir.X()
1141 ,Pt.Y() + dist*dir.Y()
1142 ,Pt.Z() + dist*dir.Z());
1144 param1 = Cyl.Radius();
1145 if(dist>RealEpsilon()) {
1146 pt2.SetCoord( Pt.X() - dist*dir.X()
1147 ,Pt.Y() - dist*dir.Y()
1148 ,Pt.Z() - dist*dir.Z());
1149 param2=Cyl.Radius();
1155 typeres=IntAna_NoGeometricSolution;
1159 //=======================================================================
1160 //function : IntAna_QuadQuadGeo
1161 //purpose : Cone - Cone
1162 //=======================================================================
1163 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con1,
1164 const gp_Cone& Con2,
1165 const Standard_Real Tol)
1166 : done(Standard_False),
1168 typeres(IntAna_Empty),
1179 myCommonGen(Standard_False),
1183 Perform(Con1,Con2,Tol);
1186 //=======================================================================
1187 //function : Perform
1189 //=======================================================================
1190 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con1,
1191 const gp_Cone& Con2,
1192 const Standard_Real Tol)
1196 Standard_Real tg1, tg2, aDA1A2, aTol2;
1197 gp_Pnt aPApex1, aPApex2;
1199 Standard_Real TOL_APEX_CONF = 1.e-10;
1202 tg1=Tan(Con1.SemiAngle());
1203 tg2=Tan(Con2.SemiAngle());
1205 if((tg1 * tg2) < 0.) {
1210 aPApex1=Con1.Apex();
1211 aPApex2=Con2.Apex();
1212 aDA1A2=aPApex1.SquareDistance(aPApex2);
1214 AxeOperator A1A2(Con1.Axis(),Con2.Axis());
1220 gp_Pnt P=Con1.Apex();
1221 gp_Dir D=Con1.Position().Direction();
1222 Standard_Real d=gp_Vec(D).Dot(gp_Vec(P,Con2.Apex()));
1224 if(Abs(tg1-tg2)>myEPSILON_ANGLE_CONE) {
1225 if (fabs(d) < TOL_APEX_CONF) {
1226 typeres = IntAna_Point;
1231 x=(d*tg2)/(tg1+tg2);
1232 pt1.SetCoord( P.X() + x*D.X()
1237 x=(d*tg2)/(tg2-tg1);
1238 pt2.SetCoord( P.X() + x*D.X()
1244 typeres=IntAna_Circle;
1247 if (fabs(d) < TOL_APEX_CONF) {
1248 typeres=IntAna_Same;
1251 typeres=IntAna_Circle;
1254 pt1.SetCoord( P.X() + x*D.X()
1257 param1 = Abs(x * tg1);
1261 } //-- fin A1A2.Same
1263 else if((Abs(tg1-tg2)<myEPSILON_ANGLE_CONE) && (A1A2.Parallel())) {
1264 //-- voir AnVer12mai98
1265 Standard_Real DistA1A2=A1A2.Distance();
1266 gp_Dir DA1=Con1.Position().Direction();
1267 gp_Vec O1O2(Con1.Apex(),Con2.Apex());
1268 gp_Dir O1O2n(O1O2); // normalization of the vector before projection
1269 Standard_Real O1O2_DA1=gp_Vec(DA1).Dot(gp_Vec(O1O2n));
1271 gp_Vec O1_Proj_A2(O1O2n.X()-O1O2_DA1*DA1.X(),
1272 O1O2n.Y()-O1O2_DA1*DA1.Y(),
1273 O1O2n.Z()-O1O2_DA1*DA1.Z());
1274 gp_Dir DB1=gp_Dir(O1_Proj_A2);
1276 Standard_Real yO1O2=O1O2.Dot(gp_Vec(DA1));
1277 Standard_Real ABSTG1 = Abs(tg1);
1278 Standard_Real X2 = (DistA1A2/ABSTG1 - yO1O2)*0.5;
1279 Standard_Real X1 = X2+yO1O2;
1281 gp_Pnt P1(Con1.Apex().X() + X1*( DA1.X() + ABSTG1*DB1.X()),
1282 Con1.Apex().Y() + X1*( DA1.Y() + ABSTG1*DB1.Y()),
1283 Con1.Apex().Z() + X1*( DA1.Z() + ABSTG1*DB1.Z()));
1285 gp_Pnt MO1O2(0.5*(Con1.Apex().X()+Con2.Apex().X()),
1286 0.5*(Con1.Apex().Y()+Con2.Apex().Y()),
1287 0.5*(Con1.Apex().Z()+Con2.Apex().Z()));
1288 gp_Vec P1MO1O2(P1,MO1O2);
1290 gp_Dir DA1_X_DB1=DA1.Crossed(DB1);
1291 gp_Dir OrthoPln = DA1_X_DB1.Crossed(gp_Dir(P1MO1O2));
1293 IntAna_QuadQuadGeo INTER_QUAD_PLN(gp_Pln(P1,OrthoPln),Con1,Tol,Tol);
1294 if(INTER_QUAD_PLN.IsDone()) {
1295 switch(INTER_QUAD_PLN.TypeInter()) {
1296 case IntAna_Ellipse: {
1297 typeres=IntAna_Ellipse;
1298 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1300 dir1 = E.Position().Direction();
1301 dir2 = E.Position().XDirection();
1302 param1 = E.MajorRadius();
1303 param1bis = E.MinorRadius();
1307 case IntAna_Circle: {
1308 typeres=IntAna_Circle;
1309 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1311 dir1 = C.Position().XDirection();
1312 dir2 = C.Position().YDirection();
1313 param1 = C.Radius();
1317 case IntAna_Hyperbola: {
1318 typeres=IntAna_Hyperbola;
1319 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1320 pt1 = pt2 = H.Location();
1321 dir1 = H.Position().Direction();
1322 dir2 = H.Position().XDirection();
1323 param1 = param2 = H.MajorRadius();
1324 param1bis = param2bis = H.MinorRadius();
1329 typeres=IntAna_Line;
1330 gp_Lin H=INTER_QUAD_PLN.Line(1);
1331 pt1 = pt2 = H.Location();
1332 dir1 = dir2 = H.Position().Direction();
1333 param1 = param2 = 0.0;
1334 param1bis = param2bis = 0.0;
1339 typeres=IntAna_NoGeometricSolution;
1342 }// else if((Abs(tg1-tg2)<EPSILON_ANGLE_CONE) && (A1A2.Parallel()))
1344 else if (aDA1A2<aTol2) {
1346 // When apices are coinsided there can be 3 possible cases
1347 // 3.1 - empty solution (iRet=0)
1348 // 3.2 - one line when cone1 touches cone2 (iRet=1)
1349 // 3.3 - two lines when cone1 intersects cone2 (iRet=2)
1351 Standard_Integer iRet;
1352 Standard_Real aGamma, aBeta1, aBeta2;
1353 Standard_Real aD1, aR1, aTgBeta1, aTgBeta2, aHalfPI;
1354 Standard_Real aCosGamma, aSinGamma, aDx, aR2, aRD2, aD2;
1355 gp_Pnt2d aP0, aPA1, aP1, aPA2;
1359 // Preliminary analysis. Determination of iRet
1364 aPA1.SetCoord(aD1, 0.);
1365 aP0.SetCoord(0., 0.);
1369 aGamma=aAx1.Angle(aAx2);
1370 if (aGamma>aHalfPI){
1373 aCosGamma=Cos(aGamma);
1374 aSinGamma=Sin(aGamma);
1376 aBeta1=Con1.SemiAngle();
1377 aTgBeta1=Tan(aBeta1);
1378 aTgBeta1=Abs(aTgBeta1);
1380 aBeta2=Con2.SemiAngle();
1381 aTgBeta2=Tan(aBeta2);
1382 aTgBeta2=Abs(aTgBeta2);
1385 aP1.SetCoord(aD1, aR1);
1388 aVAx2.SetCoord(aCosGamma, aSinGamma);
1389 gp_Dir2d aDAx2(aVAx2);
1390 gp_Lin2d aLAx2(aP0, aDAx2);
1392 gp_Vec2d aV(aP0, aP1);
1394 aPA2=aP0.Translated(aDx*aDAx2);
1397 aDx=aPA2.Distance(aP0);
1401 aRD2=aPA2.Distance(aP1);
1403 if (aRD2>(aR2+Tol)) {
1405 typeres=IntAna_Empty; //nothing
1409 iRet=1; //touch case => 1 line
1410 if (aRD2<(aR2-Tol)) {
1411 iRet=2;//intersection => couple of lines
1414 // Finding the solution in 3D
1417 gp_Pnt aQApex1, aQA1, aQA2, aQX, aQX1, aQX2;
1418 gp_Dir aD3Ax1, aD3Ax2;
1420 IntAna_QuadQuadGeo aIntr;
1422 aQApex1=Con1.Apex();
1423 aD3Ax1=aAx1.Direction();
1424 aQA1.SetCoord(aQApex1.X()+aD1*aD3Ax1.X(),
1425 aQApex1.Y()+aD1*aD3Ax1.Y(),
1426 aQApex1.Z()+aD1*aD3Ax1.Z());
1428 aDx=aD3Ax1.Dot(aAx2.Direction());
1432 aD3Ax2=aAx2.Direction();
1434 aD2=aD1*sqrt((1.+aTgBeta1*aTgBeta1)/(1.+aTgBeta2*aTgBeta2));
1436 aQA2.SetCoord(aQApex1.X()+aD2*aD3Ax2.X(),
1437 aQApex1.Y()+aD2*aD3Ax2.Y(),
1438 aQApex1.Z()+aD2*aD3Ax2.Z());
1440 gp_Pln aPln1(aQA1, aD3Ax1);
1441 gp_Pln aPln2(aQA2, aD3Ax2);
1443 aIntr.Perform(aPln1, aPln2, Tol, Tol);
1444 if (!aIntr.IsDone()) {
1445 iRet=-1; // just in case. it must not be so
1446 typeres=IntAna_NoGeometricSolution;
1451 const gp_Dir& aDLin=aLin.Direction();
1452 gp_Vec aVLin(aDLin);
1453 gp_Pnt aOrig=aLin.Location();
1454 gp_Vec aVr(aQA1, aOrig);
1456 aQX=aOrig.Translated(aDx*aVLin);
1460 typeres=IntAna_Line;
1471 gp_Vec aVX(aQApex1, aQX);
1478 aDa=aQA1.Distance(aQX);
1479 aDx=sqrt(aR1*aR1-aDa*aDa);
1480 aQX1=aQX.Translated(aDx*aVLin);
1481 aQX2=aQX.Translated(-aDx*aVLin);
1485 gp_Vec aVX1(aQApex1, aQX1);
1487 gp_Vec aVX2(aQApex1, aQX2);
1490 } //else if (aDA1A2<aTol2) {
1491 //Case when cones have common generatrix
1492 else if(A1A2.Intersect()) {
1493 //Check if apex of one cone belongs another one
1494 Standard_Real u, v, tol2 = Tol*Tol;
1495 ElSLib::Parameters(Con2, aPApex1, u, v);
1496 gp_Pnt p = ElSLib::Value(u, v, Con2);
1497 if(aPApex1.SquareDistance(p) > tol2) {
1498 typeres=IntAna_NoGeometricSolution;
1502 ElSLib::Parameters(Con1, aPApex2, u, v);
1503 p = ElSLib::Value(u, v, Con1);
1504 if(aPApex2.SquareDistance(p) > tol2) {
1505 typeres=IntAna_NoGeometricSolution;
1509 //Cones have a common generatrix passing through apexes
1510 myCommonGen = Standard_True;
1512 //common generatrix of cones
1513 gp_Lin aGen(aPApex1, gp_Dir(gp_Vec(aPApex1, aPApex2)));
1515 //Intersection point of axes
1516 gp_Pnt aPAxeInt = A1A2.PtIntersect();
1518 //Characteristic point of intersection curve
1519 u = ElCLib::Parameter(aGen, aPAxeInt);
1520 myPChar = ElCLib::Value(u, aGen);
1523 //Other generatrixes of cones laying in maximal plane
1524 gp_Lin aGen1 = aGen.Rotated(Con1.Axis(), M_PI);
1525 gp_Lin aGen2 = aGen.Rotated(Con2.Axis(), M_PI);
1527 //Intersection point of generatrixes
1528 gp_Dir aN; //solution plane normal
1529 gp_Dir aD1 = aGen1.Direction();
1531 gp_Dir aD2(aD1.Crossed(aGen.Direction()));
1533 if(aD1.IsParallel(aGen2.Direction(), Precision::Angular())) {
1534 aN = aD1.Crossed(aD2);
1536 else if(aGen1.SquareDistance(aGen2) > tol2) {
1537 //Something wrong ???
1538 typeres=IntAna_NoGeometricSolution;
1542 gp_Dir D1 = aGen1.Position().Direction();
1543 gp_Dir D2 = aGen2.Position().Direction();
1544 gp_Pnt O1 = aGen1.Location();
1545 gp_Pnt O2 = aGen2.Location();
1546 Standard_Real D1DotD2 = D1.Dot(D2);
1547 Standard_Real aSin = 1.-D1DotD2*D1DotD2;
1548 gp_Vec O1O2 (O1,O2);
1549 Standard_Real U2 = (D1.XYZ()*(O1O2.Dot(D1))-(O1O2.XYZ())).Dot(D2.XYZ());
1551 gp_Pnt aPGint(ElCLib::Value(U2, aGen2));
1553 aD1 = gp_Dir(gp_Vec(aPGint, myPChar));
1554 aN = aD1.Crossed(aD2);
1556 //Plane that must contain intersection curves
1557 gp_Pln anIntPln(myPChar, aN);
1559 IntAna_QuadQuadGeo INTER_QUAD_PLN(anIntPln,Con1,Tol,Tol);
1561 if(INTER_QUAD_PLN.IsDone()) {
1562 switch(INTER_QUAD_PLN.TypeInter()) {
1563 case IntAna_Ellipse: {
1564 typeres=IntAna_Ellipse;
1565 gp_Elips E=INTER_QUAD_PLN.Ellipse(1);
1567 dir1 = E.Position().Direction();
1568 dir2 = E.Position().XDirection();
1569 param1 = E.MajorRadius();
1570 param1bis = E.MinorRadius();
1574 case IntAna_Circle: {
1575 typeres=IntAna_Circle;
1576 gp_Circ C=INTER_QUAD_PLN.Circle(1);
1578 dir1 = C.Position().XDirection();
1579 dir2 = C.Position().YDirection();
1580 param1 = C.Radius();
1584 case IntAna_Parabola: {
1585 typeres=IntAna_Parabola;
1586 gp_Parab Prb=INTER_QUAD_PLN.Parabola(1);
1587 pt1 = Prb.Location();
1588 dir1 = Prb.Position().Direction();
1589 dir2 = Prb.Position().XDirection();
1590 param1 = Prb.Focal();
1594 case IntAna_Hyperbola: {
1595 typeres=IntAna_Hyperbola;
1596 gp_Hypr H=INTER_QUAD_PLN.Hyperbola(1);
1597 pt1 = pt2 = H.Location();
1598 dir1 = H.Position().Direction();
1599 dir2 = H.Position().XDirection();
1600 param1 = param2 = H.MajorRadius();
1601 param1bis = param2bis = H.MinorRadius();
1606 typeres=IntAna_NoGeometricSolution;
1612 typeres=IntAna_NoGeometricSolution;
1615 //=======================================================================
1616 //function : IntAna_QuadQuadGeo
1617 //purpose : Sphere - Cone
1618 //=======================================================================
1619 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
1621 const Standard_Real Tol)
1622 : done(Standard_False),
1624 typeres(IntAna_Empty),
1635 myCommonGen(Standard_False),
1639 Perform(Sph,Con,Tol);
1641 //=======================================================================
1642 //function : Perform
1644 //=======================================================================
1645 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
1647 const Standard_Real)
1653 AxeOperator A1A2(Con.Axis(),Sph.Position().Axis());
1654 gp_Pnt Pt=Sph.Location();
1656 if((A1A2.Intersect() && (Pt.Distance(A1A2.PtIntersect())==0.0))
1658 gp_Pnt ConApex= Con.Apex();
1659 Standard_Real dApexSphCenter=Pt.Distance(ConApex);
1661 if(dApexSphCenter>RealEpsilon()) {
1662 ConDir = gp_Dir(gp_Vec(ConApex,Pt));
1665 ConDir = Con.Position().Direction();
1668 Standard_Real Rad=Sph.Radius();
1669 Standard_Real tga=Tan(Con.SemiAngle());
1673 //-- x: Roots of (x**2 + y**2 = Rad**2)
1674 //-- tga = y / (x+dApexSphCenter)
1675 Standard_Real tgatga = tga * tga;
1676 math_DirectPolynomialRoots Eq( 1.0+tgatga
1677 ,2.0*tgatga*dApexSphCenter
1678 ,-Rad*Rad + dApexSphCenter*dApexSphCenter*tgatga);
1680 Standard_Integer nbsol=Eq.NbSolutions();
1682 typeres=IntAna_Empty;
1685 typeres=IntAna_Circle;
1687 Standard_Real x = Eq.Value(1);
1688 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1690 pt1.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1691 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1692 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1693 param1 = tga * dApexSphCenterpx;
1694 param1 = Abs(param1);
1696 if(param1<=myEPSILON_MINI_CIRCLE_RADIUS) {
1697 typeres=IntAna_PointAndCircle;
1702 Standard_Real x=Eq.Value(2);
1703 Standard_Real dApexSphCenterpx = dApexSphCenter+x;
1705 pt2.SetCoord( ConApex.X() + (dApexSphCenterpx) * ConDir.X()
1706 ,ConApex.Y() + (dApexSphCenterpx) * ConDir.Y()
1707 ,ConApex.Z() + (dApexSphCenterpx) * ConDir.Z());
1708 param2 = tga * dApexSphCenterpx;
1709 param2 = Abs(param2);
1711 if(param2<=myEPSILON_MINI_CIRCLE_RADIUS) {
1712 typeres=IntAna_PointAndCircle;
1719 done=Standard_False;
1723 typeres=IntAna_NoGeometricSolution;
1727 //=======================================================================
1728 //function : IntAna_QuadQuadGeo
1729 //purpose : Sphere - Sphere
1730 //=======================================================================
1731 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo( const gp_Sphere& Sph1
1732 ,const gp_Sphere& Sph2
1733 ,const Standard_Real Tol)
1734 : done(Standard_False),
1736 typeres(IntAna_Empty),
1747 myCommonGen(Standard_False),
1751 Perform(Sph1,Sph2,Tol);
1753 //=======================================================================
1754 //function : Perform
1756 //=======================================================================
1757 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph1,
1758 const gp_Sphere& Sph2,
1759 const Standard_Real Tol)
1762 gp_Pnt O1=Sph1.Location();
1763 gp_Pnt O2=Sph2.Location();
1764 Standard_Real dO1O2=O1.Distance(O2);
1765 Standard_Real R1=Sph1.Radius();
1766 Standard_Real R2=Sph2.Radius();
1767 Standard_Real Rmin,Rmax;
1768 typeres=IntAna_Empty;
1769 param2bis=0.0; //-- pour eviter param2bis not used ....
1771 if(R1>R2) { Rmin=R2; Rmax=R1; } else { Rmin=R1; Rmax=R2; }
1773 if(dO1O2<=Tol && (Abs(R1-R2) <= Tol)) {
1774 typeres = IntAna_Same;
1777 if(dO1O2<=Tol) { return; }
1778 gp_Dir Dir=gp_Dir(gp_Vec(O1,O2));
1779 Standard_Real t = Rmax - dO1O2 - Rmin;
1781 //----------------------------------------------------------------------
1782 //-- |----------------- R1 --------------------|
1783 //-- |----dO1O2-----|-----------R2----------|
1786 //-- |------ R1 ------|---------dO1O2----------|
1787 //-- |-------------------R2-----------------------|
1789 //----------------------------------------------------------------------
1790 if(t >= 0.0 && t <=Tol) {
1791 typeres = IntAna_Point;
1794 if(R1==Rmax) t2=(R1 + (R2 + dO1O2)) * 0.5;
1795 else t2=(-R1+(dO1O2-R2))*0.5;
1797 pt1.SetCoord( O1.X() + t2*Dir.X()
1798 ,O1.Y() + t2*Dir.Y()
1799 ,O1.Z() + t2*Dir.Z());
1802 //-----------------------------------------------------------------
1803 //-- |----------------- dO1O2 --------------------|
1804 //-- |----R1-----|-----------R2----------|-Tol-|
1806 //-- |----------------- Rmax --------------------|
1807 //-- |----Rmin----|-------dO1O2-------|-Tol-|
1809 //-----------------------------------------------------------------
1810 if((dO1O2 > (R1+R2+Tol)) || (Rmax > (dO1O2+Rmin+Tol))) {
1811 typeres=IntAna_Empty;
1814 //---------------------------------------------------------------
1817 //---------------------------------------------------------------
1818 Standard_Real Alpha=0.5*(R1*R1-R2*R2+dO1O2*dO1O2)/(dO1O2);
1819 Standard_Real Beta = R1*R1-Alpha*Alpha;
1820 Beta = (Beta>0.0)? Sqrt(Beta) : 0.0;
1822 if(Beta<= myEPSILON_MINI_CIRCLE_RADIUS) {
1823 typeres = IntAna_Point;
1824 Alpha = (R1 + (dO1O2 - R2)) * 0.5;
1827 typeres = IntAna_Circle;
1831 pt1.SetCoord( O1.X() + Alpha*Dir.X()
1832 ,O1.Y() + Alpha*Dir.Y()
1833 ,O1.Z() + Alpha*Dir.Z());
1841 //=======================================================================
1842 //function : IntAna_QuadQuadGeo
1843 //purpose : Plane - Torus
1844 //=======================================================================
1845 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Pln& Pln,
1846 const gp_Torus& Tor,
1847 const Standard_Real Tol)
1848 : done(Standard_False),
1850 typeres(IntAna_Empty),
1861 myCommonGen(Standard_False),
1865 Perform(Pln,Tor,Tol);
1867 //=======================================================================
1868 //function : Perform
1870 //=======================================================================
1871 void IntAna_QuadQuadGeo::Perform(const gp_Pln& Pln,
1872 const gp_Torus& Tor,
1873 const Standard_Real Tol)
1875 done = Standard_True;
1877 Standard_Real aRMin, aRMaj;
1879 aRMin = Tor.MinorRadius();
1880 aRMaj = Tor.MajorRadius();
1881 if (aRMin >= aRMaj) {
1882 typeres = IntAna_NoGeometricSolution;
1886 const gp_Ax1 aPlnAx = Pln.Axis();
1887 const gp_Ax1 aTorAx = Tor.Axis();
1889 Standard_Boolean bParallel, bNormal;
1891 bParallel = aTorAx.IsParallel(aPlnAx, myEPSILON_AXES_PARA);
1892 bNormal = !bParallel ? aTorAx.IsNormal(aPlnAx, myEPSILON_AXES_PARA) : Standard_False;
1893 if (!bNormal && !bParallel) {
1894 typeres = IntAna_NoGeometricSolution;
1898 Standard_Real aDist;
1900 gp_Pnt aTorLoc = aTorAx.Location();
1902 Standard_Real aDt, X, Y, Z, A, B, C, D;
1904 Pln.Coefficients(A,B,C,D);
1905 aTorLoc.Coord(X,Y,Z);
1906 aDist = A*X + B*Y + C*Z + D;
1908 if ((Abs(aDist) - aRMin) > Tol) {
1909 typeres=IntAna_Empty;
1913 typeres = IntAna_Circle;
1915 pt1.SetCoord(X - aDist*A, Y - aDist*B, Z - aDist*C);
1916 aDt = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
1917 param1 = aRMaj + aDt;
1918 dir1 = aTorAx.Direction();
1920 if ((Abs(aDist) < aRMin) && (aDt > Tol)) {
1922 param2 = aRMaj - aDt;
1929 aDist = Pln.Distance(aTorLoc);
1930 if (aDist > myEPSILON_DISTANCE) {
1931 typeres = IntAna_NoGeometricSolution;
1935 typeres = IntAna_Circle;
1936 param2 = param1 = aRMin;
1937 dir2 = dir1 = aPlnAx.Direction();
1940 gp_Dir aDir = aTorAx.Direction()^dir1;
1941 pt1.SetXYZ(aTorLoc.XYZ() + aRMaj*aDir.XYZ());
1942 pt2.SetXYZ(aTorLoc.XYZ() - aRMaj*aDir.XYZ());
1946 //=======================================================================
1947 //function : IntAna_QuadQuadGeo
1948 //purpose : Cylinder - Torus
1949 //=======================================================================
1950 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cylinder& Cyl,
1951 const gp_Torus& Tor,
1952 const Standard_Real Tol)
1953 : done(Standard_False),
1955 typeres(IntAna_Empty),
1966 myCommonGen(Standard_False),
1970 Perform(Cyl,Tor,Tol);
1972 //=======================================================================
1973 //function : Perform
1975 //=======================================================================
1976 void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl,
1977 const gp_Torus& Tor,
1978 const Standard_Real Tol)
1980 done = Standard_True;
1982 Standard_Real aRMin, aRMaj;
1984 aRMin = Tor.MinorRadius();
1985 aRMaj = Tor.MajorRadius();
1986 if (aRMin >= aRMaj) {
1987 typeres = IntAna_NoGeometricSolution;
1991 const gp_Ax1 aCylAx = Cyl.Axis();
1992 const gp_Ax1 aTorAx = Tor.Axis();
1994 const gp_Lin aLin(aTorAx);
1995 const gp_Pnt aLocCyl = Cyl.Location();
1997 if (!aTorAx.IsParallel(aCylAx, myEPSILON_AXES_PARA) ||
1998 (aLin.Distance(aLocCyl) > myEPSILON_DISTANCE)) {
1999 typeres = IntAna_NoGeometricSolution;
2003 Standard_Real aRCyl;
2005 aRCyl = Cyl.Radius();
2006 if (((aRCyl + Tol) < (aRMaj - aRMin)) || ((aRCyl - Tol) > (aRMaj + aRMin))) {
2007 typeres = IntAna_Empty;
2011 typeres = IntAna_Circle;
2013 Standard_Real aDist = Sqrt(Abs(aRMin*aRMin - (aRCyl-aRMaj)*(aRCyl-aRMaj)));
2014 gp_XYZ aTorLoc = aTorAx.Location().XYZ();
2016 dir1 = aTorAx.Direction();
2017 pt1.SetXYZ(aTorLoc + aDist*dir1.XYZ());
2020 if ((aDist > Tol) && (aRCyl > (aRMaj - aRMin)) &&
2021 (aRCyl < (aRMaj + aRMin))) {
2023 pt2.SetXYZ(aTorLoc - aDist*dir2.XYZ());
2029 //=======================================================================
2030 //function : IntAna_QuadQuadGeo
2031 //purpose : Cone - Torus
2032 //=======================================================================
2033 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Cone& Con,
2034 const gp_Torus& Tor,
2035 const Standard_Real Tol)
2036 : done(Standard_False),
2038 typeres(IntAna_Empty),
2049 myCommonGen(Standard_False),
2053 Perform(Con,Tor,Tol);
2055 //=======================================================================
2056 //function : Perform
2058 //=======================================================================
2059 void IntAna_QuadQuadGeo::Perform(const gp_Cone& Con,
2060 const gp_Torus& Tor,
2061 const Standard_Real Tol)
2063 done = Standard_True;
2065 Standard_Real aRMin, aRMaj;
2067 aRMin = Tor.MinorRadius();
2068 aRMaj = Tor.MajorRadius();
2069 if (aRMin >= aRMaj) {
2070 typeres = IntAna_NoGeometricSolution;
2074 const gp_Ax1 aConAx = Con.Axis();
2075 const gp_Ax1 aTorAx = Tor.Axis();
2077 const gp_Lin aLin(aTorAx);
2078 const gp_Pnt aConApex = Con.Apex();
2080 if (!aTorAx.IsParallel(aConAx, myEPSILON_AXES_PARA) ||
2081 (aLin.Distance(aConApex) > myEPSILON_DISTANCE)) {
2082 typeres = IntAna_NoGeometricSolution;
2086 Standard_Real anAngle, aDist, aParam[4];
2088 gp_Pnt aTorLoc, aPCT, aPN, aPt[4];
2091 anAngle = Con.SemiAngle();
2092 aTorLoc = aTorAx.Location();
2094 aPN.SetXYZ(aTorLoc.XYZ() + aRMaj*Tor.YAxis().Direction().XYZ());
2095 gp_Dir aDN (gp_Vec(aTorLoc, aPN));
2096 gp_Ax1 anAxCLRot(aConApex, aDN);
2097 gp_Lin aConL = aLin.Rotated(anAxCLRot, anAngle);
2098 gp_Dir aDL = aConL.Position().Direction();
2099 gp_Dir aXDir = Tor.XAxis().Direction();
2101 typeres = IntAna_Empty;
2103 for (i = 0; i < 2; ++i) {
2107 aPCT.SetXYZ(aTorLoc.XYZ() + aRMaj*aXDir.XYZ());
2109 aDist = aConL.Distance(aPCT);
2110 if (aDist > aRMin+Tol) {
2114 typeres = IntAna_Circle;
2116 gp_XYZ aPh = aPCT.XYZ() - aDist*aConL.Normal(aPCT).Direction().XYZ();
2117 aDist = Sqrt(Abs(aRMin*aRMin - aDist*aDist));
2120 gp_XYZ aDVal = aDist*aDL.XYZ();
2121 aP.SetXYZ(aPh + aDVal);
2122 aParam[nbint] = aLin.Distance(aP);
2123 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2124 aDir[nbint] = aTorAx.Direction();
2126 if ((aDist < aRMin) && (aDVal.Modulus() > Tol)) {
2127 aP.SetXYZ(aPh - aDVal);
2128 aParam[nbint] = aLin.Distance(aP);
2129 aPt[nbint].SetXYZ(aP.XYZ() - aParam[nbint]*aXDir.XYZ());
2130 aDir[nbint] = aDir[nbint-1];
2135 for (i = 0; i < nbint; ++i) {
2167 //=======================================================================
2168 //function : IntAna_QuadQuadGeo
2169 //purpose : Sphere - Torus
2170 //=======================================================================
2171 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Sphere& Sph,
2172 const gp_Torus& Tor,
2173 const Standard_Real Tol)
2174 : done(Standard_False),
2176 typeres(IntAna_Empty),
2187 myCommonGen(Standard_False),
2191 Perform(Sph,Tor,Tol);
2193 //=======================================================================
2194 //function : Perform
2196 //=======================================================================
2197 void IntAna_QuadQuadGeo::Perform(const gp_Sphere& Sph,
2198 const gp_Torus& Tor,
2199 const Standard_Real Tol)
2201 done = Standard_True;
2203 Standard_Real aRMin, aRMaj;
2205 aRMin = Tor.MinorRadius();
2206 aRMaj = Tor.MajorRadius();
2207 if (aRMin >= aRMaj) {
2208 typeres = IntAna_NoGeometricSolution;
2212 const gp_Ax1 aTorAx = Tor.Axis();
2213 const gp_Lin aLin(aTorAx);
2214 const gp_Pnt aSphLoc = Sph.Location();
2216 if (aLin.Distance(aSphLoc) > myEPSILON_DISTANCE) {
2217 typeres = IntAna_NoGeometricSolution;
2221 Standard_Real aRSph, aDist;
2224 gp_Dir aXDir = Tor.XAxis().Direction();
2225 aTorLoc.SetXYZ(aTorAx.Location().XYZ() + aRMaj*aXDir.XYZ());
2226 aRSph = Sph.Radius();
2228 gp_Vec aVec12(aTorLoc, aSphLoc);
2229 aDist = aVec12.Magnitude();
2230 if (((aDist - Tol) > (aRMin + aRSph)) ||
2231 ((aDist + Tol) < Abs(aRMin - aRSph))) {
2232 typeres = IntAna_Empty;
2236 typeres = IntAna_Circle;
2238 Standard_Real anAlpha, aBeta;
2240 anAlpha = 0.5*(aRMin*aRMin - aRSph*aRSph + aDist*aDist ) / aDist;
2241 aBeta = Sqrt(Abs(aRMin*aRMin - anAlpha*anAlpha));
2243 gp_Dir aDir12(aVec12);
2244 gp_XYZ aPh = aTorLoc.XYZ() + anAlpha*aDir12.XYZ();
2245 gp_Dir aDC = Tor.YAxis().Direction()^aDir12;
2248 gp_XYZ aDVal = aBeta*aDC.XYZ();
2249 aP.SetXYZ(aPh + aDVal);
2250 param1 = aLin.Distance(aP);
2251 pt1.SetXYZ(aP.XYZ() - param1*aXDir.XYZ());
2252 dir1 = aTorAx.Direction();
2254 if ((aDist < (aRSph + aRMin)) && (aDist > Abs(aRSph - aRMin)) &&
2255 (aDVal.Modulus() > Tol)) {
2256 aP.SetXYZ(aPh - aDVal);
2257 param2 = aLin.Distance(aP);
2258 pt2.SetXYZ(aP.XYZ() - param2*aXDir.XYZ());
2264 //=======================================================================
2265 //function : IntAna_QuadQuadGeo
2266 //purpose : Torus - Torus
2267 //=======================================================================
2268 IntAna_QuadQuadGeo::IntAna_QuadQuadGeo(const gp_Torus& Tor1,
2269 const gp_Torus& Tor2,
2270 const Standard_Real Tol)
2271 : done(Standard_False),
2273 typeres(IntAna_Empty),
2284 myCommonGen(Standard_False),
2288 Perform(Tor1,Tor2,Tol);
2290 //=======================================================================
2291 //function : Perform
2293 //=======================================================================
2294 void IntAna_QuadQuadGeo::Perform(const gp_Torus& Tor1,
2295 const gp_Torus& Tor2,
2296 const Standard_Real Tol)
2298 done = Standard_True;
2300 Standard_Real aRMin1, aRMin2, aRMaj1, aRMaj2;
2302 aRMin1 = Tor1.MinorRadius();
2303 aRMaj1 = Tor1.MajorRadius();
2304 aRMin2 = Tor2.MinorRadius();
2305 aRMaj2 = Tor2.MajorRadius();
2306 if (aRMin1 >= aRMaj1 || aRMin2 >= aRMaj2) {
2307 typeres = IntAna_NoGeometricSolution;
2311 const gp_Ax1 anAx1 = Tor1.Axis();
2312 const gp_Ax1 anAx2 = Tor2.Axis();
2315 if (!anAx1.IsParallel(anAx2, myEPSILON_AXES_PARA) ||
2316 (aL1.Distance(anAx2.Location()) > myEPSILON_DISTANCE)) {
2317 typeres = IntAna_NoGeometricSolution;
2321 gp_Pnt aLoc1, aLoc2;
2323 aLoc1 = anAx1.Location();
2324 aLoc2 = anAx2.Location();
2326 if (aLoc1.IsEqual(aLoc2, Tol) &&
2327 (Abs(aRMin1 - aRMin2) <= Tol) &&
2328 (Abs(aRMaj1 - aRMaj2) <= Tol)) {
2329 typeres = IntAna_Same;
2333 Standard_Real aDist;
2336 gp_Dir aXDir1 = Tor1.XAxis().Direction();
2337 aP1.SetXYZ(aLoc1.XYZ() + aRMaj1*aXDir1.XYZ());
2338 aP2.SetXYZ(aLoc2.XYZ() + aRMaj2*aXDir1.XYZ());
2340 gp_Vec aV12(aP1, aP2);
2341 aDist = aV12.Magnitude();
2342 if (((aDist - Tol) > (aRMin1 + aRMin2)) ||
2343 ((aDist + Tol) < Abs(aRMin1 - aRMin2))) {
2344 typeres = IntAna_Empty;
2348 typeres = IntAna_Circle;
2350 Standard_Real anAlpha, aBeta;
2352 anAlpha = 0.5*(aRMin1*aRMin1 - aRMin2*aRMin2 + aDist*aDist ) / aDist;
2353 aBeta = Sqrt(Abs(aRMin1*aRMin1 - anAlpha*anAlpha));
2355 gp_Dir aDir12(aV12);
2356 gp_XYZ aPh = aP1.XYZ() + anAlpha*aDir12.XYZ();
2357 gp_Dir aDC = Tor1.YAxis().Direction()^aDir12;
2360 gp_XYZ aDVal = aBeta*aDC.XYZ();
2361 aP.SetXYZ(aPh + aDVal);
2362 param1 = aL1.Distance(aP);
2363 pt1.SetXYZ(aP.XYZ() - param1*aXDir1.XYZ());
2364 dir1 = anAx1.Direction();
2366 if ((aDist < (aRMin1 + aRMin2)) && (aDist > Abs(aRMin1 - aRMin2)) &&
2367 aDVal.Modulus() > Tol) {
2368 aP.SetXYZ(aPh - aDVal);
2369 param2 = aL1.Distance(aP);
2370 pt2.SetXYZ(aP.XYZ() - param2*aXDir1.XYZ());
2376 //=======================================================================
2378 //purpose : Returns a Point
2379 //=======================================================================
2380 gp_Pnt IntAna_QuadQuadGeo::Point(const Standard_Integer n) const
2382 if(!done) { StdFail_NotDone::Raise(); }
2383 if(n>nbint || n<1) { Standard_DomainError::Raise(); }
2384 if(typeres==IntAna_PointAndCircle) {
2385 if(n!=1) { Standard_DomainError::Raise(); }
2386 if(param1==0.0) return(pt1);
2389 else if(typeres==IntAna_Point) {
2390 if(n==1) return(pt1);
2394 // WNT (what can you expect from MicroSoft ?)
2395 return gp_Pnt(0,0,0);
2397 //=======================================================================
2399 //purpose : Returns a Line
2400 //=======================================================================
2401 gp_Lin IntAna_QuadQuadGeo::Line(const Standard_Integer n) const
2403 if(!done) { StdFail_NotDone::Raise(); }
2404 if((n>nbint) || (n<1) || (typeres!=IntAna_Line)) {
2405 Standard_DomainError::Raise();
2407 if(n==1) { return(gp_Lin(pt1,dir1)); }
2408 else { return(gp_Lin(pt2,dir2)); }
2410 //=======================================================================
2412 //purpose : Returns a Circle
2413 //=======================================================================
2414 gp_Circ IntAna_QuadQuadGeo::Circle(const Standard_Integer n) const
2416 if(!done) { StdFail_NotDone::Raise(); }
2417 if(typeres==IntAna_PointAndCircle) {
2418 if(n!=1) { Standard_DomainError::Raise(); }
2419 if(param2==0.0) return(gp_Circ(DirToAx2(pt1,dir1),param1));
2420 return(gp_Circ(DirToAx2(pt2,dir2),param2));
2422 else if((n>nbint) || (n<1) || (typeres!=IntAna_Circle)) {
2423 Standard_DomainError::Raise();
2425 if (n==1) { return(gp_Circ(DirToAx2(pt1,dir1),param1));}
2426 else if (n==2) { return(gp_Circ(DirToAx2(pt2,dir2),param2));}
2427 else if (n==3) { return(gp_Circ(DirToAx2(pt3,dir3),param3));}
2428 else { return(gp_Circ(DirToAx2(pt4,dir4),param4));}
2431 //=======================================================================
2432 //function : Ellipse
2433 //purpose : Returns a Elips
2434 //=======================================================================
2435 gp_Elips IntAna_QuadQuadGeo::Ellipse(const Standard_Integer n) const
2437 if(!done) { StdFail_NotDone::Raise(); }
2438 if((n>nbint) || (n<1) || (typeres!=IntAna_Ellipse)) {
2439 Standard_DomainError::Raise();
2443 Standard_Real R1=param1, R2=param1bis, aTmp;
2445 aTmp=R1; R1=R2; R2=aTmp;
2447 gp_Ax2 anAx2(pt1, dir1 ,dir2);
2448 gp_Elips anElips (anAx2, R1, R2);
2452 Standard_Real R1=param2, R2=param2bis, aTmp;
2454 aTmp=R1; R1=R2; R2=aTmp;
2456 gp_Ax2 anAx2(pt2, dir2 ,dir1);
2457 gp_Elips anElips (anAx2, R1, R2);
2461 //=======================================================================
2462 //function : Parabola
2463 //purpose : Returns a Parabola
2464 //=======================================================================
2465 gp_Parab IntAna_QuadQuadGeo::Parabola(const Standard_Integer n) const
2468 StdFail_NotDone::Raise();
2470 if (typeres!=IntAna_Parabola) {
2471 Standard_DomainError::Raise();
2473 if((n>nbint) || (n!=1)) {
2474 Standard_OutOfRange::Raise();
2476 return(gp_Parab(gp_Ax2( pt1
2481 //=======================================================================
2482 //function : Hyperbola
2483 //purpose : Returns a Hyperbola
2484 //=======================================================================
2485 gp_Hypr IntAna_QuadQuadGeo::Hyperbola(const Standard_Integer n) const
2488 StdFail_NotDone::Raise();
2490 if((n>nbint) || (n<1) || (typeres!=IntAna_Hyperbola)) {
2491 Standard_DomainError::Raise();
2494 return(gp_Hypr(gp_Ax2( pt1
2497 ,param1,param1bis));
2500 return(gp_Hypr(gp_Ax2( pt2
2503 ,param2,param2bis));
2506 //=======================================================================
2507 //function : HasCommonGen
2509 //=======================================================================
2510 Standard_Boolean IntAna_QuadQuadGeo::HasCommonGen() const
2514 //=======================================================================
2517 //=======================================================================
2518 const gp_Pnt& IntAna_QuadQuadGeo::PChar() const
2522 //=======================================================================
2523 //function : RefineDir
2525 //=======================================================================
2526 void RefineDir(gp_Dir& aDir)
2528 Standard_Integer k, m, n;
2529 Standard_Real aC[3];
2531 aDir.Coord(aC[0], aC[1], aC[2]);
2535 for (k=0; k<3; ++k) {
2536 if (aC[k]==1. || aC[k]==-1.) {
2539 else if (aC[k]!=0.) {
2545 Standard_Real aEps, aR1, aR2, aNum;
2551 for (k=0; k<3; ++k) {
2553 aNum=(m)? aC[k] : -aC[k];
2554 if (aNum>aR1 && aNum<aR2) {
2567 aDir.SetCoord(aC[0], aC[1], aC[2]);