Line data    Source code

       1             : #ifndef CRPROPA_PROPAGATIONCK_H
       2             : #define CRPROPA_PROPAGATIONCK_H
       3             : 
       4             : #include "crpropa/Module.h"
       5             : #include "crpropa/Units.h"
       6             : #include "crpropa/magneticField/MagneticField.h"
       7             : #include "kiss/logger.h"
       8             : 
       9             : namespace crpropa {
      10             : /**
      11             :  * \addtogroup Propagation 
      12             :  * @{
      13             :  */
      14             : 
      15             : /**
      16             :  @class PropagationCK
      17             :  @brief Rectilinear propagation through magnetic fields using the Cash-Karp method.
      18             : 
      19             :  This module solves the equations of motion of a relativistic charged particle when propagating through a magnetic field.\n
      20             :  It uses the Runge-Kutta integration method with Cash-Karp coefficients.\n
      21             :  The step size control tries to keep the relative error close to, but smaller than the designated tolerance.
      22             :  Additionally a minimum and maximum size for the steps can be set.
      23             :  For neutral particles a rectilinear propagation is applied and a next step of the maximum step size proposed.
      24             :  */
      25             : class PropagationCK: public Module {
      26             : public:
      27         286 :         class Y {
      28             :         public:
      29             :                 Vector3d x, u; /*< phase-point: position and direction */
      30             : 
      31             :                 Y() {
      32             :                 }
      33             : 
      34          34 :                 Y(const Vector3d &x, const Vector3d &u) :
      35             :                                 x(x), u(u) {
      36             :                 }
      37             : 
      38             :                 Y(double f) :
      39             :                                 x(Vector3d(f, f, f)), u(Vector3d(f, f, f)) {
      40             :                 }
      41             : 
      42             :                 Y operator *(double f) const {
      43             :                         return Y(x * f, u * f);
      44             :                 }
      45             : 
      46             :                 Y &operator +=(const Y &y) {
      47             :                         x += y.x;
      48             :                         u += y.u;
      49             :                         return *this;
      50             :                 }
      51             :         };
      52             : 
      53             : private:
      54             :         std::vector<double> a, b, bs; /*< Cash-Karp coefficients */
      55             :         ref_ptr<MagneticField> field;
      56             :         double tolerance; /*< target relative error of the numerical integration */
      57             :         double minStep; /*< minimum step size of the propagation */
      58             :         double maxStep; /*< maximum step size of the propagation */
      59             : 
      60             : public:
      61             :         /** Constructor for the adaptive Kash Carp.
      62             :          * @param field
      63             :          * @param tolerance      tolerance is criterion for step adjustment. Step adjustment takes place only if minStep < maxStep
      64             :          * @param minStep          minStep/c_light is the minimum integration time step
      65             :          * @param maxStep          maxStep/c_light is the maximum integration time step. 
      66             :          */
      67             :     PropagationCK(ref_ptr<MagneticField> field = NULL, double tolerance = 1e-4,
      68             :                         double minStep = (0.1 * kpc), double maxStep = (1 * Gpc));
      69             : 
      70             :         void process(Candidate *candidate) const;
      71             : 
      72             :         // derivative of phase point, dY/dt = d/dt(x, u) = (v, du/dt)
      73             :         // du/dt = q*c^2/E * (u x B)
      74             :         Y dYdt(const Y &y, ParticleState &p, double z) const;
      75             : 
      76             :         void tryStep(const Y &y, Y &out, Y &error, double t,
      77             :                         ParticleState &p, double z) const;
      78             : 
      79             :         void setField(ref_ptr<MagneticField> field);
      80             :         void setTolerance(double tolerance);
      81             :         void setMinimumStep(double minStep);
      82             :         void setMaximumStep(double maxStep);
      83             : 
      84             :          /** get functions for the parameters of the class PropagationCK, similar to the set functions */
      85             :         ref_ptr<MagneticField> getField() const;
      86             :         
      87             :         /** get magnetic field vector at current candidate position
      88             :          * @param pos   current position of the candidate
      89             :          * @param z      current redshift is needed to calculate the magnetic field
      90             :          * @return        magnetic field vector at the position pos */
      91             :         Vector3d getFieldAtPosition(Vector3d pos, double z) const;
      92             : 
      93             :         double getTolerance() const;
      94             :         double getMinimumStep() const;
      95             :         double getMaximumStep() const;
      96             :         std::string getDescription() const;
      97             : };
      98             : /** @}*/
      99             : 
     100             : } // namespace crpropa
     101             : 
     102             : #endif // CRPROPA_PROPAGATIONCK_H