BFieldCond.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
*/
 
 
#include "MagFieldElements/BFieldCond.h"
#include "CxxUtils/inline_hints.h"
/*
 * Compute the Biot-Savart field due to this conductor
 * If deriv[9] is passed, also computes the field derivatives
 * The results are _added_ to B[] and deriv[].
 */
 
ATH_FLATTEN
// We compile this package with optimization, even in debug builds; otherwise,
// the heavy use of Eigen makes it too slow.  However, from here we may call
// to out-of-line Eigen code that is linked from other DSOs; in that case,
// it would not be optimized.  Avoid this by forcing all Eigen code
// to be inlined here if possible.
void
BFieldCond::addBiotSavart(double scaleFactor,
                          const double* ATH_RESTRICT xyz,
                          double* ATH_RESTRICT B_in,
                          double* ATH_RESTRICT deriv) const
{
  static const double mu04pi = 1.0e-7;    // mu_0/4pi
  static const double minvsq = 10. * 10.; // (1 cm)^2
 
  const Eigen::Map<const Eigen::Vector3d> pos(xyz);
  Eigen::Map<Eigen::Vector3d> B(B_in);
 
  const Eigen::Vector3d r1 = pos - m_p1;
  const Eigen::Vector3d r2 = pos - m_p2;
  const Eigen::Vector3d v = m_u.cross(r1);
  const double vsq = std::max(v.squaredNorm(), minvsq);
 
  const double r1u = r1.dot(m_u);
 
  const double r2u = r2.dot(m_u);
  const double r1n = r1.norm();
 
  const double r2n = r2.norm();
 
  const double f0 = mu04pi * m_curr * scaleFactor / vsq;
  const double f1 = f0 * (m_finite ? r1u / r1n - r2u / r2n : 2.0);
  const double f2 = 2.0 * f1 / vsq;
 
  B += f1 * v;
 
  if (deriv) {
    Eigen::Map<Eigen::Matrix<double, 3, 3, Eigen::RowMajor>> D(deriv);
 
    D(0, 1) -= f1 * m_u(2);
    D(0, 2) += f1 * m_u(1);
    D(1, 0) += f1 * m_u(2);
    D(1, 2) -= f1 * m_u(0);
    D(2, 0) -= f1 * m_u(1);
    D(2, 1) += f1 * m_u(0);
 
    if (vsq > minvsq) {
      Eigen::Vector3d w = f2 * (m_u * r1u - r1);
 
      if (m_finite) {
        // Finite conductor segment
        w += f0 * ((m_u - r1 * r1u / (r1n * r1n)) / r1n -
                   (m_u - r2 * r2u / (r2n * r2n)) / r2n);
      }
 
      D += v * w.transpose();
    }
  }
}