Line data Source code
1 : #include "crpropa/module/Acceleration.h"
2 : #include <crpropa/Common.h>
3 : #include <crpropa/Random.h>
4 : #include <cmath>
5 :
6 : namespace crpropa {
7 :
8 0 : AbstractAccelerationModule::AbstractAccelerationModule(double _stepLength)
9 0 : : crpropa::Module(), stepLength(_stepLength) {}
10 :
11 :
12 0 : void AbstractAccelerationModule::add(StepLengthModifier *modifier) {
13 0 : modifiers.push_back(modifier);
14 0 : }
15 :
16 :
17 0 : void AbstractAccelerationModule::scatter(
18 : crpropa::Candidate *candidate,
19 : const crpropa::Vector3d &scatter_center_velocity) const {
20 : // particle momentum in lab frame
21 0 : const double E = candidate->current.getEnergy();
22 0 : const crpropa::Vector3d p = candidate->current.getMomentum();
23 :
24 : // transform to rest frame of scatter center (p: prime)
25 0 : const double beta = scatter_center_velocity.getR() / crpropa::c_light;
26 0 : const double gamma = 1. / sqrt(1 - beta * beta);
27 0 : const double Ep = gamma * (E - scatter_center_velocity.dot(p));
28 : const crpropa::Vector3d pp = (p - scatter_center_velocity* E /
29 : (crpropa::c_light * crpropa::c_light)) * gamma;
30 :
31 : // scatter into random direction
32 0 : const crpropa::Vector3d pp_new = crpropa::Random::instance().randVector() * pp.getR();
33 :
34 : // transform back
35 0 : const double E_new = gamma * (Ep + scatter_center_velocity.dot(pp_new));
36 : const crpropa::Vector3d p_new = (pp_new + scatter_center_velocity * Ep /
37 : (crpropa::c_light * crpropa::c_light)) * gamma;
38 :
39 : // update candidate properties
40 0 : candidate->current.setEnergy(E_new);
41 0 : candidate->current.setDirection(p_new / p_new.getR());
42 0 : }
43 :
44 :
45 0 : void AbstractAccelerationModule::process(crpropa::Candidate *candidate) const {
46 0 : double currentStepLength = stepLength;
47 0 : for (auto m : modifiers) {
48 0 : currentStepLength = m->modify(currentStepLength, candidate);
49 : }
50 :
51 0 : double step = candidate->getCurrentStep();
52 0 : while (step > 0) {
53 0 : double randDistance = -1. * log(crpropa::Random::instance().rand()) * currentStepLength;
54 :
55 0 : if (step < randDistance) {
56 0 : candidate->limitNextStep(0.1 * currentStepLength);
57 0 : return;
58 : }
59 0 : scatter(candidate, scatterCenterVelocity(candidate));
60 0 : step -= randDistance;
61 : }
62 : }
63 :
64 :
65 0 : SecondOrderFermi::SecondOrderFermi(double scatterVelocity, double stepLength,
66 0 : unsigned int sizeOfPitchangleTable)
67 : : AbstractAccelerationModule(stepLength),
68 0 : scatterVelocity(scatterVelocity) {
69 0 : setDescription("SecondOrderFermi Acceleration");
70 0 : angle.resize(sizeOfPitchangleTable);
71 0 : angleCDF.resize(sizeOfPitchangleTable);
72 :
73 : // have a discretized table of beamed pitch angles
74 0 : for (unsigned int i =0; i < sizeOfPitchangleTable; i++) {
75 0 : angle[i] = i * M_PI / (sizeOfPitchangleTable-1);
76 0 : angleCDF[i] = (angle[i] +scatterVelocity / crpropa::c_light * sin(angle[i])) / M_PI;
77 : }
78 0 : }
79 :
80 :
81 0 : crpropa::Vector3d SecondOrderFermi::scatterCenterVelocity(crpropa::Candidate *candidate) const
82 : {
83 0 : size_t idx = crpropa::closestIndex(crpropa::Random::instance().rand(), angleCDF);
84 0 : crpropa::Vector3d rv = crpropa::Random::instance().randVector();
85 0 : crpropa::Vector3d rotationAxis = candidate->current.getDirection().cross(rv);
86 :
87 0 : rv = candidate->current.getDirection().getRotated(rotationAxis, M_PI - angle[idx]);
88 0 : return rv * scatterVelocity;
89 : }
90 :
91 :
92 0 : DirectedFlowOfScatterCenters::DirectedFlowOfScatterCenters(
93 0 : const Vector3d &scatterCenterVelocity)
94 0 : : __scatterVelocity(scatterCenterVelocity) {}
95 :
96 :
97 0 : double DirectedFlowOfScatterCenters::modify(double steplength, Candidate* candidate)
98 : {
99 0 : double directionModifier = (-1. * __scatterVelocity.dot(candidate->current.getDirection()) + c_light) / c_light;
100 0 : return steplength / directionModifier;
101 : }
102 :
103 :
104 0 : DirectedFlowScattering::DirectedFlowScattering(
105 0 : crpropa::Vector3d scatterCenterVelocity, double stepLength)
106 : : __scatterVelocity(scatterCenterVelocity),
107 0 : AbstractAccelerationModule(stepLength) {
108 :
109 : // In a directed field of scatter centers, the probability to encounter a
110 : // scatter center depends on the direction of the candidate.
111 0 : StepLengthModifier *mod = new DirectedFlowOfScatterCenters(__scatterVelocity);
112 0 : this->add(mod);
113 0 : }
114 :
115 :
116 0 : crpropa::Vector3d DirectedFlowScattering::scatterCenterVelocity(
117 : crpropa::Candidate *candidate) const { // does not depend on candidate here.
118 0 : return __scatterVelocity;
119 : }
120 :
121 :
122 0 : QuasiLinearTheory::QuasiLinearTheory(double referenecEnergy,
123 : double turbulenceIndex,
124 0 : double minimumRigidity)
125 0 : : __referenceEnergy(referenecEnergy), __turbulenceIndex(turbulenceIndex),
126 0 : __minimumRigidity(minimumRigidity) {}
127 :
128 :
129 0 : double QuasiLinearTheory::modify(double steplength, Candidate* candidate)
130 : {
131 0 : if (candidate->current.getRigidity() < __minimumRigidity)
132 : {
133 0 : return steplength * std::pow(__minimumRigidity /
134 0 : (__referenceEnergy / eV), 2. - __turbulenceIndex);
135 : }
136 : else
137 : {
138 0 : return steplength * std::pow(candidate->current.getRigidity() /
139 0 : (__referenceEnergy / eV), 2. - __turbulenceIndex);
140 : }
141 : }
142 :
143 :
144 0 : ParticleSplitting::ParticleSplitting(Surface *surface, int crossingThreshold,
145 0 : int numberSplits, double minWeight, std::string counterid)
146 0 : : surface(surface), crossingThreshold(crossingThreshold),
147 0 : numberSplits(numberSplits), minWeight(minWeight), counterid(counterid){};
148 :
149 0 : void ParticleSplitting::process(Candidate *candidate) const {
150 : const double currentDistance =
151 0 : surface->distance(candidate->current.getPosition());
152 : const double previousDistance =
153 0 : surface->distance(candidate->previous.getPosition());
154 :
155 0 : if (currentDistance * previousDistance > 0)
156 : // candidate remains on the same side
157 : return;
158 :
159 0 : if (candidate->getWeight() < minWeight)
160 : return;
161 :
162 : int num_crossings = 1;
163 0 : if (candidate->hasProperty(counterid))
164 0 : num_crossings = candidate->getProperty(counterid).toInt32() + 1;
165 0 : candidate->setProperty(counterid, num_crossings);
166 :
167 0 : if (num_crossings % crossingThreshold != 0)
168 : return;
169 :
170 0 : candidate->updateWeight(1. / numberSplits);
171 :
172 0 : for (size_t i = 1; i < numberSplits; i++) {
173 : // No recursive split as the weights of the secondaries created
174 : // before the split are not affected
175 0 : ref_ptr<Candidate> new_candidate = candidate->clone(false);
176 0 : new_candidate->parent = candidate;
177 0 : uint64_t snr = Candidate::getNextSerialNumber();
178 0 : Candidate::setNextSerialNumber(snr + 1);
179 0 : new_candidate->setSerialNumber(snr);
180 : candidate->addSecondary(new_candidate);
181 : }
182 : };
183 :
184 : } // namespace crpropa
|
