Line data Source code
1 : /** Unit tests for core features of CRPropa
2 : Candidate
3 : ParticleState
4 : Random
5 : Common functions
6 : */
7 :
8 : #include <complex>
9 :
10 : #include "crpropa/Candidate.h"
11 : #include "crpropa/base64.h"
12 : #include "crpropa/Common.h"
13 : #include "crpropa/Units.h"
14 : #include "crpropa/ParticleID.h"
15 : #include "crpropa/ParticleMass.h"
16 : #include "crpropa/Random.h"
17 : #include "crpropa/Grid.h"
18 : #include "crpropa/GridTools.h"
19 : #include "crpropa/Geometry.h"
20 : #include "crpropa/EmissionMap.h"
21 : #include "crpropa/Vector3.h"
22 :
23 : #include <HepPID/ParticleIDMethods.hh>
24 : #include "gtest/gtest.h"
25 :
26 : namespace crpropa {
27 :
28 2 : TEST(ParticleState, position) {
29 1 : ParticleState particle;
30 : Vector3d v(1, 3, 5);
31 1 : particle.setPosition(v * Mpc);
32 2 : EXPECT_TRUE(particle.getPosition() == v * Mpc);
33 1 : }
34 :
35 2 : TEST(ParticleState, energy) {
36 1 : ParticleState particle;
37 1 : particle.setEnergy(10 * EeV);
38 1 : EXPECT_EQ(particle.getEnergy(), 10 * EeV);
39 1 : }
40 :
41 2 : TEST(ParticleState, direction) {
42 1 : ParticleState particle;
43 : Vector3d v(1, 2, 3);
44 1 : particle.setDirection(v);
45 2 : EXPECT_TRUE(particle.getDirection() == v.getUnitVector());
46 1 : }
47 :
48 2 : TEST(ParticleState, velocity) {
49 1 : ParticleState particle;
50 : Vector3d v(1, 1, 0);
51 1 : particle.setDirection(v);
52 2 : EXPECT_TRUE(particle.getVelocity() == v.getUnitVector() * c_light);
53 1 : }
54 :
55 2 : TEST(ParticleState, momentum) {
56 1 : ParticleState particle;
57 : Vector3d v(0, 1, 0);
58 1 : particle.setDirection(v);
59 1 : particle.setEnergy(100 * EeV);
60 2 : EXPECT_TRUE(particle.getMomentum() == v * (particle.getEnergy() / c_light));
61 1 : }
62 :
63 2 : TEST(ParticleState, id) {
64 1 : ParticleState particle;
65 1 : particle.setId(nucleusId(12, 6));
66 1 : EXPECT_EQ(particle.getId(), 1000060120);
67 1 : }
68 :
69 : #ifndef CRPROPA_TESTS_SKIP_EXCEPTIONS
70 1 : TEST(ParticleState, idException) {
71 1 : EXPECT_THROW(nucleusId(5, 6), std::runtime_error);
72 1 : }
73 : #endif
74 :
75 2 : TEST(ParticleState, Charge) {
76 1 : ParticleState particle;
77 :
78 1 : particle.setId(nucleusId(56, 26)); // iron
79 1 : EXPECT_DOUBLE_EQ(26 * eplus, particle.getCharge());
80 :
81 1 : particle.setId(-nucleusId(56, 26)); // anti-iron
82 1 : EXPECT_DOUBLE_EQ(-26 * eplus, particle.getCharge());
83 :
84 1 : particle.setId(11); // electron
85 1 : EXPECT_DOUBLE_EQ(-1 * eplus, particle.getCharge());
86 :
87 1 : particle.setId(-11); // positron
88 1 : EXPECT_DOUBLE_EQ(1 * eplus, particle.getCharge());
89 :
90 1 : particle.setId(12); // electron neutrino
91 1 : EXPECT_DOUBLE_EQ(0, particle.getCharge());
92 :
93 1 : particle.setId(-12); // electron anti-neutrino
94 1 : EXPECT_DOUBLE_EQ(0, particle.getCharge());
95 1 : }
96 :
97 2 : TEST(ParticleState, Rigidity) {
98 1 : ParticleState particle;
99 :
100 1 : particle.setId(nucleusId(1, 1)); // proton
101 1 : particle.setEnergy(1 * EeV);
102 1 : EXPECT_EQ(particle.getRigidity(), 1e18);
103 1 : }
104 :
105 2 : TEST(ParticleState, Mass) {
106 1 : ParticleState particle;
107 :
108 1 : particle.setId(nucleusId(1, 1)); // proton
109 1 : EXPECT_DOUBLE_EQ(mass_proton, particle.getMass());
110 :
111 1 : particle.setId(nucleusId(1, 0)); // neutron
112 1 : EXPECT_DOUBLE_EQ(mass_neutron, particle.getMass());
113 :
114 1 : int id = nucleusId(56, 26);
115 1 : particle.setId(id); // iron
116 1 : EXPECT_DOUBLE_EQ(nuclearMass(id), particle.getMass());
117 :
118 1 : particle.setId(-id); // anti-iron
119 1 : EXPECT_DOUBLE_EQ(nuclearMass(-id), particle.getMass());
120 :
121 : // approximation for unkown nucleus A * amu - Z * mass_electron
122 : int A = 238; int Z = 92; // Uranium92
123 1 : EXPECT_DOUBLE_EQ(nuclearMass(A, Z), A*amu - Z*mass_electron);
124 1 : }
125 :
126 2 : TEST(ParticleState, lorentzFactor) {
127 1 : ParticleState particle;
128 1 : particle.setId(nucleusId(1, 1));
129 1 : particle.setEnergy(1e12 * eV);
130 1 : EXPECT_DOUBLE_EQ(particle.getLorentzFactor(),
131 : 1e12 * eV / mass_proton / c_squared);
132 1 : }
133 :
134 1 : TEST(ParticleID, nucleusId) {
135 1 : EXPECT_EQ(nucleusId(3,2), 1000020030);
136 1 : }
137 :
138 1 : TEST(ParticleID, chargeNumber) {
139 1 : EXPECT_EQ(chargeNumber(1000020030), 2);
140 1 : }
141 :
142 1 : TEST(ParticleID, massNumber) {
143 1 : EXPECT_EQ(massNumber(2112), 1);
144 1 : EXPECT_EQ(massNumber(1000020030), 3);
145 1 : }
146 :
147 1 : TEST(ParticleID, isNucleus) {
148 1 : EXPECT_TRUE(isNucleus(1000020030));
149 1 : EXPECT_FALSE(isNucleus(11));
150 1 : }
151 :
152 1 : TEST(HepPID, consistencyWithReferenceImplementation) {
153 : // Tests the performance improved version against the default one
154 1 : unsigned long testPID = rand() % 1000000000 + 1000000000;
155 8 : for(size_t i=1; i < 8; i++) {
156 7 : HepPID::location loc = (HepPID::location) i;
157 7 : unsigned short newResult = HepPID::digit(loc, testPID);
158 : //original implementation
159 7 : int numerator = (int) std::pow(10.0,(loc-1));
160 7 : EXPECT_EQ(newResult, (HepPID::abspid(testPID)/numerator)%10);
161 : }
162 1 : }
163 :
164 1 : TEST(HepPID, charge) {
165 1 : EXPECT_DOUBLE_EQ(HepPID::charge(11), -1.);
166 1 : }
167 :
168 1 : TEST(Candidate, currentStep) {
169 1 : Candidate candidate;
170 1 : candidate.setCurrentStep(1 * Mpc);
171 1 : EXPECT_DOUBLE_EQ(candidate.getCurrentStep(), 1 * Mpc);
172 1 : }
173 :
174 1 : TEST(Candidate, limitNextStep) {
175 1 : Candidate candidate;
176 1 : candidate.setNextStep(5 * Mpc);
177 1 : EXPECT_DOUBLE_EQ(candidate.getNextStep(), 5 * Mpc);
178 1 : candidate.limitNextStep(2 * Mpc);
179 1 : EXPECT_DOUBLE_EQ(candidate.getNextStep(), 2 * Mpc);
180 1 : candidate.limitNextStep(3 * Mpc);
181 1 : EXPECT_DOUBLE_EQ(candidate.getNextStep(), 2 * Mpc);
182 1 : }
183 :
184 1 : TEST(Candidate, isActive) {
185 1 : Candidate candidate;
186 1 : EXPECT_TRUE(candidate.isActive());
187 1 : candidate.setActive(false);
188 1 : EXPECT_FALSE(candidate.isActive());
189 1 : }
190 :
191 1 : TEST(Candidate, property) {
192 1 : Candidate candidate;
193 2 : candidate.setProperty("foo", "bar");
194 2 : EXPECT_TRUE(candidate.hasProperty("foo"));
195 2 : std::string value = candidate.getProperty("foo");
196 1 : EXPECT_EQ("bar", value);
197 1 : }
198 :
199 1 : TEST(Candidate, weight) {
200 1 : Candidate candidate;
201 1 : EXPECT_EQ (1., candidate.getWeight());
202 :
203 1 : candidate.setWeight(5.);
204 1 : EXPECT_EQ (5., candidate.getWeight());
205 :
206 1 : candidate.updateWeight(3.);
207 1 : EXPECT_EQ (15., candidate.getWeight());
208 1 : }
209 :
210 1 : TEST(Candidate, addSecondary) {
211 1 : Candidate c;
212 1 : c.setRedshift(5);
213 1 : c.setTrajectoryLength(23);
214 1 : c.setWeight(3.);
215 1 : c.previous.setId(nucleusId(56,26));
216 1 : c.previous.setEnergy(1000);
217 1 : c.previous.setPosition(Vector3d(1,2,3));
218 1 : c.previous.setDirection(Vector3d(0,0,1));
219 :
220 1 : c.addSecondary(nucleusId(1,1), 200);
221 2 : c.addSecondary(nucleusId(1,1), 200, 5.);
222 1 : Candidate s1 = *c.secondaries[0];
223 1 : Candidate s2 = *c.secondaries[1];
224 :
225 1 : EXPECT_EQ(nucleusId(1,1), s1.current.getId());
226 1 : EXPECT_EQ(200, s1.current.getEnergy());
227 1 : EXPECT_EQ(5, s1.getRedshift());
228 1 : EXPECT_EQ(23, s1.getTrajectoryLength());
229 1 : EXPECT_EQ(1000, s1.created.getEnergy());
230 1 : EXPECT_EQ(3., s1.getWeight());
231 2 : EXPECT_TRUE(Vector3d(1,2,3) == s1.created.getPosition());
232 2 : EXPECT_TRUE(Vector3d(0,0,1) == s1.created.getDirection());
233 2 : EXPECT_TRUE(s1.getTagOrigin() == "SEC");
234 :
235 1 : EXPECT_EQ(15., s2.getWeight());
236 1 : }
237 :
238 1 : TEST(Candidate, candidateTag) {
239 1 : Candidate c;
240 :
241 : // test default tag
242 2 : EXPECT_TRUE(c.getTagOrigin() == "PRIM");
243 :
244 : // test setting tag
245 1 : c.setTagOrigin("myTag");
246 2 : EXPECT_TRUE(c.getTagOrigin() == "myTag");
247 1 : }
248 :
249 1 : TEST(Candidate, serialNumber) {
250 1 : Candidate::setNextSerialNumber(42);
251 1 : Candidate c;
252 1 : EXPECT_EQ(43, c.getSourceSerialNumber());
253 1 : }
254 :
255 1 : TEST(common, digit) {
256 1 : EXPECT_EQ(1, digit(1234, 1000));
257 1 : EXPECT_EQ(2, digit(1234, 100));
258 1 : EXPECT_EQ(3, digit(1234, 10));
259 1 : EXPECT_EQ(4, digit(1234, 1));
260 1 : }
261 :
262 1 : TEST(common, interpolate) {
263 : // create vectors x = (0, 0.02, ... 2) and y = 2x + 3 = (3, ... 7)
264 1 : std::vector<double> xD(101), yD(101);
265 102 : for (int i = 0; i <= 100; i++) {
266 101 : xD[i] = i * 0.02;
267 101 : yD[i] = 2 * xD[i] + 3;
268 : }
269 :
270 : // interpolating tabulated values of a linear function should produce exact results
271 1 : Random &R = Random::instance();
272 : double x, ytrue, yinterp;
273 10001 : for (int i = 0; i < 10000; i++) {
274 10000 : x = R.rand() * 2; // random value between 0 and 2
275 10000 : ytrue = 2 * x + 3;
276 10000 : yinterp = interpolate(x, xD, yD);
277 10000 : EXPECT_DOUBLE_EQ(ytrue, yinterp);
278 : }
279 :
280 : // test interpolation in first bin
281 : x = 0.01;
282 : ytrue = 2 * x + 3;
283 1 : yinterp = interpolate(x, xD, yD);
284 1 : EXPECT_DOUBLE_EQ(ytrue, yinterp);
285 :
286 : // test interpolation in last bin
287 : x = 1.99;
288 : ytrue = 2 * x + 3;
289 1 : yinterp = interpolate(x, xD, yD);
290 1 : EXPECT_DOUBLE_EQ(ytrue, yinterp);
291 :
292 : // value out of range, return lower bound
293 1 : EXPECT_EQ(3, interpolate(-0.001, xD, yD));
294 :
295 : // value out of range, return upper bound
296 1 : EXPECT_EQ(7, interpolate(2.001, xD, yD));
297 1 : }
298 :
299 1 : TEST(common, interpolateEquidistant) {
300 1 : std::vector<double> yD(100);
301 101 : for (int i = 0; i < 100; i++) {
302 100 : yD[i] = pow(1 + i * 2. / 99., 2);
303 : }
304 :
305 : // interpolated value should be close to computed
306 1 : double y = interpolateEquidistant(1.5001, 1, 3, yD);
307 1 : EXPECT_NEAR(pow(1.5001, 2), y, 1e-4);
308 :
309 : // value out of range, return lower bound
310 1 : EXPECT_EQ(1, interpolateEquidistant(0.9, 1, 3, yD));
311 :
312 : // value out of range, return lower bound
313 1 : EXPECT_EQ(9, interpolateEquidistant(3.1, 1, 3, yD));
314 1 : }
315 :
316 1 : TEST(common, pow_integer) {
317 1 : EXPECT_EQ(pow_integer<0>(1.23), 1);
318 1 : EXPECT_FLOAT_EQ(pow_integer<1>(1.234), 1.234);
319 1 : EXPECT_FLOAT_EQ(pow_integer<2>(1.234), pow(1.234, 2));
320 1 : EXPECT_FLOAT_EQ(pow_integer<3>(1.234), pow(1.234, 3));
321 1 : }
322 :
323 2 : TEST(common, gaussInt) {
324 9 : EXPECT_NEAR(gaussInt(([](double x){ return x*x; }), 0, 10), 1000/3., 1e-4);
325 9 : EXPECT_NEAR(gaussInt(([](double x){ return sin(x)*sin(x); }), 0, M_PI), M_PI/2., 1e-4);
326 1 : }
327 :
328 1 : TEST(Random, seed) {
329 1 : Random &a = Random::instance();
330 1 : Random &b = Random::instance();
331 :
332 1 : a.seed(42);
333 1 : double r1 = a.rand();
334 :
335 1 : a.seed(42);
336 1 : double r2 = a.rand();
337 :
338 1 : a.seed(42);
339 1 : double r3 = b.rand();
340 :
341 : // seeding should give same random numbers
342 1 : EXPECT_EQ(r1, r2);
343 :
344 : // seeding should work for all instances
345 1 : EXPECT_EQ(r1, r3);
346 1 : }
347 :
348 1 : TEST(Random, bigSeedStorage) {
349 1 : Random a;
350 : std::vector<uint32_t> bigSeed;
351 :
352 : const size_t nComp = 42;
353 : double values[nComp];
354 43 : for (size_t i = 0; i < nComp; i++)
355 : {
356 42 : values[i] = a.rand();
357 : }
358 1 : bigSeed = a.getSeed();
359 1 : Random b;
360 : //b.load(bigSeed);
361 1 : b.seed(&bigSeed[0], bigSeed.size());
362 43 : for (size_t i = 0; i < nComp; i++)
363 : {
364 42 : EXPECT_EQ(values[i], b.rand());
365 : }
366 :
367 1 : a.seed(42);
368 1 : bigSeed = a.getSeed();
369 1 : EXPECT_EQ(bigSeed.size(), 1);
370 1 : EXPECT_EQ(bigSeed[0], 42);
371 1 : b.seed(bigSeed[0]);
372 43 : for (size_t i = 0; i < nComp; i++)
373 : {
374 42 : EXPECT_EQ(a.rand(), b.rand());
375 : }
376 :
377 1 : }
378 :
379 1 : TEST(base64, de_en_coding) {
380 1 : Random a;
381 100 : for (int N=1; N < 100; N++) {
382 : std::vector<uint32_t> data;
383 99 : data.reserve(N);
384 5049 : for (int i =0; i<N; i++)
385 4950 : data.push_back(a.randInt());
386 :
387 99 : std::string encoded_data = Base64::encode((unsigned char*)&data[0], sizeof(data[0]) * data.size() / sizeof(unsigned char));
388 :
389 99 : std::string decoded_data = Base64::decode(encoded_data);
390 99 : size_t S = decoded_data.size() * sizeof(decoded_data[0]) / sizeof(uint32_t);
391 5049 : for (int i=0; i < S; i++) {
392 4950 : EXPECT_EQ(((uint32_t*)decoded_data.c_str())[i], data[i]);
393 : }
394 : }
395 :
396 1 : }
397 :
398 1 : TEST(Random, base64Seed) {
399 :
400 1 : std::string seed = "I1+8ANzXYwAqAAAAAwAAAA==";
401 : std::vector<uint32_t> bigSeed;
402 1 : bigSeed.push_back(12345123);
403 1 : bigSeed.push_back(6543324);
404 1 : bigSeed.push_back(42);
405 1 : bigSeed.push_back(3);
406 1 : Random a, b;
407 1 : a.seed(seed);
408 1 : b.seed(&bigSeed[0], bigSeed.size());
409 :
410 : const size_t nComp = 42;
411 : double values[nComp];
412 43 : for (size_t i = 0; i < nComp; i++)
413 : {
414 42 : EXPECT_EQ(a.rand(), b.rand());
415 : }
416 1 : }
417 :
418 2 : TEST(Grid, PeriodicClamp) {
419 : // Test correct determination of lower and upper neighbor
420 : int lo, hi;
421 :
422 : periodicClamp(23.12, 8, lo, hi);
423 1 : EXPECT_EQ(7, lo);
424 1 : EXPECT_EQ(0, hi);
425 :
426 : periodicClamp(-23.12, 8, lo, hi);
427 1 : EXPECT_EQ(0, lo);
428 1 : EXPECT_EQ(1, hi);
429 1 : }
430 :
431 1 : TEST(Grid, PeriodicBoundary) {
432 : // Test correct determination of periodic continuated index
433 : // periodic indices for n=8 should repeat like ...0123456701234567...
434 : int index;
435 :
436 1 : index = periodicBoundary(0, 8);
437 1 : EXPECT_EQ(0, index);
438 1 : index = periodicBoundary(7, 8);
439 1 : EXPECT_EQ(7, index);
440 1 : index = periodicBoundary(8, 8);
441 1 : EXPECT_EQ(0, index);
442 1 : index = periodicBoundary(9, 8);
443 1 : EXPECT_EQ(1, index);
444 1 : }
445 :
446 1 : TEST(Grid, ReflectiveClamp) {
447 : // Test correct determination of lower and upper neighbor
448 : // reflective indices for n=8 should repeat like ...67765432100123456776...
449 : int lo, hi;
450 : double res;
451 :
452 1 : reflectiveClamp(23.12, 8, lo, hi, res);
453 1 : EXPECT_EQ(7, lo);
454 1 : EXPECT_EQ(7, hi);
455 1 : EXPECT_FLOAT_EQ(7.12, res);
456 :
457 1 : reflectiveClamp(-23.12, 8, lo, hi, res);
458 1 : EXPECT_EQ(6, lo);
459 1 : EXPECT_EQ(7, hi);
460 1 : EXPECT_FLOAT_EQ(6.12, res);
461 1 : }
462 :
463 2 : TEST(Grid, ReflectiveBoundary) {
464 : // Test correct determination of reflected index
465 : // reflective indices for n=8 should repeat like ...67765432100123456776...
466 : int index;
467 :
468 1 : index = reflectiveBoundary(8, 8);
469 1 : EXPECT_EQ(7, index);
470 1 : index = reflectiveBoundary(9, 8);
471 1 : EXPECT_EQ(6, index);
472 1 : index = reflectiveBoundary(0, 8);
473 1 : EXPECT_EQ(0, index);
474 1 : index = reflectiveBoundary(-1, 8);
475 1 : EXPECT_EQ(0, index);
476 1 : index = reflectiveBoundary(-8, 8);
477 1 : EXPECT_EQ(7, index);
478 1 : index = reflectiveBoundary(-9, 8);
479 1 : EXPECT_EQ(7, index);
480 1 : }
481 :
482 1 : TEST(Grid1f, SimpleTest) {
483 : // Test construction and parameters
484 1 : size_t Nx = 5;
485 1 : size_t Ny = 8;
486 1 : size_t Nz = 10;
487 : double spacing = 2.0;
488 : Vector3d origin(1., 2., 3.);
489 :
490 1 : Grid1f grid(origin, Nx, Ny, Nz, spacing);
491 :
492 1 : EXPECT_TRUE(origin == grid.getOrigin());
493 1 : EXPECT_EQ(Nx, grid.getNx());
494 1 : EXPECT_EQ(Ny, grid.getNy());
495 1 : EXPECT_EQ(Nz, grid.getNz());
496 1 : EXPECT_EQ(Vector3d(spacing), grid.getSpacing());
497 1 : EXPECT_EQ(5 * 8 * 10, grid.getGrid().size());
498 :
499 : // Test index handling: get position of grid point (2, 3, 4)
500 1 : size_t some_index = 2 * Ny * Nz + 3 * Nz + 4;
501 : Vector3d some_grid_point = origin + Vector3d(2, 3, 4) * spacing + Vector3d(spacing / 2.);
502 2 : EXPECT_EQ(some_grid_point, grid.positionFromIndex(some_index));
503 :
504 : //Test if value on gridpoint is correctly retrieved
505 1 : grid.get(2, 3, 4) = 7;
506 1 : EXPECT_FLOAT_EQ(7., grid.getGrid()[some_index]);
507 : //trilinear interpolated
508 1 : EXPECT_FLOAT_EQ(7., grid.interpolate(some_grid_point));
509 : //tricubic interpolated
510 : grid.setInterpolationType(TRICUBIC);
511 1 : EXPECT_FLOAT_EQ(7., grid.interpolate(some_grid_point));
512 : //nearest neighbour interpolated
513 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
514 1 : EXPECT_FLOAT_EQ(7., grid.interpolate(some_grid_point));
515 :
516 : //Test if grid is set to zero outside of volume for clipVolume=true
517 : grid.setClipVolume(true);
518 1 : EXPECT_FLOAT_EQ(0, grid.interpolate(Vector3d(100, 0, 12)));
519 1 : }
520 :
521 2 : TEST(Grid1f, GridPropertiesConstructor) {
522 : // Test constructor for vector spacing
523 : size_t Nx = 5;
524 : size_t Ny = 8;
525 : size_t Nz = 10;
526 : Vector3d origin = Vector3d(1., 2., 3.);
527 : Vector3d spacing = Vector3d(1., 5., 3.);
528 : GridProperties properties(origin, Nx, Ny, Nz, spacing);
529 1 : Grid1f grid(properties);
530 :
531 1 : EXPECT_EQ(spacing, grid.getSpacing());
532 :
533 : // Test index handling: get position of grid point (1, 7, 6)
534 : size_t some_index = 1 * Ny * Nz + 7 * Nz + 6;
535 : Vector3d some_grid_point = origin + Vector3d(1, 7, 6) * spacing + spacing / 2.;
536 :
537 : //Test if value on gridpoint is correctly retrieved
538 1 : grid.get(1, 7, 6) = 12;
539 1 : EXPECT_FLOAT_EQ(12., grid.getGrid()[some_index]);
540 : //trilinear interpolated
541 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
542 : //tricubic interpolated
543 : grid.setInterpolationType(TRICUBIC);
544 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
545 : //nearest neighbour interpolated
546 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
547 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
548 1 : }
549 :
550 2 : TEST(Grid1f, TestVectorSpacing) {
551 : // Test constructor for vector spacing
552 : size_t Nx = 5;
553 : size_t Ny = 8;
554 : size_t Nz = 10;
555 : Vector3d origin = Vector3d(1., 2., 3.);
556 : Vector3d spacing = Vector3d(1., 5., 3.);
557 :
558 1 : Grid1f grid(origin, Nx, Ny, Nz, spacing);
559 :
560 1 : EXPECT_EQ(spacing, grid.getSpacing());
561 :
562 : // Test index handling: get position of grid point (1, 7, 6)
563 : size_t some_index = 1 * Ny * Nz + 7 * Nz + 6;
564 : Vector3d some_grid_point = origin + Vector3d(1, 7, 6) * spacing + spacing / 2.;
565 :
566 : //Test if value on gridpoint is correctly retrieved
567 1 : grid.get(1, 7, 6) = 12;
568 1 : EXPECT_FLOAT_EQ(12., grid.getGrid()[some_index]);
569 : //trilinear interpolated
570 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
571 : //tricubic interpolated
572 : grid.setInterpolationType(TRICUBIC);
573 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
574 : //nearest neighbour interpolated
575 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
576 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
577 1 : }
578 :
579 2 : TEST(Grid1f, ClosestValue) {
580 : // Check some closest values / nearest neighbour interpolation
581 1 : Grid1f grid(Vector3d(0.), 2, 2, 2, 1.);
582 1 : grid.get(0, 0, 0) = 1;
583 1 : grid.get(0, 0, 1) = 2;
584 1 : grid.get(0, 1, 0) = 3;
585 1 : grid.get(0, 1, 1) = 4;
586 1 : grid.get(1, 0, 0) = 5;
587 1 : grid.get(1, 0, 1) = 6;
588 1 : grid.get(1, 1, 0) = 7;
589 1 : grid.get(1, 1, 1) = 8;
590 :
591 : // grid points are at 0.5 and 1.5
592 1 : EXPECT_FLOAT_EQ(1, grid.closestValue(Vector3d(0.1, 0.2, 0.6)));
593 1 : EXPECT_FLOAT_EQ(2, grid.closestValue(Vector3d(0.2, 0.1, 1.3)));
594 1 : EXPECT_FLOAT_EQ(3, grid.closestValue(Vector3d(0.3, 1.2, 0.2)));
595 1 : EXPECT_FLOAT_EQ(7, grid.closestValue(Vector3d(1.7, 1.8, 0.4)));
596 :
597 : //Test if grid is set to zero outside of volume for clipVolume=true
598 1 : EXPECT_NE(0, grid.interpolate(Vector3d(0, 0, 12)));
599 : grid.setClipVolume(true);
600 1 : double b = grid.interpolate(Vector3d(0, 0, 12));
601 1 : EXPECT_FLOAT_EQ(0, b);
602 1 : }
603 :
604 2 : TEST(Grid1f, clipVolume) {
605 : // Check volume clipping for gridproperties constructor
606 : size_t N = 2;
607 : Vector3d origin = Vector3d(0.);
608 : double spacing = 2;
609 : GridProperties properties(origin, N, spacing);
610 1 : Grid1f grid(properties);
611 1 : grid.get(0, 0, 0) = 1;
612 1 : grid.get(0, 0, 1) = 2;
613 1 : grid.get(0, 1, 0) = 3;
614 1 : grid.get(0, 1, 1) = 4;
615 1 : grid.get(1, 0, 0) = 5;
616 1 : grid.get(1, 0, 1) = 6;
617 1 : grid.get(1, 1, 0) = 7;
618 1 : grid.get(1, 1, 1) = 8;
619 :
620 : //Test if grid is set to zero outside of volume for clipVolume=true
621 1 : EXPECT_NE(0, grid.interpolate(Vector3d(0, 0, 12)));
622 : grid.setClipVolume(true);
623 1 : double b = grid.interpolate(Vector3d(0, 0, 10));
624 1 : EXPECT_FLOAT_EQ(0, b);
625 1 : }
626 :
627 2 : TEST(Grid3f, Interpolation) {
628 : // Explicitly test trilinear and tricubic interpolation
629 : double spacing = 2.793;
630 : int n = 3;
631 1 : Grid3f grid(Vector3d(0.), n, n, n, spacing);
632 : grid.get(0, 0, 1) = Vector3f(1.7, 0., 0.); // set one value
633 :
634 : Vector3d b;
635 :
636 : //trilinear
637 :
638 : // grid points are at [0.5, 1.5, ...] * spacing
639 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.5) * spacing);
640 1 : EXPECT_FLOAT_EQ(1.7, b.x);
641 :
642 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.4) * spacing);
643 1 : EXPECT_FLOAT_EQ(1.7 * 0.9, b.x);
644 :
645 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.6) * spacing);
646 1 : EXPECT_FLOAT_EQ(1.7 * 0.9, b.x);
647 :
648 1 : b = grid.interpolate(Vector3d(0.5, 0.35, 1.6) * spacing);
649 1 : EXPECT_FLOAT_EQ(1.7 * 0.9 * 0.85, b.x);
650 :
651 1 : b = grid.interpolate(Vector3d(0.5, 2.65, 1.6) * spacing); // using periodic repetition
652 1 : EXPECT_FLOAT_EQ(1.7 * 0.9 * 0.15, b.x);
653 :
654 : //tricubic
655 : #ifdef HAVE_SIMD
656 : grid.setInterpolationType(TRICUBIC);
657 :
658 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.5) * spacing);
659 1 : EXPECT_FLOAT_EQ(1.7, b.x);
660 :
661 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.4) * spacing);
662 1 : EXPECT_FLOAT_EQ(1.66005015373, b.x);
663 :
664 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.6) * spacing);
665 1 : EXPECT_FLOAT_EQ(1.66005003452, b.x);
666 :
667 1 : b = grid.interpolate(Vector3d(0.5, 0.35, 1.6) * spacing);
668 1 : EXPECT_FLOAT_EQ(1.57507634163, b.x);
669 :
670 1 : b = grid.interpolate(Vector3d(0.5, 2.65, 1.6) * spacing);
671 1 : EXPECT_FLOAT_EQ(0.190802007914, b.x);
672 : #endif // HAVE_SIMD
673 1 : }
674 :
675 2 : TEST(VectordGrid, Scale) {
676 : // Test scaling a field
677 1 : ref_ptr<Grid3f> grid = new Grid3f(Vector3d(0.), 3, 1);
678 4 : for (int ix = 0; ix < 3; ix++)
679 12 : for (int iy = 0; iy < 3; iy++)
680 36 : for (int iz = 0; iz < 3; iz++)
681 27 : grid->get(ix, iy, iz) = Vector3f(1, 0, 0);
682 :
683 1 : scaleGrid(grid, 5);
684 4 : for (int ix = 0; ix < 3; ix++)
685 12 : for (int iy = 0; iy < 3; iy++)
686 36 : for (int iz = 0; iz < 3; iz++)
687 27 : EXPECT_FLOAT_EQ(5, grid->interpolate(Vector3d(0.7, 0, 0.1)).x);
688 1 : }
689 :
690 2 : TEST(Grid3f, Periodicity) {
691 : // Test for periodic boundaries: grid(x+a*n) = grid(x)
692 : size_t n = 3;
693 : double spacing = 3;
694 : double size = n * spacing;
695 1 : Grid3f grid(Vector3d(0.), n, spacing);
696 4 : for (int ix = 0; ix < 3; ix++)
697 12 : for (int iy = 0; iy < 3; iy++)
698 36 : for (int iz = 0; iz < 3; iz++)
699 27 : grid.get(ix, iy, iz) = Vector3f(iz + ix, iy * iz, ix - iz * iy);
700 :
701 : Vector3d pos(1.2, 2.3, 0.7);
702 :
703 : //trilinear interpolated
704 1 : Vector3f b = grid.interpolate(pos);
705 1 : Vector3f b2 = grid.interpolate(pos + Vector3d(1, 0, 0) * size);
706 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
707 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
708 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
709 :
710 1 : b2 = grid.interpolate(pos + Vector3d(0, 5, 0) * size);
711 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
712 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
713 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
714 :
715 1 : b2 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
716 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
717 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
718 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
719 :
720 : //tricubic interpolated
721 : #ifdef HAVE_SIMD
722 : grid.setInterpolationType(TRICUBIC);
723 1 : b = grid.interpolate(pos);
724 1 : b2 = grid.interpolate(pos + Vector3d(1, 0, 0) * size);
725 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
726 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
727 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
728 :
729 1 : b2 = grid.interpolate(pos + Vector3d(0, 5, 0) * size);
730 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
731 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
732 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
733 :
734 1 : b2 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
735 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
736 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
737 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
738 : #endif // HAVE_SIMD
739 :
740 : //nearest neighbour interpolated
741 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
742 1 : b = grid.interpolate(pos);
743 1 : b2 = grid.interpolate(pos + Vector3d(1, 0, 0) * size);
744 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
745 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
746 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
747 :
748 1 : b2 = grid.interpolate(pos + Vector3d(0, 5, 0) * size);
749 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
750 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
751 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
752 :
753 1 : b2 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
754 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
755 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
756 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
757 :
758 : //Test if grid is set to zero outside of volume for clipVolume=true
759 : grid.setClipVolume(true);
760 1 : Vector3f b3 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
761 1 : EXPECT_FLOAT_EQ(0., b3.x);
762 1 : EXPECT_FLOAT_EQ(0., b3.y);
763 1 : EXPECT_FLOAT_EQ(0., b3.z);
764 1 : }
765 :
766 2 : TEST(Grid3f, Reflectivity) {
767 : // Test for reflective boundaries: grid(pos) = grid(x+a) = grid(-x-a)
768 : size_t n = 3;
769 : double spacing = 3;
770 : double size = n * spacing;
771 1 : Grid3f grid(Vector3d(0.), n, spacing);
772 : grid.setReflective(true); //set reflective boundary
773 4 : for (int ix = 0; ix < 3; ix++)
774 12 : for (int iy = 0; iy < 3; iy++)
775 36 : for (int iz = 0; iz < 3; iz++)
776 27 : grid.get(ix, iy, iz) = Vector3f(iz + ix, iy * iz, ix - iz * iy);
777 :
778 : Vector3d pos(1.2, 2.3, 0.7);
779 :
780 : //trilinear interpolated
781 1 : Vector3f b = grid.interpolate(pos + Vector3d(1,0,0) * spacing);
782 1 : Vector3f b2 = grid.interpolate(pos *(-1) - Vector3d(1,0,0) * spacing);
783 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
784 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
785 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
786 :
787 1 : b = grid.interpolate(pos + Vector3d(0,5,0) * spacing);
788 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,5,0) * spacing);
789 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
790 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
791 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
792 :
793 1 : b = grid.interpolate(pos + Vector3d(0,0,-2) * spacing);
794 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,0,-2) * spacing);
795 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
796 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
797 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
798 :
799 : //tricubic interpolated
800 : #ifdef HAVE_SIMD
801 : grid.setInterpolationType(TRICUBIC);
802 1 : b = grid.interpolate(pos + Vector3d(1,0,0) * spacing);
803 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(1,0,0) * spacing);
804 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
805 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
806 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
807 :
808 1 : b = grid.interpolate(pos + Vector3d(0,5,0) * spacing);
809 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,5,0) * spacing);
810 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
811 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
812 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
813 :
814 1 : b = grid.interpolate(pos + Vector3d(0,0,-2) * spacing);
815 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,0,-2) * spacing);
816 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
817 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
818 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
819 : #endif //HAVE_SIMD
820 :
821 : //nearest neighbour interpolated
822 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
823 1 : b = grid.interpolate(pos + Vector3d(1,0,0) * spacing);
824 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(1,0,0) * spacing);
825 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
826 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
827 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
828 :
829 1 : b = grid.interpolate(pos + Vector3d(0,5,0) * spacing);
830 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,5,0) * spacing);
831 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
832 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
833 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
834 :
835 1 : b = grid.interpolate(pos + Vector3d(0,0,-2) * spacing);
836 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,0,-2) * spacing);
837 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
838 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
839 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
840 :
841 : //Test if grid is set to zero outside of volume for clipVolume=true
842 : grid.setClipVolume(true);
843 1 : Vector3f b3 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
844 1 : EXPECT_FLOAT_EQ(0., b3.x);
845 1 : EXPECT_FLOAT_EQ(0., b3.y);
846 1 : EXPECT_FLOAT_EQ(0., b3.z);
847 1 : }
848 :
849 2 : TEST(Grid3f, DumpLoad) {
850 : // Dump and load a field grid
851 1 : ref_ptr<Grid3f> grid1 = new Grid3f(Vector3d(0.), 3, 1);
852 1 : ref_ptr<Grid3f> grid2 = new Grid3f(Vector3d(0.), 3, 1);
853 :
854 4 : for (int ix = 0; ix < 3; ix++)
855 12 : for (int iy = 0; iy < 3; iy++)
856 36 : for (int iz = 0; iz < 3; iz++)
857 27 : grid1->get(ix, iy, iz) = Vector3f(1, 2, 3);
858 :
859 2 : dumpGrid(grid1, "testDump.raw");
860 2 : loadGrid(grid2, "testDump.raw");
861 :
862 4 : for (int ix = 0; ix < 3; ix++) {
863 12 : for (int iy = 0; iy < 3; iy++) {
864 36 : for (int iz = 0; iz < 3; iz++) {
865 27 : Vector3f b1 = grid1->get(ix, iy, iz);
866 : Vector3f b2 = grid2->get(ix, iy, iz);
867 27 : EXPECT_FLOAT_EQ(b1.x, b2.x);
868 27 : EXPECT_FLOAT_EQ(b1.y, b2.y);
869 27 : EXPECT_FLOAT_EQ(b1.z, b2.z);
870 : }
871 : }
872 : }
873 1 : }
874 :
875 2 : TEST(Grid3f, DumpLoadTxt) {
876 : // Dump and load a field grid
877 1 : ref_ptr<Grid3f> grid1 = new Grid3f(Vector3d(0.), 3, 1);
878 1 : ref_ptr<Grid3f> grid2 = new Grid3f(Vector3d(0.), 3, 1);
879 :
880 4 : for (int ix = 0; ix < 3; ix++)
881 12 : for (int iy = 0; iy < 3; iy++)
882 36 : for (int iz = 0; iz < 3; iz++)
883 27 : grid1->get(ix, iy, iz) = Vector3f(ix, iy, iz);
884 :
885 2 : dumpGridToTxt(grid1, "testDump.txt", 1e4);
886 2 : loadGridFromTxt(grid2, "testDump.txt", 1e-4);
887 :
888 4 : for (int ix = 0; ix < 3; ix++) {
889 12 : for (int iy = 0; iy < 3; iy++) {
890 36 : for (int iz = 0; iz < 3; iz++) {
891 27 : Vector3f b1 = grid1->get(ix, iy, iz);
892 : Vector3f b2 = grid2->get(ix, iy, iz);
893 27 : EXPECT_FLOAT_EQ(b1.x, b2.x);
894 27 : EXPECT_FLOAT_EQ(b1.y, b2.y);
895 27 : EXPECT_FLOAT_EQ(b1.z, b2.z);
896 : }
897 : }
898 : }
899 1 : }
900 :
901 2 : TEST(Grid3f, Speed) {
902 : // Dump and load a field grid
903 1 : Grid3f grid(Vector3d(0.), 3, 3);
904 4 : for (int ix = 0; ix < 3; ix++)
905 12 : for (int iy = 0; iy < 3; iy++)
906 36 : for (int iz = 0; iz < 3; iz++)
907 27 : grid.get(ix, iy, iz) = Vector3f(1, 2, 3);
908 :
909 : Vector3d b;
910 100001 : for (int i = 0; i < 100000; i++)
911 100000 : b = grid.interpolate(Vector3d(i));
912 1 : }
913 :
914 2 : TEST(CylindricalProjectionMap, functions) {
915 : Vector3d v;
916 1 : v.setRThetaPhi(1.0, 1.2, 2.4);
917 1 : EXPECT_NEAR(v.getPhi(), 2.4, .00001);
918 1 : EXPECT_NEAR(v.getTheta(), 1.2, .000001);
919 :
920 :
921 :
922 2 : CylindricalProjectionMap cpm(24, 12);
923 1 : size_t bin = 50;
924 1 : Vector3d d = cpm.directionFromBin(bin);
925 1 : size_t bin2 = cpm.binFromDirection(d);
926 1 : EXPECT_EQ(bin, bin2);
927 1 : }
928 :
929 1 : TEST(EmissionMap, functions) {
930 :
931 1 : EmissionMap em(360, 180, 100, 1 * EeV, 100 * EeV);
932 1 : double e = em.energyFromBin(50);
933 1 : size_t b = em.binFromEnergy(50 * EeV);
934 :
935 : Vector3d d(1.0, 0.0, 0.0);
936 :
937 1 : em.fillMap(1, 50 * EeV, d);
938 :
939 : Vector3d d2;
940 :
941 1 : bool r = em.drawDirection(1, 50 * EeV, d2);
942 1 : EXPECT_TRUE(r);
943 1 : EXPECT_TRUE(d.getAngleTo(d2) < (2. * M_PI / 180.));
944 :
945 1 : r = em.drawDirection(1, 30 * EeV, d2);
946 1 : EXPECT_FALSE(r);
947 :
948 1 : r = em.drawDirection(2, 50 * EeV, d2);
949 1 : EXPECT_FALSE(r);
950 :
951 1 : }
952 :
953 1 : TEST(EmissionMap, merge) {
954 1 : EmissionMap em1, em2;
955 1 : em1.fillMap(1, 50 * EeV, Vector3d(1.0, 0.0, 0.0));
956 1 : em2.fillMap(1, 50 * EeV, Vector3d(0.0, 1.0, 0.0));
957 1 : em2.fillMap(2, 50 * EeV, Vector3d(0.0, 1.0, 0.0));
958 :
959 1 : em1.merge(&em2);
960 :
961 1 : EXPECT_EQ(em1.getMaps().size(), 2);
962 :
963 1 : ref_ptr<CylindricalProjectionMap> cpm = em1.getMap(1, 50 * EeV);
964 1 : size_t bin = cpm->binFromDirection(Vector3d(0.0, 1.0, 0.0));
965 1 : EXPECT_TRUE(cpm->getPdf()[bin] > 0);
966 1 : }
967 :
968 1 : TEST(Variant, copyToBuffer) {
969 1 : double a = 23.42;
970 : Variant v(a);
971 : double b;
972 1 : v.copyToBuffer(&b);
973 1 : EXPECT_EQ(a, b);
974 1 : }
975 :
976 1 : TEST(Variant, stringConversion) {
977 1 : Variant v, w;
978 : {
979 1 : int32_t a = 12;
980 1 : v = Variant::fromInt32(a);
981 1 : EXPECT_EQ(a, v.asInt32());
982 :
983 1 : w = Variant::fromString(v.toString(), v.getType());
984 1 : EXPECT_EQ(a, w.asInt32());
985 : }
986 :
987 : {
988 1 : int64_t a = 12;
989 1 : v = Variant::fromInt64(a);
990 1 : EXPECT_EQ(a, v.asInt64());
991 :
992 1 : w = Variant::fromString(v.toString(), v.getType());
993 1 : EXPECT_EQ(a, w.asInt64());
994 : }
995 :
996 : {
997 : Vector3d a(1, 2, 3);
998 : Variant v = Variant::fromVector3d(a);
999 : Vector3d u = v.asVector3d();
1000 1 : EXPECT_EQ(a.getX(), u.getX());
1001 1 : EXPECT_EQ(a.getY(), u.getY());
1002 1 : EXPECT_EQ(a.getZ(), u.getZ());
1003 1 : }
1004 :
1005 : {
1006 1 : std::complex<double> a1(1, 1);
1007 1 : std::complex<double> a2(2, 0);
1008 : Vector3<std::complex<double>> a(a1, a1, a2);
1009 : Variant v = Variant::fromVector3c(a);
1010 : Vector3<std::complex<double>> u = v.asVector3c();
1011 1 : EXPECT_EQ(a1, u.getX());
1012 1 : EXPECT_EQ(a1, u.getY());
1013 1 : EXPECT_EQ(a2, u.getZ());
1014 1 : }
1015 :
1016 : {
1017 : std::complex<double> a(1, 2);
1018 : Variant v = Variant::fromComplexDouble(a);
1019 1 : std::complex<double> u = v.asComplexDouble();
1020 1 : EXPECT_EQ(u.real(), 1);
1021 1 : EXPECT_EQ(u.imag(), 2);
1022 1 : }
1023 :
1024 : {
1025 : std::vector<Variant> a;
1026 1 : a.push_back(Variant::fromDouble(1));
1027 1 : a.push_back(Variant::fromDouble(2));
1028 1 : a.push_back(Variant::fromDouble(3));
1029 1 : a.push_back(Variant::fromDouble(4));
1030 : Variant v = Variant::fromVector(a);
1031 1 : std::vector<Variant> u = v.asVector();
1032 1 : EXPECT_EQ(a[0], Variant::fromDouble(u[0]));
1033 1 : EXPECT_EQ(a[1], Variant::fromDouble(u[1]));
1034 1 : EXPECT_EQ(a[2], Variant::fromDouble(u[2]));
1035 1 : EXPECT_EQ(a[3], Variant::fromDouble(u[3]));
1036 1 : }
1037 :
1038 :
1039 1 : }
1040 :
1041 2 : TEST(Geometry, Plane) {
1042 1 : Plane p(Vector3d(0,0,1), Vector3d(0,0,1));
1043 1 : EXPECT_DOUBLE_EQ(-1., p.distance(Vector3d(0, 0, 0)));
1044 1 : EXPECT_DOUBLE_EQ(9., p.distance(Vector3d(1, 1, 10)));
1045 1 : }
1046 :
1047 2 : TEST(Geometry, Sphere) {
1048 1 : Sphere s(Vector3d(0,0,0), 1.);
1049 1 : EXPECT_DOUBLE_EQ(-1., s.distance(Vector3d(0, 0, 0)));
1050 1 : EXPECT_DOUBLE_EQ(9., s.distance(Vector3d(10, 0, 0)));
1051 1 : }
1052 :
1053 2 : TEST(Geometry, ParaxialBox) {
1054 1 : ParaxialBox b(Vector3d(0,0,0), Vector3d(3,4,5));
1055 1 : EXPECT_NEAR(-.1, b.distance(Vector3d(0.1, 0.1, 0.1)), 1E-10);
1056 1 : EXPECT_NEAR(-.1, b.distance(Vector3d(0.1, 3.8, 0.1)), 1E-10);
1057 1 : EXPECT_NEAR(-.2, b.distance(Vector3d(0.9, 3.8, 0.9)), 1E-10);
1058 1 : EXPECT_NEAR(7., b.distance(Vector3d(10., 0., 0.)), 1E-10);
1059 1 : EXPECT_NEAR(7., b.distance(Vector3d(10., 2., 0.)), 1E-10);
1060 1 : EXPECT_NEAR(8., b.distance(Vector3d(-8., 0., 0.)), 1E-10);
1061 1 : }
1062 :
1063 0 : int main(int argc, char **argv) {
1064 0 : ::testing::InitGoogleTest(&argc, argv);
1065 0 : return RUN_ALL_TESTS();
1066 : }
1067 :
1068 : } // namespace crpropa
|
