threevector.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2015,2017-2020,2022,2024-2025
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */
 
#ifndef SRC_INCLUDE_SMASH_THREEVECTOR_H_
#define SRC_INCLUDE_SMASH_THREEVECTOR_H_
 
#include <array>
#include <cmath>
#include <ostream>
 
#include "constants.h"
 
namespace smash {
 
class ThreeVector {
 public:
  ThreeVector() : x_({0., 0., 0.}) {}
 
  ThreeVector(double y1, double y2, double y3) : x_({y1, y2, y3}) {}
 
  explicit ThreeVector(std::array<double, 3> arr)
      : x_({arr[0], arr[1], arr[2]}) {}
 
  double &operator[](std::size_t i) { return x_[i]; }
  double operator[](std::size_t i) const { return x_[i]; }
 
  double inline x1() const;
  void inline set_x1(double x);
  double inline x2() const;
  void inline set_x2(double y);
  double inline x3() const;
  void inline set_x3(double z);
  double inline sqr() const;
  double inline abs() const;
  double inline get_phi() const;
  double inline get_theta() const;
  void inline rotate(double phi, double theta, double psi);
  void inline rotate_around_y(double theta);
  void inline rotate_around_z(double theta);
  void inline rotate_z_axis_to(ThreeVector &r);
  ThreeVector inline operator-() const;
  ThreeVector inline operator+=(const ThreeVector &v);
  ThreeVector inline operator-=(const ThreeVector &v);
  ThreeVector inline operator*=(const double &a);
  ThreeVector inline operator/=(const double &a);
 
  bool operator==(const ThreeVector &rhs) const { return x_ == rhs.x_; }
  bool operator!=(const ThreeVector &rhs) const { return x_ != rhs.x_; }
 
  ThreeVector inline cross_product(const ThreeVector &b) const;
  using iterator = std::array<double, 3>::iterator;
  using const_iterator = std::array<double, 3>::const_iterator;
 
  iterator begin() { return x_.begin(); }
 
  iterator end() { return x_.end(); }
 
  const_iterator begin() const { return x_.begin(); }
  const_iterator end() const { return x_.end(); }
 
  const_iterator cbegin() const { return x_.cbegin(); }
  const_iterator cend() const { return x_.cend(); }
 
 private:
  std::array<double, 3> x_;
};
 
std::ostream &operator<<(std::ostream &, const ThreeVector &);
 
double inline ThreeVector::x1() const { return x_[0]; }
 
void inline ThreeVector::set_x1(const double x) { x_[0] = x; }
 
double inline ThreeVector::x2() const { return x_[1]; }
 
void inline ThreeVector::set_x2(const double y) { x_[1] = y; }
 
double inline ThreeVector::x3() const { return x_[2]; }
 
void inline ThreeVector::set_x3(const double z) { x_[2] = z; }
 
ThreeVector inline ThreeVector::operator-() const {
  ThreeVector neg(-x_[0], -x_[1], -x_[2]);
  return neg;
}
 
ThreeVector inline ThreeVector::operator+=(const ThreeVector &v) {
  x_[0] += v.x_[0];
  x_[1] += v.x_[1];
  x_[2] += v.x_[2];
  return *this;
}
 
ThreeVector inline operator+(ThreeVector a, const ThreeVector &b) {
  a += b;
  return a;
}
 
ThreeVector inline ThreeVector::operator-=(const ThreeVector &v) {
  x_[0] -= v.x_[0];
  x_[1] -= v.x_[1];
  x_[2] -= v.x_[2];
  return *this;
}
 
ThreeVector inline operator-(ThreeVector a, const ThreeVector &b) {
  a -= b;
  return a;
}
 
ThreeVector inline ThreeVector::operator*=(const double &a) {
  x_[0] *= a;
  x_[1] *= a;
  x_[2] *= a;
  return *this;
}
 
inline ThreeVector operator*(ThreeVector a, const double &b) {
  a *= b;
  return a;
}
 
inline ThreeVector operator*(const double &a, ThreeVector b) {
  b *= a;
  return b;
}
 
inline double operator*(ThreeVector a, const ThreeVector &b) {
  return a.x1() * b.x1() + a.x2() * b.x2() + a.x3() * b.x3();
}
 
ThreeVector inline ThreeVector::cross_product(const ThreeVector &b) const {
  return ThreeVector(x_[1] * b.x3() - x_[2] * b.x2(),
                     x_[2] * b.x1() - x_[0] * b.x3(),
                     x_[0] * b.x2() - x_[1] * b.x1());
}
 
ThreeVector inline ThreeVector::operator/=(const double &a) {
  const double a_inv = 1.0 / a;
  x_[0] *= a_inv;
  x_[1] *= a_inv;
  x_[2] *= a_inv;
  return *this;
}
 
ThreeVector inline operator/(ThreeVector a, const double &b) {
  a /= b;
  return a;
}
 
double inline ThreeVector::sqr() const { return (*this) * (*this); }
 
double inline ThreeVector::abs() const { return std::sqrt((*this) * (*this)); }
 
double inline ThreeVector::get_phi() const {
  if (std::abs(x1()) < really_small && std::abs(x2()) < really_small) {
    return 0.;
  } else {
    return std::atan2(x2(), x1());
  }
}
 
double inline ThreeVector::get_theta() const {
  double r = abs();
  return (r > 0.) ? std::acos(x3() / r) : 0.;
}
 
void inline ThreeVector::rotate(double phi, double theta, double psi) {
  // Compute the cosine and sine for each angle.
  const double cos_phi = std::cos(phi);
  const double sin_phi = std::sin(phi);
  const double cos_theta = std::cos(theta);
  const double sin_theta = std::sin(theta);
  const double cos_psi = std::cos(psi);
  const double sin_psi = std::sin(psi);
  // Get original coordinates.
  std::array<double, 3> x_old = x_;
  // Compute new coordinates.
  x_[0] = (cos_phi * cos_psi - sin_phi * cos_theta * sin_psi) * x_old[0] +
          (-cos_phi * sin_psi - sin_phi * cos_theta * cos_psi) * x_old[1] +
          (sin_phi * sin_theta) * x_old[2];
  x_[1] = (sin_phi * cos_psi + cos_phi * cos_theta * sin_psi) * x_old[0] +
          (-sin_phi * sin_psi + cos_phi * cos_theta * cos_psi) * x_old[1] +
          (-cos_phi * sin_theta) * x_old[2];
  x_[2] = (sin_theta * sin_psi) * x_old[0] + (sin_theta * cos_psi) * x_old[1] +
          (cos_theta)*x_old[2];
}
 
void inline ThreeVector::rotate_around_y(double theta) {
  const double cost = std::cos(theta);
  const double sint = std::sin(theta);
  // Get original coordinates.
  std::array<double, 3> x_old = x_;
  // Compute new coordinates.
  x_[0] = cost * x_old[0] + sint * x_old[2];
  // x_[1] is unchanged
  x_[2] = -sint * x_old[0] + cost * x_old[2];
}
 
void inline ThreeVector::rotate_around_z(double theta) {
  const double cost = std::cos(theta);
  const double sint = std::sin(theta);
  // Get original coordinates.
  std::array<double, 3> x_old = x_;
  // Compute new coordinates.
  x_[0] = cost * x_old[0] - sint * x_old[1];
  x_[1] = sint * x_old[0] + cost * x_old[1];
  // x_[2] is unchanged
}
 
void inline ThreeVector::rotate_z_axis_to(ThreeVector &r) {
  rotate_around_y(r.get_theta());
  rotate_around_z(r.get_phi());
}
 
}  // namespace smash
 
#endif  // SRC_INCLUDE_SMASH_THREEVECTOR_H_