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(ParticleMass, particleMass) {
153 : //particleMass(int id) interfaces nuclearMass for nuclei
154 1 : EXPECT_DOUBLE_EQ(nuclearMass(nucleusId(1,1)), particleMass(nucleusId(1,1)));
155 : //particleMass(int id) for electron/positron, photon and neutrino
156 1 : EXPECT_DOUBLE_EQ(mass_electron,particleMass(11));
157 1 : EXPECT_DOUBLE_EQ(mass_electron,particleMass(-11));
158 1 : EXPECT_DOUBLE_EQ(0.0,particleMass(22));
159 1 : EXPECT_DOUBLE_EQ(0.0,particleMass(14));
160 1 : }
161 :
162 1 : TEST(HepPID, consistencyWithReferenceImplementation) {
163 : // Tests the performance improved version against the default one
164 1 : unsigned long testPID = rand() % 1000000000 + 1000000000;
165 8 : for(size_t i=1; i < 8; i++) {
166 7 : HepPID::location loc = (HepPID::location) i;
167 7 : unsigned short newResult = HepPID::digit(loc, testPID);
168 : //original implementation
169 7 : int numerator = (int) std::pow(10.0,(loc-1));
170 7 : EXPECT_EQ(newResult, (HepPID::abspid(testPID)/numerator)%10);
171 : }
172 1 : }
173 :
174 1 : TEST(HepPID, charge) {
175 1 : EXPECT_DOUBLE_EQ(HepPID::charge(11), -1.);
176 1 : }
177 :
178 1 : TEST(Candidate, currentStep) {
179 1 : Candidate candidate;
180 1 : EXPECT_DOUBLE_EQ(candidate.getTrajectoryLength(), 0);
181 1 : EXPECT_DOUBLE_EQ(candidate.getTime(), 0);
182 :
183 1 : candidate.setCurrentStep(1 * Mpc);
184 :
185 1 : EXPECT_DOUBLE_EQ(candidate.getCurrentStep(), 1 * Mpc);
186 1 : EXPECT_DOUBLE_EQ(candidate.getTrajectoryLength(), 1 * Mpc);
187 1 : EXPECT_DOUBLE_EQ(candidate.getTime(), 1 * Mpc / c_light);
188 1 : }
189 :
190 1 : TEST(Candidate, limitNextStep) {
191 1 : Candidate candidate;
192 1 : candidate.setNextStep(5 * Mpc);
193 1 : EXPECT_DOUBLE_EQ(candidate.getNextStep(), 5 * Mpc);
194 1 : candidate.limitNextStep(2 * Mpc);
195 1 : EXPECT_DOUBLE_EQ(candidate.getNextStep(), 2 * Mpc);
196 1 : candidate.limitNextStep(3 * Mpc);
197 1 : EXPECT_DOUBLE_EQ(candidate.getNextStep(), 2 * Mpc);
198 1 : }
199 :
200 1 : TEST(Candidate, isActive) {
201 1 : Candidate candidate;
202 1 : EXPECT_TRUE(candidate.isActive());
203 1 : candidate.setActive(false);
204 1 : EXPECT_FALSE(candidate.isActive());
205 1 : }
206 :
207 1 : TEST(Candidate, property) {
208 1 : Candidate candidate;
209 2 : candidate.setProperty("foo", "bar");
210 2 : EXPECT_TRUE(candidate.hasProperty("foo"));
211 2 : std::string value = candidate.getProperty("foo");
212 1 : EXPECT_EQ("bar", value);
213 1 : }
214 :
215 1 : TEST(Candidate, weight) {
216 1 : Candidate candidate;
217 1 : EXPECT_EQ (1., candidate.getWeight());
218 :
219 1 : candidate.setWeight(5.);
220 1 : EXPECT_EQ (5., candidate.getWeight());
221 :
222 1 : candidate.updateWeight(3.);
223 1 : EXPECT_EQ (15., candidate.getWeight());
224 1 : }
225 :
226 1 : TEST(Candidate, addSecondary) {
227 1 : Candidate c;
228 1 : c.setRedshift(5);
229 1 : c.setTrajectoryLength(23);
230 1 : c.setTime(12);
231 1 : c.setWeight(3.);
232 1 : c.previous.setId(nucleusId(56,26));
233 1 : c.previous.setEnergy(1000);
234 1 : c.previous.setPosition(Vector3d(1,2,3));
235 1 : c.previous.setDirection(Vector3d(0,0,1));
236 :
237 1 : c.addSecondary(nucleusId(1,1), 200);
238 1 : c.addSecondary(nucleusId(1,1), 200, 5.);
239 1 : c.addSecondary(11, 200);
240 1 : c.addSecondary(14, 200);
241 2 : c.addSecondary(22, 200);
242 1 : Candidate s1 = *c.secondaries[0]; //proton
243 1 : Candidate s2 = *c.secondaries[1]; //proton
244 1 : Candidate s3 = *c.secondaries[2]; //electron
245 1 : Candidate s4 = *c.secondaries[3]; //neutrino
246 1 : Candidate s5 = *c.secondaries[4]; //photon
247 :
248 1 : EXPECT_EQ(nucleusId(1,1), s1.current.getId());
249 1 : EXPECT_EQ(200, s1.current.getEnergy());
250 1 : EXPECT_EQ(5, s1.getRedshift());
251 1 : EXPECT_EQ(23, s1.getTrajectoryLength());
252 1 : EXPECT_EQ(12, s1.getTime());
253 1 : EXPECT_EQ(1000, s1.created.getEnergy());
254 1 : EXPECT_EQ(3., s1.getWeight());
255 2 : EXPECT_TRUE(Vector3d(1,2,3) == s1.created.getPosition());
256 2 : EXPECT_TRUE(Vector3d(0,0,1) == s1.created.getDirection());
257 1 : EXPECT_TRUE(s1.getTagOrigin() == "SEC");
258 1 : EXPECT_EQ(mass_electron,s3.current.getMass());
259 1 : EXPECT_EQ(0.0,s4.current.getMass());
260 1 : EXPECT_EQ(0.0,s5.current.getMass());
261 :
262 1 : EXPECT_EQ(15., s2.getWeight());
263 1 : }
264 :
265 1 : TEST(Candidate, candidateTag) {
266 1 : Candidate c;
267 :
268 : // test default tag
269 1 : EXPECT_TRUE(c.getTagOrigin() == "PRIM");
270 :
271 : // test setting tag
272 1 : c.setTagOrigin("myTag");
273 1 : EXPECT_TRUE(c.getTagOrigin() == "myTag");
274 1 : }
275 :
276 1 : TEST(Candidate, serialNumber) {
277 1 : Candidate::setNextSerialNumber(42);
278 1 : Candidate c;
279 1 : EXPECT_EQ(43, c.getSourceSerialNumber());
280 1 : }
281 :
282 1 : TEST(common, digit) {
283 1 : EXPECT_EQ(1, digit(1234, 1000));
284 1 : EXPECT_EQ(2, digit(1234, 100));
285 1 : EXPECT_EQ(3, digit(1234, 10));
286 1 : EXPECT_EQ(4, digit(1234, 1));
287 1 : }
288 :
289 1 : TEST(common, interpolate) {
290 : // create vectors x = (0, 0.02, ... 2) and y = 2x + 3 = (3, ... 7)
291 1 : std::vector<double> xD(101), yD(101);
292 102 : for (int i = 0; i <= 100; i++) {
293 101 : xD[i] = i * 0.02;
294 101 : yD[i] = 2 * xD[i] + 3;
295 : }
296 :
297 : // interpolating tabulated values of a linear function should produce exact results
298 1 : Random &R = Random::instance();
299 : double x, ytrue, yinterp;
300 10001 : for (int i = 0; i < 10000; i++) {
301 10000 : x = R.rand() * 2; // random value between 0 and 2
302 10000 : ytrue = 2 * x + 3;
303 10000 : yinterp = interpolate(x, xD, yD);
304 10000 : EXPECT_DOUBLE_EQ(ytrue, yinterp);
305 : }
306 :
307 : // test interpolation in first bin
308 : x = 0.01;
309 : ytrue = 2 * x + 3;
310 1 : yinterp = interpolate(x, xD, yD);
311 1 : EXPECT_DOUBLE_EQ(ytrue, yinterp);
312 :
313 : // test interpolation in last bin
314 : x = 1.99;
315 : ytrue = 2 * x + 3;
316 1 : yinterp = interpolate(x, xD, yD);
317 1 : EXPECT_DOUBLE_EQ(ytrue, yinterp);
318 :
319 : // value out of range, return lower bound
320 1 : EXPECT_EQ(3, interpolate(-0.001, xD, yD));
321 :
322 : // value out of range, return upper bound
323 1 : EXPECT_EQ(7, interpolate(2.001, xD, yD));
324 1 : }
325 :
326 1 : TEST(common, interpolateEquidistant) {
327 1 : std::vector<double> yD(100);
328 101 : for (int i = 0; i < 100; i++) {
329 100 : yD[i] = pow(1 + i * 2. / 99., 2);
330 : }
331 :
332 : // interpolated value should be close to computed
333 1 : double y = interpolateEquidistant(1.5001, 1, 3, yD);
334 1 : EXPECT_NEAR(pow(1.5001, 2), y, 1e-4);
335 :
336 : // value out of range, return lower bound
337 1 : EXPECT_EQ(1, interpolateEquidistant(0.9, 1, 3, yD));
338 :
339 : // value out of range, return lower bound
340 1 : EXPECT_EQ(9, interpolateEquidistant(3.1, 1, 3, yD));
341 1 : }
342 :
343 1 : TEST(common, pow_integer) {
344 1 : EXPECT_EQ(pow_integer<0>(1.23), 1);
345 1 : EXPECT_FLOAT_EQ(pow_integer<1>(1.234), 1.234);
346 1 : EXPECT_FLOAT_EQ(pow_integer<2>(1.234), pow(1.234, 2));
347 1 : EXPECT_FLOAT_EQ(pow_integer<3>(1.234), pow(1.234, 3));
348 1 : }
349 :
350 2 : TEST(common, gaussInt) {
351 9 : EXPECT_NEAR(gaussInt(([](double x){ return x*x; }), 0, 10), 1000/3., 1e-4);
352 9 : EXPECT_NEAR(gaussInt(([](double x){ return sin(x)*sin(x); }), 0, M_PI), M_PI/2., 1e-4);
353 1 : }
354 :
355 1 : TEST(Random, seed) {
356 1 : Random &a = Random::instance();
357 1 : Random &b = Random::instance();
358 :
359 1 : a.seed(42);
360 1 : double r1 = a.rand();
361 :
362 1 : a.seed(42);
363 1 : double r2 = a.rand();
364 :
365 1 : a.seed(42);
366 1 : double r3 = b.rand();
367 :
368 : // seeding should give same random numbers
369 1 : EXPECT_EQ(r1, r2);
370 :
371 : // seeding should work for all instances
372 1 : EXPECT_EQ(r1, r3);
373 1 : }
374 :
375 1 : TEST(Random, bigSeedStorage) {
376 1 : Random a;
377 : std::vector<uint32_t> bigSeed;
378 :
379 : const size_t nComp = 42;
380 : double values[nComp];
381 43 : for (size_t i = 0; i < nComp; i++)
382 : {
383 42 : values[i] = a.rand();
384 : }
385 1 : bigSeed = a.getSeed();
386 1 : Random b;
387 : //b.load(bigSeed);
388 1 : b.seed(&bigSeed[0], bigSeed.size());
389 43 : for (size_t i = 0; i < nComp; i++)
390 : {
391 42 : EXPECT_EQ(values[i], b.rand());
392 : }
393 :
394 1 : a.seed(42);
395 1 : bigSeed = a.getSeed();
396 1 : EXPECT_EQ(bigSeed.size(), 1);
397 1 : EXPECT_EQ(bigSeed[0], 42);
398 1 : b.seed(bigSeed[0]);
399 43 : for (size_t i = 0; i < nComp; i++)
400 : {
401 42 : EXPECT_EQ(a.rand(), b.rand());
402 : }
403 :
404 2 : }
405 :
406 1 : TEST(base64, de_en_coding) {
407 1 : Random a;
408 100 : for (int N=1; N < 100; N++) {
409 : std::vector<uint32_t> data;
410 99 : data.reserve(N);
411 5049 : for (int i =0; i<N; i++)
412 4950 : data.push_back(a.randInt());
413 :
414 99 : std::string encoded_data = Base64::encode((unsigned char*)&data[0], sizeof(data[0]) * data.size() / sizeof(unsigned char));
415 :
416 99 : std::string decoded_data = Base64::decode(encoded_data);
417 99 : size_t S = decoded_data.size() * sizeof(decoded_data[0]) / sizeof(uint32_t);
418 5049 : for (int i=0; i < S; i++) {
419 4950 : EXPECT_EQ(((uint32_t*)decoded_data.c_str())[i], data[i]);
420 : }
421 99 : }
422 :
423 1 : }
424 :
425 1 : TEST(Random, base64Seed) {
426 :
427 1 : std::string seed = "I1+8ANzXYwAqAAAAAwAAAA==";
428 : std::vector<uint32_t> bigSeed;
429 1 : bigSeed.push_back(12345123);
430 1 : bigSeed.push_back(6543324);
431 1 : bigSeed.push_back(42);
432 1 : bigSeed.push_back(3);
433 1 : Random a, b;
434 1 : a.seed(seed);
435 1 : b.seed(&bigSeed[0], bigSeed.size());
436 :
437 : const size_t nComp = 42;
438 : double values[nComp];
439 43 : for (size_t i = 0; i < nComp; i++)
440 : {
441 42 : EXPECT_EQ(a.rand(), b.rand());
442 : }
443 2 : }
444 :
445 2 : TEST(Grid, PeriodicClamp) {
446 : // Test correct determination of lower and upper neighbor
447 : int lo, hi;
448 :
449 : periodicClamp(23.12, 8, lo, hi);
450 1 : EXPECT_EQ(7, lo);
451 1 : EXPECT_EQ(0, hi);
452 :
453 : periodicClamp(-23.12, 8, lo, hi);
454 1 : EXPECT_EQ(0, lo);
455 1 : EXPECT_EQ(1, hi);
456 1 : }
457 :
458 1 : TEST(Grid, PeriodicBoundary) {
459 : // Test correct determination of periodic continuated index
460 : // periodic indices for n=8 should repeat like ...0123456701234567...
461 : int index;
462 :
463 1 : index = periodicBoundary(0, 8);
464 1 : EXPECT_EQ(0, index);
465 1 : index = periodicBoundary(7, 8);
466 1 : EXPECT_EQ(7, index);
467 1 : index = periodicBoundary(8, 8);
468 1 : EXPECT_EQ(0, index);
469 1 : index = periodicBoundary(9, 8);
470 1 : EXPECT_EQ(1, index);
471 1 : }
472 :
473 1 : TEST(Grid, ReflectiveClamp) {
474 : // Test correct determination of lower and upper neighbor
475 : // reflective indices for n=8 should repeat like ...67765432100123456776...
476 : int lo, hi;
477 : double res;
478 :
479 1 : reflectiveClamp(23.12, 8, lo, hi, res);
480 1 : EXPECT_EQ(7, lo);
481 1 : EXPECT_EQ(7, hi);
482 1 : EXPECT_FLOAT_EQ(7.12, res);
483 :
484 1 : reflectiveClamp(-23.12, 8, lo, hi, res);
485 1 : EXPECT_EQ(6, lo);
486 1 : EXPECT_EQ(7, hi);
487 1 : EXPECT_FLOAT_EQ(6.12, res);
488 1 : }
489 :
490 2 : TEST(Grid, ReflectiveBoundary) {
491 : // Test correct determination of reflected index
492 : // reflective indices for n=8 should repeat like ...67765432100123456776...
493 : int index;
494 :
495 1 : index = reflectiveBoundary(8, 8);
496 1 : EXPECT_EQ(7, index);
497 1 : index = reflectiveBoundary(9, 8);
498 1 : EXPECT_EQ(6, index);
499 1 : index = reflectiveBoundary(0, 8);
500 1 : EXPECT_EQ(0, index);
501 1 : index = reflectiveBoundary(-1, 8);
502 1 : EXPECT_EQ(0, index);
503 1 : index = reflectiveBoundary(-8, 8);
504 1 : EXPECT_EQ(7, index);
505 1 : index = reflectiveBoundary(-9, 8);
506 1 : EXPECT_EQ(7, index);
507 1 : }
508 :
509 1 : TEST(Grid1f, SimpleTest) {
510 : // Test construction and parameters
511 1 : size_t Nx = 5;
512 1 : size_t Ny = 8;
513 1 : size_t Nz = 10;
514 : double spacing = 2.0;
515 : Vector3d origin(1., 2., 3.);
516 :
517 1 : Grid1f grid(origin, Nx, Ny, Nz, spacing);
518 :
519 1 : EXPECT_TRUE(origin == grid.getOrigin());
520 1 : EXPECT_EQ(Nx, grid.getNx());
521 1 : EXPECT_EQ(Ny, grid.getNy());
522 1 : EXPECT_EQ(Nz, grid.getNz());
523 1 : EXPECT_EQ(Vector3d(spacing), grid.getSpacing());
524 1 : EXPECT_EQ(5 * 8 * 10, grid.getGrid().size());
525 :
526 : // Test index handling: get position of grid point (2, 3, 4)
527 : size_t some_index = 2 * Ny * Nz + 3 * Nz + 4;
528 : Vector3d some_grid_point = origin + Vector3d(2, 3, 4) * spacing + Vector3d(spacing / 2.);
529 1 : EXPECT_EQ(some_grid_point, grid.positionFromIndex(some_index));
530 :
531 : //Test if value on gridpoint is correctly retrieved
532 1 : grid.get(2, 3, 4) = 7;
533 1 : EXPECT_FLOAT_EQ(7., grid.getGrid()[some_index]);
534 : //trilinear interpolated
535 1 : EXPECT_FLOAT_EQ(7., grid.interpolate(some_grid_point));
536 : //tricubic interpolated
537 : grid.setInterpolationType(TRICUBIC);
538 1 : EXPECT_FLOAT_EQ(7., grid.interpolate(some_grid_point));
539 : //nearest neighbour interpolated
540 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
541 1 : EXPECT_FLOAT_EQ(7., grid.interpolate(some_grid_point));
542 :
543 : //Test if grid is set to zero outside of volume for clipVolume=true
544 : grid.setClipVolume(true);
545 1 : EXPECT_FLOAT_EQ(0, grid.interpolate(Vector3d(100, 0, 12)));
546 1 : }
547 :
548 2 : TEST(Grid1f, GridPropertiesConstructor) {
549 : // Test constructor for vector spacing
550 : size_t Nx = 5;
551 : size_t Ny = 8;
552 : size_t Nz = 10;
553 : Vector3d origin = Vector3d(1., 2., 3.);
554 : Vector3d spacing = Vector3d(1., 5., 3.);
555 : GridProperties properties(origin, Nx, Ny, Nz, spacing);
556 1 : Grid1f grid(properties);
557 :
558 1 : EXPECT_EQ(spacing, grid.getSpacing());
559 :
560 : // Test index handling: get position of grid point (1, 7, 6)
561 : size_t some_index = 1 * Ny * Nz + 7 * Nz + 6;
562 : Vector3d some_grid_point = origin + Vector3d(1, 7, 6) * spacing + spacing / 2.;
563 :
564 : //Test if value on gridpoint is correctly retrieved
565 1 : grid.get(1, 7, 6) = 12;
566 1 : EXPECT_FLOAT_EQ(12., grid.getGrid()[some_index]);
567 : //trilinear interpolated
568 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
569 : //tricubic interpolated
570 : grid.setInterpolationType(TRICUBIC);
571 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
572 : //nearest neighbour interpolated
573 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
574 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
575 1 : }
576 :
577 2 : TEST(Grid1f, TestVectorSpacing) {
578 : // Test constructor for vector spacing
579 : size_t Nx = 5;
580 : size_t Ny = 8;
581 : size_t Nz = 10;
582 : Vector3d origin = Vector3d(1., 2., 3.);
583 : Vector3d spacing = Vector3d(1., 5., 3.);
584 :
585 1 : Grid1f grid(origin, Nx, Ny, Nz, spacing);
586 :
587 1 : EXPECT_EQ(spacing, grid.getSpacing());
588 :
589 : // Test index handling: get position of grid point (1, 7, 6)
590 : size_t some_index = 1 * Ny * Nz + 7 * Nz + 6;
591 : Vector3d some_grid_point = origin + Vector3d(1, 7, 6) * spacing + spacing / 2.;
592 :
593 : //Test if value on gridpoint is correctly retrieved
594 1 : grid.get(1, 7, 6) = 12;
595 1 : EXPECT_FLOAT_EQ(12., grid.getGrid()[some_index]);
596 : //trilinear interpolated
597 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
598 : //tricubic interpolated
599 : grid.setInterpolationType(TRICUBIC);
600 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
601 : //nearest neighbour interpolated
602 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
603 1 : EXPECT_FLOAT_EQ(12., grid.interpolate(some_grid_point));
604 1 : }
605 :
606 2 : TEST(Grid1f, ClosestValue) {
607 : // Check some closest values / nearest neighbour interpolation
608 1 : Grid1f grid(Vector3d(0.), 2, 2, 2, 1.);
609 1 : grid.get(0, 0, 0) = 1;
610 1 : grid.get(0, 0, 1) = 2;
611 1 : grid.get(0, 1, 0) = 3;
612 1 : grid.get(0, 1, 1) = 4;
613 1 : grid.get(1, 0, 0) = 5;
614 1 : grid.get(1, 0, 1) = 6;
615 1 : grid.get(1, 1, 0) = 7;
616 1 : grid.get(1, 1, 1) = 8;
617 :
618 : // grid points are at 0.5 and 1.5
619 1 : EXPECT_FLOAT_EQ(1, grid.closestValue(Vector3d(0.1, 0.2, 0.6)));
620 1 : EXPECT_FLOAT_EQ(2, grid.closestValue(Vector3d(0.2, 0.1, 1.3)));
621 1 : EXPECT_FLOAT_EQ(3, grid.closestValue(Vector3d(0.3, 1.2, 0.2)));
622 1 : EXPECT_FLOAT_EQ(7, grid.closestValue(Vector3d(1.7, 1.8, 0.4)));
623 :
624 : //Test if grid is set to zero outside of volume for clipVolume=true
625 1 : EXPECT_NE(0, grid.interpolate(Vector3d(0, 0, 12)));
626 : grid.setClipVolume(true);
627 1 : double b = grid.interpolate(Vector3d(0, 0, 12));
628 1 : EXPECT_FLOAT_EQ(0, b);
629 1 : }
630 :
631 2 : TEST(Grid1f, clipVolume) {
632 : // Check volume clipping for gridproperties constructor
633 : size_t N = 2;
634 : Vector3d origin = Vector3d(0.);
635 : double spacing = 2;
636 : GridProperties properties(origin, N, spacing);
637 1 : Grid1f grid(properties);
638 1 : grid.get(0, 0, 0) = 1;
639 1 : grid.get(0, 0, 1) = 2;
640 1 : grid.get(0, 1, 0) = 3;
641 1 : grid.get(0, 1, 1) = 4;
642 1 : grid.get(1, 0, 0) = 5;
643 1 : grid.get(1, 0, 1) = 6;
644 1 : grid.get(1, 1, 0) = 7;
645 1 : grid.get(1, 1, 1) = 8;
646 :
647 : //Test if grid is set to zero outside of volume for clipVolume=true
648 1 : EXPECT_NE(0, grid.interpolate(Vector3d(0, 0, 12)));
649 : grid.setClipVolume(true);
650 1 : double b = grid.interpolate(Vector3d(0, 0, 10));
651 1 : EXPECT_FLOAT_EQ(0, b);
652 1 : }
653 :
654 2 : TEST(Grid3f, Interpolation) {
655 : // Explicitly test trilinear and tricubic interpolation
656 : double spacing = 2.793;
657 : int n = 3;
658 1 : Grid3f grid(Vector3d(0.), n, n, n, spacing);
659 : grid.get(0, 0, 1) = Vector3f(1.7, 0., 0.); // set one value
660 :
661 : Vector3d b;
662 :
663 : //trilinear
664 :
665 : // grid points are at [0.5, 1.5, ...] * spacing
666 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.5) * spacing);
667 1 : EXPECT_FLOAT_EQ(1.7, b.x);
668 :
669 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.4) * spacing);
670 1 : EXPECT_FLOAT_EQ(1.7 * 0.9, b.x);
671 :
672 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.6) * spacing);
673 1 : EXPECT_FLOAT_EQ(1.7 * 0.9, b.x);
674 :
675 1 : b = grid.interpolate(Vector3d(0.5, 0.35, 1.6) * spacing);
676 1 : EXPECT_FLOAT_EQ(1.7 * 0.9 * 0.85, b.x);
677 :
678 1 : b = grid.interpolate(Vector3d(0.5, 2.65, 1.6) * spacing); // using periodic repetition
679 1 : EXPECT_FLOAT_EQ(1.7 * 0.9 * 0.15, b.x);
680 :
681 : //tricubic
682 : #ifdef HAVE_SIMD
683 : grid.setInterpolationType(TRICUBIC);
684 :
685 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.5) * spacing);
686 1 : EXPECT_FLOAT_EQ(1.7, b.x);
687 :
688 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.4) * spacing);
689 1 : EXPECT_FLOAT_EQ(1.66005015373, b.x);
690 :
691 1 : b = grid.interpolate(Vector3d(0.5, 0.5, 1.6) * spacing);
692 1 : EXPECT_FLOAT_EQ(1.66005003452, b.x);
693 :
694 1 : b = grid.interpolate(Vector3d(0.5, 0.35, 1.6) * spacing);
695 1 : EXPECT_FLOAT_EQ(1.57507634163, b.x);
696 :
697 1 : b = grid.interpolate(Vector3d(0.5, 2.65, 1.6) * spacing);
698 1 : EXPECT_FLOAT_EQ(0.190802007914, b.x);
699 : #endif // HAVE_SIMD
700 1 : }
701 :
702 1 : TEST(VectordGrid, Scale) {
703 : // Test scaling a field
704 1 : ref_ptr<Grid3f> grid = new Grid3f(Vector3d(0.), 3, 1);
705 4 : for (int ix = 0; ix < 3; ix++)
706 12 : for (int iy = 0; iy < 3; iy++)
707 36 : for (int iz = 0; iz < 3; iz++)
708 27 : grid->get(ix, iy, iz) = Vector3f(1, 0, 0);
709 :
710 1 : scaleGrid(grid, 5);
711 4 : for (int ix = 0; ix < 3; ix++)
712 12 : for (int iy = 0; iy < 3; iy++)
713 36 : for (int iz = 0; iz < 3; iz++)
714 27 : EXPECT_FLOAT_EQ(5, grid->interpolate(Vector3d(0.7, 0, 0.1)).x);
715 1 : }
716 :
717 2 : TEST(Grid3f, Periodicity) {
718 : // Test for periodic boundaries: grid(x+a*n) = grid(x)
719 : size_t n = 3;
720 : double spacing = 3;
721 : double size = n * spacing;
722 1 : Grid3f grid(Vector3d(0.), n, spacing);
723 4 : for (int ix = 0; ix < 3; ix++)
724 12 : for (int iy = 0; iy < 3; iy++)
725 36 : for (int iz = 0; iz < 3; iz++)
726 27 : grid.get(ix, iy, iz) = Vector3f(iz + ix, iy * iz, ix - iz * iy);
727 :
728 : Vector3d pos(1.2, 2.3, 0.7);
729 :
730 : //trilinear interpolated
731 1 : Vector3f b = grid.interpolate(pos);
732 1 : Vector3f b2 = grid.interpolate(pos + Vector3d(1, 0, 0) * size);
733 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
734 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
735 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
736 :
737 1 : b2 = grid.interpolate(pos + Vector3d(0, 5, 0) * size);
738 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
739 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
740 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
741 :
742 1 : b2 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
743 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
744 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
745 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
746 :
747 : //tricubic interpolated
748 : #ifdef HAVE_SIMD
749 : grid.setInterpolationType(TRICUBIC);
750 1 : b = grid.interpolate(pos);
751 1 : b2 = grid.interpolate(pos + Vector3d(1, 0, 0) * size);
752 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
753 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
754 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
755 :
756 1 : b2 = grid.interpolate(pos + Vector3d(0, 5, 0) * size);
757 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
758 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
759 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
760 :
761 1 : b2 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
762 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
763 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
764 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
765 : #endif // HAVE_SIMD
766 :
767 : //nearest neighbour interpolated
768 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
769 1 : b = grid.interpolate(pos);
770 1 : b2 = grid.interpolate(pos + Vector3d(1, 0, 0) * size);
771 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
772 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
773 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
774 :
775 1 : b2 = grid.interpolate(pos + Vector3d(0, 5, 0) * size);
776 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
777 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
778 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
779 :
780 1 : b2 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
781 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
782 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
783 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
784 :
785 : //Test if grid is set to zero outside of volume for clipVolume=true
786 : grid.setClipVolume(true);
787 1 : Vector3f b3 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
788 1 : EXPECT_FLOAT_EQ(0., b3.x);
789 1 : EXPECT_FLOAT_EQ(0., b3.y);
790 1 : EXPECT_FLOAT_EQ(0., b3.z);
791 1 : }
792 :
793 2 : TEST(Grid3f, Reflectivity) {
794 : // Test for reflective boundaries: grid(pos) = grid(x+a) = grid(-x-a)
795 : size_t n = 3;
796 : double spacing = 3;
797 : double size = n * spacing;
798 1 : Grid3f grid(Vector3d(0.), n, spacing);
799 : grid.setReflective(true); //set reflective boundary
800 4 : for (int ix = 0; ix < 3; ix++)
801 12 : for (int iy = 0; iy < 3; iy++)
802 36 : for (int iz = 0; iz < 3; iz++)
803 27 : grid.get(ix, iy, iz) = Vector3f(iz + ix, iy * iz, ix - iz * iy);
804 :
805 : Vector3d pos(1.2, 2.3, 0.7);
806 :
807 : //trilinear interpolated
808 1 : Vector3f b = grid.interpolate(pos + Vector3d(1,0,0) * spacing);
809 1 : Vector3f b2 = grid.interpolate(pos *(-1) - Vector3d(1,0,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,5,0) * spacing);
815 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,5,0) * 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 :
820 1 : b = grid.interpolate(pos + Vector3d(0,0,-2) * spacing);
821 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,0,-2) * spacing);
822 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
823 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
824 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
825 :
826 : //tricubic interpolated
827 : #ifdef HAVE_SIMD
828 : grid.setInterpolationType(TRICUBIC);
829 1 : b = grid.interpolate(pos + Vector3d(1,0,0) * spacing);
830 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(1,0,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,5,0) * spacing);
836 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,5,0) * 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 1 : b = grid.interpolate(pos + Vector3d(0,0,-2) * spacing);
842 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,0,-2) * spacing);
843 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
844 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
845 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
846 : #endif //HAVE_SIMD
847 :
848 : //nearest neighbour interpolated
849 : grid.setInterpolationType(NEAREST_NEIGHBOUR);
850 1 : b = grid.interpolate(pos + Vector3d(1,0,0) * spacing);
851 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(1,0,0) * spacing);
852 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
853 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
854 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
855 :
856 1 : b = grid.interpolate(pos + Vector3d(0,5,0) * spacing);
857 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,5,0) * spacing);
858 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
859 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
860 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
861 :
862 1 : b = grid.interpolate(pos + Vector3d(0,0,-2) * spacing);
863 1 : b2 = grid.interpolate(pos *(-1) - Vector3d(0,0,-2) * spacing);
864 1 : EXPECT_FLOAT_EQ(b.x, b2.x);
865 1 : EXPECT_FLOAT_EQ(b.y, b2.y);
866 1 : EXPECT_FLOAT_EQ(b.z, b2.z);
867 :
868 : //Test if grid is set to zero outside of volume for clipVolume=true
869 : grid.setClipVolume(true);
870 1 : Vector3f b3 = grid.interpolate(pos + Vector3d(0, 0, -2) * size);
871 1 : EXPECT_FLOAT_EQ(0., b3.x);
872 1 : EXPECT_FLOAT_EQ(0., b3.y);
873 1 : EXPECT_FLOAT_EQ(0., b3.z);
874 1 : }
875 :
876 1 : TEST(Grid3f, DumpLoad) {
877 : // Dump and load a field grid
878 1 : ref_ptr<Grid3f> grid1 = new Grid3f(Vector3d(0.), 3, 1);
879 1 : ref_ptr<Grid3f> grid2 = new Grid3f(Vector3d(0.), 3, 1);
880 :
881 4 : for (int ix = 0; ix < 3; ix++)
882 12 : for (int iy = 0; iy < 3; iy++)
883 36 : for (int iz = 0; iz < 3; iz++)
884 27 : grid1->get(ix, iy, iz) = Vector3f(1, 2, 3);
885 :
886 2 : dumpGrid(grid1, "testDump.raw");
887 2 : loadGrid(grid2, "testDump.raw");
888 :
889 4 : for (int ix = 0; ix < 3; ix++) {
890 12 : for (int iy = 0; iy < 3; iy++) {
891 36 : for (int iz = 0; iz < 3; iz++) {
892 27 : Vector3f b1 = grid1->get(ix, iy, iz);
893 : Vector3f b2 = grid2->get(ix, iy, iz);
894 27 : EXPECT_FLOAT_EQ(b1.x, b2.x);
895 27 : EXPECT_FLOAT_EQ(b1.y, b2.y);
896 27 : EXPECT_FLOAT_EQ(b1.z, b2.z);
897 : }
898 : }
899 : }
900 1 : }
901 :
902 1 : TEST(Grid3f, DumpLoadTxt) {
903 : // Dump and load a field grid
904 1 : ref_ptr<Grid3f> grid1 = new Grid3f(Vector3d(0.), 3, 1);
905 1 : ref_ptr<Grid3f> grid2 = new Grid3f(Vector3d(0.), 3, 1);
906 :
907 4 : for (int ix = 0; ix < 3; ix++)
908 12 : for (int iy = 0; iy < 3; iy++)
909 36 : for (int iz = 0; iz < 3; iz++)
910 27 : grid1->get(ix, iy, iz) = Vector3f(ix, iy, iz);
911 :
912 2 : dumpGridToTxt(grid1, "testDump.txt", 1e4);
913 2 : loadGridFromTxt(grid2, "testDump.txt", 1e-4);
914 :
915 4 : for (int ix = 0; ix < 3; ix++) {
916 12 : for (int iy = 0; iy < 3; iy++) {
917 36 : for (int iz = 0; iz < 3; iz++) {
918 27 : Vector3f b1 = grid1->get(ix, iy, iz);
919 : Vector3f b2 = grid2->get(ix, iy, iz);
920 27 : EXPECT_FLOAT_EQ(b1.x, b2.x);
921 27 : EXPECT_FLOAT_EQ(b1.y, b2.y);
922 27 : EXPECT_FLOAT_EQ(b1.z, b2.z);
923 : }
924 : }
925 : }
926 1 : }
927 :
928 1 : TEST(Grid1f, DumpLoadTxtGridProperties) {
929 : // grid to dump
930 1 : ref_ptr<Grid1f> grid = new Grid1f(Vector3d(0.5, 1.5, 2.5), 3, 2, 4, Vector3d(0.2, 1.2, 2.2));
931 : grid -> setInterpolationType(TRICUBIC);
932 : grid -> setReflective(true);
933 : grid -> setClipVolume(true);
934 :
935 : // set some values for the grid
936 4 : for (int ix = 0; ix < grid -> getNx(); ix++) {
937 9 : for (int iy = 0; iy < grid -> getNy(); iy++) {
938 30 : for (int iz = 0; iz < grid -> getNz(); iz++)
939 24 : grid -> get(ix, iy, iz) = ix + iy + iz;
940 : }
941 : }
942 :
943 : // store with properties
944 2 : dumpGridToTxt(grid, "testDump.txt", 1., true);
945 :
946 : // reload grid
947 2 : ref_ptr<Grid1f> loadedGrid = loadGrid1fFromTxt("testDump.txt");
948 :
949 : // check grid properties
950 1 : EXPECT_EQ(grid -> getNx(), loadedGrid -> getNx());
951 1 : EXPECT_EQ(grid -> getNy(), loadedGrid -> getNy());
952 1 : EXPECT_EQ(grid -> getNz(), loadedGrid -> getNz());
953 :
954 1 : EXPECT_TRUE(loadedGrid -> getClipVolume());
955 :
956 : Vector3d orig = grid -> getOrigin();
957 : Vector3d loadedOrigin = loadedGrid -> getOrigin();
958 1 : EXPECT_EQ(orig.x, loadedOrigin.x);
959 1 : EXPECT_EQ(orig.y, loadedOrigin.y);
960 1 : EXPECT_EQ(orig.z, loadedOrigin.z);
961 :
962 : Vector3d spacing = grid -> getSpacing();
963 : Vector3d loadedSpacing = loadedGrid -> getSpacing();
964 1 : EXPECT_EQ(spacing.x, loadedSpacing.x);
965 1 : EXPECT_EQ(spacing.y, loadedSpacing.y);
966 1 : EXPECT_EQ(spacing.z, loadedSpacing.z);
967 :
968 1 : EXPECT_EQ(grid -> getInterpolationType(), loadedGrid -> getInterpolationType());
969 :
970 1 : EXPECT_EQ(grid -> isReflective(), loadedGrid -> isReflective());
971 :
972 : // compare loaded values
973 4 : for (int ix = 0; ix < grid -> getNx(); ix++) {
974 9 : for (int iy = 0; iy < grid -> getNy(); iy++) {
975 30 : for (int iz = 0; iz < grid -> getNz(); iz++)
976 24 : EXPECT_EQ(grid -> get(ix, iy, iz), loadedGrid -> get(ix, iy, iz));
977 : }
978 : }
979 1 : }
980 :
981 1 : TEST(Grid3f, DumpLoadTxtGridProperties) {
982 : // grid to dump
983 1 : ref_ptr<Grid3f> grid = new Grid3f(Vector3d(0.5, 1.5, 2.5), 3, 2, 4, Vector3d(0.2, 1.2, 2.2));
984 : grid -> setInterpolationType(NEAREST_NEIGHBOUR);
985 : grid -> setReflective(true);
986 : grid -> setClipVolume(true);
987 :
988 : // set some values for the grid
989 4 : for (int ix = 0; ix < grid -> getNx(); ix++) {
990 9 : for (int iy = 0; iy < grid -> getNy(); iy++) {
991 30 : for (int iz = 0; iz < grid -> getNz(); iz++)
992 24 : grid -> get(ix, iy, iz) = Vector3f(-iy, ix, 2 * iz);
993 : }
994 : }
995 :
996 : // store with properties
997 2 : dumpGridToTxt(grid, "testDump.txt", 1., true);
998 :
999 : // reload grid
1000 2 : ref_ptr<Grid3f> loadedGrid = loadGrid3fFromTxt("testDump.txt");
1001 :
1002 : // check grid properties
1003 1 : EXPECT_EQ(grid -> getNx(), loadedGrid -> getNx());
1004 1 : EXPECT_EQ(grid -> getNy(), loadedGrid -> getNy());
1005 1 : EXPECT_EQ(grid -> getNz(), loadedGrid -> getNz());
1006 :
1007 1 : EXPECT_TRUE(loadedGrid -> getClipVolume());
1008 :
1009 : Vector3d orig = grid -> getOrigin();
1010 : Vector3d loadedOrigin = loadedGrid -> getOrigin();
1011 1 : EXPECT_EQ(orig.x, loadedOrigin.x);
1012 1 : EXPECT_EQ(orig.y, loadedOrigin.y);
1013 1 : EXPECT_EQ(orig.z, loadedOrigin.z);
1014 :
1015 : Vector3d spacing = grid -> getSpacing();
1016 : Vector3d loadedSpacing = loadedGrid -> getSpacing();
1017 1 : EXPECT_EQ(spacing.x, loadedSpacing.x);
1018 1 : EXPECT_EQ(spacing.y, loadedSpacing.y);
1019 1 : EXPECT_EQ(spacing.z, loadedSpacing.z);
1020 :
1021 1 : EXPECT_EQ(grid -> getInterpolationType(), loadedGrid -> getInterpolationType());
1022 :
1023 1 : EXPECT_EQ(grid -> isReflective(), loadedGrid -> isReflective());
1024 :
1025 : // compare loaded values
1026 4 : for (int ix = 0; ix < grid -> getNx(); ix++) {
1027 9 : for (int iy = 0; iy < grid -> getNy(); iy++) {
1028 30 : for (int iz = 0; iz < grid -> getNz(); iz++) {
1029 : Vector3f gridValue = grid -> get(ix, iy, iz);
1030 : Vector3f loadedValue = loadedGrid -> get(ix, iy, iz);
1031 24 : EXPECT_EQ(gridValue.x, loadedValue.x);
1032 24 : EXPECT_EQ(gridValue.y, loadedValue.y);
1033 24 : EXPECT_EQ(gridValue.z, loadedValue.z);
1034 : }
1035 : }
1036 : }
1037 1 : }
1038 :
1039 2 : TEST(Grid3f, Speed) {
1040 : // Dump and load a field grid
1041 1 : Grid3f grid(Vector3d(0.), 3, 3);
1042 4 : for (int ix = 0; ix < 3; ix++)
1043 12 : for (int iy = 0; iy < 3; iy++)
1044 36 : for (int iz = 0; iz < 3; iz++)
1045 27 : grid.get(ix, iy, iz) = Vector3f(1, 2, 3);
1046 :
1047 : Vector3d b;
1048 100001 : for (int i = 0; i < 100000; i++)
1049 100000 : b = grid.interpolate(Vector3d(i));
1050 1 : }
1051 :
1052 2 : TEST(CylindricalProjectionMap, functions) {
1053 : Vector3d v;
1054 1 : v.setRThetaPhi(1.0, 1.2, 2.4);
1055 1 : EXPECT_NEAR(v.getPhi(), 2.4, .00001);
1056 1 : EXPECT_NEAR(v.getTheta(), 1.2, .000001);
1057 :
1058 1 : CylindricalProjectionMap cpm(24, 12);
1059 1 : size_t bin = 50;
1060 1 : Vector3d d = cpm.directionFromBin(bin);
1061 1 : size_t bin2 = cpm.binFromDirection(d);
1062 1 : EXPECT_EQ(bin, bin2);
1063 1 : }
1064 :
1065 1 : TEST(EmissionMap, functions) {
1066 :
1067 1 : EmissionMap em(360, 180, 100, 1 * EeV, 100 * EeV);
1068 1 : double e = em.energyFromBin(50);
1069 1 : size_t b = em.binFromEnergy(50 * EeV);
1070 :
1071 : Vector3d d(1.0, 0.0, 0.0);
1072 :
1073 1 : em.fillMap(1, 50 * EeV, d);
1074 :
1075 : Vector3d d2;
1076 :
1077 1 : bool r = em.drawDirection(1, 50 * EeV, d2);
1078 1 : EXPECT_TRUE(r);
1079 1 : EXPECT_TRUE(d.getAngleTo(d2) < (2. * M_PI / 180.));
1080 :
1081 1 : r = em.drawDirection(1, 30 * EeV, d2);
1082 1 : EXPECT_FALSE(r);
1083 :
1084 1 : r = em.drawDirection(2, 50 * EeV, d2);
1085 1 : EXPECT_FALSE(r);
1086 1 : }
1087 :
1088 1 : TEST(EmissionMap, merge) {
1089 1 : EmissionMap em1, em2;
1090 1 : em1.fillMap(1, 50 * EeV, Vector3d(1.0, 0.0, 0.0));
1091 1 : em2.fillMap(1, 50 * EeV, Vector3d(0.0, 1.0, 0.0));
1092 1 : em2.fillMap(2, 50 * EeV, Vector3d(0.0, 1.0, 0.0));
1093 :
1094 1 : em1.merge(&em2);
1095 :
1096 1 : EXPECT_EQ(em1.getMaps().size(), 2);
1097 :
1098 1 : ref_ptr<CylindricalProjectionMap> cpm = em1.getMap(1, 50 * EeV);
1099 1 : size_t bin = cpm->binFromDirection(Vector3d(0.0, 1.0, 0.0));
1100 1 : EXPECT_TRUE(cpm->getPdf()[bin] > 0);
1101 1 : }
1102 :
1103 1 : TEST(Variant, copyToBuffer) {
1104 1 : double a = 23.42;
1105 : Variant v(a);
1106 : double b;
1107 1 : v.copyToBuffer(&b);
1108 1 : EXPECT_EQ(a, b);
1109 1 : }
1110 :
1111 1 : TEST(Variant, stringConversion) {
1112 1 : Variant v, w;
1113 : {
1114 1 : int32_t a = 12;
1115 1 : v = Variant::fromInt32(a);
1116 1 : EXPECT_EQ(a, v.asInt32());
1117 :
1118 2 : w = Variant::fromString(v.toString(), v.getType());
1119 1 : EXPECT_EQ(a, w.asInt32());
1120 : }
1121 :
1122 : {
1123 1 : int64_t a = 12;
1124 1 : v = Variant::fromInt64(a);
1125 1 : EXPECT_EQ(a, v.asInt64());
1126 :
1127 2 : w = Variant::fromString(v.toString(), v.getType());
1128 1 : EXPECT_EQ(a, w.asInt64());
1129 : }
1130 :
1131 : {
1132 : Vector3d a(1, 2, 3);
1133 : Variant v = Variant::fromVector3d(a);
1134 : Vector3d u = v.asVector3d();
1135 1 : EXPECT_EQ(a.getX(), u.getX());
1136 1 : EXPECT_EQ(a.getY(), u.getY());
1137 1 : EXPECT_EQ(a.getZ(), u.getZ());
1138 1 : }
1139 :
1140 : {
1141 1 : std::complex<double> a1(1, 1);
1142 1 : std::complex<double> a2(2, 0);
1143 : Vector3<std::complex<double>> a(a1, a1, a2);
1144 : Variant v = Variant::fromVector3c(a);
1145 : Vector3<std::complex<double>> u = v.asVector3c();
1146 1 : EXPECT_EQ(a1, u.getX());
1147 1 : EXPECT_EQ(a1, u.getY());
1148 1 : EXPECT_EQ(a2, u.getZ());
1149 1 : }
1150 :
1151 : {
1152 : std::complex<double> a(1, 2);
1153 : Variant v = Variant::fromComplexDouble(a);
1154 1 : std::complex<double> u = v.asComplexDouble();
1155 1 : EXPECT_EQ(u.real(), 1);
1156 1 : EXPECT_EQ(u.imag(), 2);
1157 1 : }
1158 :
1159 : {
1160 : std::vector<Variant> a;
1161 1 : a.push_back(Variant::fromDouble(1));
1162 1 : a.push_back(Variant::fromDouble(2));
1163 1 : a.push_back(Variant::fromDouble(3));
1164 1 : a.push_back(Variant::fromDouble(4));
1165 : Variant v = Variant::fromVector(a);
1166 1 : std::vector<Variant> u = v.asVector();
1167 1 : EXPECT_EQ(a[0], Variant::fromDouble(u[0]));
1168 1 : EXPECT_EQ(a[1], Variant::fromDouble(u[1]));
1169 1 : EXPECT_EQ(a[2], Variant::fromDouble(u[2]));
1170 1 : EXPECT_EQ(a[3], Variant::fromDouble(u[3]));
1171 1 : }
1172 :
1173 :
1174 1 : }
1175 :
1176 2 : TEST(Geometry, Plane) {
1177 1 : Plane p(Vector3d(0,0,1), Vector3d(0,0,1));
1178 1 : EXPECT_DOUBLE_EQ(-1., p.distance(Vector3d(0, 0, 0)));
1179 1 : EXPECT_DOUBLE_EQ(9., p.distance(Vector3d(1, 1, 10)));
1180 1 : }
1181 :
1182 2 : TEST(Geometry, Sphere) {
1183 1 : Sphere s(Vector3d(0,0,0), 1.);
1184 1 : EXPECT_DOUBLE_EQ(-1., s.distance(Vector3d(0, 0, 0)));
1185 1 : EXPECT_DOUBLE_EQ(9., s.distance(Vector3d(10, 0, 0)));
1186 1 : }
1187 :
1188 2 : TEST(Geometry, ParaxialBox) {
1189 1 : ParaxialBox b(Vector3d(0,0,0), Vector3d(3,4,5));
1190 1 : EXPECT_NEAR(-.1, b.distance(Vector3d(0.1, 0.1, 0.1)), 1E-10);
1191 1 : EXPECT_NEAR(-.1, b.distance(Vector3d(0.1, 3.8, 0.1)), 1E-10);
1192 1 : EXPECT_NEAR(-.2, b.distance(Vector3d(0.9, 3.8, 0.9)), 1E-10);
1193 1 : EXPECT_NEAR(7., b.distance(Vector3d(10., 0., 0.)), 1E-10);
1194 1 : EXPECT_NEAR(7., b.distance(Vector3d(10., 2., 0.)), 1E-10);
1195 1 : EXPECT_NEAR(8., b.distance(Vector3d(-8., 0., 0.)), 1E-10);
1196 1 : }
1197 :
1198 0 : int main(int argc, char **argv) {
1199 0 : ::testing::InitGoogleTest(&argc, argv);
1200 0 : return RUN_ALL_TESTS();
1201 : }
1202 :
1203 : } // namespace crpropa
|
