//
// special.hpp
//
// LGPL Version 2.1 HEADER START
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
//
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
// MA 02110-1301  USA
//
// LGPL Version 2.1 HEADER END
//

//
// Copyright (c) 2020-2021, Regents of the University of Minnesota.
// All rights reserved.
//
// Contributors:
//    Yaser Afshar
//

#ifndef SPECIAL_HPP
#define SPECIAL_HPP

#include <cmath>

#include "helper.hpp"

/*!
 * \file special.hpp
 *
 * \brief This file contains special helper functions
 *
 * The special helper functions are:
 * \c IsNotOne return true if the input is not (close to) one within a tolerance
 * \c Dot3 reurn dot product of 2 vectors
 * \c Cross3 compute the cross product of 2 vectors
 * \c Len3 return the length of vector v
 * \c CosTheta return the computed cos theta
 */

namespace special {
static inline bool IsNotOne(double const &f, double const tol = 1e-20) {
  return std::abs(f - 1.0) >= tol;
}

static inline constexpr double Dot3(VectorOfSizeDIM const v1,
                                    VectorOfSizeDIM const v2) {
  return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}

static inline void Cross3(VectorOfSizeDIM const v1, VectorOfSizeDIM const v2,
                          VectorOfSizeDIM ans) {
  ans[0] = v1[1] * v2[2] - v1[2] * v2[1];
  ans[1] = v1[2] * v2[0] - v1[0] * v2[2];
  ans[2] = v1[0] * v2[1] - v1[1] * v2[0];
}

static inline double Len3(VectorOfSizeDIM const v) {
  return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
}

static inline void CosTheta(VectorOfSizeDIM const dr_ij, double const &r_ij,
                            VectorOfSizeDIM const dr_ik, double const &r_ik,
                            VectorOfSizeDIM dri, VectorOfSizeDIM drj,
                            VectorOfSizeDIM drk) {
  double const r_ij_inv = 1.0 / r_ij;
  double const r_ik_inv = 1.0 / r_ik;
  double const s1 = r_ij_inv * r_ik_inv;
  double const costheta = Dot3(dr_ij, dr_ik) * s1;
  double const s2 = -costheta / r_ij / r_ij;
  drj[0] = s1 * dr_ik[0] + s2 * dr_ij[0];
  drj[1] = s1 * dr_ik[1] + s2 * dr_ij[1];
  drj[2] = s1 * dr_ik[2] + s2 * dr_ij[2];
  double const s3 = -costheta * r_ik_inv * r_ik_inv;
  drk[0] = s1 * dr_ij[0] + s3 * dr_ik[0];
  drk[1] = s1 * dr_ij[1] + s3 * dr_ik[1];
  drk[2] = s1 * dr_ij[2] + s3 * dr_ik[2];
  dri[0] = -(drj[0] + drk[0]);
  dri[1] = -(drj[1] + drk[1]);
  dri[2] = -(drj[2] + drk[2]);
}
}  // namespace special

#endif  // SPECIAL