angles.h

Go to the documentation of this file.

/*
 *
 *    Copyright (c) 2014-2018
 *      SMASH Team
 *
 *    GNU General Public License (GPLv3 or later)
 *
 */

#ifndef SRC_INCLUDE_ANGLES_H_
#define SRC_INCLUDE_ANGLES_H_

#include <iostream>
#include <stdexcept>

#include "constants.h"
#include "iomanipulators.h"
#include "random.h"
#include "threevector.h"

namespace smash {

/*
 * Internals
 * ---------
 *
 * The object internally stores the azimuthal angle \f$\varphi\f$ and
 * the cosine of the polar angle \f$\cos\vartheta\f$. Nobody should rely
 * on this never changing, though; the interface user should be totally
 * oblivious to this.
 *
 * Possible future improvements
 * ----------------------------
 *
 * More distributions need to be implemented once there is a physics
 * case to use them.
 */
class Angles {
 public:
  Angles() : phi_(0), costheta_(0) {}
  Angles(double ph, double cost) {
    set_phi(ph);
    set_costheta(cost);
  }
  void distribute_isotropically();
  void set_phi(const double phi);
  void set_costheta(const double cos);
  void set_theta(const double theta);
  bool add_to_theta(const double delta);
  bool add_to_theta(const double delta, const bool reverse);
  double phi() const;
  double costheta() const;
  double sintheta() const;
  double x() const;
  double y() const;
  double z() const;
  ThreeVector inline threevec() const;
  double theta() const;

  struct InvalidTheta : public std::invalid_argument {
    using std::invalid_argument::invalid_argument;
  };

 private:
  double phi_;
  double costheta_;
};

inline std::ostream &operator<<(std::ostream &out, const Angles &a) {
  return out << "φ:" << field << a.phi() << ", cos ϑ:" << field << a.costheta();
}

void inline Angles::distribute_isotropically() {
  /* Isotropic distribution: phi in [0, 2pi) and cos(theta) in [-1,1]. */
  phi_ = random::uniform(0.0, twopi);
  costheta_ = random::uniform(-1.0, 1.0);
}

void inline Angles::set_phi(const double newphi) {
  /* Make sure that phi is in the range [0,2pi).  */
  phi_ = newphi;
  if (newphi < 0 || newphi >= twopi) {
    phi_ -= twopi * std::floor(newphi / twopi);
  }
}

void inline Angles::set_costheta(const double newcos) {
  costheta_ = newcos;
  /* check if costheta_ is in -1..1. If not, well. Error handling here
   * is a lot harder than in the above. Still, I will silently do the
   * same as above. Note, though, that costheta = 1 is allowed, even if
   * it cannot be generated by distribute_isotropically(). */
  if ((costheta_ < -1. - really_small) || (costheta_ > 1. + really_small)) {
    throw InvalidTheta("Wrong value for costheta (must be in [-1,1]): " +
                       std::to_string(costheta_));
  }
  if (costheta_ < -1.) {
    costheta_ = -1.;
  } else if (costheta_ > 1.) {
    costheta_ = 1.;
  }
}
void inline Angles::set_theta(const double newtheta) {
  /* no error handling necessary, because this gives a sensible answer
   * for every real number. */
  set_costheta(std::cos(newtheta));
}

bool inline Angles::add_to_theta(const double delta) {
  if (delta < -M_PI || delta > M_PI) {
    throw InvalidTheta("Cannot advance polar angle by " +
                       std::to_string(delta));
  }
  double theta_plus_delta = delta + theta();
  /* if sum is not in [0, PI], force it to be there:
   * "upper" overflow:
   * theta + delta + the_new_angle = 2*M_PI */
  if (theta_plus_delta > M_PI) {
    set_theta(twopi - theta_plus_delta);
    // set_phi takes care that phi_ is in [0 .. 2*M_PI]
    set_phi(phi() + M_PI);
    return true;  // meaning "we did change phi"
    // "lower" overflow: theta + delta switches sign
  } else if (theta_plus_delta < 0) {
    set_theta(-theta_plus_delta);
    set_phi(phi() + M_PI);
    return true;  // meaning "we did change phi"
    // no overflow: set theta, do not touch phi:
  } else {
    set_theta(theta_plus_delta);
  }
  return false;  // meaning "we did NOT change phi"
}
bool inline Angles::add_to_theta(const double delta, const bool reverse) {
  double plusminus_one = reverse ? -1.0 : +1.0;
  bool this_reverse = add_to_theta(plusminus_one * delta);
  /* if we had to reverse first time and now reverse again OR if we
   * didn't reverse in either part, we do not reverse in total.
   * else: if we reverse in one, but not the other part, we reverse in
   * total. */
  return this_reverse ^ reverse;
}

double inline Angles::costheta() const { return costheta_; }
double inline Angles::phi() const { return phi_; }
double inline Angles::sintheta() const {
  return std::sqrt(1.0 - costheta_ * costheta_);
}
double inline Angles::x() const { return sintheta() * std::cos(phi_); }
double inline Angles::y() const { return sintheta() * std::sin(phi_); }
double inline Angles::z() const { return costheta_; }

ThreeVector inline Angles::threevec() const {
  return ThreeVector(x(), y(), z());
}

double inline Angles::theta() const { return std::acos(costheta_); }

}  // namespace smash

#endif  // SRC_INCLUDE_ANGLES_H_