/**
 * Returns selected information for a "peak" at QLabFrame.
 *
 * @param qFrame      An arbitrary position in Q-space.  This does not have to
 *be the
 *                    position of a peak.
 * @param labCoords  Set true if the position is in the lab coordinate system,
 *false if
 *                    it is in the sample coordinate system.
 * @return a vector whose elements contain different information about the
 *"peak" at that position.
 *         each element is a pair of description of information and the string
 *form for the corresponding
 *         value.
 */
int PeaksWorkspace::peakInfoNumber(Kernel::V3D qFrame, bool labCoords) const {
  std::vector<std::pair<std::string, std::string>> Result;
  std::ostringstream oss;
  oss << std::setw(12) << std::fixed << std::setprecision(3) << (qFrame.norm());
  std::pair<std::string, std::string> QMag("|Q|", oss.str());
  Result.push_back(QMag);

  oss.str("");
  oss.clear();
  oss << std::setw(12) << std::fixed << std::setprecision(3)
      << (2.0 * M_PI / qFrame.norm());

  std::pair<std::string, std::string> dspc("d-spacing", oss.str());
  oss.str("");
  oss.clear();
  Result.push_back(dspc);

  int seqNum = -1;
  double minDist = 10000000;

  for (int i = 0; i < getNumberPeaks(); i++) {
    Peak pk = getPeak(i);
    V3D Q = pk.getQLabFrame();
    if (!labCoords)
      Q = pk.getQSampleFrame();
    double D = qFrame.distance(Q);
    if (D < minDist) {
      minDist = D;
      seqNum = i + 1;
    }
  }
  return seqNum;
}

/**
 * Determines whether point is within the object or on the surface
 * @param point :: Point to be tested
 * @returns true if point is within object or on surface
 */
bool MeshObject::isValid(const Kernel::V3D &point) const {

  BoundingBox bb = getBoundingBox();
  if (!bb.isPointInside(point)) {
    return false;
  }

  Kernel::V3D direction(0.0, 0.0, 1.0); // direction to look for intersections
  std::vector<Kernel::V3D> intersectionPoints;
  std::vector<TrackDirection> entryExitFlags;

  getIntersections(point, direction, intersectionPoints, entryExitFlags);

  if (intersectionPoints.empty()) {
    return false;
  }

  // True if point is on surface
  for (const auto &intersectionPoint : intersectionPoints) {
    if (point.distance(intersectionPoint) < M_TOLERANCE) {
      return true;
    }
  }

  // Look for nearest point then check its entry-exit flag
  double nearestPointDistance = point.distance(intersectionPoints[0]);
  size_t nearestPointIndex = 0;
  for (size_t i = 1; i < intersectionPoints.size(); ++i) {
    if (point.distance(intersectionPoints[i]) < nearestPointDistance) {
      nearestPointDistance = point.distance(intersectionPoints[i]);
      nearestPointIndex = i;
    }
  }
  return (entryExitFlags[nearestPointIndex] == TrackDirection::LEAVING);
}

std::string makeAxisName(const Kernel::V3D &Dir,
                         const std::vector<std::string> &QNames) {
  double eps(1.e-3);
  Kernel::V3D absDir(fabs(Dir.X()), fabs(Dir.Y()), fabs(Dir.Z()));
  std::string mainName;

  if ((absDir[0] >= absDir[1]) && (absDir[0] >= absDir[2])) {
    mainName = QNames[0];
  } else if (absDir[1] >= absDir[2]) {
    mainName = QNames[1];
  } else {
    mainName = QNames[2];
  }

  std::string name("["), separator = ",";
  for (size_t i = 0; i < 3; i++) {

    if (i == 2)
      separator = "]";
    if (absDir[i] < eps) {
      name += "0" + separator;
      continue;
    }
    if (Dir[i] < 0) {
      name += "-";
    }
    if (std::fabs(absDir[i] - 1) < eps) {
      name.append(mainName).append(separator);
      continue;
    }
    name.append(sprintfd(absDir[i], eps)).append(mainName).append(separator);
  }

  return name;
}

/**
 * Return the XML for a sphere.
 */
std::string sphereXML(double radius, const Kernel::V3D &centre, const std::string &id) {
    std::ostringstream xml;
    xml << "<sphere id=\"" << id << "\">"
        << "<centre x=\"" << centre.X() << "\"  y=\"" << centre.Y() << "\" z=\""
        << centre.Z() << "\" />"
        << "<radius val=\"" << radius << "\" />"
        << "</sphere>";
    return xml.str();
}

/**
 Calculate if the point R is on
 the cone (Note: have to be careful here
 since angle calcuation calcuates an angle.
 We need a distance for tolerance!)
 @param R :: Point to check
 @return 1 if on surface and 0 if not not on surface
 */
int 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 (fabs(eqn - cangle) > Tolerance) ? 0 : 1;
}

/**
 * Set the new detector position given the r,theta and phi.
 * @param detectorInfo :: Reference to the DetectorInfo
 * @param index :: Index into detectorInfo
 * @param l2 :: A single l2
 * @param theta :: A single theta
 * @param phi :: A single phi
 */
void UpdateInstrumentFromFile::setDetectorPosition(
    Geometry::DetectorInfo &detectorInfo, const size_t index, const float l2,
    const float theta, const float phi) {
  if (m_ignoreMonitors && detectorInfo.isMonitor(index))
    return;

  Kernel::V3D pos;
  pos.spherical(l2, theta, phi);
  detectorInfo.setPosition(index, pos);
}

double Line::distance(const Kernel::V3D &A) const
/**
Distance of a point from the line
@param A :: test Point
@return absolute distance (not signed)
*/
{
    const double lambda = Direct.scalar_prod(A - Origin);
    Kernel::V3D L = getPoint(lambda);
    L -= A;
    return L.norm();
}

int Line::lambdaPair(
    const int ix,
    const std::pair<std::complex<double>, std::complex<double>> &SQ,
    std::list<Kernel::V3D> &PntOut) const
/**
Helper function to decide which roots to take.
The assumption is that lambda has been solved by quadratic
equation and we require the points that correspond to these
values.
Note: have changed this so that only positive roots are returned.
This makes the quadratic solutions consistent with the ones returned
when asking if a line hits a plane. It is not clear if some other use
cases exist.
@param ix : number of solutions in SQ (0,1,2)
@param SQ : solutions to lambda (the distance along the line
@param PntOut : Output vector of points (added to)
@return Number of real unique points found.
*/
{
    // check results
    if (ix < 1)
        return 0;

    int nCnt(0); // number of good points

    Kernel::V3D Ans;
    if (SQ.first.imag() == 0.0 && SQ.first.real() >= 0.0) // +ve roots only
    {
        const double lambda = SQ.first.real();
        Kernel::V3D Ans = getPoint(lambda);
        PntOut.push_back(Ans);
        if (ix < 2) // only one unique root.
            return 1;
        nCnt = 1;
    }
    if (SQ.second.imag() == 0.0 && SQ.second.real() >= 0.0) // +ve roots only
    {
        const double lambda = SQ.second.real();
        if (!nCnt) // first point wasn't good.
        {
            PntOut.push_back(getPoint(lambda));
            return 1;
        }
        Kernel::V3D Ans2 = getPoint(lambda);
        // If points too close return only 1 item.
        if (Ans.distance(Ans2) < Tolerance)
            return 1;

        PntOut.push_back(Ans2);
        return 2;
    }
    return 0; // both point imaginary
}

/**
 * Estalish if points are coplanar.
 * @param vertices : All vertices to consider
 * @return : Return True only if all coplanar
 */
bool MeshObject2D::pointsCoplanar(const std::vector<Kernel::V3D> &vertices) {
  if (!CoplanarChecks::sufficientPoints(vertices))
    return false;

  Kernel::V3D normal = CoplanarChecks::surfaceNormal(vertices);
  // Check that a valid normal was found amongst collection of vertices
  if (normal.norm2() == 0) {
    // If all points are colinear. Not a plane.
    return false;
  }

  return CoplanarChecks::allCoplanar(vertices, normal);
}

/**
 * Find maximum angular half width of the bounding box from the observer, that
 * is
 * the greatest angle between the centre point and any corner point
 * @param observer :: Viewing point
 * @returns The value of the angular half-width
*/
double BoundingBox::angularWidth(const Kernel::V3D &observer) const {
  Kernel::V3D centre = centrePoint() - observer;
  std::vector<Kernel::V3D> pts;
  this->getFullBox(pts, observer);

  std::vector<Kernel::V3D>::const_iterator ip;
  double centre_norm_inv = 1.0 / centre.norm();
  double thetaMax(-1.0);
  for (ip = pts.begin(); ip != pts.end(); ++ip) {
    double theta = acos(ip->scalar_prod(centre) * centre_norm_inv / ip->norm());
    if (theta > thetaMax)
      thetaMax = theta;
  }
  return thetaMax;
}

/**
 * Determines wither point is on the surface.
 * @param point :: Point to check
 * @returns true if the point is on the surface
 */
bool MeshObject::isOnSide(const Kernel::V3D &point) const {

  BoundingBox bb = getBoundingBox();
  if (!bb.isPointInside(point)) {
    return false;
  }

  const std::vector<Kernel::V3D> directions = {
      Kernel::V3D{0, 0, 1}, Kernel::V3D{0, 1, 0},
      Kernel::V3D{1, 0, 0}}; // directions to look for intersections
  // We have to look in several directions in case a point is on a face
  // or edge parallel to the first direction or also the second direction.
  for (const auto &direction : directions) {
    std::vector<Kernel::V3D> intersectionPoints;
    std::vector<TrackDirection> entryExitFlags;

    getIntersections(point, direction, intersectionPoints, entryExitFlags);

    if (intersectionPoints.empty()) {
      return false;
    }

    for (const auto &intersectionPoint : intersectionPoints) {
      if (point.distance(intersectionPoint) < M_TOLERANCE) {
        return true;
      }
    }
  }
  return false;
}

double Cylinder::distance(const Kernel::V3D &A) const
/**
 Calculates the distance of point A from
 the surface of the  cylinder.
 @param A :: Point to calculate distance from
 @return :: +ve distance to the surface.

 \todo INCOMPLETE AS Does not deal with the case of
 non axis centre line  (has been updated?? )
 */
{
  // First find the normal going to the point
  const Kernel::V3D Amov = A - Centre;
  double lambda = Amov.scalar_prod(Normal);
  const Kernel::V3D Ccut = Normal * lambda;
  // The distance is from the centre line to the
  return fabs(Ccut.distance(Amov) - Radius);
}

    /**
     * Set the new detector position given the r,theta and phi.
     * @param det :: A pointer to the detector
     * @param l2 :: A single l2
     * @param theta :: A single theta
     * @param phi :: A single phi
     */
    void UpdateInstrumentFromFile::setDetectorPosition(const Geometry::IDetector_const_sptr & det, const float l2,
                                                       const float theta, const float phi)
    {
      if( m_ignoreMonitors && det->isMonitor() ) return;

      Geometry::ParameterMap & pmap = m_workspace->instrumentParameters();
      Kernel::V3D pos;
      if (!m_ignorePhi)
      {
        pos.spherical(l2, theta, phi);
      }
      else
      {
        double r,t,p;
        det->getPos().getSpherical(r,t,p);
        pos.spherical(l2, theta, p);
      }
      Geometry::ComponentHelper::moveComponent(*det, pmap, pos, Geometry::ComponentHelper::Absolute);
    }

std::pair<double, double> LoadILLSANS::calculateQMaxQMin() {
  double min = std::numeric_limits<double>::max(),
         max = std::numeric_limits<double>::min();
  g_log.debug("Calculating Qmin Qmax...");
  std::size_t nHist = m_localWorkspace->getNumberHistograms();
  for (std::size_t i = 0; i < nHist; ++i) {
    Geometry::IDetector_const_sptr det = m_localWorkspace->getDetector(i);
    if (!det->isMonitor()) {
      const MantidVec &lambdaBinning = m_localWorkspace->readX(i);
      Kernel::V3D detPos = det->getPos();
      double r, theta, phi;
      detPos.getSpherical(r, theta, phi);
      double v1 = calculateQ(*(lambdaBinning.begin()), theta);
      double v2 = calculateQ(*(lambdaBinning.end() - 1), theta);
      // std::cout << "i=" << i << " theta="<<theta << " lambda_i=" <<
      // *(lambdaBinning.begin()) << " lambda_f=" << *(lambdaBinning.end()-1) <<
      // " v1=" << v1 << " v2=" << v2 << '\n';
      if (i == 0) {
        min = v1;
        max = v1;
      }
      if (v1 < min) {
        min = v1;
      }
      if (v2 < min) {
        min = v2;
      }
      if (v1 > max) {
        max = v1;
      }
      if (v2 > max) {
        max = v2;
      }
    } else
      g_log.debug() << "Detector " << i << " is a Monitor : " << det->getID()
                    << '\n';
  }

  g_log.debug() << "Calculating Qmin Qmax. Done : [" << min << "," << max
                << "]\n";

  return std::pair<double, double>(min, max);
}

double
Torus::distance(const Kernel::V3D& Pt) const
  /**
    Calculates the distance from the point to the Torus
    does not calculate the point on the Torusthat is closest
    @param Pt :: Point to calcuate from

    - normalise to a cone vertex at the origin
    - calculate the angle between the axis and the Point
    - Calculate the distance to P
    @return distance to Pt
  */
{
  const Kernel::V3D Px=Pt-Centre;
  // test is the centre to point distance is zero
  if(Px.norm()<Tolerance)
    return Px.norm();
  return Px.norm();
}

double Plane::distance(const Kernel::V3D &A) const
/**
  Determine the distance of point A from the plane
  returns a value relative to the normal
  @param A :: point to get distance from
  @return singed distance from point
*/
{
  return A.scalar_prod(NormV) - Dist;
}

int Line::intersect(std::list<Kernel::V3D> &PntOut, const Cylinder &Cyl) const
/**
For the line that intersects the cylinder generate
add the point to the VecOut, return number of points
added. It does not check the points for validity.

@param PntOut :: Vector of points found by the line/cylinder intersection
@param Cyl :: Cylinder to intersect line with
@return Number of points found by intersection
*/
{
    const Kernel::V3D Cent = Cyl.getCentre();
    const Kernel::V3D Ax = Origin - Cent;
    const Kernel::V3D N = Cyl.getNormal();
    const double R = Cyl.getRadius();
    const double vDn = N.scalar_prod(Direct);
    const double vDA = N.scalar_prod(Ax);
    // First solve the equation of intersection
    double C[3];
    C[0] = 1.0 - (vDn * vDn);
    C[1] = 2.0 * (Ax.scalar_prod(Direct) - vDA * vDn);
    C[2] = Ax.scalar_prod(Ax) - (R * R + vDA * vDA);
    std::pair<std::complex<double>, std::complex<double>> SQ;
    const int ix = solveQuadratic(C, SQ);
    // This takes the centre displacement into account:
    return lambdaPair(ix, SQ, PntOut);
}

    double Cone::distance(const Kernel::V3D& Pt) const
    /**
     Calculates the distance from the point to the Cone
     does not calculate the point on the cone that is closest
     @param Pt :: Point to calcuate from

     - normalise to a cone vertex at the origin
     - calculate the angle between the axis and the Point
     - Calculate the distance to P
     @return distance to Pt
     */
    {
      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);
    }

/**
 * This function calculates the exponential contribution to the He3 tube
 * efficiency.
 * @param spectraIndex :: the current index to calculate
 * @param idet :: the current detector pointer
 * @throw out_of_range if twice tube thickness is greater than tube diameter
 * @return the exponential contribution for the given detector
 */
double He3TubeEfficiency::calculateExponential(
    std::size_t spectraIndex,
    boost::shared_ptr<const Geometry::IDetector> idet) {
  // Get the parameters for the current associated tube
  double pressure =
      this->getParameter("TubePressure", spectraIndex, "tube_pressure", idet);
  double tubethickness =
      this->getParameter("TubeThickness", spectraIndex, "tube_thickness", idet);
  double temperature = this->getParameter("TubeTemperature", spectraIndex,
                                          "tube_temperature", idet);

  double detRadius(0.0);
  Kernel::V3D detAxis;
  this->getDetectorGeometry(idet, detRadius, detAxis);
  double detDiameter = 2.0 * detRadius;
  double twiceTubeThickness = 2.0 * tubethickness;

  // now get the sin of the angle, it's the magnitude of the cross product of
  // unit vector along the detector tube axis and a unit vector directed from
  // the sample to the detector center
  Kernel::V3D vectorFromSample = idet->getPos() - this->samplePos;
  vectorFromSample.normalize();
  Kernel::Quat rot = idet->getRotation();
  // rotate the original cylinder object axis to get the detector axis in the
  // actual instrument
  rot.rotate(detAxis);
  detAxis.normalize();
  // Scalar product is quicker than cross product
  double cosTheta = detAxis.scalar_prod(vectorFromSample);
  double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);

  const double straight_path = detDiameter - twiceTubeThickness;
  if (std::fabs(straight_path - 0.0) < TOL) {
    throw std::out_of_range("Twice tube thickness cannot be greater than "
                            "or equal to the tube diameter");
  }

  const double pathlength = straight_path / sinTheta;
  return EXP_SCALAR_CONST * (pressure / temperature) * pathlength;
}

void Cone::setNorm(const Kernel::V3D &A)
/**
 Sets the Normal and the Base Equation
 @param A :: New Normal direction
 */
{
    if (A.norm() > Tolerance) {
        Normal = A;
        Normal.normalize();
        setBaseEqn();
    }
    return;
}