include/IRT2/ToricSurface.h

0001 #pragma once
0002
0003 #include "ParametricSurface.h"
0004
0005 //
0006 // Parameterization: R - major radius, r - donut tube radius;
0007 //
0008 // NB: this class is actually disfunctional (a crossing point with a straight line
0009 // in 3D ends up with a quartic equation);
0010 //
0011
0012 namespace IRT2 {
0013
0014 class ToricSurface: public ParametricSurface {
0015  public:
0016  ToricSurface(): m_Concave(true), m_R(0.0), m_r(0.0), m_Alfa(0.0) {};
0017  ToricSurface(const TVector3 &x0, const TVector3 &nz, double R, double r,
0018           double umin = 0.0, double umax = 2*M_PI,
0019           double vmin = 0.0, double vmax = 2*M_PI):
0020   ParametricSurface(x0, umin, umax, vmin, vmax), m_Concave(true), m_R(R), m_r(r) {
0021     if (!m_Nz.x() && !m_Nz.y()) {
0022       m_Alfa = 0.0;
0023       m_Nr = TVector3(0,0,1);
0024     } else {
0025       auto axis = TVector3(0,0,1).Cross(m_Nz);
0026       m_Alfa = axis.Mag();
0027       m_Nr = axis.Unit();
0028     } //if
0029
0030   };
0031   ~ToricSurface() {};
0032
0033   // FIXME: no range check?;
0034   TVector3 GetSpacePoint(double theta, double phi) const {
0035     // Here theta is rotation angle around the tube, phi is rotation around the torus axis;
0036     double x = (m_R + m_r*cos(theta))*cos(phi), y = (m_R + m_r*cos(theta))*sin(phi), z = m_r*sin(theta);
0037     TVector3 local(x, y, z);
0038     // FIXME: sign check;
0039     local.Rotate(m_Alfa, m_Nr);
0040
0041     return GetCenter() + local;
0042   };
0043
0044   TVector3 GetNormal(const TVector3 &xx) const {
0045     //assert(0);
0046     auto dx = xx - GetCenter();
0047     // Rotate to align torus axis with (0,0,1);
0048     dx.Rotate(-m_Alfa, m_Nr);
0049     TVector3 pc(m_R, 0.0, 0.0);
0050     pc.Rotate(dx.Phi(), TVector3(0,0,1));
0051     dx -= pc;
0052
0053     // Rotate back;
0054     dx.Rotate( m_Alfa, m_Nr);
0055
0056     return (m_Concave ? -1.0 : 1.0)*dx.Unit();
0057   };
0058   //double GetRadius( void ) const { return m_Radius; };
0059
0060   bool GetCrossing(const TVector3 &x0, const TVector3 &n0, TVector3 *crs, bool check_normal = true) const;
0061   double GetDistance(const TVector3 &xx) const;
0062
0063   void SetConvex( void ) { m_Concave = false; };
0064   ParametricSurface *_Clone(double angle, const TVector3 &axis) const {
0065     auto copy = new ToricSurface(*this);
0066
0067     copy->m_Center.Rotate(angle, axis);
0068
0069     copy->m_Nz.Rotate(angle, axis);
0070     copy->m_Nr.Rotate(angle, axis);
0071
0072     return copy;
0073   };
0074
0075   // Theta;
0076   double GetU(const TVector3 &xx) const {
0077     auto dx = xx - GetCenter();
0078     // FIXME: sign check;
0079     dx.Rotate(-m_Alfa, m_Nr);
0080     dx.Rotate(-dx.Phi(), TVector3(0,0,1));
0081     TVector3 pc(m_R, 0.0, 0.0);
0082     dx -= pc;
0083     TVector2 qc(dx.X(), dx.Z());
0084
0085     // Is this phi definition the same as in GetSpacePoint()?;
0086     return qc.Phi();
0087   };
0088   // Phi;
0089   double GetV(const TVector3 &xx) const {
0090     auto dx = xx - GetCenter();
0091     // FIXME: sign check;
0092     dx.Rotate(-m_Alfa, m_Nr);
0093
0094     return dx.Phi();
0095   };
0096
0097  private:
0098   bool m_Concave;
0099   TVector3 m_Nz, m_Nr;
0100   double m_R, m_r, m_Alfa;
0101
0102 #ifndef DISABLE_ROOT_IO
0103   ClassDef(ToricSurface, 1);
0104 #endif
0105 };
0106
0107 } // namespace IRT2