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