src/IntAna/IntAna_IntQuadQuad.cxx

   1 // Created on: 1992-06-29
   2 // Created by: Laurent BUCHARD
   3 // Copyright (c) 1992-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   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.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <stdio.h>
  18
  19 #include <Standard_Stream.hxx>
  20
  21 #ifndef OCCT_DEBUG
  22 #define No_Standard_RangeError
  23 #define No_Standard_OutOfRange
  24 #endif
  25
  26 //======================================================================
  27 //== I n t e r s e c t i o n   C O N E           Q U A D R I Q U E
  28 //==                           C Y L I N D R E   Q U A D R I Q U E
  29 //======================================================================
  30
  31 #include <ElSLib.hxx>
  32 #include <gp_Ax2.hxx>
  33 #include <gp_Ax3.hxx>
  34 #include <gp_Cone.hxx>
  35 #include <gp_Cylinder.hxx>
  36 #include <gp_Pnt.hxx>
  37 #include <IntAna_Curve.hxx>
  38 #include <IntAna_IntQuadQuad.hxx>
  39 #include <IntAna_Quadric.hxx>
  40 #include <math_TrigonometricFunctionRoots.hxx>
  41 #include <Standard_DomainError.hxx>
  42 #include <Standard_OutOfRange.hxx>
  43 #include <StdFail_NotDone.hxx>
  44
  45 //=======================================================================
  46 //function : AddSpecialPoints
  47 //purpose  : Sometimes the boundaries theTheta1 and theTheta2 are
  48 //            computed with some inaccuracy. At that, some special points
  49 //            (cone apex or sphere pole(s)), which are true intersection
  50 //            points lie out of the domain [theTheta1, theTheta2] of the ALine.
  51 //           This function corrects these boundaries to make them be included
  52 //            in the domain of the ALine.
  53 //           Parameters Theta1 and Theta2 must be initialized
  54 //            before calling this function.
  55 //=======================================================================
  56 template <class gpSmth>
  57 static void AddSpecialPoints(const IntAna_Quadric& theQuad,
  58                              const gpSmth& theGpObj,
  59                              Standard_Real& theTheta1,
  60                              Standard_Real& theTheta2)
  61 {
  62   const Standard_Real aPeriod = M_PI + M_PI;
  63   const NCollection_List<gp_Pnt> &aLSP = theQuad.SpecialPoints();
  64
  65   if (aLSP.IsEmpty())
  66     return;
  67
  68   Standard_Real aU = 0.0, aV = 0.0;
  69   Standard_Real aMaxDelta = 0.0;
  70   for (NCollection_List<gp_Pnt>::Iterator anItr(aLSP); anItr.More(); anItr.Next())
  71   {
  72     const gp_Pnt &aPt = anItr.Value();
  73     ElSLib::Parameters(theGpObj, aPt, aU, aV);
  74     const gp_Pnt aPProj(ElSLib::Value(aU, aV, theGpObj));
  75
  76     if (aPt.SquareDistance(aPProj) > Precision::SquareConfusion())
  77     {
  78       // aPt is not an intersection point
  79       continue;
  80     }
  81
  82     Standard_Real aDelta1 = Min(aU - theTheta1, 0.0),
  83                   aDelta2 = Max(aU - theTheta2, 0.0);
  84
  85     if (aDelta1 < -M_PI)
  86     {
  87       // Must be aDelta1 = Min(aU - theTheta1 + aPeriod, 0.0).
  88       // But aU - theTheta1 + aPeriod >= 0 always.
  89       aDelta1 = 0.0;
  90     }
  91
  92     if (aDelta2 > M_PI)
  93     {
  94       // Must be aDelta2 = Max(aU - theTheta2 - aPeriod, 0.0).
  95       // But aU - theTheta2 - aPeriod <= 0 always.
  96       aDelta2 = 0.0;
  97     }
  98
  99     const Standard_Real aDelta = Max(-aDelta1, aDelta2);
 100     aMaxDelta = Max(aMaxDelta, aDelta);
 101   }
 102
 103   if(aMaxDelta != 0.0)
 104   {
 105     theTheta1 -= aMaxDelta;
 106     theTheta2 += aMaxDelta;
 107     if ((theTheta2 - theTheta1) > aPeriod)
 108     {
 109       theTheta2 = theTheta1 + aPeriod;
 110     }
 111   }
 112 }
 113
 114 //=======================================================================
 115 //class : TrigonometricRoots
 116 //purpose: Classe Interne (Donne des racines classees d un polynome trigo)
 117 //======================================================================
 118 class TrigonometricRoots {
 119
 120  private:
 121   Standard_Real Roots[4];
 122   Standard_Boolean done;
 123   Standard_Integer NbRoots;
 124   Standard_Boolean infinite_roots;
 125
 126  public:
 127   TrigonometricRoots(const Standard_Real CC,
 128                      const Standard_Real SC,
 129                      const Standard_Real C,
 130                      const Standard_Real S,
 131                      const Standard_Real Cte,
 132                      const Standard_Real Binf,
 133                      const Standard_Real Bsup);
 134   //IsDone
 135   Standard_Boolean IsDone() {
 136     return done;
 137   }
 138   //IsARoot
 139   Standard_Boolean IsARoot(Standard_Real u) {
 140     Standard_Integer i;
 141     Standard_Real aEps=RealEpsilon();
 142     Standard_Real PIpPI = M_PI + M_PI;
 143     //
 144     for(i=0 ; i<NbRoots; ++i) {
 145       if(Abs(u - Roots[i])<=aEps) {
 146         return Standard_True;
 147       }
 148       if(Abs(u - Roots[i]-PIpPI)<=aEps) {
 149         return Standard_True;
 150       }
 151     }
 152     return Standard_False;
 153   }
 154   //NbSolutions
 155   Standard_Integer NbSolutions() {
 156     if(!done) {
 157       throw StdFail_NotDone();
 158     }
 159     return NbRoots;
 160   }
 161   //InfiniteRoots
 162   Standard_Boolean InfiniteRoots() {
 163     if(!done) {
 164       throw StdFail_NotDone();
 165     }
 166     return infinite_roots;
 167   }
 168   //Value
 169   Standard_Real Value(const Standard_Integer n) {
 170     if((!done)||(n>NbRoots)) {
 171       throw StdFail_NotDone();
 172     }
 173     return Roots[n-1];
 174   }
 175 };
 176 //=======================================================================
 177 //function : TrigonometricRoots::TrigonometricRoots
 178 //purpose  :
 179 //=======================================================================
 180 TrigonometricRoots::TrigonometricRoots(const Standard_Real CC,
 181                                        const Standard_Real SC,
 182                                        const Standard_Real C,
 183                                        const Standard_Real S,
 184                                        const Standard_Real Cte,
 185                                        const Standard_Real Binf,
 186                                        const Standard_Real Bsup)
 187 : infinite_roots(Standard_False)
 188 {
 189   Standard_Integer i, j, SvNbRoots;
 190   Standard_Boolean Triee;
 191   Standard_Real PIpPI = M_PI + M_PI;
 192   //
 193   done=Standard_False;
 194   //
 195   //-- F= AA*CN*CN+2*BB*CN*SN+CC*CN+DD*SN+EE;
 196   math_TrigonometricFunctionRoots MTFR(CC, SC, C, S, Cte, Binf, Bsup);
 197   if(!MTFR.IsDone()) {
 198     return;
 199   }
 200   //
 201   done=Standard_True;
 202   if(MTFR.InfiniteRoots()) {
 203     infinite_roots=Standard_True;
 204     return;
 205   }
 206   //
 207   NbRoots=MTFR.NbSolutions();
 208   //
 209   for(i=0; i<NbRoots; ++i) {
 210     Roots[i]=MTFR.Value(i+1);
 211     if(Roots[i]<0.){
 212       Roots[i]+=PIpPI;
 213     }
 214     if(Roots[i]>PIpPI) {
 215       Roots[i]-=PIpPI;
 216     }
 217   }
 218   //
 219   //-- La recherche directe donne n importe quoi.
 220   SvNbRoots=NbRoots;
 221   for(i=0; i<SvNbRoots; ++i) {
 222     Standard_Real y;
 223     Standard_Real co=cos(Roots[i]);
 224     Standard_Real si=sin(Roots[i]);
 225     y=co*(CC*co + (SC+SC)*si + C) + S*si + Cte;
 226     if(Abs(y)>1e-8) {
 227       done=Standard_False;
 228       return; //-- le 1er avril 98
 229     }
 230   }
 231   //
 232   do {
 233     Triee=Standard_True;
 234     for(i=1,j=0; i<SvNbRoots; ++i,++j) {
 235       if(Roots[i]<Roots[j]) {
 236         Triee=Standard_False;
 237         Standard_Real t=Roots[i];
 238         Roots[i]=Roots[j];
 239         Roots[j]=t;
 240       }
 241     }
 242   }
 243   while(!Triee);
 244   //
 245   infinite_roots=Standard_False;
 246   //
 247   if(!NbRoots) {//--!!!!! Detection du cas Pol = Cte ( 1e-50 ) !!!!
 248     if((Abs(CC) + Abs(SC) + Abs(C) + Abs(S)) < 1e-10) {
 249       if(Abs(Cte) < 1e-10)  {
 250         infinite_roots=Standard_True;
 251       }
 252     }
 253   }
 254 }
 255 //=======================================================================
 256 //class    : MyTrigonometricFunction
 257 //purpose  :
 258 //  Classe interne : Donne Value et Derivative sur un polynome
 259 //                   trigonometrique
 260 //======================================================================
 261 class MyTrigonometricFunction {
 262
 263  private:
 264   Standard_Real CC,SS,SC,S,C,Cte;
 265
 266  public:
 267   //
 268   MyTrigonometricFunction(const Standard_Real xCC,
 269                           const Standard_Real xSS,
 270                           const Standard_Real xSC,
 271                           const Standard_Real xC,
 272                           const Standard_Real xS,
 273                           const Standard_Real xCte) {
 274     CC=xCC;
 275     SS=xSS;
 276     SC=xSC;
 277     S=xS;
 278     C=xC;
 279     Cte=xCte;
 280   }
 281
 282   Standard_Real Value(const Standard_Real& U) {
 283     Standard_Real sinus, cosinus, aRet;
 284     //
 285     sinus=sin(U);
 286     cosinus=cos(U);
 287     aRet= CC*cosinus*cosinus +
 288           SS*sinus*sinus +
 289           2.0*(sinus*(SC*cosinus+S)+cosinus*C)+
 290           Cte;
 291     //
 292     return aRet;
 293   }
 294   //
 295   Standard_Real Derivative(const Standard_Real& U) {
 296     Standard_Real sinus, cosinus;
 297     //
 298     sinus=sin(U);
 299     cosinus=cos(U);
 300     //
 301     return(2.0*((sinus*cosinus)*(SS-CC)
 302                 +S*cosinus
 303                 -C*sinus
 304                 +SC*(cosinus*cosinus-sinus*sinus)));
 305   }
 306 };
 307
 308 //////////
 309 //=======================================================================
 310 //function : IntAna_IntQuadQuad::IntAna_IntQuadQuad
 311 //purpose  : C o n s t r u c t e u r    v i d e
 312 //=======================================================================
 313 IntAna_IntQuadQuad::IntAna_IntQuadQuad(void) {
 314   done=Standard_False;
 315   identical = Standard_False;
 316   NbCurves=0;
 317   Nbpoints=0;
 318   myNbMaxCurves=12;
 319   myEpsilon=0.00000001;
 320   myEpsilonCoeffPolyNull=0.00000001;
 321   memset (nextcurve, 0, sizeof (nextcurve));
 322   memset (previouscurve, 0, sizeof (previouscurve));
 323 }
 324 //=======================================================================
 325 //function : IntAna_IntQuadQuad::IntAna_IntQuadQuad
 326 //purpose  : I n t e r s e c t i o n   C y l i n d r e   Q u a d r i q u e
 327 //=======================================================================
 328 IntAna_IntQuadQuad::IntAna_IntQuadQuad(const gp_Cylinder& Cyl,
 329                                        const IntAna_Quadric& Quad,
 330                                        const Standard_Real Tol) {
 331   myNbMaxCurves=12;
 332   myEpsilon=0.00000001;
 333   myEpsilonCoeffPolyNull=0.00000001;
 334   Perform(Cyl,Quad,Tol);
 335 }
 336 //=======================================================================
 337 //function : Perform
 338 //purpose  : I n t e r s e c t i o n   C y l i n d r e   Q u a d r i q u e
 339 //=======================================================================
 340 void IntAna_IntQuadQuad::Perform(const gp_Cylinder& Cyl,
 341                                  const IntAna_Quadric& Quad,
 342                                  const Standard_Real)
 343 {
 344   done = Standard_True;
 345   identical= Standard_False;
 346   NbCurves=0;
 347   Nbpoints=0;
 348   //
 349   Standard_Boolean UN_SEUL_Z_PAR_THETA, DEUX_Z_PAR_THETA,
 350                    Z_POSITIF, Z_INDIFFERENT, Z_NEGATIF;
 351   //
 352   UN_SEUL_Z_PAR_THETA=Standard_False;
 353   DEUX_Z_PAR_THETA=Standard_True;
 354   Z_POSITIF=Standard_True;
 355   Z_INDIFFERENT=Standard_True;
 356   Z_NEGATIF=Standard_False;
 357   //
 358   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, aRealEpsilon, RCyl, R2;
 359   Standard_Real PIpPI = M_PI + M_PI;
 360   //
 361   for(Standard_Integer raz = 0 ; raz < myNbMaxCurves ; raz++) {
 362     previouscurve[raz] = nextcurve[raz] = 0;
 363   }
 364   //
 365   RCyl=Cyl.Radius();
 366   aRealEpsilon=RealEpsilon();
 367   //----------------------------------------------------------------------
 368   //-- Coefficients de la quadrique :
 369   //--      2       2       2
 370   //-- Qxx x + Qyy y + Qzz z + 2 ( Qxy x y + Qxz x z + Qyz y z )
 371   //-- + 2 ( Qx x + Qy y + Qz z ) + Q1
 372   //----------------------------------------------------------------------
 373   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1);
 374
 375   //----------------------------------------------------------------------
 376   //-- Ecriture des Coefficients de la Quadrique dans le repere liee
 377   //-- au Cylindre
 378   //--                 Cyl.Position() -> Ax2
 379   //----------------------------------------------------------------------
 380   Quad.NewCoefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,Cyl.Position());
 381
 382   //----------------------------------------------------------------------
 383   //-- Parametrage du Cylindre Cyl :
 384   //--     X = Rcyl * Cos(Theta)
 385   //--     Y = Rcyl * Sin(Theta)
 386   //--     Z = Z
 387   //----------------------------------------------------------------------
 388   //-- Donne une Equation de la forme :
 389   //--   F(Sin(Theta),Cos(Theta),ZCyl) = 0
 390   //--   (Equation du second degre en ZCyl)
 391   //--    ZCyl**2  CoeffZ2(Theta) + ZCyl CoeffZ1(Theta) + CoeffZ0(Theta)
 392   //----------------------------------------------------------------------
 393   //-- CoeffZ0 = Q1 + 2*Qx*RCyl*Cos[Theta] + Qxx*RCyl^2*Cos[Theta]^2
 394   //--           2*RCyl*Sin[Theta]* ( Qy + Qxy*RCyl*Cos[Theta]);
 395   //--           Qyy*RCyl^2*Sin[Theta]^2;
 396   //-- CoeffZ1 =2.0 * ( Qz + RCyl*(Qxz*Cos[Theta] + Sin[Theta]*Qyz)) ;
 397   //-- CoeffZ2 =  Qzz;
 398   //-- !!!! Attention , si on norme sur Qzz pour detecter le cas 1er degre
 399   //----------------------------------------------------------------------
 400   //-- On Cherche Les Racines en Theta du discriminant de cette equation :
 401   //-- DIS(Theta) = C_1 + C_SS Sin()**2 + C_CC Cos()**2 + 2 C_SC Sin() Cos()
 402   //--                  + 2 C_S  Sin()    + 2 C_C Cos()
 403   //--
 404   //-- Si Qzz = 0   Alors  On Resout Z=Fct(Theta)  sur le domaine de Theta
 405   //----------------------------------------------------------------------
 406
 407   if(Abs(Qzz)<myEpsilonCoeffPolyNull) {
 408     done=Standard_False; //-- le 12 mars 98
 409     return;
 410   }
 411   else { //#0
 412     //----------------------------------------------------------------------
 413     //--- Cas  ou  F(Z,Theta) est du second degre en Z                  ----
 414     //----------------------------------------------------------------------
 415     R2   = RCyl*RCyl;
 416
 417     Standard_Real C_1  = Qz*Qz - Qzz*Q1;
 418     Standard_Real C_SS = R2    * ( Qyz*Qyz - Qyy*Qzz );
 419     Standard_Real C_CC = R2    * ( Qxz*Qxz - Qxx*Qzz );
 420     Standard_Real C_S  = RCyl  * ( Qyz*Qz  - Qy*Qzz  );
 421     Standard_Real C_C  = RCyl  * ( Qxz*Qz  - Qx*Qzz  );
 422     Standard_Real C_SC = R2    * ( Qxz*Qyz - Qxy*Qzz );
 423     //
 424     MyTrigonometricFunction MTF(C_CC,C_SS,C_SC,C_C,C_S,C_1);
 425     TrigonometricRoots PolDIS(C_CC-C_SS,
 426                               C_SC,
 427                               C_C+C_C,
 428                               C_S+C_S,
 429                               C_1+C_SS, 0., PIpPI);
 430     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 431     if(PolDIS.IsDone()==Standard_False) {
 432       done=Standard_False;
 433       return;
 434     }
 435     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 436     if(PolDIS.InfiniteRoots()) {
 437       TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 438                                         myEpsilon,0.,PIpPI,
 439                                         UN_SEUL_Z_PAR_THETA,
 440                                         Z_POSITIF);
 441       TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 442                                         myEpsilon,0.,PIpPI,
 443                                         UN_SEUL_Z_PAR_THETA,
 444                                         Z_NEGATIF);
 445       NbCurves=2;
 446     }
 447
 448     else { //#1
 449       //---------------------------------------------------------------
 450       //-- La Recherche des Zero de DIS a reussi
 451       //---------------------------------------------------------------
 452       Standard_Integer nbsolDIS=PolDIS.NbSolutions();
 453       if(nbsolDIS==0) {
 454         //--------------------------------------------------------------
 455         //-- Le Discriminant a un signe constant :
 456         //--
 457         //-- Si Positif  ---> 2 Courbes
 458         //-- Sinon       ---> Pas de solution
 459         //--------------------------------------------------------------
 460         if(MTF.Value(M_PI) >= -aRealEpsilon) {
 461
 462           TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 463                                             myEpsilon,0.0,PIpPI,
 464                                             UN_SEUL_Z_PAR_THETA,
 465                                             Z_POSITIF);
 466           TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 467                                             myEpsilon,0.0,PIpPI,
 468                                             UN_SEUL_Z_PAR_THETA,
 469                                             Z_NEGATIF);
 470
 471           NbCurves=2;
 472         }
 473         else {
 474           //------------------------------------------------------------
 475           //-- Le Discriminant est toujours Negatif
 476           //------------------------------------------------------------
 477           NbCurves=0;
 478         }
 479       }
 480       else { //#2
 481         //------------------------------------------------------------
 482         //-- Le Discriminant a des racines
 483         //-- (Le Discriminant est une fonction periodique donc ... )
 484         //------------------------------------------------------------
 485         if( nbsolDIS==1 ) {
 486           //------------------------------------------------------------
 487           //-- Point de Tangence pour ce Theta Solution
 488           //-- Si Signe de Discriminant >=0 pour tout Theta
 489           //--     Alors
 490           //--       Courbe Solution pour la partie Positive
 491           //--       entre les 2 racines ( Ici Tout le domaine )
 492           //-- Sinon Seulement un point Tangent
 493           //------------------------------------------------------------
 494           if(MTF.Value(PolDIS.Value(1)+M_PI) >= -aRealEpsilon ) {
 495             //------------------------------------------------------------
 496             //-- On a Un Point de Tangence + Une Courbe Solution
 497             //------------------------------------------------------------
 498             TheCurve[0].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 499                                               myEpsilon,0.0,PIpPI,
 500                                               UN_SEUL_Z_PAR_THETA,
 501                                               Z_POSITIF);
 502             TheCurve[1].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 503                                               myEpsilon,0.0,PIpPI,
 504                                               UN_SEUL_Z_PAR_THETA,
 505                                               Z_NEGATIF);
 506
 507             NbCurves=2;
 508           }
 509           else {
 510             //------------------------------------------------------------
 511             //-- On a simplement un Point de tangence
 512             //------------------------------------------------------------
 513             //--Standard_Real Theta = PolDIS(1);
 514             //--Standard_Real SecPar= -0.5 * MTFZ1.Value(Theta) / MTFZ2.Value(Theta);
 515             //--Thepoints[Nbpoints] = ElSLib::CylinderValue(Theta,SecPar,Cyl.Position(),Cyl.Radius());
 516             //--Nbpoints++;
 517             NbCurves=0;
 518           }
 519         }
 520         else {  // #3
 521           //------------------------------------------------------------
 522           //-- On detecte   :  Les racines double
 523           //--              :  Les Zones Positives [Modulo 2 PI]
 524           //-- Les Courbes solutions seront obtenues pour les valeurs
 525           //-- de Theta ou Discriminant(Theta) > 0  (>=0? en limite)
 526           //-- Par la resolution de l equation implicite avec Theta fixe
 527           //------------------------------------------------------------
 528           //-- Si tout est solution, Alors on est sur une iso
 529           //-- Ce cas devrait etre traite en amont
 530           //------------------------------------------------------------
 531           //-- On construit la fonction permettant connaissant un Theta
 532           //-- de calculer les Z+ et Z- Correspondants.
 533           //-------------------------------------------------------------
 534
 535           //-------------------------------------------------------------
 536           //-- Calcul des Intervalles ou Discriminant >=0
 537           //-- On a au plus nbsol intervalles ( en fait 2 )
 538           //--  (1,2) (2,3) .. (nbsol,1+2PI)
 539           //-------------------------------------------------------------
 540           Standard_Integer i;
 541           Standard_Real Theta1, Theta2, Theta3, autrepar, qwet;
 542           Standard_Boolean UnPtTg = Standard_False;
 543           //
 544           NbCurves=0;
 545           if(nbsolDIS == 2) {
 546             for(i=1; i<=nbsolDIS ; i++) {
 547               Theta1=PolDIS.Value(i);
 548               Theta2=(i<nbsolDIS)? PolDIS.Value(i+1)  :  (PolDIS.Value(1)+PIpPI);
 549               //----------------------------------------------------------------
 550               //-- On detecte les racines doubles
 551               //-- Si il n y a que 2 racines alors on prend tout l interval
 552               //----------------------------------------------------------------
 553               if(Abs(Theta2-Theta1)<=aRealEpsilon) {
 554                 UnPtTg = Standard_True;
 555                 autrepar=Theta1-0.1;
 556                 if(autrepar<0.) {
 557                   autrepar=Theta1+0.1;
 558                 }
 559                 //
 560                 qwet=MTF.Value(autrepar);
 561                 if(qwet>=0.) {
 562                   Standard_Real aParam = Theta1 + PIpPI;
 563                   AddSpecialPoints(Quad, Cyl, Theta1, aParam);
 564                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 565                                                            myEpsilon,Theta1,aParam,
 566                                                            UN_SEUL_Z_PAR_THETA,
 567                                                            Z_POSITIF);
 568                   NbCurves++;
 569                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 570                                                            myEpsilon,Theta1,aParam,
 571                                                            UN_SEUL_Z_PAR_THETA,
 572                                                            Z_NEGATIF);
 573                   NbCurves++;
 574                 }
 575               }
 576             }
 577           }
 578
 579           for(i=1; UnPtTg == (Standard_False) && (i<=nbsolDIS) ; i++) {
 580             Theta1=PolDIS.Value(i);
 581             Theta2=(i<nbsolDIS)? PolDIS.Value(i+1) : (PolDIS.Value(1)+PIpPI);
 582             //----------------------------------------------------------------
 583             //-- On detecte les racines doubles
 584             //-- Si il n y a que 2 racines alors on prend tout l interval
 585             //----------------------------------------------------------------
 586             if(Abs(Theta2-Theta1)<=1e-12) {
 587               //-- std::cout<<"\n####### Un Point de Tangence en Theta = "<<Theta1<<std::endl;
 588               //-- i++;
 589             }
 590             else {  //-- On evite les pbs numeriques (Tout est du meme signe entre les racines)
 591               qwet=MTF.Value(0.5*(Theta1+Theta2))
 592                     +MTF.Value(0.4*Theta1+0.6*Theta2)
 593                       +MTF.Value(0.6*Theta1+0.4*Theta2);
 594               if(qwet >= 0.) {
 595                 //------------------------------------------------------------
 596                 //-- On est positif a partir de Theta1
 597                 //-- L intervalle Theta1,Theta2 est retenu
 598                 //------------------------------------------------------------
 599
 600                 //-- le 8 octobre 1997 :
 601                 //-- PB: Racine en 0    pi/2 pi/2  et 2pi
 602                 //-- On cree 2 courbes : 0    -> pi/2    2zpartheta
 603                 //--                     pi/2 -> 2pi     2zpartheta
 604                 //--
 605                 //-- la courbe 0 pi/2   passe par le double pt de tangence pi/2
 606                 //-- il faut donc couper cette courbe en 2
 607                 //--
 608                 Theta3 = ((i+1)<nbsolDIS)? PolDIS.Value(i+2)  :  (PolDIS.Value(1)+PIpPI);
 609                 //ft
 610                 if((Theta3-Theta2)<5.e-8) {
 611                 //
 612                   AddSpecialPoints(Quad, Cyl, Theta1, Theta2);
 613                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 614                                                            myEpsilon,Theta1,Theta2,
 615                                                            UN_SEUL_Z_PAR_THETA,
 616                                                            Z_POSITIF);
 617                   NbCurves++;
 618                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 619                                                            myEpsilon,Theta1,Theta2,
 620                                                            UN_SEUL_Z_PAR_THETA,
 621                                                            Z_NEGATIF);
 622                   NbCurves++;
 623                 }
 624                 else {
 625                   AddSpecialPoints(Quad, Cyl, Theta1, Theta2);
 626                   TheCurve[NbCurves].SetCylinderQuadValues(Cyl,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 627                                                            myEpsilon,Theta1,Theta2,
 628                                                            DEUX_Z_PAR_THETA,
 629                                                            Z_INDIFFERENT);
 630                   NbCurves++;
 631                 }
 632               }//if(qwet >= 0.)
 633             }//else {  //-- On evite les pbs numerique ...
 634           } //for(i=1; UnPtTg == (Standard_False) && (i<=nbsolDIS) ; i++) {
 635         }//else { // #3
 636       }//else { //#2
 637     }//else { //#1
 638
 639   }//else { //#0
 640 }
 641
 642 //=======================================================================
 643 //function :IntAna_IntQuadQuad::IntAna_IntQuadQuad
 644 //purpose  :
 645 //=======================================================================
 646 IntAna_IntQuadQuad::IntAna_IntQuadQuad(const gp_Cone& Cone,
 647                                        const IntAna_Quadric& Quad,
 648                                        const Standard_Real Tol)
 649 {
 650   myNbMaxCurves=12;
 651   myEpsilon=0.00000001;
 652   myEpsilonCoeffPolyNull=0.00000001;
 653   Perform(Cone,Quad,Tol);
 654 }
 655 //=======================================================================
 656 //function :Perform
 657 //purpose  :
 658 //=======================================================================
 659 void IntAna_IntQuadQuad::Perform(const gp_Cone& Cone,
 660                                  const IntAna_Quadric& Quad,
 661                                  const Standard_Real)
 662 {
 663   //
 664   Standard_Boolean UN_SEUL_Z_PAR_THETA,
 665                    Z_POSITIF, Z_NEGATIF;
 666   //
 667   UN_SEUL_Z_PAR_THETA=Standard_False;
 668   Z_POSITIF=Standard_True;
 669   Z_NEGATIF=Standard_False;
 670   //
 671   Standard_Integer i;
 672   Standard_Real Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1;
 673   Standard_Real Theta1, Theta2, TgAngle;
 674   Standard_Real PIpPI = M_PI + M_PI;
 675   //
 676   done=Standard_True;
 677   identical = Standard_False;
 678   NbCurves=0;
 679   Nbpoints=0;
 680   //
 681   for(i=0 ; i<myNbMaxCurves ; ++i) {
 682     previouscurve[i]=0;
 683     nextcurve[i]=0;
 684   }
 685   //
 686   Quad.Coefficients(Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1);
 687   //
 688   gp_Ax3 tAx3(Cone.Position());
 689   tAx3.SetLocation(Cone.Apex());
 690   Quad.NewCoefficients(Qxx,Qyy,Qzz,
 691                        Qxy,Qxz,Qyz,
 692                        Qx,Qy,Qz,Q1,
 693                        tAx3);
 694   //
 695   TgAngle=1./(Tan(Cone.SemiAngle()));
 696   //
 697   // The parametrization of the Cone
 698   //
 699   // x= z*tan(beta)*cos(t)
 700   // y= z*tan(beta)*sin(t)
 701   // z=z
 702   //
 703   // The intersection curves are defined by the equation
 704   //
 705   //                        2
 706   //          f(z,t)= A(t)*z + B(t)*z + C(t)=0
 707   //
 708   //
 709   // 1. Try to solve A(t)=0 -> PolZ2
 710   //
 711   Standard_Integer nbsol, nbsol1, nbsolZ2;
 712   Standard_Real Z2CC, Z2SS, Z2Cte, Z2SC, Z2C, Z2S;
 713   Standard_Real Z1CC, Z1SS, Z1Cte, Z1SC, Z1C, Z1S;
 714   Standard_Real C_1,C_SS, C_CC, C_S, C_C, C_SC;
 715   //
 716   Z2CC = Qxx;
 717   Z2SS = Qyy;
 718   Z2Cte= Qzz * TgAngle*TgAngle;
 719   Z2SC = Qxy;
 720   Z2C  = (TgAngle)*Qxz;
 721   Z2S  = (TgAngle)*Qyz;
 722   //
 723   TrigonometricRoots PolZ2(Z2CC-Z2SS,Z2SC,Z2C+Z2C,Z2S+Z2S,Z2Cte+Z2SS,0.,PIpPI);
 724   if(!PolZ2.IsDone()) {
 725     done=!done;
 726     return;
 727   }
 728   //
 729   //MyTrigonometricFunction MTF2(Z2CC,Z2SS,Z2SC,Z2C,Z2S,Z2Cte);
 730   nbsolZ2 = PolZ2.NbSolutions();
 731   //
 732   // 2. Try to solve B(t)=0  -> PolZ1
 733   //
 734   Z1Cte = 2.*(TgAngle)*Qz;
 735   Z1SS  = 0.;
 736   Z1CC  = 0.;
 737   Z1S   = Qy;
 738   Z1C   = Qx;
 739   Z1SC  = 0.;
 740   //
 741   TrigonometricRoots PolZ1(Z1CC-Z1SS,Z1SC, Z1C+Z1C,Z1S+Z1S, Z1Cte+Z1SS,0.,PIpPI);
 742   if(!PolZ1.IsDone()) {
 743     done=!done;
 744     return;
 745   }
 746   MyTrigonometricFunction MTFZ1(Z1CC,Z1SS,Z1SC,Z1C,Z1S,Z1Cte);
 747   //
 748   nbsol1=PolZ1.NbSolutions();
 749   if(PolZ2.InfiniteRoots()) { // i.e A(t)=0 for any t
 750     if(!PolZ1.InfiniteRoots()) {
 751       if(!nbsol1) {
 752         //-- B(t)*z + C(t) = 0    avec C(t) != 0
 753         TheCurve[0].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,
 754                                       myEpsilon,0.,PIpPI,
 755                                       UN_SEUL_Z_PAR_THETA,
 756                                       Z_POSITIF);
 757         NbCurves=1;
 758       }
 759       else {
 760         /*
 761         Standard_Integer ii;
 762         for(ii=1; ii<= nbsol1 ; ++ii) {
 763           Standard_Real Theta=PolZ1.Value(ii);
 764           if(Abs(MTFZ1.Value(Theta))<=myEpsilon) {
 765             //-- Une droite Solution  Z=  -INF -> +INF  pour Theta
 766             //-- std::cout<<"######## Droite Solution Pour Theta = "<<Theta<<std::endl;
 767           }
 768           else {
 769             //-- std::cout<<"\n#### _+_+_+_+_+ CAS A(t) Z + B = 0   avec A et B ==0 "<<std::endl;
 770           }
 771         }
 772         */
 773       }
 774     }
 775     else {
 776       if(Abs(Q1)<=myEpsilon) {
 777         done=!done;
 778         return;
 779       }
 780       else {
 781         //-- Pas de Solutions
 782       }
 783     }
 784     return;
 785   }
 786   //
 787   //else { //#5
 788   //
 789   //                2
 790   //-- f(z,t)=A(t)*z + B(t)*z + C(t)=0   avec A n est pas toujours nul
 791   //
 792   //                                              2
 793   // 3. Try to explore s.c. Discriminant:  D(t)=B(t)-4*A(t)*C(t) => Pol
 794   //
 795   // The function D(t) is just a formula that has sense for quadratic
 796   // equation above.
 797   // For cases when A(t)=0 (say at t=ti, t=tj. etc) the equation
 798   // will be
 799   //
 800   //     f(z,t)=B(t)*z + C(t)=0, where B(t)!=0,
 801   //
 802   // and the D(t) is nonsense for it.
 803   //
 804   C_1  = TgAngle*TgAngle * ( Qz * Qz - Qzz * Q1);
 805   C_SS = Qy * Qy - Qyy * Q1;
 806   C_CC = Qx * Qx - Qxx * Q1;
 807   C_S  = TgAngle*( Qy * Qz - Qyz * Q1);
 808   C_C  = TgAngle*( Qx * Qz - Qxz * Q1);
 809   C_SC = Qx * Qy - Qxy * Q1;
 810   //
 811   TrigonometricRoots Pol(C_CC-C_SS, C_SC, C_C+C_C,C_S+C_S, C_1+C_SS,0.,PIpPI);
 812   if(!Pol.IsDone()) {
 813     done=!done;
 814     return;
 815   }
 816   //
 817   nbsol=Pol.NbSolutions();
 818   MyTrigonometricFunction MTF(C_CC,C_SS,C_SC,C_C,C_S,C_1);
 819   //
 820   if(Pol.InfiniteRoots()) {
 821     //                             2
 822     //         f(z,t)= (z(t)-zo(t)) = 0       Discriminant(t)=0 pour tout t
 823     //
 824     TheCurve[0].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
 825                                   0.,PIpPI,
 826                                   UN_SEUL_Z_PAR_THETA, Z_POSITIF);
 827     TheCurve[1].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
 828                                   0.,PIpPI,
 829                                   UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
 830     NbCurves=2;
 831     return;
 832   }
 833   //else {//#4
 834   //                     2
 835   //        f(z,t)=A(t)*z + B(t)*z + C(t)      Discriminant(t) != 0
 836   //
 837   if(!nbsol && (MTF.Value(M_PI)<0.) ) {
 838     //-- Discriminant signe constant negatif
 839     return;
 840   }
 841   //else {//# 3
 842   //
 843   //-- On Traite le cas : Discriminant(t) > 0 pour tout t en simulant 1
 844   //   racine double en 0
 845   Standard_Boolean DiscriminantConstantPositif = Standard_False;
 846   if(!nbsol) {
 847     nbsol=1;
 848     DiscriminantConstantPositif = Standard_True;
 849   }
 850   //----------------------------------------------------------------------
 851   //-- Le discriminant admet au moins une racine ( -> Point de Tangence )
 852   //-- ou est constant positif.
 853   //----------------------------------------------------------------------
 854   for(i=1; i<=nbsol; ++i) {
 855     if(DiscriminantConstantPositif) {
 856       Theta1 = 0.;
 857       Theta2 = PIpPI-myEpsilon;
 858     }
 859     else {
 860       Theta1=Pol.Value(i);
 861       Theta2=(i<nbsol)? Pol.Value(i+1) : (Pol.Value(1)+PIpPI);
 862     }
 863     //
 864     if(Abs(Theta2-Theta1)<=myEpsilon){
 865       done=Standard_False;
 866       return;// !!! pkv
 867     }
 868     //else {// #2
 869     Standard_Real qwet =MTF.Value(0.5*(Theta1+Theta2))+
 870                         MTF.Value(0.4*Theta1+0.6*Theta2)+
 871                         MTF.Value(0.6*Theta1+0.4*Theta2);
 872     if(qwet < 0.) {
 873       continue;
 874     }
 875     //---------------------------------------------------------------------
 876     //-- On a des Solutions entre Theta1 et Theta 2
 877     //---------------------------------------------------------------------
 878
 879     //---------------------------------------------------------------------
 880     //-- On Subdivise si necessaire Theta1-->Theta2
 881     //-- en Theta1-->t1   t1--->t2   ...  tn--->Theta2
 882     //--  ou t1,t2,...tn sont des racines de A(t)
 883     //--
 884     //-- Seule la courbe Z_NEGATIF est affectee
 885     //----------------------------------------------------------------------
 886     Standard_Boolean RacinesdePolZ2DansTheta1Theta2;
 887     Standard_Integer i2;
 888     Standard_Real r;
 889     //
 890     //nbsolZ2=PolZ2.NbSolutions();
 891     RacinesdePolZ2DansTheta1Theta2=Standard_False;
 892     for(i2=1; i2<=nbsolZ2 && !RacinesdePolZ2DansTheta1Theta2; ++i2) {
 893       r=PolZ2.Value(i2);
 894       if(r>Theta1 && r<Theta2) {
 895         RacinesdePolZ2DansTheta1Theta2=Standard_True;
 896       }
 897       else {
 898         r+=PIpPI;
 899         if(r>Theta1 && r<Theta2){
 900           RacinesdePolZ2DansTheta1Theta2=Standard_True;
 901         }
 902       }
 903     }
 904     //
 905     if(!RacinesdePolZ2DansTheta1Theta2) {
 906       //--------------------------------------------------------------------
 907       //-- Pas de Branches Infinies
 908       TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,myEpsilon,
 909                                            Theta1,Theta2,
 910                                            UN_SEUL_Z_PAR_THETA,Z_POSITIF);
 911       NbCurves++;
 912       TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,myEpsilon,
 913                                            Theta1,Theta2,
 914                                            UN_SEUL_Z_PAR_THETA,
 915                                            Z_NEGATIF);
 916       NbCurves++;
 917     }
 918
 919     else { //#1
 920       Standard_Boolean NoChanges;
 921       Standard_Real NewMin, NewMax, to;
 922       //
 923       NewMin=Theta1;
 924       NewMax=Theta2;
 925       NoChanges=Standard_True;
 926       //
 927       for(i2=1; i2<= (nbsolZ2+nbsolZ2) ; ++i2) {
 928         if(i2>nbsolZ2) {
 929           to=PolZ2.Value(i2-nbsolZ2) + PIpPI;
 930         }
 931         else {
 932           to=PolZ2.Value(i2);
 933         }
 934         //
 935         // to is value of the parameter t where A(t)=0, i.e.
 936         // f(z,to)=B(to)*z + C(to)=0, B(to)!=0.
 937         //
 938         // z=-C(to)/B(to) is one and only
 939         // solution in spite of the fact that D(t)>0 for any value of t.
 940         //
 941         if(to<NewMax && to>NewMin) {
 942           //-----------------------------------------------------------------
 943           //-- On coupe au moins une fois le domaine Theta1 Theta2
 944           //-----------------------------------------------------------------
 945           NoChanges=Standard_False;
 946           TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
 947                                                NewMin,to,
 948                                                UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
 949           //
 950           NbCurves++;
 951           TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
 952                                                NewMin,to,
 953                                                UN_SEUL_Z_PAR_THETA, Z_POSITIF);
 954           //------------------------------------------------------------
 955           //-- A == 0
 956           //-- Si B < 0    Alors Infini sur   Z+
 957           //-- Sinon             Infini sur   Z-
 958           //------------------------------------------------------------
 959           if(PolZ2.IsARoot(NewMin)) {
 960             if(MTFZ1.Value(NewMin) < 0.) {
 961               TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
 962             }
 963             else {
 964               TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
 965             }
 966           }
 967           if(MTFZ1.Value(to) < 0.) {
 968             TheCurve[NbCurves].SetIsLastOpen(Standard_True);
 969           }
 970           else {
 971             TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
 972           }
 973           //------------------------------------------------------------
 974           NbCurves++;
 975           NewMin=to;
 976         }//if(to<NewMax && to>NewMin)
 977       }// for(i2=1; i2<= (nbsolZ2+nbsolZ2) ; ++i2)
 978       //
 979       if(NoChanges) {
 980         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1,  myEpsilon,
 981                                              Theta1,Theta2,
 982                                              UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
 983         NbCurves++;
 984         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
 985                                              Theta1,Theta2,
 986                                              UN_SEUL_Z_PAR_THETA, Z_POSITIF);
 987         //------------------------------------------------------------
 988         //-- A == 0
 989         //-- Si B < 0    Alors Infini sur   Z+
 990         //-- Sinon             Infini sur   Z-
 991         //------------------------------------------------------------
 992         if(PolZ2.IsARoot(Theta1)) {
 993           if(MTFZ1.Value(Theta1) < 0.) {
 994             TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
 995           }
 996           else {
 997             TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
 998           }
 999         }
1000         if(PolZ2.IsARoot(Theta2)) {
1001           if(MTFZ1.Value(Theta2) < 0.) {
1002             TheCurve[NbCurves].SetIsLastOpen(Standard_True);
1003           }
1004           else {
1005             TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
1006           }
1007         }
1008         //------------------------------------------------------------
1009         NbCurves++;
1010       }//if(NoChanges) {
1011       else {// #0
1012         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
1013                                              NewMin,Theta2,
1014                                              UN_SEUL_Z_PAR_THETA, Z_NEGATIF);
1015         NbCurves++;
1016         TheCurve[NbCurves].SetConeQuadValues(Cone,Qxx,Qyy,Qzz,Qxy,Qxz,Qyz,Qx,Qy,Qz,Q1, myEpsilon,
1017                                              NewMin,Theta2,
1018                                              UN_SEUL_Z_PAR_THETA, Z_POSITIF);
1019         //------------------------------------------------------------
1020         //-- A == 0
1021         //-- Si B < 0    Alors Infini sur   Z+
1022         //-- Sinon             Infini sur   Z-
1023         //------------------------------------------------------------
1024         if(PolZ2.IsARoot(NewMin)) {
1025           if(MTFZ1.Value(NewMin) < 0.) {
1026             TheCurve[NbCurves].SetIsFirstOpen(Standard_True);
1027           }
1028           else {
1029             TheCurve[NbCurves-1].SetIsFirstOpen(Standard_True);
1030           }
1031         }
1032         if(PolZ2.IsARoot(Theta2)) {
1033           if(MTFZ1.Value(Theta2) < 0.) {
1034             TheCurve[NbCurves].SetIsLastOpen(Standard_True);
1035           }
1036           else {
1037             TheCurve[NbCurves-1].SetIsLastOpen(Standard_True);
1038           }
1039         }
1040         //------------------------------------------------------------
1041
1042         NbCurves++;
1043       }//else {// #0
1044     }//else { //#1
1045   }//for(i=1; i<=nbsol ; i++) {
1046   //}//else { //#5
1047   InternalSetNextAndPrevious();
1048 }
1049 //=======================================================================
1050 //function :InternalSetNextAndPrevious
1051 //purpose  :
1052 //-- Raccordement des courbes bout a bout
1053 //--    - Utilisation du champ : IsFirstOpen
1054 //--    -                        IsLastOpen
1055 //-- Pas de verification si ces champs retournent Vrai.
1056 //--
1057 //--
1058 //-- 1: Test de      Fin(C1)    = Debut(C2)   ->N(C1) et P(C2)
1059 //-- 2:              Debut(C1)  = Fin(C2)     ->P(C1) et N(C2)
1060 //=======================================================================
1061 void IntAna_IntQuadQuad::InternalSetNextAndPrevious()
1062 {
1063   Standard_Boolean NotLastOpenC2, NotFirstOpenC2;
1064   Standard_Integer c1,c2;
1065   Standard_Real aEps, aEPSILON_DISTANCE;
1066   Standard_Real DInfC1,DSupC1, DInfC2,DSupC2;
1067   //
1068   aEps=0.0000001;
1069   aEPSILON_DISTANCE=0.0000000001;
1070   //
1071   for(c1=0; c1<NbCurves; c1++) {
1072     nextcurve[c1] =0;
1073     previouscurve[c1] = 0;
1074   }
1075   //
1076   for(c1=0; c1 < NbCurves; c1++) {
1077     TheCurve[c1].Domain(DInfC1,DSupC1);
1078
1079     for(c2 = 0; (c2 < NbCurves) && (c2!=c1) ; c2++) {
1080
1081       NotLastOpenC2  = ! TheCurve[c2].IsLastOpen();
1082       NotFirstOpenC2 = ! TheCurve[c2].IsFirstOpen();
1083       TheCurve[c2].Domain(DInfC2,DSupC2);
1084       if(TheCurve[c1].IsFirstOpen() == Standard_False) {
1085         if(NotLastOpenC2) {
1086           if(Abs(DInfC1-DSupC2)<=aEps &&
1087              (TheCurve[c1].Value(DInfC1)
1088               .Distance(TheCurve[c2].Value(DSupC2))<aEPSILON_DISTANCE)) {
1089             previouscurve[c1] = c2+1;
1090             nextcurve[c2]     = c1+1;
1091           }
1092         }
1093         if(NotFirstOpenC2) {
1094           if(Abs(DInfC1-DInfC2)<=aEps
1095              && (TheCurve[c1].Value(DInfC1)
1096                  .Distance(TheCurve[c2].Value(DInfC2))<aEPSILON_DISTANCE)) {
1097             previouscurve[c1] = -(c2+1);
1098             previouscurve[c2] = -(c1+1);
1099           }
1100         }
1101       }
1102       if(TheCurve[c1].IsLastOpen() == Standard_False) {
1103         if(NotLastOpenC2) {
1104           if(Abs(DSupC1-DSupC2)<=aEps
1105              && (TheCurve[c1].Value(DSupC1)
1106                  .Distance(TheCurve[c2].Value(DSupC2))<aEPSILON_DISTANCE)) {
1107
1108             nextcurve[c1] = -(c2+1);
1109             nextcurve[c2] = -(c1+1);
1110           }
1111         }
1112         if(NotFirstOpenC2) {
1113           if(Abs(DSupC1-DInfC2)<=aEps
1114              && (TheCurve[c1].Value(DSupC1)
1115                  .Distance(TheCurve[c2].Value(DInfC2))<aEPSILON_DISTANCE)) {
1116             nextcurve[c1]     = c2+1;
1117             previouscurve[c2] = c1+1;
1118           }
1119         }
1120       }
1121     }
1122   }
1123 }
1124 //=======================================================================
1125 //function :HasPreviousCurve
1126 //purpose  :
1127 //=======================================================================
1128 Standard_Boolean IntAna_IntQuadQuad::HasPreviousCurve(const Standard_Integer I) const
1129 {
1130   if(!done) {
1131     throw StdFail_NotDone("IntQuadQuad Not done");
1132   }
1133   if (identical) {
1134     throw Standard_DomainError("IntQuadQuad identical");
1135   }
1136   if((I>NbCurves)||(I<=0)) {
1137     throw Standard_OutOfRange("Incorrect Curve Number 'HasPrevious Curve'");
1138   }
1139   if(previouscurve[I-1]) {
1140     return Standard_True;
1141   }
1142   return Standard_False;
1143 }
1144 //=======================================================================
1145 //function :HasNextCurve
1146 //purpose  :
1147 //=======================================================================
1148 Standard_Boolean IntAna_IntQuadQuad::HasNextCurve(const Standard_Integer I) const
1149 {
1150   if(!done) {
1151     throw StdFail_NotDone("IntQuadQuad Not done");
1152   }
1153   if (identical) {
1154     throw Standard_DomainError("IntQuadQuad identical");
1155   }
1156   if((I>NbCurves)||(I<=0)) {
1157     throw Standard_OutOfRange("Incorrect Curve Number 'HasNextCurve'");
1158   }
1159   if(nextcurve[I-1]) {
1160     return Standard_True;
1161   }
1162   return(Standard_False);
1163 }
1164 //=======================================================================
1165 //function :PreviousCurve
1166 //purpose  :
1167 //=======================================================================
1168 Standard_Integer IntAna_IntQuadQuad::PreviousCurve  (const Standard_Integer I,
1169                                                      Standard_Boolean& theOpposite) const
1170 {
1171   if(HasPreviousCurve(I)) {
1172     if(previouscurve[I-1]>0) {
1173       theOpposite = Standard_False;
1174       return(previouscurve[I-1]);
1175     }
1176     else {
1177       theOpposite = Standard_True;
1178       return( - previouscurve[I-1]);
1179     }
1180   }
1181   else {
1182     throw Standard_DomainError("Incorrect Curve Number 'PreviousCurve'");
1183   }
1184 }
1185 //=======================================================================
1186 //function :NextCurve
1187 //purpose  :
1188 //=======================================================================
1189 Standard_Integer IntAna_IntQuadQuad::NextCurve (const Standard_Integer I,
1190                                                 Standard_Boolean& theOpposite) const
1191 {
1192   if(HasNextCurve(I)) {
1193     if(nextcurve[I]>0) {
1194       theOpposite = Standard_False;
1195       return(nextcurve[I-1]);
1196     }
1197     else {
1198       theOpposite = Standard_True;
1199       return( - nextcurve[I-1]);
1200     }
1201   }
1202   else {
1203     throw Standard_DomainError("Incorrect Curve Number 'NextCurve'");
1204   }
1205 }
1206 //=======================================================================
1207 //function :Curve
1208 //purpose  :
1209 //=======================================================================
1210 const IntAna_Curve& IntAna_IntQuadQuad::Curve(const Standard_Integer i) const
1211 {
1212   if(!done) {
1213     throw StdFail_NotDone("IntQuadQuad Not done");
1214   }
1215   if (identical) {
1216     throw Standard_DomainError("IntQuadQuad identical");
1217   }
1218   if((i <= 0) || (i>NbCurves)) {
1219     throw Standard_OutOfRange("Incorrect Curve Number");
1220   }
1221   return TheCurve[i-1];
1222 }
1223 //=======================================================================
1224 //function :Point
1225 //purpose  :
1226 //=======================================================================
1227 const gp_Pnt& IntAna_IntQuadQuad::Point (const Standard_Integer i) const
1228 {
1229   if(!done) {
1230     throw StdFail_NotDone("IntQuadQuad Not done");
1231   }
1232   if (identical) {
1233     throw Standard_DomainError("IntQuadQuad identical");
1234   }
1235   if((i <= 0) || (i>Nbpoints)) {
1236     throw Standard_OutOfRange("Incorrect Point Number");
1237   }
1238   return Thepoints[i-1];
1239 }
1240 //=======================================================================
1241 //function :Parameters
1242 //purpose  :
1243 //=======================================================================
1244 void IntAna_IntQuadQuad::Parameters (const Standard_Integer, //i,
1245                                      Standard_Real& ,
1246                                      Standard_Real& ) const
1247 {
1248   std::cout << "IntAna_IntQuadQuad::Parameters(...) is not yet implemented" << std::endl;
1249 }
1250
1251 /*********************************************************************************
1252
1253 Mathematica Pour Cone Quadrique
1254 In[6]:= y[r_,t_]:=r*Sin[t]
1255
1256 In[7]:= x[r_,t_]:=r*Cos[t]
1257
1258 In[8]:= z[r_,t_]:=r*d
1259
1260 In[14]:= Quad[x_,y_,z_]:=Qxx x x + Qyy y y + Qzz z z + 2 (Qxy x y + Qxz x z + Qyz y z + Qx x + Qy y + Qz z ) + Q1
1261
1262 In[15]:= Quad[x[r,t],y[r,t],z[r,t]]
1263
1264                2      2        2       2        2       2
1265 Out[15]= Q1 + d  Qzz r  + Qxx r  Cos[t]  + Qyy r  Sin[t]  +
1266
1267                                       2
1268 >    2 (d Qz r + Qx r Cos[t] + d Qxz r  Cos[t] + Qy r Sin[t] +
1269
1270                2               2
1271 >       d Qyz r  Sin[t] + Qxy r  Cos[t] Sin[t])
1272
1273 In[16]:= QQ=%
1274
1275
1276
1277 In[17]:= Collect[QQ,r]
1278 Collect[QQ,r]
1279
1280 Out[17]= Q1 + r (2 d Qz + 2 Qx Cos[t] + 2 Qy Sin[t]) +
1281
1282       2   2                                  2
1283 >    r  (d  Qzz + 2 d Qxz Cos[t] + Qxx Cos[t]  + 2 d Qyz Sin[t] +
1284
1285                                         2
1286 >       2 Qxy Cos[t] Sin[t] + Qyy Sin[t] )
1287 ********************************************************************************/
1288
1289
1290   //**********************************************************************
1291   //***                   C O N E   - Q U A D R I Q U E                ***
1292   //***   2    2                                  2                    ***
1293   //***  R  ( d  Qzz + 2 d Qxz Cos[t] + Qxx Cos[t]  + 2 d Qyz Sin[t] + ***
1294   //***                                                                ***
1295   //***       2 Qxy Cos[t] Sin[t] + Qyy Sin[t] )                       ***
1296   //***                                                                ***
1297   //*** + R  ( 2 d Qz + 2 Qx Cos[t] + 2 Qy Sin[t] )                    ***
1298   //***                                                                ***
1299   //*** + Q1                                                           ***
1300   //**********************************************************************
1301   //FortranForm= ( DIS = QQ1 QQ1 - 4 QQ0 QQ2  ) / 4
1302   //   -  d**2*Qz**2 - d**2*Qzz*Q1 + (Qx**2 - Qxx*Q1)*Cos(t)**2 +
1303   //   -   (2*d*Qy*Qz - 2*d*Qyz*Q1)*Sin(t) + (Qy**2 - Qyy*Q1)*Sin(t)**2 +
1304   //   -   Cos(t)*(2*d*Qx*Qz - 2*d*Qxz*Q1 + (2*Qx*Qy - 2*Qxy*Q1)*Sin(t))
1305   //**********************************************************************
1306 //modified by NIZNHY-PKV Fri Dec  2 10:56:03 2005f
1307 /*
1308 static
1309   void DumpCurve(const Standard_Integer aIndex,
1310                  IntAna_Curve& aC);
1311 //=======================================================================
1312 //function : DumpCurve
1313 //purpose  :
1314 //=======================================================================
1315 void DumpCurve(const Standard_Integer aIndex,
1316                IntAna_Curve& aC)
1317 {
1318   Standard_Boolean bIsOpen, bIsConstant, bIsFirstOpen, bIsLastOpen;
1319   Standard_Integer i, aNb;
1320   Standard_Real aT1, aT2, aT, dT;
1321   gp_Pnt aP;
1322   //
1323   aC.Domain(aT1, aT2);
1324   bIsOpen=aC.IsOpen();
1325   bIsConstant=aC.IsConstant();
1326   bIsFirstOpen=aC.IsFirstOpen();
1327   bIsLastOpen=aC.IsLastOpen();
1328   //
1329   printf("\n");
1330   printf(" * IntAna_Curve #%d*\n", aIndex);
1331   printf(" Domain: [ %lf, %lf ]\n", aT1, aT2);
1332   printf(" IsOpen=%d\n", bIsOpen);
1333   printf(" IsConstant=%d\n", bIsConstant);
1334   printf(" IsFirstOpen =%d\n", bIsFirstOpen);
1335   printf(" IsLastOpen =%d\n", bIsLastOpen);
1336   //
1337   aNb=11;
1338   dT=(aT2-aT1)/(aNb-1);
1339   for (i=0; i<aNb; ++i) {
1340     aT=aT1+i*dT;
1341     if (i==(aNb-1)){
1342       aT=aT2;
1343     }
1344     //
1345     aP=aC.Value(aT);
1346     printf("point p%d_%d %lf %lf %lf\n",
1347            aIndex, i, aP.X(), aP.Y(), aP.Z());
1348   }
1349   printf(" ---- end of curve ----\n");
1350 }
1351 */
1352 //modified by NIZNHY-PKV Fri Dec  2 10:42:16 2005t