Line data Source code
1 : #ifndef CRPROPA_ACCELERATION
2 : #define CRPROPA_ACCELERATION
3 :
4 : #include <crpropa/Candidate.h>
5 : #include <crpropa/Geometry.h>
6 : #include <crpropa/Module.h>
7 : #include <crpropa/Units.h>
8 : #include <crpropa/Vector3.h>
9 :
10 : #include <string>
11 :
12 : namespace crpropa {
13 : /** @addtogroup Acceleration
14 : * @{
15 : */
16 :
17 : /// @class StepLengthModifier
18 : /// @brief Modifies the steplength of an acceleration module.
19 : class StepLengthModifier : public Referenced {
20 : public:
21 : /// Returns an update of the steplength
22 : /// @param steplength Modifies step length, e.g., based on scattering
23 : /// model.
24 : /// @param candidate Additional candidate properties are usually
25 : /// included in the calculation of the updated
26 : /// step length.
27 : virtual double modify(double steplength, Candidate *candidate) = 0;
28 : };
29 :
30 :
31 : /// @class AbstractAccelerationModule
32 : /// @brief Core functionallity for acceleration by scattering with scatter
33 : /// centers moving in a velocity field.
34 : /// @details The velocity field is implicity implemented in the derived classes
35 : /// for performance reasons. Models for the dependence of the step length of
36 : /// the scatter process are set via modifiers.
37 0 : class AbstractAccelerationModule : public Module {
38 : double stepLength;
39 : std::vector<ref_ptr<StepLengthModifier>> modifiers;
40 :
41 : public:
42 : /// The parent's constructor need to be called on initialization!
43 : AbstractAccelerationModule(double _stepLength = 1. * parsec);
44 : // add a step length modifier to the model
45 : void add(StepLengthModifier *modifier);
46 : // update the candidate
47 : void process(Candidate *candidate) const;
48 :
49 : /// Returns the velocity vector of the scatter centers in the rest frame of
50 : /// the candidate. Needs to be implemented in inheriting classes.
51 : virtual Vector3d scatterCenterVelocity(Candidate *candidate) const = 0;
52 :
53 : /// Scatter the candidate with a center with given scatter center
54 : /// velocity into a random direction. Assumes that the
55 : /// candidate is ultra-relativistic (m = 0).
56 : void scatter(Candidate *candidate,
57 : const Vector3d &scatter_center_velocity) const;
58 : };
59 :
60 :
61 : /// @class SecondOrderFermi
62 : /// @brief Implements scattering with centers moving in isotropic directions.
63 : /// All scatter centers have the same velocity.
64 : class SecondOrderFermi : public AbstractAccelerationModule {
65 : double scatterVelocity;
66 : std::vector<double> angle;
67 : std::vector<double> angleCDF;
68 :
69 : public:
70 : /** Constructor
71 : @param scatterVelocity velocity of scattering centers
72 : @param stepLength average mean free path
73 : @param sizeOfPitchangleTable number of precalculated pitch angles
74 : */
75 : SecondOrderFermi(double scatterVelocity = .1 * crpropa::c_light,
76 : double stepLength = 1. * crpropa::parsec,
77 : unsigned int sizeOfPitchangleTable = 10000);
78 : virtual crpropa::Vector3d
79 : scatterCenterVelocity(crpropa::Candidate *candidate) const;
80 : };
81 :
82 :
83 : /// @class DirectedFlowScattering
84 : /// @brief Scattering in a directed flow of scatter centers.
85 : /// @details Two of these region with different
86 : /// velocities can be used to create first order Fermi scenario.
87 : /// Thanks to Aritra Ghosh, Groningn University, for first work in 2017 on
88 : /// the shock acceleration in CRPropa leading to this module.
89 : class DirectedFlowScattering : public AbstractAccelerationModule {
90 : private:
91 : crpropa::Vector3d __scatterVelocity;
92 :
93 : public:
94 : /** Constructor
95 : * @param scatterCenterVelocity velocity of scattering centers
96 : * @param stepLength average mean free path
97 : */
98 : DirectedFlowScattering(crpropa::Vector3d scatterCenterVelocity,
99 : double stepLength = 1. * parsec);
100 : virtual crpropa::Vector3d
101 : scatterCenterVelocity(crpropa::Candidate *candidate) const;
102 : };
103 :
104 :
105 : /// @class DirectedFlowOfScatterCenters
106 : /// @brief In a directed flow, the step length depend on the direction of the
107 : /// particles as headon collisions are more likely than tail=on collisions -
108 : /// propagating against the flow is harder.
109 : class DirectedFlowOfScatterCenters : public StepLengthModifier {
110 : private:
111 : Vector3d __scatterVelocity;
112 :
113 : public:
114 : /** Constructor
115 : * @param scatterCenterVelocity velocity of scattering centers
116 : */
117 : DirectedFlowOfScatterCenters(const Vector3d &scatterCenterVelocity);
118 : double modify(double steplength, Candidate *candidate);
119 : };
120 :
121 :
122 : /// @class QuasiLinearTheory
123 : /// @brief Scales the steplength according to quasi linear theory.
124 : /// @details Following quasi-linear theory [Schlickeiser1989], the mean free
125 : /// path \\f$\\lambda\\f$ of a
126 : /// particle with energy \\f$E\\f$ and charge \\f$Z\\f$ in a field with turbulence
127 : /// spectrum \\f$\\frac{k}{k_{\\min}}^{-q}\\f$ is
128 : /// \\f$ \\lambda = {\\left(\\frac{B}{\\delta B}\\right)}^2 {\\left(R_G\\;
129 : /// k_{\\min}\\right)}^{1-q} R_G \\equiv \\lambda_0 {\\left( \\frac{E}{1
130 : /// EeV}\\frac{1}{Z} \\right)}^{2-q} \\f$
131 : /// where \\f$R_G = \\frac{E}{B Z}\\f$ is the gyro-radius of the
132 : /// particles.
133 : /// This class implements the rigidity dependent scaling factor used to modify
134 : /// the base step length.
135 : /// @par
136 : /// \b [Schlickeiser1989] R. Schlickeiser, Cosmic-Ray Transport and
137 : /// Acceleration. II. Cosmic Rays in Moving Cold Media with Application to
138 : /// Diffusive Shock Wave Acceleration,
139 : /// The Astrophysical Journal 336 (1989) 264. doi:10.1086/167010.
140 : class QuasiLinearTheory : public StepLengthModifier {
141 : private:
142 : double __referenceEnergy;
143 : double __turbulenceIndex;
144 : double __minimumRigidity;
145 :
146 : public:
147 : /** Constructor
148 : * @param referenecEnergy reference energy - break of power spectrum
149 : * @param turbulenceIndex power law index of the isotropic magnetic
150 : * turbulence power spectrum; default is set
151 : * to Kolmogorov turbulence.
152 : * @param minimumRigidity minimal rigidity
153 : */
154 : QuasiLinearTheory(double referenecEnergy = 1. * EeV,
155 : double turbulenceIndex = 5. / 3,
156 : double minimumRigidity = 0);
157 : double modify(double steplength, Candidate *candidate);
158 : };
159 :
160 :
161 : /// @class ParticleSplitting
162 : /// @brief Implements particle splitting, i.e. inverse thinning, to speed up
163 : /// the simulation.
164 : /// @details After crossing a surface a given number of times, the particle is
165 : /// split to N partilces with weight 1/N. This eases performance constraints in
166 : /// acceleration simulations due to the power law nature of many acceleration
167 : /// mechanisms.
168 : /// Thanks to Matthew Weiss, Penn State University for the first work on this
169 : /// feature in 2017.
170 : class ParticleSplitting : public Module {
171 : int numberSplits;
172 : int crossingThreshold;
173 : double minWeight;
174 : ref_ptr<Surface> surface;
175 : std::string counterid;
176 :
177 : public:
178 : /** Constructor
179 : @param surface The surface to monitor
180 : @param crossingThreshold Number of crossings after which a particle is split
181 : @param numberSplits Number of particles the candidate is split into
182 : @param minWeight Minimum weight to consider. Particles with
183 : a lower weight are not split again.
184 : @param counterid An unique string to identify the particle
185 : property used for counting. Useful if
186 : multiple splitting modules are present.
187 : */
188 : ParticleSplitting(Surface *surface, int crossingThreshold = 50,
189 : int numberSplits = 5, double minWeight = 0.01,
190 : std::string counterid = "ParticleSplittingCounter");
191 :
192 : // update the candidate
193 : void process(Candidate *candidate) const;
194 : };
195 :
196 :
197 : /** @} */ // end of group Acceleration
198 :
199 : } // namespace crpropa
200 :
201 : #endif
|
