[occt.git] / src / IntSurf / IntSurf_Quadric.cxx

// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.


#include <ElCLib.hxx>
#include <ElSLib.hxx>
#include <gp.hxx>
#include <gp_Cone.hxx>
#include <gp_Cylinder.hxx>
#include <gp_Pln.hxx>
#include <gp_Pnt.hxx>
#include <gp_Sphere.hxx>
#include <gp_Torus.hxx>
#include <gp_Vec.hxx>
#include <IntSurf_Quadric.hxx>
#include <StdFail_NotDone.hxx>

// ============================================================
IntSurf_Quadric::IntSurf_Quadric ():typ(GeomAbs_OtherSurface),
   prm1(0.), prm2(0.), prm3(0.), prm4(0.)
{}
// ============================================================
IntSurf_Quadric::IntSurf_Quadric (const gp_Pln& P):
      ax3(P.Position()),typ(GeomAbs_Plane)
{
  ax3direc = ax3.Direct();
  P.Coefficients(prm1,prm2,prm3,prm4); 
}
// ============================================================
IntSurf_Quadric::IntSurf_Quadric (const gp_Cylinder& C):

       ax3(C.Position()),lin(ax3.Axis()),typ(GeomAbs_Cylinder)
{
  prm2=prm3=prm4=0.0;
  ax3direc=ax3.Direct();
  prm1=C.Radius();
}
// ============================================================
IntSurf_Quadric::IntSurf_Quadric (const gp_Sphere& S):

       ax3(S.Position()),lin(ax3.Axis()),typ(GeomAbs_Sphere)
{
  prm2=prm3=prm4=0.0;
  ax3direc = ax3.Direct();
  prm1=S.Radius();
}
// ============================================================
IntSurf_Quadric::IntSurf_Quadric (const gp_Cone& C):

       ax3(C.Position()),typ(GeomAbs_Cone)
{
  ax3direc = ax3.Direct();
  lin.SetPosition(ax3.Axis());
  prm1 = C.RefRadius();
  prm2 = C.SemiAngle();
  prm3 = Cos(prm2);
  prm4 = 0.0;
}
// ============================================================
IntSurf_Quadric::IntSurf_Quadric (const gp_Torus& T):

       ax3(T.Position()),typ(GeomAbs_Torus)
{
  ax3direc = ax3.Direct();
  lin.SetPosition(ax3.Axis());
  prm1 = T.MajorRadius();
  prm2 = T.MinorRadius();
  prm3 = 0.0;
  prm4 = 0.0;
}
// ============================================================
void IntSurf_Quadric::SetValue (const gp_Pln& P)
{
  typ = GeomAbs_Plane;
  ax3 = P.Position();
  ax3direc = ax3.Direct();
  P.Coefficients(prm1,prm2,prm3,prm4); 
}
// ============================================================
void IntSurf_Quadric::SetValue (const gp_Cylinder& C)
{
  typ = GeomAbs_Cylinder;
  ax3 = C.Position();
  ax3direc = ax3.Direct();
  lin.SetPosition(ax3.Axis());
  prm1 = C.Radius();
  prm2=prm3=prm4=0.0;
}
// ============================================================
void IntSurf_Quadric::SetValue (const gp_Sphere& S)
{
  typ = GeomAbs_Sphere;
  ax3 = S.Position();
  ax3direc = ax3.Direct();
  lin.SetPosition(ax3.Axis());
  prm1 = S.Radius();
  prm2=prm3=prm4=0.0;
}
// ============================================================
void IntSurf_Quadric::SetValue (const gp_Cone& C)
{
  typ = GeomAbs_Cone;
  ax3 = C.Position();
  ax3direc = ax3.Direct();
  lin.SetPosition(ax3.Axis());
  prm1 = C.RefRadius();
  prm2 = C.SemiAngle();
  prm3 = Cos(prm2);
  prm4 = 0.0;
}
// ============================================================
void IntSurf_Quadric::SetValue (const gp_Torus& T)
{
  typ = GeomAbs_Torus;
  ax3 = T.Position();
  ax3direc = ax3.Direct();
  lin.SetPosition(ax3.Axis());
  prm1 = T.MajorRadius();
  prm2 = T.MinorRadius();
  prm3 = 0.0;
  prm4 = 0.0;
}
// ============================================================
Standard_Real IntSurf_Quadric::Distance (const gp_Pnt& P) const {
  switch (typ) {
  case GeomAbs_Plane:   // plan
    return prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
  case GeomAbs_Cylinder:   // cylindre
    return (lin.Distance(P) - prm1);
  case GeomAbs_Sphere:   // sphere
    return (lin.Location().Distance(P) - prm1);
  case GeomAbs_Cone:   // cone
    {
      Standard_Real dist = lin.Distance(P);
      Standard_Real U,V;
      ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
      gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
      Standard_Real distp = lin.Distance(Pp);
      dist = (dist-distp)/prm3;
      return(dist);
    }
  case GeomAbs_Torus:   // torus
    {
      gp_Pnt O, Pp, PT;
      //
      O = ax3.Location();
      gp_Vec OZ (ax3.Direction());
      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
      //
      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
      //
      Standard_Real dist = P.Distance(PT) - prm2;
      return dist;
    }
  default:
    {
    }
    break;
  }
  return(0.0);
}
// ============================================================
gp_Vec IntSurf_Quadric::Gradient (const gp_Pnt& P) const {
  gp_Vec grad;
  switch (typ) {
  case GeomAbs_Plane:   // plan
    grad.SetCoord(prm1,prm2,prm3);
    break;
  case GeomAbs_Cylinder:   // cylindre
    {
      gp_XYZ PP(lin.Location().XYZ());
      PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
      grad.SetXYZ(P.XYZ()-PP);
      Standard_Real N = grad.Magnitude();
      if(N>1e-14) { grad.Divide(N); } 
      else { grad.SetCoord(0.0,0.0,0.0); } 
    }
    break;
  case GeomAbs_Sphere:   // sphere
    {
      gp_XYZ PP(P.XYZ());
      grad.SetXYZ((PP-lin.Location().XYZ()));
      Standard_Real N = grad.Magnitude();
      if(N>1e-14) { grad.Divide(N); } 
      else { grad.SetCoord(0.0,0.0,0.0); }
    }
    break;
  case GeomAbs_Cone:   // cone
    {
      Standard_Real U,V;
      ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
      gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
      gp_Vec D1u,D1v;
      ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
      grad=D1u.Crossed(D1v);
      if(ax3direc==Standard_False) { 
        grad.Reverse();
      }
      grad.Normalize();
    }
    break;
  case GeomAbs_Torus:   // torus
    {
      gp_Pnt O, Pp, PT;
      //
      O = ax3.Location();
      gp_Vec OZ (ax3.Direction());
      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
      //
      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
      //
      grad.SetXYZ(P.XYZ() - PT.XYZ());
      Standard_Real N = grad.Magnitude();
      if(N>1e-14) { grad.Divide(N); } 
      else { grad.SetCoord(0., 0., 0.); }
    }
    break;
  default:
    {}
    break;
  }
  return grad;
}
// ============================================================
void IntSurf_Quadric::ValAndGrad (const gp_Pnt& P,
                                  Standard_Real& Dist,
                                  gp_Vec& Grad) const
{

  switch (typ) {
  case GeomAbs_Plane:   
    {
      Dist = prm1*P.X() + prm2*P.Y() + prm3*P.Z() + prm4;
      Grad.SetCoord(prm1,prm2,prm3);
    }
    break;
  case GeomAbs_Cylinder:   
    {
      Dist = lin.Distance(P) - prm1;
      gp_XYZ PP(lin.Location().XYZ());
      PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
      Grad.SetXYZ((P.XYZ()-PP));
      Standard_Real N = Grad.Magnitude();
      if(N>1e-14) { Grad.Divide(N); } 
      else { Grad.SetCoord(0.0,0.0,0.0); }
    }
    break;
  case GeomAbs_Sphere:  
    {
      Dist = lin.Location().Distance(P) - prm1;
      gp_XYZ PP(P.XYZ());
      Grad.SetXYZ((PP-lin.Location().XYZ()));
      Standard_Real N = Grad.Magnitude();
      if(N>1e-14) { Grad.Divide(N); } 
      else { Grad.SetCoord(0.0,0.0,0.0); }
    }
    break;
  case GeomAbs_Cone:  
    {
      Standard_Real dist = lin.Distance(P);
      Standard_Real U,V;
      gp_Vec D1u,D1v;
      gp_Pnt Pp;
      ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
      ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
      Standard_Real distp = lin.Distance(Pp);
      dist = (dist-distp)/prm3;
      Dist = dist;
      Grad=D1u.Crossed(D1v);
      if(ax3direc==Standard_False) { 
        Grad.Reverse();
      }
      //-- lbr le 7 mars 96 
      //-- Si le gardient est nul, on est sur l axe 
      //-- et dans ce cas dist vaut 0 
      //-- On peut donc renvoyer une valeur quelconque. 
      if(   Grad.X() > 1e-13 
         || Grad.Y() > 1e-13 
         || Grad.Z() > 1e-13) { 
        Grad.Normalize();
      }
    }
    break;
  case GeomAbs_Torus:   
    {
      gp_Pnt O, Pp, PT;
      //
      O = ax3.Location();
      gp_Vec OZ (ax3.Direction());
      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
      //
      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
      //
      Dist = P.Distance(PT) - prm2;
      //
      Grad.SetXYZ(P.XYZ()-PT.XYZ());
      Standard_Real N = Grad.Magnitude();
      if(N>1e-14) { Grad.Divide(N); } 
      else { Grad.SetCoord(0., 0., 0.); }
    }
    break;
  default:
    {}
    break;
  }
}
// ============================================================
gp_Pnt IntSurf_Quadric::Value(const Standard_Real U,
                              const Standard_Real V) const
{
  switch (typ) {

  case GeomAbs_Plane:  
    return ElSLib::PlaneValue(U,V,ax3);
  case GeomAbs_Cylinder:
    return ElSLib::CylinderValue(U,V,ax3,prm1);
  case GeomAbs_Sphere:
    return ElSLib::SphereValue(U,V,ax3,prm1);
  case GeomAbs_Cone:  
    return ElSLib::ConeValue(U,V,ax3,prm1,prm2);
  case GeomAbs_Torus:  
    return ElSLib::TorusValue(U,V,ax3,prm1,prm2);
  default:
    {
      gp_Pnt p(0,0,0);
      return(p);
    }
    //break;
  }
// pop : pour NT
//  return gp_Pnt(0,0,0);
}
// ============================================================
void IntSurf_Quadric::D1(const Standard_Real U,
                         const Standard_Real V,
                         gp_Pnt& P,
                         gp_Vec& D1U,
                         gp_Vec& D1V) const
{
  switch (typ) {
  case GeomAbs_Plane:  
    ElSLib::PlaneD1(U,V,ax3,P,D1U,D1V);
    break;
  case GeomAbs_Cylinder:
    ElSLib::CylinderD1(U,V,ax3,prm1,P,D1U,D1V);
    break;
  case GeomAbs_Sphere: 
    ElSLib::SphereD1(U,V,ax3,prm1,P,D1U,D1V);
    break;
  case GeomAbs_Cone:  
    ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
    break;
  case GeomAbs_Torus:  
    ElSLib::TorusD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
    break;
  default:
    {
    }
    break;
  }
}
// ============================================================
gp_Vec IntSurf_Quadric::DN(const Standard_Real U,
                           const Standard_Real V,
                           const Standard_Integer Nu,
                           const Standard_Integer Nv) const
{
  switch (typ) {
  case GeomAbs_Plane: 
    return ElSLib::PlaneDN(U,V,ax3,Nu,Nv);
  case GeomAbs_Cylinder:
    return ElSLib::CylinderDN(U,V,ax3,prm1,Nu,Nv);
  case GeomAbs_Sphere:
    return ElSLib::SphereDN(U,V,ax3,prm1,Nu,Nv);
  case GeomAbs_Cone:  
    return ElSLib::ConeDN(U,V,ax3,prm1,prm2,Nu,Nv);
  case GeomAbs_Torus:  
    return ElSLib::TorusDN(U,V,ax3,prm1,prm2,Nu,Nv);
  default:
    {
      gp_Vec v(0,0,0);
      return(v);
    }
    //break;
  }
// pop : pour NT
//  return gp_Vec(0,0,0);
}
// ============================================================
gp_Vec IntSurf_Quadric::Normale(const Standard_Real U,
                                const Standard_Real V) const
{
  switch (typ) {
  case GeomAbs_Plane:  
    if(ax3direc) 
      return ax3.Direction();
    else 
      return ax3.Direction().Reversed();
  case GeomAbs_Cylinder:  
    return Normale(Value(U,V));
  case GeomAbs_Sphere:  
    return Normale(Value(U,V));
  case GeomAbs_Cone:   
    {
      gp_Pnt P;
      gp_Vec D1u,D1v;
      ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1u,D1v);
      if(D1u.Magnitude()<0.0000001) { 
        gp_Vec Vn(0.0,0.0,0.0); 
        return(Vn);
      }
      return(D1u.Crossed(D1v));
    }
  case GeomAbs_Torus:
    return Normale(Value(U,V));
  default:     
    {
      gp_Vec v(0,0,0);
      return(v);
    }     
  //  break;
  }
// pop : pour NT
//  return gp_Vec(0,0,0);
}
// ============================================================
gp_Vec IntSurf_Quadric::Normale (const gp_Pnt& P) const
{
  switch (typ) {
  case GeomAbs_Plane:   
    if(ax3direc) 
      return ax3.Direction();
    else 
      return ax3.Direction().Reversed();
  case GeomAbs_Cylinder:   
    {
      if(ax3direc) { 
        return lin.Normal(P).Direction();
      }
      else { 
        gp_Dir D(lin.Normal(P).Direction());
        D.Reverse();
        return(D);
      }
    }
  case GeomAbs_Sphere:   
    {
      if(ax3direc) { 
        gp_Vec ax3P(ax3.Location(),P);
        return gp_Dir(ax3P);
      }
      else { 
        gp_Vec Pax3(P,ax3.Location());
        return gp_Dir(Pax3);
      }
    }
  case GeomAbs_Cone:   
    {
      Standard_Real U,V;
      ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
      return Normale(U,V);
    }
  case GeomAbs_Torus:
    {
      gp_Pnt O, Pp, PT;
      //
      O = ax3.Location();
      gp_Vec OZ (ax3.Direction());
      Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
      //
      gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ? 
        ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
      PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
      if (PT.SquareDistance(P) < 1e-14) {
        return gp_Dir(OZ);
      }
      gp_Dir aD(ax3direc ? gp_Vec(PT, P) : gp_Vec(P, PT));
      return aD;
    }
  default:      
    {
      gp_Vec v(0,0,0);
      return(v);
    } //    break;
  }
}
// ============================================================
void IntSurf_Quadric::Parameters (const gp_Pnt& P,
                                  Standard_Real& U,
                                  Standard_Real& V) const
{
  switch (typ) {
  case GeomAbs_Plane:   
    ElSLib::PlaneParameters(ax3,P,U,V);
    break;
  case GeomAbs_Cylinder: 
    ElSLib::CylinderParameters(ax3,prm1,P,U,V);
    break;
  case GeomAbs_Sphere:   
    ElSLib::SphereParameters(ax3,prm1,P,U,V);
    break;
  case GeomAbs_Cone: 
    ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
    break;
  case GeomAbs_Torus: 
    ElSLib::TorusParameters(ax3,prm1,prm2,P,U,V);
    break;
  default:
    break;
  }
}
// ============================================================
Commit	Line	Data
b311480e	1	// Copyright (c) 1995-1999 Matra Datavision
973c2be1	2	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	3	//
973c2be1	4	// This file is part of Open CASCADE Technology software library.
b311480e	5	//
d5f74e42	6	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	7	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	8	// by the Free Software Foundation, with special exception defined in the file
	9	// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
	10	// distribution for complete text of the license and disclaimer of any warranty.
b311480e	11	//
973c2be1	12	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	13	// commercial license or contractual agreement.
b311480e	14
7fd59977	15
	16	#include <ElCLib.hxx>
	17	#include <ElSLib.hxx>
	18	#include <gp.hxx>
42cf5bc1	19	#include <gp_Cone.hxx>
	20	#include <gp_Cylinder.hxx>
	21	#include <gp_Pln.hxx>
	22	#include <gp_Pnt.hxx>
	23	#include <gp_Sphere.hxx>
	24	#include <gp_Torus.hxx>
	25	#include <gp_Vec.hxx>
	26	#include <IntSurf_Quadric.hxx>
	27	#include <StdFail_NotDone.hxx>
7fd59977	28
7fd59977	29	// ============================================================
	30	IntSurf_Quadric::IntSurf_Quadric ():typ(GeomAbs_OtherSurface),
	31	prm1(0.), prm2(0.), prm3(0.), prm4(0.)
	32	{}
	33	// ============================================================
	34	IntSurf_Quadric::IntSurf_Quadric (const gp_Pln& P):
	35	ax3(P.Position()),typ(GeomAbs_Plane)
	36	{
	37	ax3direc = ax3.Direct();
	38	P.Coefficients(prm1,prm2,prm3,prm4);
	39	}
	40	// ============================================================
	41	IntSurf_Quadric::IntSurf_Quadric (const gp_Cylinder& C):
	42
	43	ax3(C.Position()),lin(ax3.Axis()),typ(GeomAbs_Cylinder)
	44	{
	45	prm2=prm3=prm4=0.0;
	46	ax3direc=ax3.Direct();
	47	prm1=C.Radius();
	48	}
	49	// ============================================================
	50	IntSurf_Quadric::IntSurf_Quadric (const gp_Sphere& S):
	51
	52	ax3(S.Position()),lin(ax3.Axis()),typ(GeomAbs_Sphere)
	53	{
	54	prm2=prm3=prm4=0.0;
	55	ax3direc = ax3.Direct();
	56	prm1=S.Radius();
	57	}
	58	// ============================================================
	59	IntSurf_Quadric::IntSurf_Quadric (const gp_Cone& C):
	60
	61	ax3(C.Position()),typ(GeomAbs_Cone)
	62	{
	63	ax3direc = ax3.Direct();
	64	lin.SetPosition(ax3.Axis());
	65	prm1 = C.RefRadius();
	66	prm2 = C.SemiAngle();
	67	prm3 = Cos(prm2);
	68	prm4 = 0.0;
	69	}
	70	// ============================================================
7eed5d29	71	IntSurf_Quadric::IntSurf_Quadric (const gp_Torus& T):
	72
	73	ax3(T.Position()),typ(GeomAbs_Torus)
	74	{
	75	ax3direc = ax3.Direct();
	76	lin.SetPosition(ax3.Axis());
	77	prm1 = T.MajorRadius();
	78	prm2 = T.MinorRadius();
	79	prm3 = 0.0;
	80	prm4 = 0.0;
	81	}
	82	// ============================================================
7fd59977	83	void IntSurf_Quadric::SetValue (const gp_Pln& P)
	84	{
	85	typ = GeomAbs_Plane;
	86	ax3 = P.Position();
	87	ax3direc = ax3.Direct();
	88	P.Coefficients(prm1,prm2,prm3,prm4);
	89	}
	90	// ============================================================
	91	void IntSurf_Quadric::SetValue (const gp_Cylinder& C)
	92	{
	93	typ = GeomAbs_Cylinder;
	94	ax3 = C.Position();
	95	ax3direc = ax3.Direct();
	96	lin.SetPosition(ax3.Axis());
	97	prm1 = C.Radius();
	98	prm2=prm3=prm4=0.0;
	99	}
	100	// ============================================================
	101	void IntSurf_Quadric::SetValue (const gp_Sphere& S)
	102	{
	103	typ = GeomAbs_Sphere;
	104	ax3 = S.Position();
	105	ax3direc = ax3.Direct();
	106	lin.SetPosition(ax3.Axis());
	107	prm1 = S.Radius();
	108	prm2=prm3=prm4=0.0;
	109	}
	110	// ============================================================
	111	void IntSurf_Quadric::SetValue (const gp_Cone& C)
	112	{
	113	typ = GeomAbs_Cone;
	114	ax3 = C.Position();
	115	ax3direc = ax3.Direct();
	116	lin.SetPosition(ax3.Axis());
	117	prm1 = C.RefRadius();
	118	prm2 = C.SemiAngle();
	119	prm3 = Cos(prm2);
	120	prm4 = 0.0;
	121	}
	122	// ============================================================
7eed5d29	123	void IntSurf_Quadric::SetValue (const gp_Torus& T)
	124	{
	125	typ = GeomAbs_Torus;
	126	ax3 = T.Position();
	127	ax3direc = ax3.Direct();
	128	lin.SetPosition(ax3.Axis());
	129	prm1 = T.MajorRadius();
	130	prm2 = T.MinorRadius();
	131	prm3 = 0.0;
	132	prm4 = 0.0;
	133	}
	134	// ============================================================
7fd59977	135	Standard_Real IntSurf_Quadric::Distance (const gp_Pnt& P) const {
	136	switch (typ) {
	137	case GeomAbs_Plane: // plan
	138	return prm1P.X() + prm2P.Y() + prm3*P.Z() + prm4;
	139	case GeomAbs_Cylinder: // cylindre
	140	return (lin.Distance(P) - prm1);
	141	case GeomAbs_Sphere: // sphere
	142	return (lin.Location().Distance(P) - prm1);
	143	case GeomAbs_Cone: // cone
	144	{
	145	Standard_Real dist = lin.Distance(P);
	146	Standard_Real U,V;
	147	ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
	148	gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
	149	Standard_Real distp = lin.Distance(Pp);
	150	dist = (dist-distp)/prm3;
	151	return(dist);
	152	}
7eed5d29	153	case GeomAbs_Torus: // torus
	154	{
	155	gp_Pnt O, Pp, PT;
	156	//
	157	O = ax3.Location();
	158	gp_Vec OZ (ax3.Direction());
	159	Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
	160	//
	161	gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
	162	ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
	163	PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
	164	//
	165	Standard_Real dist = P.Distance(PT) - prm2;
	166	return dist;
	167	}
7fd59977	168	default:
	169	{
	170	}
	171	break;
	172	}
	173	return(0.0);
	174	}
	175	// ============================================================
	176	gp_Vec IntSurf_Quadric::Gradient (const gp_Pnt& P) const {
	177	gp_Vec grad;
	178	switch (typ) {
	179	case GeomAbs_Plane: // plan
	180	grad.SetCoord(prm1,prm2,prm3);
	181	break;
	182	case GeomAbs_Cylinder: // cylindre
	183	{
	184	gp_XYZ PP(lin.Location().XYZ());
	185	PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
	186	grad.SetXYZ(P.XYZ()-PP);
	187	Standard_Real N = grad.Magnitude();
	188	if(N>1e-14) { grad.Divide(N); }
	189	else { grad.SetCoord(0.0,0.0,0.0); }
	190	}
	191	break;
	192	case GeomAbs_Sphere: // sphere
	193	{
	194	gp_XYZ PP(P.XYZ());
	195	grad.SetXYZ((PP-lin.Location().XYZ()));
	196	Standard_Real N = grad.Magnitude();
	197	if(N>1e-14) { grad.Divide(N); }
	198	else { grad.SetCoord(0.0,0.0,0.0); }
	199	}
	200	break;
	201	case GeomAbs_Cone: // cone
	202	{
	203	Standard_Real U,V;
	204	ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
	205	gp_Pnt Pp = ElSLib::ConeValue(U,V,ax3,prm1,prm2);
	206	gp_Vec D1u,D1v;
	207	ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
	208	grad=D1u.Crossed(D1v);
	209	if(ax3direc==Standard_False) {
7eed5d29	210	grad.Reverse();
7fd59977	211	}
	212	grad.Normalize();
	213	}
	214	break;
7eed5d29	215	case GeomAbs_Torus: // torus
	216	{
	217	gp_Pnt O, Pp, PT;
	218	//
	219	O = ax3.Location();
	220	gp_Vec OZ (ax3.Direction());
	221	Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
	222	//
	223	gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
	224	ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
	225	PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
	226	//
	227	grad.SetXYZ(P.XYZ() - PT.XYZ());
	228	Standard_Real N = grad.Magnitude();
	229	if(N>1e-14) { grad.Divide(N); }
	230	else { grad.SetCoord(0., 0., 0.); }
	231	}
	232	break;
7fd59977	233	default:
	234	{}
	235	break;
	236	}
	237	return grad;
	238	}
	239	// ============================================================
	240	void IntSurf_Quadric::ValAndGrad (const gp_Pnt& P,
7eed5d29	241	Standard_Real& Dist,
7eed5d29	242	gp_Vec& Grad) const
7fd59977	243	{
	244
	245	switch (typ) {
	246	case GeomAbs_Plane:
	247	{
	248	Dist = prm1P.X() + prm2P.Y() + prm3*P.Z() + prm4;
	249	Grad.SetCoord(prm1,prm2,prm3);
	250	}
	251	break;
	252	case GeomAbs_Cylinder:
	253	{
	254	Dist = lin.Distance(P) - prm1;
	255	gp_XYZ PP(lin.Location().XYZ());
	256	PP.Add(ElCLib::Parameter(lin,P)*lin.Direction().XYZ());
	257	Grad.SetXYZ((P.XYZ()-PP));
	258	Standard_Real N = Grad.Magnitude();
	259	if(N>1e-14) { Grad.Divide(N); }
	260	else { Grad.SetCoord(0.0,0.0,0.0); }
	261	}
	262	break;
	263	case GeomAbs_Sphere:
	264	{
	265	Dist = lin.Location().Distance(P) - prm1;
	266	gp_XYZ PP(P.XYZ());
	267	Grad.SetXYZ((PP-lin.Location().XYZ()));
	268	Standard_Real N = Grad.Magnitude();
	269	if(N>1e-14) { Grad.Divide(N); }
	270	else { Grad.SetCoord(0.0,0.0,0.0); }
	271	}
	272	break;
	273	case GeomAbs_Cone:
	274	{
	275	Standard_Real dist = lin.Distance(P);
	276	Standard_Real U,V;
	277	gp_Vec D1u,D1v;
	278	gp_Pnt Pp;
	279	ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
	280	ElSLib::ConeD1(U,V,ax3,prm1,prm2,Pp,D1u,D1v);
	281	Standard_Real distp = lin.Distance(Pp);
	282	dist = (dist-distp)/prm3;
	283	Dist = dist;
	284	Grad=D1u.Crossed(D1v);
	285	if(ax3direc==Standard_False) {
7eed5d29	286	Grad.Reverse();
7fd59977	287	}
	288	//-- lbr le 7 mars 96
	289	//-- Si le gardient est nul, on est sur l axe
	290	//-- et dans ce cas dist vaut 0
	291	//-- On peut donc renvoyer une valeur quelconque.
	292	if( Grad.X() > 1e-13
7eed5d29	293	\|\| Grad.Y() > 1e-13
	294	\|\| Grad.Z() > 1e-13) {
	295	Grad.Normalize();
7fd59977	296	}
	297	}
	298	break;
7eed5d29	299	case GeomAbs_Torus:
	300	{
	301	gp_Pnt O, Pp, PT;
	302	//
	303	O = ax3.Location();
	304	gp_Vec OZ (ax3.Direction());
	305	Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
	306	//
	307	gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
	308	ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
	309	PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
	310	//
	311	Dist = P.Distance(PT) - prm2;
	312	//
	313	Grad.SetXYZ(P.XYZ()-PT.XYZ());
	314	Standard_Real N = Grad.Magnitude();
	315	if(N>1e-14) { Grad.Divide(N); }
	316	else { Grad.SetCoord(0., 0., 0.); }
	317	}
	318	break;
7fd59977	319	default:
	320	{}
	321	break;
	322	}
	323	}
	324	// ============================================================
	325	gp_Pnt IntSurf_Quadric::Value(const Standard_Real U,
7eed5d29	326	const Standard_Real V) const
7fd59977	327	{
	328	switch (typ) {
	329
	330	case GeomAbs_Plane:
	331	return ElSLib::PlaneValue(U,V,ax3);
	332	case GeomAbs_Cylinder:
	333	return ElSLib::CylinderValue(U,V,ax3,prm1);
	334	case GeomAbs_Sphere:
	335	return ElSLib::SphereValue(U,V,ax3,prm1);
	336	case GeomAbs_Cone:
	337	return ElSLib::ConeValue(U,V,ax3,prm1,prm2);
7eed5d29	338	case GeomAbs_Torus:
7eed5d29	339	return ElSLib::TorusValue(U,V,ax3,prm1,prm2);
7fd59977	340	default:
	341	{
	342	gp_Pnt p(0,0,0);
	343	return(p);
	344	}
d3f26155	345	//break;
7fd59977	346	}
7fd59977	347	// pop : pour NT
d3f26155	348	// return gp_Pnt(0,0,0);
7fd59977	349	}
	350	// ============================================================
	351	void IntSurf_Quadric::D1(const Standard_Real U,
7eed5d29	352	const Standard_Real V,
	353	gp_Pnt& P,
	354	gp_Vec& D1U,
	355	gp_Vec& D1V) const
7fd59977	356	{
	357	switch (typ) {
	358	case GeomAbs_Plane:
	359	ElSLib::PlaneD1(U,V,ax3,P,D1U,D1V);
	360	break;
	361	case GeomAbs_Cylinder:
	362	ElSLib::CylinderD1(U,V,ax3,prm1,P,D1U,D1V);
	363	break;
	364	case GeomAbs_Sphere:
	365	ElSLib::SphereD1(U,V,ax3,prm1,P,D1U,D1V);
	366	break;
	367	case GeomAbs_Cone:
	368	ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
	369	break;
7eed5d29	370	case GeomAbs_Torus:
	371	ElSLib::TorusD1(U,V,ax3,prm1,prm2,P,D1U,D1V);
	372	break;
7fd59977	373	default:
	374	{
	375	}
	376	break;
	377	}
	378	}
	379	// ============================================================
	380	gp_Vec IntSurf_Quadric::DN(const Standard_Real U,
7eed5d29	381	const Standard_Real V,
	382	const Standard_Integer Nu,
	383	const Standard_Integer Nv) const
7fd59977	384	{
	385	switch (typ) {
	386	case GeomAbs_Plane:
	387	return ElSLib::PlaneDN(U,V,ax3,Nu,Nv);
	388	case GeomAbs_Cylinder:
	389	return ElSLib::CylinderDN(U,V,ax3,prm1,Nu,Nv);
	390	case GeomAbs_Sphere:
	391	return ElSLib::SphereDN(U,V,ax3,prm1,Nu,Nv);
	392	case GeomAbs_Cone:
	393	return ElSLib::ConeDN(U,V,ax3,prm1,prm2,Nu,Nv);
7eed5d29	394	case GeomAbs_Torus:
7eed5d29	395	return ElSLib::TorusDN(U,V,ax3,prm1,prm2,Nu,Nv);
7fd59977	396	default:
	397	{
	398	gp_Vec v(0,0,0);
	399	return(v);
	400	}
d3f26155	401	//break;
7fd59977	402	}
7fd59977	403	// pop : pour NT
d3f26155	404	// return gp_Vec(0,0,0);
7fd59977	405	}
	406	// ============================================================
	407	gp_Vec IntSurf_Quadric::Normale(const Standard_Real U,
7eed5d29	408	const Standard_Real V) const
7fd59977	409	{
	410	switch (typ) {
	411	case GeomAbs_Plane:
	412	if(ax3direc)
	413	return ax3.Direction();
	414	else
	415	return ax3.Direction().Reversed();
	416	case GeomAbs_Cylinder:
	417	return Normale(Value(U,V));
	418	case GeomAbs_Sphere:
	419	return Normale(Value(U,V));
	420	case GeomAbs_Cone:
	421	{
	422	gp_Pnt P;
	423	gp_Vec D1u,D1v;
	424	ElSLib::ConeD1(U,V,ax3,prm1,prm2,P,D1u,D1v);
	425	if(D1u.Magnitude()<0.0000001) {
7eed5d29	426	gp_Vec Vn(0.0,0.0,0.0);
7eed5d29	427	return(Vn);
7fd59977	428	}
	429	return(D1u.Crossed(D1v));
	430	}
7eed5d29	431	case GeomAbs_Torus:
7eed5d29	432	return Normale(Value(U,V));
7fd59977	433	default:
	434	{
	435	gp_Vec v(0,0,0);
	436	return(v);
	437	}
d3f26155	438	// break;
7fd59977	439	}
7fd59977	440	// pop : pour NT
d3f26155	441	// return gp_Vec(0,0,0);
7fd59977	442	}
	443	// ============================================================
	444	gp_Vec IntSurf_Quadric::Normale (const gp_Pnt& P) const
	445	{
	446	switch (typ) {
	447	case GeomAbs_Plane:
	448	if(ax3direc)
	449	return ax3.Direction();
	450	else
	451	return ax3.Direction().Reversed();
	452	case GeomAbs_Cylinder:
	453	{
	454	if(ax3direc) {
7eed5d29	455	return lin.Normal(P).Direction();
7fd59977	456	}
7fd59977	457	else {
7eed5d29	458	gp_Dir D(lin.Normal(P).Direction());
	459	D.Reverse();
	460	return(D);
7fd59977	461	}
	462	}
	463	case GeomAbs_Sphere:
	464	{
	465	if(ax3direc) {
7eed5d29	466	gp_Vec ax3P(ax3.Location(),P);
7eed5d29	467	return gp_Dir(ax3P);
7fd59977	468	}
7fd59977	469	else {
7eed5d29	470	gp_Vec Pax3(P,ax3.Location());
7eed5d29	471	return gp_Dir(Pax3);
7fd59977	472	}
	473	}
	474	case GeomAbs_Cone:
	475	{
	476	Standard_Real U,V;
	477	ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
7eed5d29	478	return Normale(U,V);
	479	}
	480	case GeomAbs_Torus:
	481	{
	482	gp_Pnt O, Pp, PT;
	483	//
	484	O = ax3.Location();
	485	gp_Vec OZ (ax3.Direction());
	486	Pp = P.Translated(OZ.Multiplied(-(gp_Vec(O,P).Dot(ax3.Direction()))));
	487	//
	488	gp_Dir DOPp = (O.SquareDistance(Pp) < 1e-14) ?
	489	ax3.XDirection() : gp_Dir(gp_Vec(O, Pp));
	490	PT.SetXYZ(O.XYZ() + DOPp.XYZ()*prm1);
	491	if (PT.SquareDistance(P) < 1e-14) {
	492	return gp_Dir(OZ);
	493	}
	494	gp_Dir aD(ax3direc ? gp_Vec(PT, P) : gp_Vec(P, PT));
	495	return aD;
7fd59977	496	}
	497	default:
	498	{
	499	gp_Vec v(0,0,0);
	500	return(v);
d3f26155	501	} // break;
7fd59977	502	}
7fd59977	503	}
	504	// ============================================================
	505	void IntSurf_Quadric::Parameters (const gp_Pnt& P,
7eed5d29	506	Standard_Real& U,
7eed5d29	507	Standard_Real& V) const
7fd59977	508	{
	509	switch (typ) {
	510	case GeomAbs_Plane:
	511	ElSLib::PlaneParameters(ax3,P,U,V);
	512	break;
	513	case GeomAbs_Cylinder:
	514	ElSLib::CylinderParameters(ax3,prm1,P,U,V);
	515	break;
	516	case GeomAbs_Sphere:
	517	ElSLib::SphereParameters(ax3,prm1,P,U,V);
	518	break;
	519	case GeomAbs_Cone:
7eed5d29	520	ElSLib::ConeParameters(ax3,prm1,prm2,P,U,V);
	521	break;
	522	case GeomAbs_Torus:
	523	ElSLib::TorusParameters(ax3,prm1,prm2,P,U,V);
7fd59977	524	break;
7fd59977	525	default:
7fd59977	526	break;
	527	}
	528	}
	529	// ============================================================
	530
	531
	532
	533
	534