src/Contap/Contap_ContAna.cxx

   1 // File:      Contap_ContAna.cxx
   2 // Created:   Thu Mar  4 10:51:23 1993
   3 // Author:    Jacques GOUSSARD
   4 // Copyright: OPEN CASCADE 2000
   5
   6 #include <Contap_ContAna.ixx>
   7
   8 #include <gp_XYZ.hxx>
   9 #include <gp.hxx>
  10
  11 #define Tolpetit 1.e-8
  12
  13 Contap_ContAna::Contap_ContAna (): done(Standard_False) {}
  14
  15 void Contap_ContAna::Perform (const gp_Sphere& S,
  16                               const gp_Dir& D)
  17 {
  18   done  = Standard_False;
  19   typL  = GeomAbs_Circle;
  20   pt1   = S.Location();
  21   dir1  = D;
  22   if (Abs(D.Dot(S.XAxis().Direction())) < 0.9999999999999) {
  23     dir2 = D.Crossed(S.XAxis().Direction());
  24   }
  25   else {
  26     dir2 = D.Crossed(S.YAxis().Direction());
  27   }
  28   prm   = S.Radius();
  29   nbSol = 1;
  30   done  = Standard_True;
  31 }
  32
  33 void Contap_ContAna::Perform (const gp_Sphere& S,
  34                               const gp_Dir& D,
  35                               const Standard_Real Angle)
  36 {
  37   done  = Standard_False;
  38   typL  = GeomAbs_Circle;
  39
  40   dir1  = D;
  41   if (Abs(D.Dot(S.XAxis().Direction())) < 0.9999999999999) {
  42     dir2 = D.Crossed(S.XAxis().Direction());
  43   }
  44   else {
  45     dir2 = D.Crossed(S.YAxis().Direction());
  46   }
  47   Standard_Real alpha = (S.Direct() ? Angle : -Angle);
  48   pt1.SetXYZ(S.Location().XYZ() - S.Radius()*sin(alpha)*D.XYZ()) ;
  49   prm   = S.Radius()*cos(alpha);
  50   nbSol = 1;
  51   done  = Standard_True;
  52 }
  53
  54 void Contap_ContAna::Perform (const gp_Sphere& S,
  55                               const gp_Pnt& Eye)
  56 {
  57   done = Standard_False;
  58
  59   Standard_Real radius = S.Radius();
  60   Standard_Real dist = Eye.Distance(S.Location());
  61   if (dist <= radius) {
  62     nbSol = 0;
  63   }
  64   else {
  65     prm = radius*sqrt(1.-radius*radius/(dist*dist));
  66     if (prm < Tolpetit) {
  67       nbSol = 0;
  68     }
  69     else {
  70       gp_XYZ locxyz(S.Location().XYZ());
  71       dir1.SetXYZ(Eye.XYZ()-locxyz);
  72       pt1.SetXYZ(locxyz + (radius*radius/dist)*dir1.XYZ());
  73       if (Abs(dir1.Dot(S.XAxis().Direction())) < 0.9999999999999) {
  74         dir2 = dir1.Crossed(S.XAxis().Direction());
  75       }
  76       else {
  77         dir2 = dir1.Crossed(S.YAxis().Direction());
  78       }
  79       nbSol = 1;
  80       typL = GeomAbs_Circle;
  81     }
  82   }
  83   done = Standard_True;
  84 }
  85
  86 void Contap_ContAna::Perform (const gp_Cylinder& C,
  87                               const gp_Dir& D)
  88 {
  89   done = Standard_False;
  90
  91   gp_XYZ normale(C.Position().Direction().XYZ());
  92   normale.Cross(D.XYZ());
  93   if (normale.Modulus() <= 1e-15) {
  94     nbSol = 0;
  95   }
  96   else {
  97     normale.Normalize();
  98     typL = GeomAbs_Line;
  99     dir1 = C.Position().Direction();
 100     dir2 = dir1;
 101     pt1.SetXYZ(C.Location().XYZ() + C.Radius()*normale);
 102     pt2.SetXYZ(C.Location().XYZ() - C.Radius()*normale);
 103     nbSol = 2;
 104   }
 105
 106   done = Standard_True;
 107 }
 108
 109 void Contap_ContAna::Perform (const gp_Cylinder& C,
 110                               const gp_Dir& D,
 111                               const Standard_Real Angle)
 112 {
 113   done = Standard_False;
 114
 115   Standard_Real Coefcos = D.Dot(C.Position().XDirection());
 116   Standard_Real Coefsin = D.Dot(C.Position().YDirection());
 117   Standard_Real Coefcst = cos(M_PI*0.5 + Angle);
 118
 119   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 120   Standard_Real norm2 = sqrt(norm1);
 121
 122   if (Abs(Coefcst) < norm2) {
 123     typL = GeomAbs_Line;
 124     nbSol = 2;
 125     dir1 = dir2 = C.Position().Direction();
 126
 127     if (!C.Direct()) { // The normal is inverted.
 128       Coefcos = -Coefcos;
 129       Coefsin = -Coefsin;
 130     }
 131
 132     // Necessary to solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 133     // and the origins of solution are in the reference of the
 134     // cylinder in (R*cost0, R*sint0,0) and (R*cost1,R*sint1,0)
 135     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 136     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 137     // and by using trigonometric relations the values of cosinus
 138     // and sinus to the solutions are obtained.
 139
 140     prm = Sqrt(norm1 - Coefcst*Coefcst);
 141     Standard_Real cost0,sint0,cost1,sint1;
 142
 143     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 144     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 145
 146     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 147     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 148
 149     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 150     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 151     gp_XYZ dirxyz;
 152
 153     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir);
 154     dirxyz.Multiply(C.Radius());
 155     pt1.SetXYZ(C.Location().XYZ().Added(dirxyz));
 156
 157     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir);
 158     dirxyz.Multiply(C.Radius());
 159     pt2.SetXYZ(C.Location().XYZ().Added(dirxyz));
 160   }
 161   else {
 162     nbSol = 0;
 163   }
 164
 165   done = Standard_True;
 166 }
 167
 168 void Contap_ContAna::Perform (const gp_Cylinder& C,
 169                               const gp_Pnt& Eye)
 170 {
 171   done = Standard_False;
 172
 173   Standard_Real radius = C.Radius();
 174   gp_Lin theaxis(C.Axis());
 175   Standard_Real dist = theaxis.Distance(Eye);
 176   if (dist <= radius) {
 177     nbSol = 0;
 178   }
 179   else {
 180     typL = GeomAbs_Line;
 181     prm = radius*sqrt(1.-radius*radius/(dist*dist));
 182     dir1 = C.Axis().Direction();
 183     dir2 = dir1;
 184     gp_XYZ axeye(theaxis.Normal(Eye).Direction().XYZ()); // orientate the axis to the outside
 185     gp_XYZ normale((theaxis.Direction().Crossed(axeye)).XYZ());
 186 //      normale.Normalize();
 187     pt1.SetXYZ(C.Location().XYZ() + (radius*radius/dist)*axeye);
 188     pt2.SetXYZ(pt1.XYZ() - prm*normale);
 189     pt1.SetXYZ(pt1.XYZ() + prm*normale);
 190     nbSol = 2;
 191   }
 192   done = Standard_True;
 193 }
 194
 195 void Contap_ContAna::Perform (const gp_Cone& C,
 196                               const gp_Dir& D)
 197 {
 198   done = Standard_False;
 199
 200   Standard_Real Tgtalpha = Tan(C.SemiAngle());
 201
 202   Standard_Real Coefcos = D.Dot(C.Position().XDirection());
 203   Standard_Real Coefsin = D.Dot(C.Position().YDirection());
 204   Standard_Real Coefcst = D.Dot(C.Axis().Direction())*Tgtalpha;
 205
 206   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 207   Standard_Real norm2 = Sqrt(norm1);
 208 //  if (Abs(Abs(Coefcst)-norm2) <= Tolpetit) { // tol angulaire 1.e-8
 209 //    typL = GeomAbs_Line;
 210 //    nbSol = 1;
 211 //    pt1 = C.Apex();
 212 //    dir1 = D;
 213 //  }
 214 //  else if (Abs(Coefcst) < norm2) {
 215
 216   if (Abs(Coefcst) < norm2) {
 217     typL = GeomAbs_Line;
 218     nbSol = 2;
 219     pt1 = C.Apex();
 220     pt2 = pt1;
 221     // Necessary to solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 222     // and director vectors of solutions are
 223     // cos(t0) * XDirection + sin(t0) * YDirection + ZDirection/Tgtalpha
 224     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 225     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 226     // and by using trigonometric relations the values of cosinus
 227     // and sinus to the solutions are obtained.
 228
 229     prm = Sqrt(norm1 - Coefcst*Coefcst);
 230     Standard_Real cost0,sint0,cost1,sint1;
 231
 232     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 233     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 234
 235     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 236     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 237
 238     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 239     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 240     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 241     gp_XYZ dirxyz;
 242     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,1./Tgtalpha,Zdir);
 243     dir1.SetXYZ(dirxyz);
 244     pt1.SetXYZ(pt1.XYZ()+dirxyz);
 245     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,1./Tgtalpha,Zdir);
 246     dir2.SetXYZ(dirxyz);
 247     pt2.SetXYZ(pt2.XYZ()+dirxyz);
 248   }
 249   else {
 250     nbSol = 0;
 251   }
 252   done = Standard_True;
 253 }
 254
 255 void Contap_ContAna::Perform (const gp_Cone& C,
 256                               const gp_Dir& D,
 257                               const Standard_Real Angle)
 258 {
 259   done = Standard_False;
 260   nbSol = 0;
 261
 262   Standard_Real Ang = C.SemiAngle();
 263   Standard_Real Cosa = cos(Ang);
 264   Standard_Real Sina = sin(Ang);
 265
 266   Standard_Real Coefcos = D.Dot(C.Position().XDirection());
 267   Standard_Real Coefsin = D.Dot(C.Position().YDirection());
 268
 269   Standard_Real Coefcst1 = cos(M_PI*0.5 + Angle);
 270
 271   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 272   Standard_Real norm2 = Sqrt(norm1);
 273
 274   Standard_Real Coefnz = D.Dot(C.Axis().Direction())*Sina;
 275   Standard_Real Coefcst = (Coefcst1 + Coefnz)/Cosa;
 276
 277   if (Abs(Coefcst) < norm2) {
 278     typL = GeomAbs_Line;
 279     nbSol+= 2;
 280     pt1 = C.Apex();
 281     pt2 = pt1;
 282
 283     // It is requiredto solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 284     // and the director vectors of solutions are
 285     // cos(t0) * XDirection + sin(t0) * YDirection + ZDirection/Tgtalpha
 286     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 287     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 288     // and by using trigonometric relations the values of cosinus
 289     // and sinus to the solutions are obtained.
 290
 291     prm = Sqrt(norm1 - Coefcst*Coefcst);
 292     Standard_Real cost0,sint0,cost1,sint1;
 293
 294     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 295     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 296
 297     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 298     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 299
 300     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 301     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 302     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 303     if (!C.Direct()) {
 304       Zdir.Reverse();
 305     }
 306     gp_XYZ dirxyz;
 307     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,Cosa/Sina,Zdir);
 308     dir1.SetXYZ(dirxyz);
 309     pt1.SetXYZ(pt1.XYZ()+dirxyz);
 310     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,Cosa/Sina,Zdir);
 311     dir2.SetXYZ(dirxyz);
 312     pt2.SetXYZ(pt2.XYZ()+dirxyz);
 313   }
 314
 315   Coefcst = (Coefcst1 - Coefnz)/Cosa;
 316
 317   if (Abs(Coefcst) < norm2) {
 318     typL = GeomAbs_Line;
 319     nbSol+= 2;
 320     pt3 = C.Apex();
 321     pt4 = pt3;
 322
 323     prm = Sqrt(norm1 - Coefcst*Coefcst);
 324     Standard_Real cost0,sint0,cost1,sint1;
 325
 326     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 327     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 328
 329     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 330     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 331
 332     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 333     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 334     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 335     if (!C.Direct()) {
 336       Zdir.Reverse();
 337     }
 338     gp_XYZ dirxyz;
 339     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,-Cosa/Sina,Zdir);
 340     dir3.SetXYZ(dirxyz);
 341     pt3.SetXYZ(pt3.XYZ()+dirxyz);
 342     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,-Cosa/Sina,Zdir);
 343     dir4.SetXYZ(dirxyz);
 344     pt4.SetXYZ(pt4.XYZ()+dirxyz);
 345     if (nbSol == 2) {
 346       pt1 = pt3;
 347       pt2 = pt4;
 348       dir1 = dir3;
 349       dir2 = dir4;
 350     }
 351   }
 352
 353   done = Standard_True;
 354 }
 355
 356 void Contap_ContAna::Perform (const gp_Cone& C,
 357                               const gp_Pnt& Eye)
 358 {
 359   done = Standard_False;
 360
 361   Standard_Real Tgtalpha = Tan(C.SemiAngle());
 362
 363   gp_XYZ apexeye(Eye.XYZ());
 364   apexeye.Subtract(C.Apex().XYZ());
 365
 366   Standard_Real Coefcos = apexeye.Dot(C.Position().XDirection().XYZ());
 367   Standard_Real Coefsin = apexeye.Dot(C.Position().YDirection().XYZ());
 368   Standard_Real Coefcst = apexeye.Dot(C.Axis().Direction().XYZ())*Tgtalpha;
 369
 370   Standard_Real norm1 = Coefcos*Coefcos + Coefsin*Coefsin;
 371   Standard_Real norm2 = Sqrt(Coefcos*Coefcos + Coefsin*Coefsin);
 372 //  if (Abs(Abs(Coefcst)-norm2) <= Tolpetit) { // tol angulaire 1.e-8
 373 //    typL = GeomAbs_Line;
 374 //    nbSol = 1;
 375 //    pt1 = C.Apex();
 376 //    dir1.SetXYZ(apexeye);
 377 //  }
 378 //  else if (Abs(Coefcst) < norm2) {
 379
 380   if (Abs(Coefcst) < norm2) {
 381     typL = GeomAbs_Line;
 382     nbSol = 2;
 383     pt1 = C.Apex();
 384     pt2 = pt1;
 385     // It is required to solve Coefcos*cos(t) + Coefsin*sin(t) = Coefcst
 386     // and the director vectors of solutions are
 387     // cos(t0) * XDirection + sin(t0) * YDirection + ZDirection/Tgtalpha
 388     // By setting cos(phi) = Coefcos/Sqrt(Coefcos**2 + Coefsin**2) and
 389     //           sin(phi) = Coefsin/Sqrt(Coefcos**2 + Coefsin**2)
 390     // and by using trigonometric relations the values of cosinus
 391     // and sinus to the solutions are obtained.
 392
 393     prm = Sqrt(norm1 - Coefcst*Coefcst);
 394     Standard_Real cost0,sint0,cost1,sint1;
 395
 396     cost0 = (Coefcos*Coefcst - Coefsin*prm)/norm1;
 397     cost1 = (Coefcos*Coefcst + Coefsin*prm)/norm1;
 398
 399     sint0 = ( Coefcos*prm + Coefsin*Coefcst)/norm1;
 400     sint1 = (-Coefcos*prm + Coefsin*Coefcst)/norm1;
 401
 402     gp_XYZ Xdir(C.Position().XDirection().XYZ());
 403     gp_XYZ Ydir(C.Position().YDirection().XYZ());
 404     gp_XYZ Zdir(C.Axis().Direction().XYZ());
 405     gp_XYZ dirxyz;
 406     dirxyz.SetLinearForm(cost0,Xdir,sint0,Ydir,1./Tgtalpha,Zdir);
 407     dir1.SetXYZ(dirxyz);
 408     pt1.SetXYZ(pt1.XYZ()+dirxyz);
 409     dirxyz.SetLinearForm(cost1,Xdir,sint1,Ydir,1./Tgtalpha,Zdir);
 410     dir2.SetXYZ(dirxyz);
 411     pt2.SetXYZ(pt2.XYZ()+dirxyz);
 412   }
 413   else {
 414     nbSol = 0;
 415   }
 416   done = Standard_True;
 417 }
 418
 419 gp_Lin Contap_ContAna::Line (const Standard_Integer Index) const
 420 {
 421   if (!done) {StdFail_NotDone::Raise();}
 422   if (typL != GeomAbs_Line || nbSol == 0) {Standard_DomainError::Raise();}
 423   if (Index <=0 || Index > nbSol) {Standard_OutOfRange::Raise();}
 424   switch (Index) {
 425   case 1:
 426     return gp_Lin(pt1,dir1);
 427   case 2:
 428     return gp_Lin(pt2,dir2);
 429   case 3:
 430     return gp_Lin(pt3,dir3);
 431   case 4:
 432     return gp_Lin(pt4,dir4);
 433   }
 434   Standard_OutOfRange::Raise("Program error in Contap_ContAna");
 435   return gp_Lin();
 436 }