Line data Source code
1 : #include "crpropa/Candidate.h"
2 : #include "crpropa/ParticleID.h"
3 : #include "crpropa/module/SimplePropagation.h"
4 : #include "crpropa/module/PropagationBP.h"
5 : #include "crpropa/module/PropagationCK.h"
6 : #include "crpropa/magneticField/turbulentField/PlaneWaveTurbulence.h"
7 :
8 : #include "gtest/gtest.h"
9 :
10 : #include <string>
11 : #include <iostream>
12 :
13 : namespace crpropa {
14 :
15 1 : TEST(testSimplePropagation, step) {
16 1 : double minStep = 20;
17 1 : double maxStep = 100;
18 1 : SimplePropagation propa(minStep, maxStep);
19 :
20 1 : ParticleState p;
21 1 : p.setPosition(Vector3d(0, 0, 0));
22 1 : p.setDirection(Vector3d(0, 1, 0));
23 :
24 1 : Candidate c(p);
25 1 : c.setNextStep(10);
26 :
27 1 : propa.process(&c);
28 :
29 1 : EXPECT_EQ(minStep, c.getCurrentStep());
30 1 : EXPECT_EQ(maxStep, c.getNextStep());
31 2 : EXPECT_EQ(Vector3d(0, 0, 0), c.created.getPosition());
32 2 : EXPECT_EQ(Vector3d(0, 1, 0), c.created.getDirection());
33 2 : EXPECT_EQ(Vector3d(0, 0, 0), c.previous.getPosition());
34 2 : EXPECT_EQ(Vector3d(0, 1, 0), c.previous.getDirection());
35 2 : EXPECT_EQ(Vector3d(0, 20, 0), c.current.getPosition());
36 2 : EXPECT_EQ(Vector3d(0, 1, 0), c.current.getDirection());
37 2 : }
38 :
39 :
40 1 : TEST(testPropagationCK, zeroField) {
41 1 : PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 0)));
42 :
43 : double minStep = 0.1 * kpc;
44 1 : propa.setMinimumStep(minStep);
45 :
46 1 : ParticleState p;
47 1 : p.setId(nucleusId(1, 1));
48 1 : p.setEnergy(100 * EeV);
49 1 : p.setPosition(Vector3d(0, 0, 0));
50 1 : p.setDirection(Vector3d(0, 1, 0));
51 1 : Candidate c(p);
52 1 : c.setNextStep(0);
53 :
54 1 : propa.process(&c);
55 :
56 1 : EXPECT_DOUBLE_EQ(minStep, c.getCurrentStep()); // perform minimum step
57 1 : EXPECT_DOUBLE_EQ(5 * minStep, c.getNextStep()); // acceleration by factor 5
58 1 : }
59 :
60 : #ifndef CRPROPA_TESTS_SKIP_EXCEPTIONS
61 1 : TEST(testPropagationCK, exceptions) {
62 : // minStep should be smaller than maxStep
63 2 : EXPECT_THROW(PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)), 0.42, 10 , 0), std::runtime_error);
64 : // Too large tolerance: tolerance should be between 0 and 1
65 2 : EXPECT_THROW(PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)), 42., 10 * kpc , 20 * kpc), std::runtime_error);
66 :
67 2 : PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
68 :
69 : // set maximum step, so that it can be tested what happens if a larger minStep is set.
70 1 : propa.setMaximumStep(1 * Mpc);
71 :
72 : // this tests _that_ the expected exception is thrown
73 1 : EXPECT_THROW(propa.setTolerance(2.), std::runtime_error);
74 1 : EXPECT_THROW(propa.setMinimumStep(-1.), std::runtime_error);
75 1 : EXPECT_THROW(propa.setMinimumStep(2 * Mpc), std::runtime_error);
76 :
77 : // set minimum step, so that it can be tested what happens if a smaller maxStep is set.
78 1 : propa.setMinimumStep(0.5 * Mpc);
79 :
80 1 : EXPECT_THROW(propa.setMaximumStep(0.1 * Mpc), std::runtime_error);
81 1 : }
82 : #endif
83 :
84 1 : TEST(testPropagationCK, constructor) {
85 : // Test construction and parameters
86 1 : ref_ptr<MagneticField> bField = new UniformMagneticField(Vector3d(0, 0, 1 * nG));
87 :
88 1 : double minStep = 1.;
89 1 : double maxStep = 100.;
90 1 : double tolerance = 0.01;
91 :
92 1 : PropagationCK propa(bField, tolerance, minStep, maxStep);
93 :
94 1 : EXPECT_EQ(minStep, propa.getMinimumStep());
95 1 : EXPECT_EQ(maxStep, propa.getMaximumStep());
96 1 : EXPECT_EQ(tolerance, propa.getTolerance());
97 2 : EXPECT_EQ(bField, propa.getField());
98 :
99 : // Update parameters
100 1 : minStep = 10.;
101 1 : maxStep = 10.;
102 1 : propa.setTolerance(0.0001);
103 1 : bField = new UniformMagneticField(Vector3d(10 * nG, 0, 1 * nG));
104 :
105 1 : propa.setTolerance(tolerance);
106 1 : propa.setMinimumStep(minStep);
107 1 : propa.setMaximumStep(maxStep);
108 1 : propa.setField(bField);
109 :
110 1 : EXPECT_EQ(minStep, propa.getMinimumStep());
111 1 : EXPECT_EQ(maxStep, propa.getMaximumStep());
112 1 : EXPECT_EQ(tolerance, propa.getTolerance());
113 2 : EXPECT_EQ(bField, propa.getField());
114 :
115 : // The propagation should be initialized with the default constructor
116 2 : PropagationCK propaCKField(bField);
117 1 : EXPECT_EQ(propaCKField.getMaximumStep(), 1 * Gpc);
118 2 : }
119 :
120 :
121 : // Test if the step size is reduced correctly if the error is too large with respect to the tolerance: r > 1
122 1 : TEST(testPropagationCK, reduceStep) {
123 1 : PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 100 * nG)));
124 :
125 : double minStep = 0.1 * kpc;
126 : double maxStep = 1 * Gpc;
127 1 : propa.setMinimumStep(minStep);
128 1 : propa.setMaximumStep(maxStep);
129 : // small tolerance leads to large values of r
130 1 : propa.setTolerance(1e-15);
131 :
132 1 : ParticleState p;
133 1 : p.setId(nucleusId(1, 1));
134 1 : p.setEnergy(100 * TeV);
135 1 : p.setPosition(Vector3d(0, 0, 0));
136 1 : p.setDirection(Vector3d(0, 1, 0));
137 1 : Candidate c(p);
138 : // large step leads to large errors and thus in combination with the low tolerance to high values of r
139 1 : c.setNextStep(maxStep);
140 :
141 1 : propa.process(&c);
142 :
143 : // adaptive algorithm should propagate particle with minimum step size due to the low value for the tolerance
144 1 : EXPECT_DOUBLE_EQ(minStep, c.getCurrentStep()); // perform minimum step because of large r due to small tolerance
145 1 : EXPECT_DOUBLE_EQ(minStep, c.getNextStep()); // stay at minimum step because of large r due to small tolerance
146 1 : }
147 :
148 :
149 : // Test if the step size is increased correctly if the error is small with respect to the tolerance: r < 1
150 1 : TEST(testPropagationCK, increaseStep) {
151 1 : PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
152 :
153 : double minStep = 0.001 * pc;
154 : double maxStep = 3.125 * pc;
155 1 : propa.setMinimumStep(minStep);
156 1 : propa.setMaximumStep(maxStep);
157 : // large tolerance leads to small values of r. Consequently, the step size can be increased.
158 1 : propa.setTolerance(0.9);
159 :
160 1 : ParticleState p;
161 1 : p.setId(nucleusId(1, 1));
162 1 : p.setEnergy(100 * EeV);
163 1 : p.setPosition(Vector3d(0, 0, 0));
164 1 : p.setDirection(Vector3d(0, 1, 0));
165 1 : Candidate c(p);
166 :
167 : // each step the step size can be increased by a factor of 5.
168 6 : for (int i = 1; i < 6; i++){
169 5 : propa.process(&c);
170 5 : EXPECT_DOUBLE_EQ(minStep*pow(5, i) / pc, c.getNextStep()/pc);
171 : }
172 : // after 5 steps the maxStep is reached. The current step is, however, less.
173 1 : EXPECT_DOUBLE_EQ(maxStep/pc/5., c.getCurrentStep()/pc);
174 1 : EXPECT_DOUBLE_EQ(maxStep/pc, c.getNextStep()/pc);
175 1 : }
176 :
177 :
178 1 : TEST(testPropagationCK, proton) {
179 1 : PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
180 :
181 : double minStep = 0.1 * kpc;
182 1 : propa.setMinimumStep(minStep);
183 :
184 1 : ParticleState p;
185 1 : p.setId(nucleusId(1, 1));
186 1 : p.setEnergy(100 * EeV);
187 1 : p.setPosition(Vector3d(0, 0, 0));
188 1 : p.setDirection(Vector3d(0, 1, 0));
189 1 : Candidate c(p);
190 1 : c.setNextStep(0);
191 :
192 1 : propa.process(&c);
193 :
194 1 : EXPECT_DOUBLE_EQ(minStep, c.getCurrentStep()); // perform minimum step
195 1 : EXPECT_DOUBLE_EQ(5 * minStep, c.getNextStep()); // acceleration by factor 5
196 1 : }
197 :
198 :
199 : // Test the numerical results for parallel magnetic field lines along the z-axis
200 1 : TEST(testPropagationCK, gyration) {
201 1 : PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
202 :
203 : double step = 10. * Mpc; // gyroradius is 108.1 Mpc
204 1 : propa.setMaximumStep(step);
205 1 : propa.setMinimumStep(step);
206 :
207 :
208 1 : ParticleState p;
209 1 : p.setId(nucleusId(1, 1));
210 1 : p.setEnergy(100 * EeV);
211 1 : p.setPosition(Vector3d(0, 0, 0));
212 1 : p.setDirection(Vector3d(1, 1, 1));
213 1 : Candidate c(p);
214 1 : c.setNextStep(0);
215 1 : propa.process(&c);
216 :
217 1 : double dirX = c.current.getDirection().x;
218 1 : double dirY = c.current.getDirection().y;
219 1 : double dirZ = c.current.getDirection().z;
220 1 : double posZ = c.current.getPosition().z;
221 :
222 : // Test if the analytical solution is achieved for the components of the momentum with the CK method as expected in
223 : // the background magnetic field.
224 : double precision = 1e-7;
225 : double expected = 2 / 3.;
226 1 : EXPECT_NEAR(expected, dirX * dirX + dirY * dirY, expected * precision); // constant momentum in the plane perpendicular to background magnetic field field
227 : expected = 1 / 3.;
228 1 : EXPECT_NEAR(expected, dirZ * dirZ, expected * precision); // constant momentum parallel to the background magnetic field
229 : expected = step * step / 3.;
230 1 : EXPECT_NEAR(expected, posZ * posZ, expected * precision); // constant velocity parallel to the background magnetic field
231 :
232 : // Nine new steps to have finally propagated the particle ten times
233 10 : for (int i = 0; i < 9; i++){
234 9 : propa.process(&c);
235 : }
236 :
237 1 : dirX = c.current.getDirection().x;
238 1 : dirY = c.current.getDirection().y;
239 1 : dirZ = c.current.getDirection().z;
240 1 : posZ = c.current.getPosition().z;
241 :
242 : // Compare the numerical solutions after ten steps with the analytical solution of the trajectories
243 : expected = 2 / 3.;
244 1 : EXPECT_NEAR(expected, dirX * dirX + dirY * dirY, expected * precision); // constant momentum in the plane perpendicular to background magnetic field field
245 : expected = 1 / 3.;
246 1 : EXPECT_NEAR(expected, dirZ * dirZ, expected * precision); // constant momentum parallel to the background magnetic field
247 : expected = 100 * step * step / 3.;
248 1 : EXPECT_NEAR(expected, posZ * posZ, expected * precision); // constant velocity parallel to the background magnetic field
249 1 : }
250 :
251 :
252 1 : TEST(testPropagationCK, neutron) {
253 1 : PropagationCK propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
254 1 : propa.setMinimumStep(1 * kpc);
255 1 : propa.setMaximumStep(42 * Mpc);
256 :
257 1 : ParticleState p;
258 1 : p.setId(nucleusId(1, 0));
259 1 : p.setEnergy(100 * EeV);
260 1 : p.setPosition(Vector3d(0, 0, 0));
261 1 : p.setDirection(Vector3d(0, 1, 0));
262 1 : Candidate c(p);
263 :
264 1 : propa.process(&c);
265 :
266 1 : EXPECT_DOUBLE_EQ(1 * kpc, c.getCurrentStep());
267 1 : EXPECT_DOUBLE_EQ(42 * Mpc, c.getNextStep());
268 2 : EXPECT_EQ(Vector3d(0, 1 * kpc, 0), c.current.getPosition());
269 2 : EXPECT_EQ(Vector3d(0, 1, 0), c.current.getDirection());
270 1 : }
271 :
272 :
273 1 : TEST(testPropagationBP, zeroField) {
274 1 : PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 0)), 1 * kpc);
275 :
276 : double minStep = 0.1 * kpc;
277 1 : propa.setMinimumStep(minStep);
278 1 : propa.setTolerance(0.42);
279 :
280 1 : ParticleState p;
281 1 : p.setId(nucleusId(1, 1));
282 1 : p.setEnergy(100 * EeV);
283 1 : p.setPosition(Vector3d(0, 0, 0));
284 1 : p.setDirection(Vector3d(0, 1, 0));
285 1 : Candidate c(p);
286 1 : c.setNextStep(0);
287 :
288 1 : propa.process(&c);
289 :
290 1 : EXPECT_DOUBLE_EQ(minStep, c.getCurrentStep()); // perform minimum step
291 1 : EXPECT_DOUBLE_EQ(5 * minStep, c.getNextStep()); // acceleration by factor 5
292 1 : }
293 :
294 : #ifndef CRPROPA_TESTS_SKIP_EXCEPTIONS
295 1 : TEST(testPropagationBP, exceptions) {
296 : // minStep should be smaller than maxStep
297 2 : EXPECT_THROW(PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)), 0.42, 10 , 0), std::runtime_error);
298 : // Too large tolerance: tolerance should be between 0 and 1
299 2 : EXPECT_THROW(PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)), 42., 10 * kpc , 20 * kpc), std::runtime_error);
300 :
301 2 : PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
302 :
303 : // set maximum step, so that it can be tested what happens if a larger minStep is set.
304 1 : propa.setMaximumStep(1 * Mpc);
305 :
306 : // this tests _that_ the expected exception is thrown
307 1 : EXPECT_THROW(propa.setTolerance(2.), std::runtime_error);
308 1 : EXPECT_THROW(propa.setMinimumStep(-1.), std::runtime_error);
309 1 : EXPECT_THROW(propa.setMinimumStep(2 * Mpc), std::runtime_error);
310 :
311 : // set minimum step, so that it can be tested what happens if a smaller maxStep is set.
312 1 : propa.setMinimumStep(0.5 * Mpc);
313 :
314 1 : EXPECT_THROW(propa.setMaximumStep(0.1 * Mpc), std::runtime_error);
315 1 : }
316 : #endif
317 :
318 :
319 1 : TEST(testPropagationBP, constructor) {
320 : // Test construction and parameters
321 1 : ref_ptr<MagneticField> bField = new UniformMagneticField(Vector3d(0, 0, 1 * nG));
322 :
323 1 : double minStep = 1.;
324 1 : double maxStep = 100.;
325 1 : double tolerance = 0.01;
326 :
327 1 : PropagationBP propa(bField, tolerance, minStep, maxStep);
328 :
329 1 : EXPECT_EQ(minStep, propa.getMinimumStep());
330 1 : EXPECT_EQ(maxStep, propa.getMaximumStep());
331 1 : EXPECT_EQ(tolerance, propa.getTolerance());
332 2 : EXPECT_EQ(bField, propa.getField());
333 :
334 : // Update parameters
335 1 : minStep = 10.;
336 1 : maxStep = 10.;
337 1 : propa.setTolerance(0.0001);
338 1 : bField = new UniformMagneticField(Vector3d(10 * nG, 0, 1 * nG));
339 :
340 1 : propa.setTolerance(tolerance);
341 1 : propa.setMinimumStep(minStep);
342 1 : propa.setMaximumStep(maxStep);
343 1 : propa.setField(bField);
344 :
345 1 : EXPECT_EQ(minStep, propa.getMinimumStep());
346 1 : EXPECT_EQ(maxStep, propa.getMaximumStep());
347 1 : EXPECT_EQ(tolerance, propa.getTolerance());
348 2 : EXPECT_EQ(bField, propa.getField());
349 :
350 : // Test the fixed step size version of the Boris push
351 1 : minStep = 10. * kpc;
352 2 : PropagationBP propaBP(bField, minStep);
353 1 : EXPECT_EQ(propaBP.getMaximumStep(), propaBP.getMaximumStep());
354 :
355 : // The propagation should be initialized with the default constructor
356 2 : PropagationBP propaBPField(bField);
357 1 : EXPECT_EQ(propaBPField.getMaximumStep(), 1 * kpc);
358 2 : }
359 :
360 :
361 : // Test if the step size is reduced correctly if the error is too large with respect to the tolerance: r > 1
362 1 : TEST(testPropagationBP, reduceStep) {
363 1 : PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 100 * nG)), 1 * kpc);
364 :
365 : double minStep = 0.1 * kpc;
366 : double maxStep = 1 * Gpc;
367 1 : propa.setMinimumStep(minStep);
368 1 : propa.setMaximumStep(maxStep);
369 : // small tolerance leads to large values of r
370 1 : propa.setTolerance(1e-15);
371 :
372 1 : ParticleState p;
373 1 : p.setId(nucleusId(1, 1));
374 1 : p.setEnergy(100 * TeV);
375 1 : p.setPosition(Vector3d(0, 0, 0));
376 1 : p.setDirection(Vector3d(0, 1, 0));
377 1 : Candidate c(p);
378 : // large step leads to large errors and thus in combination with the low tolerance to high values of r
379 1 : c.setNextStep(maxStep);
380 :
381 1 : propa.process(&c);
382 :
383 : // adaptive algorithm should propagate particle with minimum step size due to the low value for the tolerance
384 1 : EXPECT_DOUBLE_EQ(minStep, c.getCurrentStep()); // perform minimum step because of large r due to small tolerance
385 1 : EXPECT_DOUBLE_EQ(minStep, c.getNextStep()); // stay at minimum step because of large r due to small tolerance
386 1 : }
387 :
388 :
389 : // Test if the step size is increased correctly if the error is small with respect to the tolerance: r < 1
390 1 : TEST(testPropagationBP, increaseStep) {
391 1 : PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)), 1 * kpc);
392 :
393 : double minStep = 0.001 * pc;
394 : double maxStep = 3.125 * pc;
395 1 : propa.setMinimumStep(minStep);
396 1 : propa.setMaximumStep(maxStep);
397 : // large tolerance leads to small values of r. Consequently, the step size can be increased.
398 1 : propa.setTolerance(0.9);
399 :
400 1 : ParticleState p;
401 1 : p.setId(nucleusId(1, 1));
402 1 : p.setEnergy(100 * EeV);
403 1 : p.setPosition(Vector3d(0, 0, 0));
404 1 : p.setDirection(Vector3d(0, 1, 0));
405 1 : Candidate c(p);
406 :
407 : // each step the step size can be increased by a factor of 5.
408 6 : for (int i = 1; i < 6; i++){
409 5 : propa.process(&c);
410 5 : EXPECT_DOUBLE_EQ(minStep*pow(5, i) / pc, c.getNextStep()/pc);
411 : }
412 : // after 5 steps the maxStep is reached. The current step is, however, less.
413 1 : EXPECT_DOUBLE_EQ(maxStep/pc/5., c.getCurrentStep()/pc);
414 1 : EXPECT_DOUBLE_EQ(maxStep/pc, c.getNextStep()/pc);
415 1 : }
416 :
417 :
418 1 : TEST(testPropagationBP, proton) {
419 1 : PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
420 :
421 : double step = 0.01 * kpc;
422 1 : propa.setMinimumStep(step);
423 1 : propa.setMaximumStep(10*step);
424 1 : propa.setTolerance(0.00001);
425 :
426 1 : ParticleState p;
427 1 : p.setId(nucleusId(1, 1));
428 1 : p.setEnergy(100 * EeV);
429 1 : p.setPosition(Vector3d(0, 0, 0));
430 1 : p.setDirection(Vector3d(0, 1, 0));
431 1 : Candidate c(p);
432 1 : c.setNextStep(0);
433 :
434 1 : propa.process(&c);
435 :
436 1 : EXPECT_DOUBLE_EQ(step, c.getCurrentStep()); // perform step
437 1 : EXPECT_DOUBLE_EQ(5 * step, c.getNextStep()); // acceleration by factor 5
438 1 : }
439 :
440 :
441 : // Test the numerical results for parallel magnetic field lines along the z-axis
442 1 : TEST(testPropagationBP, gyration) {
443 1 : PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
444 :
445 : double step = 10. * Mpc; // gyroradius is 108.1 Mpc
446 1 : propa.setMaximumStep(step);
447 1 : propa.setMinimumStep(step);
448 :
449 :
450 1 : ParticleState p;
451 1 : p.setId(nucleusId(1, 1));
452 1 : p.setEnergy(100 * EeV);
453 1 : p.setPosition(Vector3d(0, 0, 0));
454 1 : p.setDirection(Vector3d(1, 1, 1));
455 1 : Candidate c(p);
456 1 : c.setNextStep(0);
457 1 : propa.process(&c);
458 :
459 1 : double dirX = c.current.getDirection().x;
460 1 : double dirY = c.current.getDirection().y;
461 1 : double dirZ = c.current.getDirection().z;
462 1 : double posZ = c.current.getPosition().z;
463 :
464 : // Test if the analytical solution is achieved for the components of the momentum with the Boris push as expected in
465 : // the background magnetic field.
466 1 : EXPECT_DOUBLE_EQ(2 / 3., dirX * dirX + dirY * dirY); // constant momentum in the perpendicular plane to background magnetic field field
467 1 : EXPECT_DOUBLE_EQ(1 / 3., dirZ * dirZ); // constant momentum parallel to the background magnetic field
468 1 : EXPECT_DOUBLE_EQ( step * step / 3., posZ * posZ); // constant velocity parallel to the background magnetic field
469 :
470 : // Nine new steps to have finally propagated the particle ten times
471 10 : for (int i = 0; i < 9; i++){
472 9 : propa.process(&c);
473 : }
474 :
475 1 : dirX = c.current.getDirection().x;
476 1 : dirY = c.current.getDirection().y;
477 1 : dirZ = c.current.getDirection().z;
478 1 : posZ = c.current.getPosition().z;
479 :
480 : // Compare the numerical solutions after ten steps with the analytical solution of the trajectories
481 1 : EXPECT_DOUBLE_EQ(2 / 3., dirX * dirX + dirY * dirY); // constant momentum in the perpendicular plane to background magnetic field field
482 1 : EXPECT_DOUBLE_EQ(1 / 3., dirZ * dirZ); // constant momentum parallel to the background magnetic field
483 1 : EXPECT_DOUBLE_EQ(100 * step * step / 3., posZ * posZ); // constant velocity parallel to the background magnetic field
484 1 : }
485 :
486 :
487 : // Test the that the optimization for fixed step sizes works
488 1 : TEST(testPropagationBP, fixedStepOptimization) {
489 : // particle 1 with fixed step sizes
490 : double fixed_step = pc;
491 2 : PropagationBP propa1(new PlaneWaveTurbulence(TurbulenceSpectrum(gauss, pc, 100*pc), 10, 1), fixed_step);
492 1 : ParticleState p1;
493 1 : p1.setId(nucleusId(1, 1));
494 1 : p1.setEnergy(100 * EeV);
495 1 : p1.setPosition(Vector3d(0, 0, 0));
496 1 : p1.setDirection(Vector3d(1, 1, 1));
497 1 : Candidate c1(p1);
498 1 : c1.setNextStep(0);
499 : // Nine new steps to have finally propagated the particle ten times
500 10 : for (int i = 0; i < 9; i++){
501 9 : propa1.process(&c1);
502 : }
503 :
504 : // particle 2 with different min and max steps. The tolerance is chosen such that particle 2 will be
505 : // propagated with the same step as particle 1, however not using the optimization for fixed step sizes
506 : double tolerance = 1;
507 2 : PropagationBP propa2(new PlaneWaveTurbulence(TurbulenceSpectrum(gauss, pc, 100*pc), 10, 1), tolerance, fixed_step, 1.1*fixed_step);
508 1 : ParticleState p2;
509 1 : p2.setId(nucleusId(1, 1));
510 1 : p2.setEnergy(100 * EeV);
511 1 : p2.setPosition(Vector3d(0, 0, 0));
512 1 : p2.setDirection(Vector3d(1, 1, 1));
513 1 : Candidate c2(p2);
514 1 : c1.setNextStep(0);
515 : // Nine new steps to have finally propagated the particle ten times
516 10 : for (int i = 0; i < 9; i++){
517 9 : propa2.process(&c2);
518 : }
519 :
520 1 : EXPECT_DOUBLE_EQ(c1.current.getDirection().x, c2.current.getDirection().x);
521 1 : EXPECT_DOUBLE_EQ(c1.current.getDirection().y, c2.current.getDirection().y);
522 1 : EXPECT_DOUBLE_EQ(c1.current.getDirection().z, c2.current.getDirection().z);
523 1 : EXPECT_DOUBLE_EQ(c1.current.getPosition().x, c2.current.getPosition().x);
524 1 : EXPECT_DOUBLE_EQ(c1.current.getPosition().y, c2.current.getPosition().y);
525 1 : EXPECT_DOUBLE_EQ(c1.current.getPosition().z, c2.current.getPosition().z);
526 1 : }
527 :
528 :
529 1 : TEST(testPropagationBP, neutron) {
530 1 : PropagationBP propa(new UniformMagneticField(Vector3d(0, 0, 1 * nG)));
531 :
532 1 : propa.setMinimumStep(1 * kpc);
533 1 : propa.setMaximumStep(1 * kpc);
534 :
535 1 : ParticleState p;
536 1 : p.setId(nucleusId(1, 0));
537 1 : p.setEnergy(100 * EeV);
538 1 : p.setPosition(Vector3d(0, 0, 0));
539 1 : p.setDirection(Vector3d(0, 1, 0));
540 1 : Candidate c(p);
541 :
542 1 : propa.process(&c);
543 :
544 1 : EXPECT_DOUBLE_EQ(1 * kpc, c.getCurrentStep());
545 1 : EXPECT_DOUBLE_EQ(1 * kpc, c.getNextStep());
546 2 : EXPECT_EQ(Vector3d(0, 1 * kpc, 0), c.current.getPosition());
547 2 : EXPECT_EQ(Vector3d(0, 1, 0), c.current.getDirection());
548 1 : }
549 :
550 :
551 0 : int main(int argc, char **argv) {
552 0 : ::testing::InitGoogleTest(&argc, argv);
553 0 : return RUN_ALL_TESTS();
554 : }
555 :
556 : } // namespace crpropa
|
