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