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
|
