Line data Source code
1 : #include "crpropa/module/PropagationCK.h"
2 :
3 : #include <limits>
4 : #include <sstream>
5 : #include <stdexcept>
6 : #include <vector>
7 :
8 : namespace crpropa {
9 :
10 : // Cash-Karp coefficients
11 : const double cash_karp_a[] = {
12 : 0., 0., 0., 0., 0., 0.,
13 : 1. / 5., 0., 0., 0., 0., 0.,
14 : 3. / 40., 9. / 40., 0., 0., 0., 0.,
15 : 3. / 10., -9. / 10., 6. / 5., 0., 0., 0.,
16 : -11. / 54., 5. / 2., -70. / 27., 35. / 27., 0., 0.,
17 : 1631. / 55296., 175. / 512., 575. / 13824., 44275. / 110592., 253. / 4096., 0.
18 : };
19 :
20 : const double cash_karp_b[] = {
21 : 37. / 378., 0, 250. / 621., 125. / 594., 0., 512. / 1771.
22 : };
23 :
24 : const double cash_karp_bs[] = {
25 : 2825. / 27648., 0., 18575. / 48384., 13525. / 55296., 277. / 14336., 1. / 4.
26 : };
27 :
28 42 : void PropagationCK::tryStep(const Y &y, Y &out, Y &error, double h,
29 : ParticleState &particle, double z) const {
30 : std::vector<Y> k;
31 42 : k.reserve(6);
32 :
33 : out = y;
34 : error = Y(0);
35 :
36 : // calculate the sum of b_i * k_i
37 294 : for (size_t i = 0; i < 6; i++) {
38 :
39 : Y y_n = y;
40 882 : for (size_t j = 0; j < i; j++)
41 630 : y_n += k[j] * a[i * 6 + j] * h;
42 :
43 : // update k_i
44 252 : k[i] = dYdt(y_n, particle, z);
45 :
46 252 : out += k[i] * b[i] * h;
47 252 : error += k[i] * (b[i] - bs[i]) * h;
48 : }
49 42 : }
50 :
51 252 : PropagationCK::Y PropagationCK::dYdt(const Y &y, ParticleState &p, double z) const {
52 : // normalize direction vector to prevent numerical losses
53 252 : Vector3d velocity = y.u.getUnitVector() * c_light;
54 :
55 : // get B field at particle position
56 252 : Vector3d B = getFieldAtPosition(y.x, z);
57 :
58 : // Lorentz force: du/dt = q*c/E * (v x B)
59 252 : Vector3d dudt = p.getCharge() * c_light / p.getEnergy() * velocity.cross(B);
60 252 : return Y(velocity, dudt);
61 : }
62 :
63 13 : PropagationCK::PropagationCK(ref_ptr<MagneticField> field, double tolerance,
64 13 : double minStep, double maxStep) :
65 13 : minStep(0) {
66 15 : setField(field);
67 13 : setTolerance(tolerance);
68 12 : setMaximumStep(maxStep);
69 12 : setMinimumStep(minStep);
70 :
71 : // load Cash-Karp coefficients
72 : a.assign(cash_karp_a, cash_karp_a + 36);
73 : b.assign(cash_karp_b, cash_karp_b + 6);
74 : bs.assign(cash_karp_bs, cash_karp_bs + 6);
75 11 : }
76 :
77 34 : void PropagationCK::process(Candidate *candidate) const {
78 : // save the new previous particle state
79 34 : ParticleState ¤t = candidate->current;
80 : candidate->previous = current;
81 :
82 34 : Y yIn(current.getPosition(), current.getDirection());
83 34 : double step = maxStep;
84 :
85 : // rectilinear propagation for neutral particles
86 34 : if (current.getCharge() == 0) {
87 2 : step = clip(candidate->getNextStep(), minStep, maxStep);
88 1 : current.setPosition(yIn.x + yIn.u * step);
89 1 : candidate->setCurrentStep(step);
90 1 : candidate->setNextStep(maxStep);
91 : return;
92 : }
93 :
94 : Y yOut, yErr;
95 33 : double newStep = step;
96 33 : double z = candidate->getRedshift();
97 :
98 :
99 : // if minStep is the same as maxStep the adaptive algorithm with its error
100 : // estimation is not needed and the computation time can be saved:
101 33 : if (minStep == maxStep){
102 10 : tryStep(yIn, yOut, yErr, step / c_light, current, z);
103 : } else {
104 23 : step = clip(candidate->getNextStep(), minStep, maxStep);
105 23 : newStep = step;
106 : double r = 42; // arbitrary value
107 :
108 : // try performing step until the target error (tolerance) or the minimum/maximum step size has been reached
109 : while (true) {
110 32 : tryStep(yIn, yOut, yErr, step / c_light, current, z);
111 32 : r = yErr.u.getR() / tolerance; // ratio of absolute direction error and tolerance
112 32 : if (r > 1) { // large direction error relative to tolerance, try to decrease step size
113 10 : if (step == minStep) // already minimum step size
114 : break;
115 : else {
116 9 : newStep = step * 0.95 * pow(r, -0.2);
117 17 : newStep = std::max(newStep, 0.1 * step); // limit step size decrease
118 9 : newStep = std::max(newStep, minStep); // limit step size to minStep
119 : step = newStep;
120 : }
121 : } else { // small direction error relative to tolerance, try to increase step size
122 22 : if (step != maxStep) { // only update once if maximum step size yet not reached
123 22 : newStep = step * 0.95 * pow(r, -0.2);
124 35 : newStep = std::min(newStep, 5 * step); // limit step size increase
125 22 : newStep = std::min(newStep, maxStep); // limit step size to maxStep
126 : }
127 : break;
128 : }
129 : }
130 : }
131 :
132 33 : current.setPosition(yOut.x);
133 33 : current.setDirection(yOut.u.getUnitVector());
134 33 : candidate->setCurrentStep(step);
135 33 : candidate->setNextStep(newStep);
136 : }
137 :
138 14 : void PropagationCK::setField(ref_ptr<MagneticField> f) {
139 14 : field = f;
140 14 : }
141 :
142 2 : ref_ptr<MagneticField> PropagationCK::getField() const {
143 2 : return field;
144 : }
145 :
146 253 : Vector3d PropagationCK::getFieldAtPosition(Vector3d pos, double z) const {
147 : Vector3d B(0, 0, 0);
148 : try {
149 : // check if field is valid and use the field vector at the
150 : // position pos with the redshift z
151 253 : if (field.valid())
152 253 : B = field->getField(pos, z);
153 0 : } catch (std::exception &e) {
154 0 : KISS_LOG_ERROR << "PropagationCK: Exception in PropagationCK::getFieldAtPosition.\n"
155 0 : << e.what();
156 0 : }
157 253 : return B;
158 : }
159 :
160 18 : void PropagationCK::setTolerance(double tol) {
161 18 : if ((tol > 1) or (tol < 0))
162 2 : throw std::runtime_error(
163 4 : "PropagationCK: target error not in range 0-1");
164 16 : tolerance = tol;
165 16 : }
166 :
167 22 : void PropagationCK::setMinimumStep(double min) {
168 22 : if (min < 0)
169 1 : throw std::runtime_error("PropagationCK: minStep < 0 ");
170 21 : if (min > maxStep)
171 2 : throw std::runtime_error("PropagationCK: minStep > maxStep");
172 19 : minStep = min;
173 19 : }
174 :
175 19 : void PropagationCK::setMaximumStep(double max) {
176 19 : if (max < minStep)
177 1 : throw std::runtime_error("PropagationCK: maxStep < minStep");
178 18 : maxStep = max;
179 18 : }
180 :
181 2 : double PropagationCK::getTolerance() const {
182 2 : return tolerance;
183 : }
184 :
185 2 : double PropagationCK::getMinimumStep() const {
186 2 : return minStep;
187 : }
188 :
189 3 : double PropagationCK::getMaximumStep() const {
190 3 : return maxStep;
191 : }
192 :
193 0 : std::string PropagationCK::getDescription() const {
194 0 : std::stringstream s;
195 0 : s << "Propagation in magnetic fields using the Cash-Karp method.";
196 0 : s << " Target error: " << tolerance;
197 0 : s << ", Minimum Step: " << minStep / kpc << " kpc";
198 0 : s << ", Maximum Step: " << maxStep / kpc << " kpc";
199 0 : return s.str();
200 0 : }
201 :
202 : } // namespace crpropa
|
