Line data    Source code

       1             : #ifndef CRPROPA_PROPAGATIONBP_H
       2             : #define CRPROPA_PROPAGATIONBP_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 PropagationBP
      17             :  @brief Propagation through magnetic fields using the Boris 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 Boris push integration method.\n
      21             :  It can be used with a fixed step size or an adaptive version which supports the step size control.
      22             :  The step size control tries to keep the relative error close to, but smaller than the designated tolerance.
      23             :  Additionally a minimum and maximum size for the steps can be set.
      24             :  For neutral particles a rectilinear propagation is applied and a next step of the maximum step size proposed.
      25             :  */
      26             : class PropagationBP: public Module {
      27             : 
      28             : public:
      29          84 :         class Y {
      30             :         public:
      31             :                 Vector3d x, u; /*< phase-point: position and direction */
      32             : 
      33             :                 Y() {
      34             :                 }
      35             : 
      36          37 :                 Y(const Vector3d &x, const Vector3d &u) :
      37             :                                 x(x), u(u) {
      38             :                 }
      39             : 
      40             :                 Y(double f) :
      41             :                                 x(Vector3d(f, f, f)), u(Vector3d(f, f, f)) {
      42             :                 }
      43             : 
      44             :                 Y operator *(double f) const {
      45             :                         return Y(x * f, u * f);
      46             :                 }
      47             : 
      48             :                 Y &operator +=(const Y &y) {
      49             :                         x += y.x;
      50             :                         u += y.u;
      51             :                         return *this;
      52             :                 }
      53             :         };
      54             : 
      55             : private:
      56             :         ref_ptr<MagneticField> field;
      57             :         double tolerance; /** target relative error of the numerical integration */
      58             :         double minStep; /** minimum step size of the propagation */
      59             :         double maxStep; /** maximum step size of the propagation */
      60             : 
      61             : public:
      62             :         /** Default constructor for the Boris push. It is constructed with a fixed step size.
      63             :          * @param field
      64             :          * @param fixedStep 
      65             :          */
      66             :         PropagationBP(ref_ptr<MagneticField> field = NULL, double fixedStep = 1. * kpc);
      67             : 
      68             :         /** Constructor for the adaptive Boris push.
      69             :          * @param field
      70             :          * @param tolerance      tolerance is criterion for step adjustment. Step adjustment takes place only if minStep < maxStep
      71             :          * @param minStep          minStep/c_light is the minimum integration time step
      72             :          * @param maxStep          maxStep/c_light is the maximum integration time step. 
      73             :          */
      74             :     PropagationBP(ref_ptr<MagneticField> field, double tolerance, double minStep, double maxStep);
      75             : 
      76             :         /** Propagates the particle. Is called once per iteration.
      77             :          * @param candidate      The Candidate is a passive object, that holds the information about the state of the cosmic ray and the simulation itself. */
      78             :         void process(Candidate *candidate) const;
      79             : 
      80             :         /** Calculates the new position and direction of the particle based on the solution of the Lorentz force
      81             :          * @param pos   current position of the candidate
      82             :          * @param dir   current direction of the candidate
      83             :          * @param step  current step size of the candidate
      84             :          * @param z             current redshift is needed to calculate the magnetic field
      85             :          * @param q             current charge of the candidate
      86             :          * @param m             current mass of the candidate
      87             :          * @return        return the new calculated position and direction of the candidate 
      88             :          */
      89             :         Y dY(Vector3d  pos, Vector3d  dir, double step, double z, double q, double m) const;
      90             : 
      91             :         /** comparison of the position after one step with the position after two steps with step/2.
      92             :          * @param x1    position after one step of size step
      93             :          * @param x2    position after two steps of size step/2
      94             :          * @param step  current step size
      95             :          * @return        measurement of the error of the step 
      96             :          */
      97             :         double errorEstimation(const Vector3d x1, const Vector3d x2, double step) const;
      98             : 
      99             :         /** Get magnetic field vector at current candidate position
     100             :          * @param pos   current position of the candidate
     101             :          * @param z      current redshift is needed to calculate the magnetic field
     102             :          * @return        magnetic field vector at the position pos 
     103             :          */
     104             :         Vector3d getFieldAtPosition(Vector3d pos, double z) const;
     105             : 
     106             :         /** Adapt step size if required and calculates the new position and direction of the particle with the usage of the function dY
     107             :          * @param y              current position and direction of candidate
     108             :          * @param out      position and direction of candidate after the step
     109             :          * @param error  error for the current step
     110             :          * @param h              current step size
     111             :          * @param p              current particle state
     112             :          * @param z              current red shift
     113             :          * @param q              current charge of the candidate 
     114             :          * @param m              current mass of the candidate
     115             :          */
     116             :         void tryStep(const Y &y, Y &out, Y &error, double h, ParticleState &p, double z, double q, double m) const;
     117             : 
     118             :         /** Set functions for the parameters of the class PropagationBP */
     119             : 
     120             :         /** Set a specific magnetic field
     121             :          * @param field  specific magnetic field 
     122             :          */
     123             :         void setField(ref_ptr<MagneticField> field);
     124             :         /** Set a specific tolerance for the step size adaption
     125             :          * @param tolerance      tolerance is criterion for step adjustment. Step adjustment takes place only if minStep < maxStep. 
     126             :          */
     127             :         void setTolerance(double tolerance);
     128             :         /** Set the minimum step for the Boris push
     129             :          * @param minStep          minStep/c_light is the minimum integration time step 
     130             :          */
     131             :         void setMinimumStep(double minStep);
     132             :         /** Set the maximum step for the Boris push
     133             :          * @param maxStep          maxStep/c_light is the maximum integration time step 
     134             :          */
     135             :         void setMaximumStep(double maxStep);
     136             : 
     137             :         /** Get functions for the parameters of the class PropagationBP, similar to the set functions */
     138             : 
     139             :         ref_ptr<MagneticField> getField() const;
     140             :         double getTolerance() const;
     141             :         double getMinimumStep() const;
     142             :         double getMaximumStep() const;
     143             :         std::string getDescription() const;
     144             : };
     145             : /** @}*/
     146             : 
     147             : } // namespace crpropa
     148             : 
     149             : #endif // PROPAGATIONBP_H