array_vector.hpp

// This file is part of the ACTS project.
//
// Copyright (C) 2016 CERN for the benefit of the ACTS project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at https://mozilla.org/MPL/2.0/.
 
#pragma once
 
// Project include(s).
#include "algebra/impl/array_operators.hpp"
#include "detray/algebra/common/math.hpp"
#include "detray/algebra/concepts.hpp"
#include "detray/algebra/generic/generic.hpp"
#include "detray/definitions/detail/qualifiers.hpp"
 
namespace detray::algebra::array {

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 2)
DETRAY_HOST_DEVICE constexpr scalar_t phi(const array_t<scalar_t, N> &v) {
  return algebra::generic::math::phi(v);
}

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 2)
DETRAY_HOST_DEVICE constexpr scalar_t perp(const array_t<scalar_t, N> &v) {
  return algebra::generic::math::perp(v);
}

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 2)
DETRAY_HOST_DEVICE constexpr scalar_t theta(const array_t<scalar_t, N> &v) {
  return algebra::generic::math::theta(v);
}

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 3)
DETRAY_HOST_DEVICE constexpr array_t<scalar_t, N> cross(
    const array_t<scalar_t, N> &a, const array_t<scalar_t, N> &b) {
  return algebra::generic::math::cross(a, b);
}
 
template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 3)
DETRAY_HOST_DEVICE constexpr array_t<scalar_t, N> cross(
    const array_t<scalar_t, N> &a, const array_t<array_t<scalar_t, N>, 1> &b) {
  return algebra::generic::math::cross(a, b);
}
 
template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 3)
DETRAY_HOST_DEVICE constexpr array_t<scalar_t, N> cross(
    const array_t<array_t<scalar_t, N>, 1> &a, const array_t<scalar_t, N> &b) {
  return algebra::generic::math::cross(a, b);
}
 
template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 3)
DETRAY_HOST_DEVICE constexpr array_t<scalar_t, N> cross(
    const array_t<array_t<scalar_t, N>, 1> &a,
    const array_t<array_t<scalar_t, N>, 1> &b) {
  return algebra::generic::math::cross(a, b);
}

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
DETRAY_HOST_DEVICE constexpr scalar_t dot(const array_t<scalar_t, N> &a,
                                          const array_t<scalar_t, N> &b) {
  array_t<scalar_t, N> tmp{};
  for (index_t i = 0; i < N; i++) {
    tmp[i] = a[i] * b[i];
  }
  scalar_t ret{0.f};
  for (index_t i = 0; i < N; i++) {
    ret += tmp[i];
  }
  return ret;
}
 
template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
DETRAY_HOST_DEVICE constexpr scalar_t dot(
    const array_t<scalar_t, N> &a, const array_t<array_t<scalar_t, N>, 1> &b) {
  array_t<scalar_t, N> tmp{};
  for (index_t i = 0; i < N; i++) {
    tmp[i] = a[i] * b[0][i];
  }
  scalar_t ret{0.f};
  for (index_t i = 0; i < N; i++) {
    ret += tmp[i];
  }
  return ret;
}
 
template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
DETRAY_HOST_DEVICE constexpr scalar_t dot(
    const array_t<array_t<scalar_t, N>, 1> &a, const array_t<scalar_t, N> &b) {
  array_t<scalar_t, N> tmp{};
  for (index_t i = 0; i < N; i++) {
    tmp[i] = a[0][i] * b[i];
  }
  scalar_t ret{0.f};
  for (index_t i = 0; i < N; i++) {
    ret += tmp[i];
  }
  return ret;
}
 
template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
DETRAY_HOST_DEVICE constexpr scalar_t dot(
    const array_t<array_t<scalar_t, N>, 1> &a,
    const array_t<array_t<scalar_t, N>, 1> &b) {
  array_t<scalar_t, N> tmp{};
  for (index_t i = 0; i < N; i++) {
    tmp[i] = a[0][i] * b[0][i];
  }
  scalar_t ret{0.f};
  for (index_t i = 0; i < N; i++) {
    ret += tmp[i];
  }
  return ret;
}

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 2)
DETRAY_HOST_DEVICE constexpr scalar_t norm(const array_t<scalar_t, N> &v) {
  return algebra::math::sqrt(dot(v, v));
}

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
  requires(N >= 3)
DETRAY_HOST_DEVICE constexpr scalar_t eta(
    const array_t<scalar_t, N> &v) noexcept {
  return algebra::math::atanh(v[2] / norm(v));
}

template <concepts::index index_t, template <typename, index_t> class array_t,
          concepts::scalar scalar_t, index_t N>
DETRAY_HOST_DEVICE constexpr array_t<scalar_t, N> normalize(
    const array_t<scalar_t, N> &v) {
  return (static_cast<scalar_t>(1.) / norm(v)) * v;
}
 
}  // namespace detray::algebra::array