Vec3D.h

Go to the documentation of this file.

/*
 * Vec3D.h
 *
 *  Created on: 18 Mar 2019
 *      Author: Valerio Formato
 */
 
#ifndef HERD_VEC3D_H_
#define HERD_VEC3D_H_
 
#include "herd/analysis/common/CooArray.h"
#include "herd/analysis/common/RefFrame.h"
 
#include <array>
#include <cmath>
#include <iostream>
#include <limits>
#include <type_traits>
 
namespace Herd {
 
class Vec3D;
inline float operator*(const Vec3D left, const Vec3D right);
 
class Vec3D : public CooArray<float> {
public:
  friend inline std::ostream &operator<<(std::ostream &os, const Vec3D &vec);
 
  constexpr Vec3D()
      : CooArray<float>(std::numeric_limits<float>::max(), std::numeric_limits<float>::max(),
                        std::numeric_limits<float>::max()) {}
 
  constexpr Vec3D(float x, float y, float z) : CooArray<float>(x, y, z) {}
 
  constexpr Vec3D(std::array<float, 3> arr) : CooArray(std::forward<std::array<float, 3>>(arr)) {}
 
  constexpr Vec3D(RefFrame::Direction direction) : CooArray<float>(0, 0, 0) {
    switch (direction) {
    case RefFrame::Direction::Xpos:
      (*this)[RefFrame::Coo::X] = 1;
      break;
    case RefFrame::Direction::Xneg:
      (*this)[RefFrame::Coo::X] = -1;
      break;
    case RefFrame::Direction::Ypos:
      (*this)[RefFrame::Coo::Y] = 1;
      break;
    case RefFrame::Direction::Yneg:
      (*this)[RefFrame::Coo::Y] = -1;
      break;
    case RefFrame::Direction::Zpos:
      (*this)[RefFrame::Coo::Z] = 1;
      break;
    case RefFrame::Direction::Zneg:
      (*this)[RefFrame::Coo::Z] = -1;
      break;
    case RefFrame::Direction::XposYpos:
      (*this)[RefFrame::Coo::X] = 1;
      (*this)[RefFrame::Coo::Y] = 1;
      (*this) /= this->Mag();
      break;
    case RefFrame::Direction::XposYneg:
      (*this)[RefFrame::Coo::X] = 1;
      (*this)[RefFrame::Coo::Y] = -1;
      (*this) /= this->Mag();
      break;
    case RefFrame::Direction::XnegYpos:
      (*this)[RefFrame::Coo::X] = -1;
      (*this)[RefFrame::Coo::Y] = 1;
      (*this) /= this->Mag();
      break;
    case RefFrame::Direction::XnegYneg:
      (*this)[RefFrame::Coo::X] = -1;
      (*this)[RefFrame::Coo::Y] = -1;
      (*this) /= this->Mag();
      break;
    default:
      break;
    }
  }
 
  void Dump() const { std::cout << *this << std::endl; }
 
  [[nodiscard]] float Mag() const { return std::sqrt((*this) * (*this)); }
 
  [[nodiscard]] Vec3D UnitVector() const {
    Vec3D result = *this;
    return result /= Mag();
  }
 
  [[nodiscard]] constexpr Vec3D CrossProduct(const Vec3D &rhs) const {
    float a = (*this)[RefFrame::Coo::Y] * rhs[RefFrame::Coo::Z] - (*this)[RefFrame::Coo::Z] * rhs[RefFrame::Coo::Y];
    float b = (*this)[RefFrame::Coo::Z] * rhs[RefFrame::Coo::X] - (*this)[RefFrame::Coo::X] * rhs[RefFrame::Coo::Z];
    float c = (*this)[RefFrame::Coo::X] * rhs[RefFrame::Coo::Y] - (*this)[RefFrame::Coo::Y] * rhs[RefFrame::Coo::X];
    Vec3D product(a, b, c);
    return product;
  }
 
  class UndefinedAngle : public Exception {
  public:
    using Exception::Exception;
  };
 
  [[nodiscard]] float PolarAngle() const {
    // polar is always defined in the [0, M_PI) range
    return M_PI - std::acos(static_cast<double>((*this)[RefFrame::Coo::Z] / Mag()));
  }
 
  [[nodiscard]] float AzimuthAngle() const {
    if (PolarAngle() == 0.0f || std::fabs(PolarAngle() - M_PI) < std::numeric_limits<float>::epsilon())
      throw UndefinedAngle("Vec3D::AzimuthAngle is undefined when polar == 0 or pi");
    return std::atan2((*this)[RefFrame::Coo::Y], (*this)[RefFrame::Coo::X]);
  }
 
  Vec3D &operator+=(Vec3D const &right) {
    (*this)[RefFrame::Coo::X] += right[RefFrame::Coo::X];
    (*this)[RefFrame::Coo::Y] += right[RefFrame::Coo::Y];
    (*this)[RefFrame::Coo::Z] += right[RefFrame::Coo::Z];
 
    return *this;
  };
  Vec3D &operator+=(Vec3D &&right) {
    (*this)[RefFrame::Coo::X] += right[RefFrame::Coo::X];
    (*this)[RefFrame::Coo::Y] += right[RefFrame::Coo::Y];
    (*this)[RefFrame::Coo::Z] += right[RefFrame::Coo::Z];
 
    return *this;
  };
 
  Vec3D operator-(Vec3D const &right) const {
    Vec3D retValue;
    retValue[RefFrame::Coo::X] = (*this)[RefFrame::Coo::X] - right[RefFrame::Coo::X];
    retValue[RefFrame::Coo::Y] = (*this)[RefFrame::Coo::Y] - right[RefFrame::Coo::Y];
    retValue[RefFrame::Coo::Z] = (*this)[RefFrame::Coo::Z] - right[RefFrame::Coo::Z];
 
    return retValue;
  };
 
  Vec3D &operator-=(Vec3D const &right) {
    (*this)[RefFrame::Coo::X] -= right[RefFrame::Coo::X];
    (*this)[RefFrame::Coo::Y] -= right[RefFrame::Coo::Y];
    (*this)[RefFrame::Coo::Z] -= right[RefFrame::Coo::Z];
 
    return *this;
  };
  Vec3D &operator-=(Vec3D &&right) {
    (*this)[RefFrame::Coo::X] -= right[RefFrame::Coo::X];
    (*this)[RefFrame::Coo::Y] -= right[RefFrame::Coo::Y];
    (*this)[RefFrame::Coo::Z] -= right[RefFrame::Coo::Z];
 
    return *this;
  };
 
  // This neat trick restricts the rhs type to be a number (bool is also viable
  // but nothing is perfect, right?)
  template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
  Vec3D &operator*=(T const &right) {
    (*this)[RefFrame::Coo::X] *= right;
    (*this)[RefFrame::Coo::Y] *= right;
    (*this)[RefFrame::Coo::Z] *= right;
 
    return *this;
  }
 
  template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
  Vec3D &operator*=(T &&right) {
    (*this)[RefFrame::Coo::X] *= right;
    (*this)[RefFrame::Coo::Y] *= right;
    (*this)[RefFrame::Coo::Z] *= right;
 
    return *this;
  }
 
  // This neat trick restricts the rhs type to be a number (bool is also viable
  // but nothing is perfect, right?)
  template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
  Vec3D &operator/=(T const &right) {
    (*this)[RefFrame::Coo::X] /= right;
    (*this)[RefFrame::Coo::Y] /= right;
    (*this)[RefFrame::Coo::Z] /= right;
 
    return *this;
  }
 
  template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
  Vec3D &operator/=(T &&right) {
    (*this)[RefFrame::Coo::X] /= right;
    (*this)[RefFrame::Coo::Y] /= right;
    (*this)[RefFrame::Coo::Z] /= right;
 
    return *this;
  }
 
  bool operator==(const Vec3D &right) const {
    return (((*this)[RefFrame::Coo::X] == right[RefFrame::Coo::X]) &&
            ((*this)[RefFrame::Coo::Y] == right[RefFrame::Coo::Y]) &&
            ((*this)[RefFrame::Coo::Z] == right[RefFrame::Coo::Z]));
  }
};
 
// VF: Check here for this neat trick to enable move semantics:
// https://stackoverflow.com/questions/9851188/does-it-make-sense-to-use-move-semantics-for-operator-and-or-operator
inline Vec3D operator+(Vec3D left, Vec3D const &right) {
  left += right;
  return left;
}
inline Vec3D operator+(const Vec3D &left, Vec3D &&right) {
  right += left;
  return std::move(right);
}
 
// This neat trick restricts the rhs type to be a number (bool is also viable
// but nothing is perfect, right?)
template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
inline Vec3D operator*(Vec3D left, T const &right) {
  left *= right;
  return left;
}
 
// symmetrization
template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
inline Vec3D operator*(T const &right, Vec3D left) {
  return left * right;
}
 
// This neat trick restricts the rhs type to be a number (bool is also viable
// but nothing is perfect, right?)
template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
inline Vec3D operator/(Vec3D left, T const &right) {
  left /= right;
  return left;
}
 
// symmetrization
template <typename T, typename std::enable_if<std::is_arithmetic<T>::value>::type * = nullptr>
inline Vec3D operator/(T const &right, Vec3D left) {
  return left / right;
}
 
inline float operator*(const Vec3D left, const Vec3D right) {
  return (left[RefFrame::Coo::X] * right[RefFrame::Coo::X] + left[RefFrame::Coo::Y] * right[RefFrame::Coo::Y] +
          left[RefFrame::Coo::Z] * right[RefFrame::Coo::Z]);
}
 
std::ostream &operator<<(std::ostream &os, const Vec3D &vec) {
  os << "( " << vec[RefFrame::Coo::X] << ", " << vec[RefFrame::Coo::Y] << ", " << vec[RefFrame::Coo::Z] << " )";
  return os;
}
} // namespace Herd
 
#endif