RelativisticChargeInMagneticFieldODE.h

Go to the documentation of this file.

#ifndef _gal_RelativisticChargeInMagneticFieldODE_h_
#define _gal_RelativisticChargeInMagneticFieldODE_h_

#include <GalacticUnits.h>


namespace galactic {

  template<class MagneticFieldFunctor>
  class RelativisticChargeInMagneticFieldODE {
  public:
    RelativisticChargeInMagneticFieldODE(const double charge,
                                         const double energy,
                                         MagneticFieldFunctor& mag)
      : fQcE(charge * lightspeed / energy), fMagField(mag) { }

    template<typename Vector6>
    void
    operator()(const double /*x*/, const Vector6& y, Vector6& dYdX)
      const
    {
      double b[3];
      fMagField(y, b);
      dYdX[0] = y[3];  // u_x
      dYdX[1] = y[4];  // u_y
      dYdX[2] = y[5];  // u_z
      dYdX[3] = fQcE * (y[4]*b[2] - y[5]*b[1]);
      dYdX[4] = fQcE * (y[5]*b[0] - y[3]*b[2]);
      dYdX[5] = fQcE * (y[3]*b[1] - y[4]*b[0]);
    }

    operator unsigned int() const
    { return 6; }

  private:
    const double fQcE;
    MagneticFieldFunctor& fMagField;
  };


}


#endif