Cone.cpp

Go to the documentation of this file.

// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
//   NScD Oak Ridge National Laboratory, European Spallation Source,
//   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include <algorithm>
#include <cmath>
#include <complex>
#include <fstream>
#include <list>
#include <map>
#include <sstream>
#include <stack>
#include <vector>
 
#include "MantidGeometry/Surfaces/BaseVisit.h"
#include "MantidGeometry/Surfaces/Cone.h"
#include "MantidGeometry/Surfaces/Line.h"
#include "MantidGeometry/Surfaces/Quadratic.h"
#include "MantidGeometry/Surfaces/Surface.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/Strings.h"
#include "MantidKernel/Tolerance.h"
#include "MantidKernel/V3D.h"
 
#ifdef ENABLE_OPENCASCADE
// Opencascade defines _USE_MATH_DEFINES without checking whether it is already
// used.
// Undefine it here before we include the headers to avoid a warning
#ifdef _MSC_VER
#undef _USE_MATH_DEFINES
#ifdef M_SQRT1_2
#undef M_SQRT1_2
#endif
#endif
 
#include "MantidKernel/WarningSuppressions.h"
GNU_DIAG_OFF("conversion")
GNU_DIAG_OFF("cast-qual")
#include <BRepPrimAPI_MakeCone.hxx>
GNU_DIAG_ON("conversion")
GNU_DIAG_ON("cast-qual")
#endif
 
namespace Mantid::Geometry {
using Kernel::Tolerance;
using Kernel::V3D;
 
Cone::Cone()
    : Quadratic(), Centre(), Normal(1, 0, 0), alpha(0.0), cangle(1.0)
{
  setBaseEqn();
}
 
Cone *Cone::doClone() const
{
  return new Cone(*this);
}
 
std::unique_ptr<Cone> Cone::clone() const { return std::unique_ptr<Cone>(doClone()); }
 
int Cone::setSurface(const std::string &Pstr)
{
  std::string Line = Pstr;
  std::string item;
  if (!Mantid::Kernel::Strings::section(Line, item) || tolower(item[0]) != 'k' || item.length() < 2 ||
      item.length() > 3)
    return -1;
 
  // Cones on X/Y/Z axis
  const std::size_t itemPt((item[1] == '/' && item.length() == 3) ? 2 : 1);
  const auto ptype = static_cast<std::size_t>(tolower(item[itemPt]) - 'x');
  if (ptype >= 3)
    return -2;
  std::vector<double> norm(3, 0.0);
  std::vector<double> cent(3, 0.0);
  norm[ptype] = 1.0;
 
  if (itemPt == 1) // kx type cone
  {
    if (!Mantid::Kernel::Strings::section(Line, cent[ptype]))
      return -3;
  } else {
    std::size_t index;
    for (index = 0; index < 3 && Mantid::Kernel::Strings::section(Line, cent[index]); index++)
      ;
    if (index != 3)
      return -4;
  }
  // The user must enter t^2 which is tan^2(angle) for MCNPX
  double tanAng;
  if (!Mantid::Kernel::Strings::section(Line, tanAng))
    return -5;
 
  Centre = Kernel::V3D(cent[0], cent[1], cent[2]);
  Normal = Kernel::V3D(norm[0], norm[1], norm[2]);
  setTanAngle(sqrt(tanAng));
  setBaseEqn();
  return 0;
}
 
int Cone::operator==(const Cone &A) const
{
  if (this == &A)
    return 1;
  if (fabs(cangle - A.cangle) > Tolerance)
    return 0;
  if (Centre.distance(A.Centre) > Tolerance)
    return 0;
  if (Normal.distance(A.Normal) > Tolerance)
    return 0;
 
  return 1;
}
 
void Cone::setBaseEqn()
{
  const double c2(cangle * cangle);
  const double CdotN(Centre.scalar_prod(Normal));
  BaseEqn[0] = c2 - Normal[0] * Normal[0];                      // A x^2
  BaseEqn[1] = c2 - Normal[1] * Normal[1];                      // B y^2
  BaseEqn[2] = c2 - Normal[2] * Normal[2];                      // C z^2
  BaseEqn[3] = -2 * Normal[0] * Normal[1];                      // D xy
  BaseEqn[4] = -2 * Normal[0] * Normal[2];                      // E xz
  BaseEqn[5] = -2 * Normal[1] * Normal[2];                      // F yz
  BaseEqn[6] = 2.0 * (Normal[0] * CdotN - Centre[0] * c2);      // G x
  BaseEqn[7] = 2.0 * (Normal[1] * CdotN - Centre[1] * c2);      // H y
  BaseEqn[8] = 2.0 * (Normal[2] * CdotN - Centre[2] * c2);      // J z
  BaseEqn[9] = c2 * Centre.scalar_prod(Centre) - CdotN * CdotN; // K const
}
 
void Cone::rotate(const Kernel::Matrix<double> &R)
{
  Centre.rotate(R);
  Normal.rotate(R);
  setBaseEqn();
}
 
void Cone::displace(const Kernel::V3D &A)
{
  Centre += A;
  setBaseEqn();
}
 
void Cone::setCentre(const Kernel::V3D &A)
{
  Centre = A;
  setBaseEqn();
}
 
void Cone::setNorm(const Kernel::V3D &A)
{
  const double norm = A.norm();
  if (norm > Tolerance) {
    Normal = A / norm;
    setBaseEqn();
  }
}
 
void Cone::setAngle(const double A)
{
  alpha = A;
  cangle = cos(M_PI * alpha / 180.0);
  setBaseEqn();
}
 
void Cone::setTanAngle(const double A)
{
  cangle = 1.0 / sqrt(A * A + 1.0); // convert tan(theta) to cos(theta)
  alpha = acos(cangle) * 180.0 / M_PI;
  setBaseEqn();
}
 
double Cone::distance(const Kernel::V3D &Pt) const
{
  const Kernel::V3D Px = Pt - Centre;
  // test is the centre to point distance is zero
  if (Px.norm() < Tolerance)
    return Px.norm();
  double Pangle = Px.scalar_prod(Normal) / Px.norm();
  if (Pangle < 0.0)
    Pangle = acos(-Pangle);
  else
    Pangle = acos(Pangle);
 
  Pangle -= M_PI * alpha / 180.0;
  return Px.norm() * sin(Pangle);
}
 
int Cone::side(const Kernel::V3D &R) const {
 
  const Kernel::V3D cR = R - Centre;
  double rptAngle = cR.scalar_prod(Normal);
  rptAngle *= rptAngle / cR.scalar_prod(cR);
  const double eqn(sqrt(rptAngle));
  if (fabs(eqn - cangle) < Tolerance)
    return 0;
  return (eqn > cangle) ? 1 : -1;
}
 
bool Cone::onSurface(const Kernel::V3D &R) const {
 
  const Kernel::V3D cR = R - Centre;
  double rptAngle = cR.scalar_prod(Normal);
  rptAngle *= rptAngle / cR.scalar_prod(cR);
  const double eqn(sqrt(rptAngle));
 
  return (std::abs(eqn - cangle) <= Tolerance);
}
 
void Cone::write(std::ostream &OX) const
{
  //               -3 -2 -1 0 1 2 3
  const char Tailends[] = "zyx xyz";
  const int Ndir = Normal.masterDir(Tolerance);
  if (Ndir == 0) {
    Quadratic::write(OX);
    return;
  }
  std::ostringstream cx;
  Quadratic::writeHeader(cx);
 
  const int Cdir = Centre.masterDir(Tolerance);
  cx.precision(Surface::Nprecision);
  // Name and transform
 
  if (Cdir || Centre.nullVector(Tolerance)) {
    cx << " k";
    cx << Tailends[Ndir + 3] << " "; // set x,y,z based on Ndir
    cx << ((Cdir > 0) ? Centre[static_cast<std::size_t>(Cdir - 1)] : Centre[static_cast<std::size_t>(-Cdir - 1)]);
    cx << " ";
  } else {
    cx << " k/";
    cx << Tailends[Ndir + 3] << " "; // set x,y,z based on Ndir
    for (std::size_t i = 0; i < 3; i++)
      cx << Centre[i] << " ";
  }
  const double TA = tan((M_PI * alpha) / 180.0); // tan^2(angle)
  cx << TA * TA;
  Mantid::Kernel::Strings::writeMCNPX(cx.str(), OX);
}
 
void Cone::getBoundingBox(double &xmax, double &ymax, double &zmax, double &xmin, double &ymin, double &zmin) {
 
  if (Normal.X() != 0 && Normal.Y() == 0 && Normal.Z() == 0) {
    if (Normal.X() < 0) {
      xmin = std::max(Centre.X(), ymin);
    } else {
      xmax = std::min(Centre.X(), ymax);
    }
    double radius = fabs(xmax - xmin) * sin((M_PI * alpha) / 180.0);
    ymin = Centre.Y() - radius;
    ymax = Centre.Y() + radius;
    zmin = Centre.Z() - radius;
    zmax = Centre.Z() + radius;
  } else if (Normal.X() == 0 && Normal.Y() != 0 && Normal.Z() == 0) {
    if (Normal.Y() < 0) {
      ymin = std::max(Centre.Y(), ymin);
    } else {
      ymax = std::min(Centre.Y(), ymax);
    }
    double radius = fabs(ymax - ymin) * sin((M_PI * alpha) / 180.0);
    xmin = Centre.X() - radius;
    xmax = Centre.X() + radius;
    zmin = Centre.Z() - radius;
    zmax = Centre.Z() + radius;
  } else if (Normal.X() == 0 && Normal.Y() == 0 && Normal.Z() != 0) {
    if (Normal.Z() < 0) {
      zmin = std::max(Centre.Z(), ymin);
    } else {
      zmax = std::min(Centre.Z(), ymax);
    }
    double radius = fabs(zmax - zmin) * sin((M_PI * alpha) / 180.0);
    xmin = Centre.X() - radius;
    xmax = Centre.X() + radius;
    ymin = Centre.Y() - radius;
    ymax = Centre.Y() + radius;
  }
}
 
#ifdef ENABLE_OPENCASCADE
TopoDS_Shape Cone::createShape() {
  gp_Ax2 gpA(gp_Pnt(Centre[0], Centre[1], Centre[2]), gp_Dir(Normal[0], Normal[1], Normal[2]));
  return BRepPrimAPI_MakeCone(gpA, 0.0, 1000.0 / tan(acos(cangle * M_PI / 180.0)), 1000.0, 2.0 * M_PI).Shape();
}
#endif
} // namespace Mantid::Geometry