Helix.h

#ifndef _include_HELIX_H_
#define _include_HELIX_H_

#include <corsika/geometry/BaseTrajectory.h>
#include <corsika/geometry/Point.h>
#include <corsika/geometry/Vector.h>
#include <corsika/units/PhysicalUnits.h>
#include <cmath>

namespace corsika::geometry {
  /*!
   * A Helix is defined by the cyclotron frequency \f$ \omega_c \f$, the initial
   * Point r0 and
   * the velocity vectors \f$ \vec{v}_{\parallel} \f$ and \f$ \vec{v}_{\perp} \f$
   * denoting the projections of the initial velocity \f$ \vec{v}_0 \f$ parallel
   * and perpendicular to the axis \f$ \vec{B} \f$, respectively, i.e.
   * \f{align*}{
        \vec{v}_{\parallel} &= \frac{\vec{v}_0 \cdot \vec{B}}{\vec{B}^2} \vec{B} \\
        \vec{v}_{\perp} &= \vec{v}_0 - \vec{v}_{\parallel}
     \f}
   */

  class Helix : public BaseTrajectory {
    using VelocityVec = Vector<corsika::units::si::SpeedType::dimension_type>;

    Point const r0;
    corsika::units::si::FrequencyType const omegaC;
    VelocityVec const vPar;
    VelocityVec const vPerp, uPerp;

    corsika::units::si::LengthType const radius;

  public:
    Helix(Point const& pR0, corsika::units::si::FrequencyType pOmegaC,
          VelocityVec const& pvPar, VelocityVec const& pvPerp)
        : r0(pR0)
        , omegaC(pOmegaC)
        , vPar(pvPar)
        , vPerp(pvPerp)
        , uPerp(vPerp.cross(vPar.normalized()))
        , radius(pvPar.norm() / abs(pOmegaC)) {}

    Point GetPosition(corsika::units::si::TimeType t) const override {
      return r0 + vPar * t +
             (vPerp * (cos(omegaC * t) - 1) + uPerp * sin(omegaC * t)) / omegaC;
    }

    auto GetRadius() const { return radius; }

    LengthType DistanceBetween(corsika::units::si::TimeType t1,
                               corsika::units::si::TimeType t2) const override {
      return (vPar + vPerp).norm() * (t2 - t1);
    }
  };

} // namespace corsika::geometry

#endif