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