src/IntAna/IntAna_IntQuadQuad.cxx

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