Line data Source code
1 : #include "crpropa/massDistribution/Massdistribution.h"
2 : #include "crpropa/massDistribution/Cordes.h"
3 : #include "crpropa/massDistribution/Ferriere.h"
4 : #include "crpropa/massDistribution/Nakanishi.h"
5 : #include "crpropa/massDistribution/ConstantDensity.h"
6 : #include "crpropa/Units.h"
7 : #include "crpropa/Grid.h"
8 :
9 : #include "gtest/gtest.h"
10 :
11 : #include <stdexcept>
12 : #include <cmath>
13 : #include <string>
14 :
15 : namespace crpropa {
16 :
17 1 : TEST(testConstantDensity, SimpleTest) {
18 : //test ConstantDensity in all types and in total density (output)
19 1 : ConstantDensity n(2/ccm,3/ccm, 2/ccm);
20 : Vector3d p(1*pc,2*pc,1*kpc); // random position for testing density
21 1 : EXPECT_DOUBLE_EQ(n.getHIDensity(p), 2e6); // density output in m^-3
22 1 : EXPECT_DOUBLE_EQ(n.getHIIDensity(p), 3e6);
23 1 : EXPECT_DOUBLE_EQ( n.getH2Density(p), 2e6);
24 1 : EXPECT_DOUBLE_EQ(n.getDensity(p), 7e6); // total density 2+3+2 = 7 (/ccm)
25 1 : EXPECT_DOUBLE_EQ(n.getNucleonDensity(p),9e6); // nucleon density 2+3+2*2 = 9 (/ccm) factor 2 for molecular hydrogen
26 :
27 : //test set/get function for used type
28 1 : bool useHI = n.getIsForHI();
29 1 : bool useHII= n.getIsForHII();
30 1 : bool useH2 = n.getIsForH2();
31 : //check if all types are activited
32 1 : EXPECT_TRUE(useHI);
33 1 : EXPECT_TRUE(useHII);
34 1 : EXPECT_TRUE(useH2);
35 :
36 : //set density number to 500
37 1 : n.setHI(500.);
38 1 : n.setHII(500.);
39 1 : n.setH2(500.);
40 :
41 : //check if output is changed to 500 in types and is 0 (all deactivated) for Hges
42 1 : EXPECT_DOUBLE_EQ(n.getHIDensity(p), 500.);
43 1 : EXPECT_DOUBLE_EQ(n.getHIIDensity(p), 500.);
44 1 : EXPECT_DOUBLE_EQ(n.getH2Density(p), 500.);
45 :
46 : //deactivate all
47 1 : n.setHI(false);
48 1 : n.setHII(false);
49 1 : n.setH2(false);
50 :
51 : //check if all isfor are set to false
52 1 : useHI = n.getIsForHI();
53 1 : useHII= n.getIsForHII();
54 1 : useH2 = n.getIsForH2();
55 1 : EXPECT_FALSE(useHI);
56 1 : EXPECT_FALSE(useHII);
57 1 : EXPECT_FALSE(useH2);
58 :
59 : //get<type>Density is independent of type activation, getDensity is not independent
60 : //check if getDensity returns 0. (should give a error message in log file)
61 1 : ::testing::internal::CaptureStderr();
62 1 : EXPECT_DOUBLE_EQ(n.getDensity(p), 0);
63 1 : EXPECT_DOUBLE_EQ(n.getNucleonDensity(p),0);
64 1 : std::string captured = testing::internal::GetCapturedStderr();
65 1 : EXPECT_NE(captured.find("WARNING"), std::string::npos);
66 1 : }
67 :
68 :
69 2 : TEST(testDensityList, SimpleTest) {
70 :
71 : DensityList MS;
72 1 : MS.addDensity(new ConstantDensity(1,1,2)); //sum 4
73 2 : MS.addDensity(new ConstantDensity(2,3,1)); //sum 6
74 :
75 : Vector3d p(50*pc,10*pc,-30*pc); //random position for testing density
76 1 : EXPECT_DOUBLE_EQ(MS.getHIDensity(p),3);
77 1 : EXPECT_DOUBLE_EQ(MS.getHIIDensity(p),4);
78 1 : EXPECT_DOUBLE_EQ(MS.getH2Density(p),3);
79 1 : EXPECT_DOUBLE_EQ(MS.getDensity(p),10); //sum of sums
80 1 : EXPECT_DOUBLE_EQ(MS.getNucleonDensity(p),13); // 3+4+2*3 factor 2 for molecular hydrogen
81 :
82 1 : }
83 :
84 2 : TEST(testCordes, checkValueAtCertainPoints) {
85 :
86 : Cordes n;
87 :
88 : //check type Information
89 1 : EXPECT_FALSE(n.getIsForHI());
90 1 : EXPECT_TRUE(n.getIsForHII());
91 1 : EXPECT_FALSE(n.getIsForH2());
92 :
93 : Vector3d p(3.1*kpc,2.9*kpc,-30*pc); //position for testing density
94 :
95 1 : EXPECT_NEAR(n.getHIIDensity(p), 184500.,1); // output in m^-3 ; uncertainty of 1e-6 cm^-3
96 1 : EXPECT_NEAR(n.getDensity(p), 184500.,1);
97 1 : EXPECT_NEAR(n.getNucleonDensity(p), 184500,1); // only HII component -> no differenz between density and nucleon density
98 1 : p.z=30*pc; // invariant density for +/- z
99 1 : EXPECT_NEAR(n.getDensity(p),184500,1);
100 :
101 1 : }
102 :
103 2 : TEST(testNakanishi, checkValueAtCertainPoints) {
104 :
105 : Nakanishi n;
106 :
107 : //check type Information
108 1 : EXPECT_TRUE(n.getIsForHI());
109 1 : EXPECT_FALSE(n.getIsForHII());
110 1 : EXPECT_TRUE(n.getIsForH2());
111 :
112 : //first position for testing density
113 : Vector3d p(4*kpc,-2.5*kpc,-0.85*kpc);
114 :
115 : //testing HI component
116 1 : EXPECT_NEAR(n.getHIPlanedensity(p),162597,1); //uncertaincy of 1e-6 cm^-3
117 1 : EXPECT_NEAR(n.getHIScaleheight(p),0.3109*kpc,0.1*pc);
118 1 : EXPECT_NEAR(n.getHIDensity(p),914,1); //uncertainc 1e-6 cm^-3
119 :
120 : //testing H2 compontent
121 1 : EXPECT_NEAR(n.getH2Planedensity(p),741999,1); //uncertaincy of 1e-6 cm^-3
122 1 : EXPECT_NEAR(n.getH2Scaleheight(p),88.2*pc,0.1*pc);
123 1 : EXPECT_NEAR(n.getH2Density(p),0,1);
124 :
125 : //testing total Density
126 1 : EXPECT_NEAR(n.getDensity(p),914,2); //double uncertaincy for both type รก 1cm^-3
127 1 : EXPECT_NEAR(n.getNucleonDensity(p),914,2); // 914 + 0*2 factor 2 for molecular hydrogen
128 :
129 : //second position for testing density
130 : p = Vector3d(50*pc,100*pc,10*pc);
131 :
132 : //testing HI component
133 1 : EXPECT_NEAR(n.getHIPlanedensity(p),543249,1);
134 1 : EXPECT_NEAR(n.getHIScaleheight(p),125.6*pc,0.1*pc);
135 1 : EXPECT_NEAR(n.getHIDensity(p),540867,1);
136 :
137 : //testing H2 component
138 1 : EXPECT_NEAR(n.getH2Planedensity(p),10556748,1);
139 1 : EXPECT_NEAR(n.getH2Scaleheight(p),57.2*pc,0.1*pc);
140 1 : EXPECT_NEAR(n.getH2Density(p),10335137,1);
141 :
142 : //testing total Density
143 1 : EXPECT_NEAR(n.getDensity(p),10876004,2);
144 1 : EXPECT_NEAR(n.getNucleonDensity(p),21211141,2); // factor 2 in molecular hydrogen
145 :
146 : //test set type function
147 1 : n.setIsForHI(false);
148 1 : EXPECT_FALSE(n.getIsForHI());
149 1 : EXPECT_TRUE(n.getIsForH2());
150 :
151 1 : n.setIsForH2(false);
152 1 : EXPECT_FALSE(n.getIsForHI());
153 1 : EXPECT_FALSE(n.getIsForH2());
154 :
155 : //check if density output is zero if all density-types are deaktivated (should give warning in log-file)
156 1 : ::testing::internal::CaptureStderr();
157 1 : EXPECT_DOUBLE_EQ(n.getDensity(p),0);
158 1 : EXPECT_DOUBLE_EQ(n.getNucleonDensity(p),0);
159 1 : std::string captured = testing::internal::GetCapturedStderr();
160 1 : EXPECT_NE(captured.find("WARNING"), std::string::npos);
161 :
162 1 : }
163 :
164 1 : TEST(testFerriere, checkValueAtCertainPoints) {
165 1 : ::testing::internal::CaptureStderr();
166 : Ferriere n;
167 :
168 : //check type information
169 :
170 1 : EXPECT_TRUE(n.getIsForHI());
171 1 : EXPECT_TRUE(n.getIsForHII());
172 1 : EXPECT_TRUE(n.getIsForH2());
173 :
174 : //testing density in inner Ring (R <= 3*kpc)
175 : Vector3d p(60*pc,-60*pc,-20*pc); //testing position in region of CMZ
176 :
177 : //test CMZ Trafo
178 : Vector3d Trafo;
179 1 : Trafo = n.CMZTransformation(p);
180 1 : EXPECT_NEAR(Trafo.x,5.9767*pc,1e-4*pc);
181 1 : EXPECT_NEAR(Trafo.y,12.8171*pc,1e-4*pc);
182 1 : EXPECT_DOUBLE_EQ(Trafo.z,p.z); //no transformation in z component
183 :
184 : //test DISK Trafo
185 1 : Trafo = n.DiskTransformation(p);
186 1 : EXPECT_NEAR(Trafo.x,11.0660*pc,1e-4*pc);
187 1 : EXPECT_NEAR(Trafo.y,82.5860*pc,1e-4*pc);
188 1 : EXPECT_NEAR(Trafo.z,-25.6338*pc,1e-4*pc);
189 :
190 : //testing density
191 1 : EXPECT_NEAR(n.getHIDensity(p),6237723,1); //uncertaincy 1e-6 cm^-3
192 1 : EXPECT_NEAR(n.getH2Density(p),35484825,1);
193 1 : EXPECT_NEAR(n.getHIIDensity(p),6243793,1);
194 1 : EXPECT_NEAR(n.getDensity(p),47966341,1);
195 1 : EXPECT_NEAR(n.getNucleonDensity(p),83451166,2); //factor 2 in molecular hydrogen; double uncertaincy
196 :
197 : Vector3d p2(-500*pc,-900*pc,35*pc); //testing position in region of the DISK
198 1 : EXPECT_NEAR(n.getHIIDensity(p2),48190,1);
199 1 : EXPECT_NEAR(n.getHIDensity(p2),5,1);
200 1 : EXPECT_NEAR(n.getH2Density(p2),0,1);
201 1 : EXPECT_NEAR(n.getDensity(p2),48195,1);
202 1 : EXPECT_NEAR(n.getNucleonDensity(p2),48195,1); // no H2 component -> no difference between density and nucleon-density
203 :
204 : //testing the outer region R>3kpc
205 : Vector3d p3(5*kpc,4*kpc,-29*pc); //testing position with 3kpc < R < R_sun
206 1 : EXPECT_NEAR(n.getHIDensity(p3),540607,1);
207 1 : EXPECT_NEAR(n.getHIIDensity(p3),66495 ,1);
208 1 : EXPECT_NEAR(n.getH2Density(p3),2492685,1);
209 1 : EXPECT_NEAR(n.getDensity(p3),3099787,1);
210 1 : EXPECT_NEAR(n.getNucleonDensity(p3), 5592472,1);
211 :
212 : Vector3d p4(10*kpc,2*kpc,50*pc); //testing position with R > R_sun
213 1 : EXPECT_NEAR(n.getHIDensity(p4),431294,1);
214 1 : EXPECT_NEAR(n.getHIIDensity(p4),22109,1);
215 1 : EXPECT_NEAR(n.getH2Density(p4),54099,1);
216 1 : EXPECT_NEAR(n.getDensity(p4),507502,1);
217 1 : EXPECT_NEAR(n.getNucleonDensity(p4),561601,1);
218 :
219 : //test get/set type function
220 1 : n.setIsForHI(false);
221 1 : EXPECT_FALSE(n.getIsForHI());
222 1 : EXPECT_TRUE(n.getIsForHII());
223 1 : EXPECT_TRUE(n.getIsForH2());
224 :
225 1 : n.setIsForHII(false);
226 1 : EXPECT_FALSE(n.getIsForHI());
227 1 : EXPECT_FALSE(n.getIsForHII());
228 1 : EXPECT_TRUE(n.getIsForH2());
229 :
230 1 : n.setIsForH2(false);
231 1 : EXPECT_FALSE(n.getIsForHI());
232 1 : EXPECT_FALSE(n.getIsForHII());
233 1 : EXPECT_FALSE(n.getIsForH2());
234 :
235 : //check if density is set to zero if all types are deactivated (should give warning in log-file)
236 1 : EXPECT_DOUBLE_EQ(n.getDensity(p),0);
237 1 : EXPECT_DOUBLE_EQ(n.getNucleonDensity(p),0);
238 1 : std::string captured = testing::internal::GetCapturedStderr();
239 1 : EXPECT_NE(captured.find("WARNING"), std::string::npos);
240 1 : }
241 :
242 1 : TEST(testGridDensity, SimpleTest) {
243 1 : ref_ptr<Grid1f> grid = new Grid1f(Vector3d(0.), 1, 1, 1, 1.);
244 1 : DensityGrid dens = DensityGrid(grid, true, false, false);
245 :
246 : // check active types
247 1 : EXPECT_TRUE(dens.getIsForHI());
248 1 : EXPECT_FALSE(dens.getIsForHII());
249 1 : EXPECT_FALSE(dens.getIsForH2());
250 :
251 : // check set function
252 1 : dens.setIsForH2(true);
253 1 : EXPECT_TRUE(dens.getIsForH2());
254 :
255 1 : dens.setIsForHII(true);
256 1 : EXPECT_TRUE(dens.getIsForHII());
257 :
258 1 : dens.setIsForHI(false);
259 1 : EXPECT_FALSE(dens.getIsForHI());
260 1 : }
261 :
262 1 : TEST(testGridDensity, testRetrunValue) {
263 : size_t Nx = 5;
264 : size_t Ny = 8;
265 : size_t Nz = 10;
266 : double spacing = 2.0;
267 : Vector3d origin(1., 2., 3.);
268 :
269 1 : ref_ptr<Grid1f> grid = new Grid1f(origin, Nx, Ny, Nz, spacing);
270 :
271 : // set some values for the grid
272 1 : grid->get(3, 2, 4) = 5;
273 1 : grid->get(3, 2, 5) = 12;
274 1 : grid->get(2, 3, 4) = 6;
275 :
276 2 : DensityGrid dens = DensityGrid(grid, true, false, false);
277 :
278 : // a point in the region where values are defined for the grid.
279 : Vector3d position = origin + Vector3d(2.2, 2.8, 4.1) * spacing;
280 1 : double valueFromGrid = grid->interpolate(position);
281 1 : double nHI = dens.getHIDensity(position);
282 1 : double nHII = dens.getHIIDensity(position);
283 1 : double nH2 = dens.getH2Density(position);
284 1 : double nNucleon = dens.getNucleonDensity(position);
285 1 : double nTotal = dens.getDensity(position);
286 :
287 : // Check for values
288 1 : EXPECT_DOUBLE_EQ(valueFromGrid, nHI);
289 1 : EXPECT_DOUBLE_EQ(nHII, 0); // HII is set to false and should be 0
290 1 : EXPECT_DOUBLE_EQ(nH2, 0); // H2 is set to false and should be 0
291 1 : EXPECT_DOUBLE_EQ(nNucleon, valueFromGrid);
292 1 : EXPECT_DOUBLE_EQ(nTotal, valueFromGrid);
293 1 : }
294 :
295 :
296 : } //namespace crpropa
|
