Line data Source code
1 : #include <stdexcept>
2 :
3 : #include "crpropa/magneticField/MagneticFieldGrid.h"
4 : #include "crpropa/magneticField/CMZField.h"
5 : #include "crpropa/magneticField/PolarizedSingleModeMagneticField.h"
6 : #include "crpropa/magneticField/GalacticMagneticField.h"
7 : #include "crpropa/magneticField/UF23Field.h"
8 : #include "crpropa/Grid.h"
9 : #include "crpropa/Units.h"
10 : #include "crpropa/Common.h"
11 :
12 : #include "gtest/gtest.h"
13 :
14 : using namespace crpropa;
15 :
16 1 : TEST(testUniformMagneticField, SimpleTest) {
17 : UniformMagneticField B(Vector3d(-1, 5, 3));
18 : Vector3d b = B.getField(Vector3d(1, 0, 0));
19 1 : EXPECT_DOUBLE_EQ(b.getX(), -1);
20 1 : EXPECT_DOUBLE_EQ(b.getY(), 5);
21 1 : EXPECT_DOUBLE_EQ(b.getZ(), 3);
22 1 : }
23 :
24 2 : TEST(testMagneticDipoleField, SimpleTest) {
25 : // Test magnetic dipole
26 : // mu0 / (4*M_PI) * m / r^3 (2*cos(theta)*e_r + sin(theta)*e_theta)
27 : MagneticDipoleField B(Vector3d(0,0,0), Vector3d(0,0,1), 1);
28 1 : Vector3d b1 = B.getField(Vector3d(0, 0, 1)); // theta = 0
29 1 : Vector3d b2 = B.getField(Vector3d(1, 0, 0)); // theta = 0
30 1 : EXPECT_NEAR(b1.getX(), 0, 1E-8);
31 1 : EXPECT_NEAR(b1.getY(), 0, 1E-8);
32 1 : EXPECT_NEAR(b1.getZ(), mu0 / (4*M_PI) * 2, 1E-8);
33 1 : EXPECT_NEAR(b2.getX(), 0, 1E-8);
34 1 : EXPECT_NEAR(b2.getY(), 0, 1E-8);
35 1 : EXPECT_NEAR(b2.getZ(), -1 * mu0 / (4*M_PI), 1E-8);
36 1 : }
37 :
38 : #ifdef CRPROPA_HAVE_MUPARSER
39 1 : TEST(testRenormalizeMagneticField, simpleTest) {
40 1 : ref_ptr<UniformMagneticField> field = new UniformMagneticField(Vector3d(2*nG, 0, 0));
41 2 : RenormalizeMagneticField modField(field, "B^2-1*nG");
42 1 : Vector3d b = modField.getField(Vector3d(5, 5, 5));
43 1 : EXPECT_NEAR(b.getR(), 3*nG, 0.001);
44 2 : }
45 : #endif
46 :
47 2 : TEST(testMagneticFieldList, SimpleTest) {
48 : // Test a list of three magnetic fields
49 : MagneticFieldList B;
50 1 : B.addField(new UniformMagneticField(Vector3d(1, 0, 0)));
51 1 : B.addField(new UniformMagneticField(Vector3d(0, 2, 0)));
52 2 : B.addField(new UniformMagneticField(Vector3d(0, 0, 3)));
53 1 : Vector3d b = B.getField(Vector3d(0.));
54 1 : EXPECT_DOUBLE_EQ(b.x, 1);
55 1 : EXPECT_DOUBLE_EQ(b.y, 2);
56 1 : EXPECT_DOUBLE_EQ(b.z, 3);
57 1 : }
58 :
59 1 : TEST(testMagneticFieldEvolution, SimpleTest) {
60 : // Test if this decorator scales the underlying field as (1+z)^m
61 1 : ref_ptr<UniformMagneticField> B = new UniformMagneticField(Vector3d(1,0,0));
62 : double z = 1.2;
63 : double m = 3;
64 2 : MagneticFieldEvolution Bz(B, m);
65 :
66 : // scaled field
67 1 : Vector3d b = Bz.getField(Vector3d(0,0,0), z);
68 1 : EXPECT_DOUBLE_EQ(b.x, pow(1+z, m));
69 :
70 : // unscaled field
71 1 : b = Bz.getField(Vector3d(0,0,0));
72 1 : EXPECT_DOUBLE_EQ(b.x, 1);
73 1 : }
74 :
75 1 : class EchoMagneticField: public MagneticField {
76 : public:
77 3 : Vector3d getField(const Vector3d &position) const {
78 3 : return position;
79 : }
80 : };
81 :
82 1 : TEST(testPeriodicMagneticField, Exceptions) {
83 1 : ref_ptr<EchoMagneticField> f = new EchoMagneticField();
84 : ref_ptr<PeriodicMagneticField> p = new PeriodicMagneticField(f,
85 2 : Vector3d(10000, 10000, 10000), Vector3d(1000, 1000, 1000), true);
86 :
87 : // box 0, 0, 0
88 1 : Vector3d v = p->getField(Vector3d(1000, 2000, 3000));
89 1 : EXPECT_DOUBLE_EQ(0, v.x);
90 1 : EXPECT_DOUBLE_EQ(1000, v.y);
91 1 : EXPECT_DOUBLE_EQ(2000, v.z);
92 :
93 : // box 1, 2, 3
94 1 : v = p->getField(Vector3d(12000, 23000, 34000));
95 1 : EXPECT_DOUBLE_EQ(9000, v.x);
96 1 : EXPECT_DOUBLE_EQ(2000, v.y);
97 1 : EXPECT_DOUBLE_EQ(7000, v.z);
98 :
99 : // box -1, -2, -3
100 1 : v = p->getField(Vector3d(0, -10000, -20000));
101 1 : EXPECT_DOUBLE_EQ(1000, v.x);
102 1 : EXPECT_DOUBLE_EQ(9000, v.y);
103 1 : EXPECT_DOUBLE_EQ(1000, v.z);
104 :
105 1 : }
106 :
107 1 : TEST(testCMZMagneticField, SimpleTest) {
108 1 : ref_ptr<CMZField> field = new CMZField();
109 :
110 : // check use-Values
111 1 : EXPECT_FALSE(field->getUseMCField());
112 1 : EXPECT_TRUE(field->getUseICField());
113 1 : EXPECT_FALSE(field->getUseNTFField());
114 1 : EXPECT_FALSE(field->getUseRadioArc());
115 :
116 : // check set function
117 1 : field->setUseMCField(true);
118 1 : EXPECT_TRUE(field->getUseMCField());
119 1 : field->setUseICField(false);
120 1 : EXPECT_FALSE(field->getUseICField());
121 1 : field->setUseNTFField(true);
122 1 : EXPECT_TRUE(field->getUseNTFField());
123 1 : field->setUseRadioArc(true);
124 1 : EXPECT_TRUE(field->getUseRadioArc());
125 1 : }
126 :
127 1 : TEST(testCMZMagneticField, TestICComponent) {
128 1 : ref_ptr<CMZField> field = new CMZField();
129 : Vector3d pos(10*pc,15*pc,-5*pc);
130 :
131 : // check IC field at given position
132 1 : Vector3d bVec = field->getField(pos);
133 1 : EXPECT_NEAR(bVec.getR(), 10.501*muG, 1E-3*muG);
134 1 : EXPECT_NEAR(bVec.x, 0.225*muG, 1E-3*muG);
135 1 : EXPECT_NEAR(bVec.y, 0.524*muG, 1E-3*muG);
136 1 : EXPECT_NEAR(bVec.z, 10.486*muG, 1E-3*muG);
137 1 : }
138 1 : TEST(testCMZMagneticField, TestNTFField){
139 1 : ref_ptr<CMZField> field = new CMZField();
140 : Vector3d pos(10*pc,15*pc,-5*pc);
141 :
142 : // check NFTField at given position
143 1 : Vector3d bVec = field->getNTFField(pos);
144 1 : EXPECT_NEAR(bVec.getR(),1.692*muG, 1e-3*muG);
145 1 : EXPECT_NEAR(bVec.x, -0.584*muG, 1e-3*muG);
146 1 : EXPECT_NEAR(bVec.y, -1.185*muG, 1e-3*muG);
147 1 : EXPECT_NEAR(bVec.z, 1.057*muG, 1e-3*muG);
148 1 : }
149 1 : TEST(testCMZMagneticField, TestRaidoArcField){
150 1 : ref_ptr<CMZField> field = new CMZField();
151 : Vector3d pos(10*pc,15*pc,-5*pc);
152 :
153 : // check RadioArcField at given position
154 1 : Vector3d bVec = field->getRadioArcField(pos);
155 1 : EXPECT_NEAR(bVec.getR(), 31.616*muG, 1e-3*muG);
156 1 : EXPECT_NEAR(bVec.x, -4.671*muG, 1e-3*muG);
157 1 : EXPECT_NEAR(bVec.y, 5.465*muG, 1e-3*muG);
158 1 : EXPECT_NEAR(bVec.z, 30.788*muG, 1e-3*muG);
159 1 : }
160 :
161 1 : TEST(testCMZMagneticField, TestAzimutalComponent){
162 1 : ref_ptr<CMZField> field = new CMZField();
163 : Vector3d mid(12*pc, 9*pc, 20*pc);
164 : Vector3d pos(9*pc, 10*pc, 25*pc);
165 :
166 : // simple Test for inner part
167 1 : Vector3d bVec = field->BAz(pos, mid, 100, 0.2, 60*pc);
168 1 : EXPECT_NEAR(bVec.x, 3939.782, 1e-3);
169 1 : EXPECT_NEAR(bVec.y, 14347.304, 1e-3);
170 1 : EXPECT_DOUBLE_EQ(bVec.z, 0);
171 :
172 : // simple Test for outer part
173 1 : bVec = field->BAz(pos, mid, 100, 0.2, 10*pc);
174 1 : EXPECT_NEAR(bVec.x, -164.659, 1e-3);
175 1 : EXPECT_NEAR(bVec.y, -1317.270, 1e-3);
176 1 : EXPECT_DOUBLE_EQ(bVec.z, 0);
177 :
178 : // test for molecular Clouds
179 1 : bVec = field->getMCField(pos);
180 1 : EXPECT_NEAR(bVec.x, -8.339*muG, 1e-3*muG);
181 1 : EXPECT_NEAR(bVec.y, -0.850*muG, 1e-3*muG);
182 1 : EXPECT_DOUBLE_EQ(bVec.z, 0);
183 1 : }
184 :
185 1 : TEST(testToroidalHaloField, SimpleTest) {
186 1 : ref_ptr<ToroidalHaloField> field = new ToroidalHaloField();
187 1 : Vector3d b = field->getField(Vector3d(0.));
188 1 : EXPECT_DOUBLE_EQ(b.x, 0);
189 1 : EXPECT_DOUBLE_EQ(b.y, 0);
190 1 : EXPECT_DOUBLE_EQ(b.z, 0);
191 :
192 1 : b = field->getField(Vector3d(1,0,0));
193 1 : EXPECT_DOUBLE_EQ(b.x, 0.5);
194 1 : EXPECT_DOUBLE_EQ(b.y, 0);
195 1 : EXPECT_DOUBLE_EQ(b.z, 0);
196 1 : }
197 :
198 1 : TEST(testLogarithmicSpiralField, SimpleTest) {
199 1 : ref_ptr<LogarithmicSpiralField> field = new LogarithmicSpiralField();
200 1 : Vector3d b = field->getField(Vector3d(8.5, 0, 0)*kpc);
201 1 : EXPECT_NEAR(b.x, -1., 1E-2);
202 1 : EXPECT_NEAR(b.y, 0, 1E-10);
203 1 : EXPECT_NEAR(b.z, 0, 1E-10);
204 1 : }
205 :
206 :
207 : // UF23 reference values
208 : std::vector<std::vector<Vector3d>>
209 1 : getUF23ReferenceValues()
210 : {
211 : using V = Vector3d;
212 : std::vector<std::vector<Vector3d>> retVal;
213 2 : retVal.push_back({
214 : V{-1.46860417e+00, 1.64555489e+00, 8.35702311e-01},
215 : V{0.00000000e+00, -4.22433497e-01, 1.83232985e-01},
216 : V{-3.00239177e-01, -2.97045767e-01, 1.83232985e-01},
217 : V{8.58382023e-04, 7.71891049e-03, 9.71705527e-01},
218 : V{-1.17276875e+00, -2.33013590e-01, 4.10798494e-01},
219 : V{2.63883569e-01, -8.79631081e-01, 5.54229961e-06},
220 : V{-1.71971907e-02, -2.00291358e-02, 4.16024875e-02},
221 : V{6.18960813e-01, -2.14437026e+00, 8.35702315e-01},
222 : V{-1.50883911e+00, 2.63874909e+00, 1.16578928e-02}
223 : });
224 2 : retVal.push_back({
225 : V{-1.39447625e+00, 1.57135437e+00, 8.65163870e-01},
226 : V{0.00000000e+00, -4.44318027e-01, 1.54430718e-01},
227 : V{-3.15695875e-01, -3.12652066e-01, 1.54430718e-01},
228 : V{2.09678202e-03, 2.63832916e-03, 9.81077691e-01},
229 : V{-1.08603376e+00, -1.86359136e-01, 3.91730436e-01},
230 : V{2.11095801e-01, -7.04082667e-01, 8.86591950e-05},
231 : V{-1.38087843e-02, -1.52011002e-02, 3.06842086e-02},
232 : V{5.76543516e-01, -2.03544911e+00, 8.65163870e-01},
233 : V{-3.77778745e-01, 6.89644787e-01, 5.85347897e-02},
234 : });
235 2 : retVal.push_back({
236 : V{-1.04170166e+00, 1.93212739e+00, 2.95362499e+00},
237 : V{0.00000000e+00, -5.87995638e-01, 4.66925506e-01},
238 : V{-4.20802701e-01, -4.10291633e-01, 4.59557272e-01},
239 : V{7.79982996e-03, 1.50482351e-02, 5.50337531e+00},
240 : V{-1.34639246e+00, 6.80103301e-02, 8.60002649e-01},
241 : V{1.37242004e-01, -8.67085225e-01, 1.25106474e-01},
242 : V{-1.69325979e-02, -2.95454738e-02, 4.72616404e-02},
243 : V{4.53873303e-01, -2.03292501e+00, 9.95205454e-01},
244 : V{-1.30599234e+00, 2.25151057e+00, 2.56704152e-01},
245 : });
246 2 : retVal.push_back({
247 : V{-1.45840221e+00, 1.64227817e+00, 8.41586777e-01},
248 : V{0.00000000e+00, -6.98341816e-01, 1.84323895e-01},
249 : V{-4.95342500e-01, -4.92154113e-01, 1.84323895e-01},
250 : V{-5.27533653e-03, 1.48965596e-02, 9.86107885e-01},
251 : V{-1.47472798e+00, -3.33905274e-01, 4.11265918e-01},
252 : V{2.31202234e-01, -7.70722755e-01, 1.12129423e-05},
253 : V{-1.54200071e-02, -2.22012150e-02, 4.22733003e-02},
254 : V{5.92201401e-01, -2.10378256e+00, 8.41586782e-01},
255 : V{4.05057122e-03, -9.76673905e-03, 8.24321504e-03},
256 : });
257 2 : retVal.push_back({
258 : V{-1.50595160e+00, 1.67869437e+00, 8.29806168e-01},
259 : V{0.00000000e+00, -2.27289811e-01, 1.86869632e-01},
260 : V{-1.62291024e-01, -1.59076795e-01, 1.86869632e-01},
261 : V{6.71536463e-04, 7.88990042e-03, 9.63291786e-01},
262 : V{-9.35561196e-01, -1.55840707e-01, 4.13609050e-01},
263 : V{2.68120686e-01, -8.93743434e-01, 3.48927149e-06},
264 : V{-1.88291055e-02, -1.94245715e-02, 4.29970299e-02},
265 : V{6.50804557e-01, -2.18817590e+00, 8.29806172e-01},
266 : V{-1.58450595e+00, 2.77817014e+00, 1.10336683e-02},
267 : });
268 2 : retVal.push_back({
269 : V{-1.33098823e+00, 1.49556466e+00, 6.88642986e-01},
270 : V{0.00000000e+00, -4.99342453e-01, 1.16143813e-01},
271 : V{-3.54225500e-01, -3.51947350e-01, 1.16143813e-01},
272 : V{2.14352986e-03, 3.57971370e-03, 8.05219599e-01},
273 : V{-1.15136627e+00, -2.48536546e-01, 2.95713475e-01},
274 : V{1.18057505e-01, -3.98308301e-01, 9.11299352e-04},
275 : V{-1.06899006e-02, -1.12335952e-02, 2.33358311e-02},
276 : V{5.61264936e-01, -1.94601963e+00, 6.88642987e-01},
277 : V{-1.18200258e+00, 2.01179003e+00, 8.02453884e-02},
278 : });
279 2 : retVal.push_back({
280 : V{-3.65508814e-01, 4.74456797e-01, 5.76165841e-01},
281 : V{0.00000000e+00, 0.00000000e+00, 6.36668108e-02},
282 : V{-3.68784336e-03, -2.20689056e-03, 6.36668108e-02},
283 : V{-5.03208784e-02, 5.09887480e-02, 6.27665021e-01},
284 : V{-4.04821314e-01, -1.01426396e-02, 2.06022291e-01},
285 : V{-2.34287239e-02, -9.72536117e-02, 2.80624948e-02},
286 : V{-4.42698169e-03, -6.23429869e-03, 1.07555956e-02},
287 : V{3.25022555e-01, -1.18950549e+00, 5.76163734e-01},
288 : V{-8.65217092e-01, 1.85270104e+00, 3.73913202e-01},
289 : });
290 2 : retVal.push_back({
291 : V{-1.92580231e+00, 2.17134243e+00, 1.13723929e+00},
292 : V{0.00000000e+00, -4.73261099e-01, 2.65974345e-01},
293 : V{-3.36780079e-01, -3.32368900e-01, 2.65974345e-01},
294 : V{-5.24558924e-04, 1.51886238e-02, 1.33757258e+00},
295 : V{-1.47213976e+00, -2.82203258e-01, 5.71920142e-01},
296 : V{3.54206333e-01, -1.18108962e+00, 1.00077423e-04},
297 : V{-2.67268842e-02, -2.88588605e-02, 6.38364561e-02},
298 : V{7.90570955e-01, -2.84039611e+00, 1.13723931e+00},
299 : V{-1.97316856e+00, 3.44153600e+00, 3.20266742e-02},
300 : });
301 1 : return retVal;
302 0 : }
303 :
304 :
305 1 : TEST(testUF23Field, SimpleTest) {
306 : const std::vector<UF23Field::ModelType> models =
307 : {
308 : UF23Field::base, UF23Field::neCL, UF23Field::expX, UF23Field::spur,
309 : UF23Field::cre10, UF23Field::synCG, UF23Field::twistX, UF23Field::nebCor
310 1 : };
311 :
312 : using V = Vector3d;
313 : const std::vector<Vector3d> testPositions =
314 : {
315 : V{1, 1, 1}, V{0, 0, -8}, V{0.1, -0.1, -8}, V{0.1, 0.1, 0.1}, V{-1, 3, 4},
316 : V{-10, -3, 2}, V{-10, -10, 20}, V{-4, -2, 1}, V{6, 5, -0.1}
317 1 : };
318 :
319 : const std::vector<std::vector<Vector3d>> referenceValues =
320 1 : getUF23ReferenceValues();
321 :
322 : const double precision = 1e-6;
323 :
324 9 : for (unsigned int i = 0; i < models.size(); ++i) {
325 : const auto& model = models[i];
326 8 : const UF23Field uf23Field(model);
327 80 : for (unsigned int j = 0; j < testPositions.size(); ++j) {
328 : const auto& position = testPositions[j];
329 72 : const auto val = uf23Field.getField(position*kpc);
330 : const auto& refVal = referenceValues[i][j] * microgauss;
331 72 : EXPECT_NEAR(val.x, refVal.x, precision);
332 72 : EXPECT_NEAR(val.y, refVal.y, precision);
333 72 : EXPECT_NEAR(val.z, refVal.z, precision);
334 : }
335 : }
336 1 : }
337 :
338 1 : TEST(testPolarizedSingleModeMagneticField, SimpleTest) {
339 1 : PolarizedSingleModeMagneticField B(2, 4, 0.5, Vector3d(1,1,1), Vector3d(0,1,0), Vector3d(1,0,0), "amplitude", "polarization", "elliptical");
340 1 : Vector3d b = B.getField(Vector3d(1,1,2));
341 1 : EXPECT_DOUBLE_EQ(b.x, 1);
342 1 : EXPECT_NEAR(b.y, 0, 1E-10);
343 1 : EXPECT_NEAR(b.z, 0, 1E-10);
344 1 : }
345 :
346 1 : int main(int argc, char **argv) {
347 1 : ::testing::InitGoogleTest(&argc, argv);
348 : return RUN_ALL_TESTS();
349 : }
|
