Line data Source code
1 : #include "crpropa/GridTools.h"
2 : #include "crpropa/magneticField/MagneticField.h"
3 :
4 : #include <fstream>
5 : #include <sstream>
6 :
7 : namespace crpropa {
8 :
9 0 : void scaleGrid(ref_ptr<Grid1f> grid, double a) {
10 0 : for (int ix = 0; ix < grid->getNx(); ix++)
11 0 : for (int iy = 0; iy < grid->getNy(); iy++)
12 0 : for (int iz = 0; iz < grid->getNz(); iz++)
13 0 : grid->get(ix, iy, iz) *= a;
14 0 : }
15 :
16 11 : void scaleGrid(ref_ptr<Grid3f> grid, double a) {
17 654 : for (int ix = 0; ix < grid->getNx(); ix++)
18 41612 : for (int iy = 0; iy < grid->getNy(); iy++)
19 2662436 : for (int iz = 0; iz < grid->getNz(); iz++)
20 2621467 : grid->get(ix, iy, iz) *= a;
21 11 : }
22 :
23 1 : Vector3f meanFieldVector(ref_ptr<Grid3f> grid) {
24 : size_t Nx = grid->getNx();
25 : size_t Ny = grid->getNy();
26 : size_t Nz = grid->getNz();
27 : Vector3f mean(0.);
28 65 : for (int ix = 0; ix < Nx; ix++)
29 4160 : for (int iy = 0; iy < Ny; iy++)
30 266240 : for (int iz = 0; iz < Nz; iz++)
31 : mean += grid->get(ix, iy, iz);
32 1 : return mean / Nx / Ny / Nz;
33 : }
34 :
35 0 : double meanFieldStrength(ref_ptr<Grid3f> grid) {
36 : size_t Nx = grid->getNx();
37 : size_t Ny = grid->getNy();
38 : size_t Nz = grid->getNz();
39 : double mean = 0;
40 0 : for (int ix = 0; ix < Nx; ix++)
41 0 : for (int iy = 0; iy < Ny; iy++)
42 0 : for (int iz = 0; iz < Nz; iz++)
43 0 : mean += grid->get(ix, iy, iz).getR();
44 0 : return mean / Nx / Ny / Nz;
45 : }
46 :
47 0 : double meanFieldStrength(ref_ptr<Grid1f> grid) {
48 : size_t Nx = grid->getNx();
49 : size_t Ny = grid->getNy();
50 : size_t Nz = grid->getNz();
51 : double mean = 0;
52 0 : for (int ix = 0; ix < Nx; ix++)
53 0 : for (int iy = 0; iy < Ny; iy++)
54 0 : for (int iz = 0; iz < Nz; iz++)
55 0 : mean += grid->get(ix, iy, iz);
56 0 : return mean / Nx / Ny / Nz;
57 : }
58 :
59 11 : double rmsFieldStrength(ref_ptr<Grid3f> grid) {
60 : size_t Nx = grid->getNx();
61 : size_t Ny = grid->getNy();
62 : size_t Nz = grid->getNz();
63 : double sumV2 = 0;
64 715 : for (int ix = 0; ix < Nx; ix++)
65 45760 : for (int iy = 0; iy < Ny; iy++)
66 2928640 : for (int iz = 0; iz < Nz; iz++)
67 2883584 : sumV2 += grid->get(ix, iy, iz).getR2();
68 11 : return std::sqrt(sumV2 / Nx / Ny / Nz);
69 : }
70 :
71 0 : double rmsFieldStrength(ref_ptr<Grid1f> grid) {
72 : size_t Nx = grid->getNx();
73 : size_t Ny = grid->getNy();
74 : size_t Nz = grid->getNz();
75 : double sumV2 = 0;
76 0 : for (int ix = 0; ix < Nx; ix++)
77 0 : for (int iy = 0; iy < Ny; iy++)
78 0 : for (int iz = 0; iz < Nz; iz++)
79 0 : sumV2 += pow(grid->get(ix, iy, iz), 2);
80 0 : return std::sqrt(sumV2 / Nx / Ny / Nz);
81 : }
82 :
83 0 : std::array<float, 3> rmsFieldStrengthPerAxis(ref_ptr<Grid3f> grid) {
84 : size_t Nx = grid->getNx();
85 : size_t Ny = grid->getNy();
86 : size_t Nz = grid->getNz();
87 : float sumV2_x = 0;
88 : float sumV2_y = 0;
89 : float sumV2_z = 0;
90 0 : for (int ix = 0; ix < Nx; ix++)
91 0 : for (int iy = 0; iy < Ny; iy++)
92 0 : for (int iz = 0; iz < Nz; iz++) {
93 0 : sumV2_x += pow(grid->get(ix, iy, iz).x, 2);
94 0 : sumV2_y += pow(grid->get(ix, iy, iz).y, 2);
95 0 : sumV2_z += pow(grid->get(ix, iy, iz).z, 2);
96 : }
97 : return {
98 0 : std::sqrt(sumV2_x / Nx / Ny / Nz),
99 0 : std::sqrt(sumV2_y / Nx / Ny / Nz),
100 0 : std::sqrt(sumV2_z / Nx / Ny / Nz)
101 0 : };
102 : }
103 :
104 0 : void fromMagneticField(ref_ptr<Grid3f> grid, ref_ptr<MagneticField> field) {
105 : Vector3d origin = grid->getOrigin();
106 : Vector3d spacing = grid->getSpacing();
107 : size_t Nx = grid->getNx();
108 : size_t Ny = grid->getNy();
109 : size_t Nz = grid->getNz();
110 0 : for (size_t ix = 0; ix < Nx; ix++)
111 0 : for (size_t iy = 0; iy < Ny; iy++)
112 0 : for (size_t iz = 0; iz < Nz; iz++) {
113 0 : Vector3d pos = Vector3d(double(ix) + 0.5, double(iy) + 0.5, double(iz) + 0.5) * spacing + origin;
114 0 : Vector3d B = field->getField(pos);
115 : grid->get(ix, iy, iz) = B;
116 : }
117 0 : }
118 :
119 0 : void fromMagneticFieldStrength(ref_ptr<Grid1f> grid, ref_ptr<MagneticField> field) {
120 : Vector3d origin = grid->getOrigin();
121 : Vector3d spacing = grid->getSpacing();
122 : size_t Nx = grid->getNx();
123 : size_t Ny = grid->getNy();
124 : size_t Nz = grid->getNz();
125 0 : for (size_t ix = 0; ix < Nx; ix++)
126 0 : for (size_t iy = 0; iy < Ny; iy++)
127 0 : for (size_t iz = 0; iz < Nz; iz++) {
128 0 : Vector3d pos = Vector3d(double(ix) + 0.5, double(iy) + 0.5, double(iz) + 0.5) * spacing + origin;
129 0 : double s = field->getField(pos).getR();
130 0 : grid->get(ix, iy, iz) = s;
131 : }
132 0 : }
133 :
134 1 : void loadGrid(ref_ptr<Grid3f> grid, std::string filename, double c) {
135 1 : std::ifstream fin(filename.c_str(), std::ios::binary);
136 1 : if (!fin) {
137 0 : std::stringstream ss;
138 0 : ss << "load Grid3f: " << filename << " not found";
139 0 : throw std::runtime_error(ss.str());
140 0 : }
141 :
142 : // get length of file and compare to size of grid
143 1 : fin.seekg(0, fin.end);
144 1 : size_t length = fin.tellg() / sizeof(float);
145 1 : fin.seekg (0, fin.beg);
146 :
147 : size_t nx = grid->getNx();
148 : size_t ny = grid->getNy();
149 : size_t nz = grid->getNz();
150 :
151 1 : if (length != (3 * nx * ny * nz))
152 0 : throw std::runtime_error("loadGrid: file and grid size do not match");
153 :
154 4 : for (int ix = 0; ix < grid->getNx(); ix++) {
155 12 : for (int iy = 0; iy < grid->getNy(); iy++) {
156 36 : for (int iz = 0; iz < grid->getNz(); iz++) {
157 : Vector3f &b = grid->get(ix, iy, iz);
158 27 : fin.read((char*) &(b.x), sizeof(float));
159 27 : fin.read((char*) &(b.y), sizeof(float));
160 27 : fin.read((char*) &(b.z), sizeof(float));
161 27 : b *= c;
162 : }
163 : }
164 : }
165 1 : fin.close();
166 1 : }
167 :
168 0 : void loadGrid(ref_ptr<Grid1f> grid, std::string filename, double c) {
169 0 : std::ifstream fin(filename.c_str(), std::ios::binary);
170 0 : if (!fin) {
171 0 : std::stringstream ss;
172 0 : ss << "load Grid1f: " << filename << " not found";
173 0 : throw std::runtime_error(ss.str());
174 0 : }
175 :
176 : // get length of file and compare to size of grid
177 0 : fin.seekg(0, fin.end);
178 0 : size_t length = fin.tellg() / sizeof(float);
179 0 : fin.seekg (0, fin.beg);
180 :
181 : size_t nx = grid->getNx();
182 : size_t ny = grid->getNy();
183 : size_t nz = grid->getNz();
184 :
185 0 : if (length != (nx * ny * nz))
186 0 : throw std::runtime_error("loadGrid: file and grid size do not match");
187 :
188 0 : for (int ix = 0; ix < nx; ix++) {
189 0 : for (int iy = 0; iy < ny; iy++) {
190 0 : for (int iz = 0; iz < nz; iz++) {
191 : float &b = grid->get(ix, iy, iz);
192 0 : fin.read((char*) &b, sizeof(float));
193 0 : b *= c;
194 : }
195 : }
196 : }
197 0 : fin.close();
198 0 : }
199 :
200 1 : void dumpGrid(ref_ptr<Grid3f> grid, std::string filename, double c) {
201 1 : std::ofstream fout(filename.c_str(), std::ios::binary);
202 1 : if (!fout) {
203 0 : std::stringstream ss;
204 0 : ss << "dump Grid3f: " << filename << " not found";
205 0 : throw std::runtime_error(ss.str());
206 0 : }
207 4 : for (int ix = 0; ix < grid->getNx(); ix++) {
208 12 : for (int iy = 0; iy < grid->getNy(); iy++) {
209 36 : for (int iz = 0; iz < grid->getNz(); iz++) {
210 27 : Vector3f b = grid->get(ix, iy, iz) * c;
211 27 : fout.write((char*) &(b.x), sizeof(float));
212 27 : fout.write((char*) &(b.y), sizeof(float));
213 27 : fout.write((char*) &(b.z), sizeof(float));
214 : }
215 : }
216 : }
217 1 : fout.close();
218 1 : }
219 :
220 0 : void dumpGrid(ref_ptr<Grid1f> grid, std::string filename, double c) {
221 0 : std::ofstream fout(filename.c_str(), std::ios::binary);
222 0 : if (!fout) {
223 0 : std::stringstream ss;
224 0 : ss << "dump Grid1f: " << filename << " not found";
225 0 : throw std::runtime_error(ss.str());
226 0 : }
227 0 : for (int ix = 0; ix < grid->getNx(); ix++) {
228 0 : for (int iy = 0; iy < grid->getNy(); iy++) {
229 0 : for (int iz = 0; iz < grid->getNz(); iz++) {
230 0 : float b = grid->get(ix, iy, iz) * c;
231 0 : fout.write((char*) &b, sizeof(float));
232 : }
233 : }
234 : }
235 0 : fout.close();
236 0 : }
237 :
238 1 : void loadGridFromTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
239 1 : std::ifstream fin(filename.c_str());
240 1 : if (!fin) {
241 0 : std::stringstream ss;
242 0 : ss << "load Grid3f: " << filename << " not found";
243 0 : throw std::runtime_error(ss.str());
244 0 : }
245 : // skip header lines
246 1 : while (fin.peek() == '#')
247 0 : fin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
248 :
249 4 : for (int ix = 0; ix < grid->getNx(); ix++) {
250 12 : for (int iy = 0; iy < grid->getNy(); iy++) {
251 36 : for (int iz = 0; iz < grid->getNz(); iz++) {
252 : Vector3f &b = grid->get(ix, iy, iz);
253 27 : fin >> b.x >> b.y >> b.z;
254 27 : b *= c;
255 27 : if (fin.eof())
256 0 : throw std::runtime_error("load Grid3f: file too short");
257 : }
258 : }
259 : }
260 1 : fin.close();
261 1 : }
262 :
263 0 : void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
264 0 : std::ifstream fin(filename.c_str());
265 0 : if (!fin) {
266 0 : std::stringstream ss;
267 0 : ss << "load Grid1f: " << filename << " not found";
268 0 : throw std::runtime_error(ss.str());
269 0 : }
270 : // skip header lines
271 0 : while (fin.peek() == '#')
272 0 : fin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
273 :
274 0 : for (int ix = 0; ix < grid->getNx(); ix++) {
275 0 : for (int iy = 0; iy < grid->getNy(); iy++) {
276 0 : for (int iz = 0; iz < grid->getNz(); iz++) {
277 : float &b = grid->get(ix, iy, iz);
278 : fin >> b;
279 0 : b *= c;
280 0 : if (fin.eof())
281 0 : throw std::runtime_error("load Grid1f: file too short");
282 : }
283 : }
284 : }
285 0 : fin.close();
286 0 : }
287 :
288 1 : void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
289 1 : std::ofstream fout(filename.c_str());
290 1 : if (!fout) {
291 0 : std::stringstream ss;
292 0 : ss << "dump Grid3f: " << filename << " not found";
293 0 : throw std::runtime_error(ss.str());
294 0 : }
295 4 : for (int ix = 0; ix < grid->getNx(); ix++) {
296 12 : for (int iy = 0; iy < grid->getNy(); iy++) {
297 36 : for (int iz = 0; iz < grid->getNz(); iz++) {
298 27 : Vector3f b = grid->get(ix, iy, iz) * c;
299 27 : fout << b << "\n";
300 : }
301 : }
302 : }
303 1 : fout.close();
304 1 : }
305 :
306 0 : void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
307 0 : std::ofstream fout(filename.c_str());
308 0 : if (!fout) {
309 0 : std::stringstream ss;
310 0 : ss << "dump Grid1f: " << filename << " not found";
311 0 : throw std::runtime_error(ss.str());
312 0 : }
313 0 : for (int ix = 0; ix < grid->getNx(); ix++) {
314 0 : for (int iy = 0; iy < grid->getNy(); iy++) {
315 0 : for (int iz = 0; iz < grid->getNz(); iz++) {
316 0 : float b = grid->get(ix, iy, iz) * c;
317 0 : fout << b << "\n";
318 : }
319 : }
320 : }
321 0 : fout.close();
322 0 : }
323 :
324 : #ifdef CRPROPA_HAVE_FFTW3F
325 :
326 0 : std::vector<std::pair<int, float>> gridPowerSpectrum(ref_ptr<Grid3f> grid) {
327 :
328 0 : double rms = rmsFieldStrength(grid);
329 : size_t n = grid->getNx(); // size of array
330 :
331 : // arrays to hold the complex vector components of the B(k)-field
332 : fftwf_complex *Bkx, *Bky, *Bkz;
333 0 : Bkx = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n);
334 0 : Bky = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n);
335 0 : Bkz = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n);
336 :
337 : fftwf_complex *Bx = (fftwf_complex *)Bkx;
338 : fftwf_complex *By = (fftwf_complex *)Bky;
339 : fftwf_complex *Bz = (fftwf_complex *)Bkz;
340 :
341 : // save to temp
342 : int i;
343 0 : for (size_t ix = 0; ix < n; ix++) {
344 0 : for (size_t iy = 0; iy < n; iy++) {
345 0 : for (size_t iz = 0; iz < n; iz++) {
346 0 : i = ix * n * n + iy * n + iz;
347 : Vector3<float> &b = grid->get(ix, iy, iz);
348 0 : Bx[i][0] = b.x / rms;
349 0 : By[i][0] = b.y / rms;
350 0 : Bz[i][0] = b.z / rms;
351 : }
352 : }
353 : }
354 :
355 : // in-place, real to complex, inverse Fourier transformation on each component
356 : // note that the last elements of B(x) are unused now
357 : fftwf_plan plan_x =
358 0 : fftwf_plan_dft_3d(n, n, n, Bx, Bkx, FFTW_FORWARD, FFTW_ESTIMATE);
359 0 : fftwf_execute(plan_x);
360 0 : fftwf_destroy_plan(plan_x);
361 :
362 : fftwf_plan plan_y =
363 0 : fftwf_plan_dft_3d(n, n, n, By, Bky, FFTW_FORWARD, FFTW_ESTIMATE);
364 0 : fftwf_execute(plan_y);
365 0 : fftwf_destroy_plan(plan_y);
366 :
367 : fftwf_plan plan_z =
368 0 : fftwf_plan_dft_3d(n, n, n, Bz, Bkz, FFTW_FORWARD, FFTW_ESTIMATE);
369 0 : fftwf_execute(plan_z);
370 0 : fftwf_destroy_plan(plan_z);
371 :
372 : float power;
373 : std::map<size_t, std::pair<float, int>> spectrum;
374 : int k;
375 :
376 0 : for (size_t ix = 0; ix < n; ix++) {
377 0 : for (size_t iy = 0; iy < n; iy++) {
378 0 : for (size_t iz = 0; iz < n; iz++) {
379 0 : i = ix * n * n + iy * n + iz;
380 0 : k = static_cast<int>(
381 0 : std::floor(std::sqrt(ix * ix + iy * iy + iz * iz)));
382 0 : if (k > n / 2. || k == 0)
383 0 : continue;
384 0 : power = ((Bkx[i][0] * Bkx[i][0] + Bkx[i][1] * Bkx[i][1]) +
385 0 : (Bky[i][0] * Bky[i][0] + Bky[i][1] * Bky[i][1]) +
386 0 : (Bkz[i][0] * Bkz[i][0] + Bkz[i][1] * Bkz[i][1]));
387 0 : if (spectrum.find(k) == spectrum.end()) {
388 0 : spectrum[k].first = power;
389 0 : spectrum[k].second = 1;
390 : } else {
391 0 : spectrum[k].first += power;
392 0 : spectrum[k].second += 1;
393 : }
394 : }
395 : }
396 : }
397 :
398 0 : fftwf_free(Bkx);
399 0 : fftwf_free(Bky);
400 0 : fftwf_free(Bkz);
401 :
402 : std::vector<std::pair<int, float>> points;
403 0 : for (std::map<size_t, std::pair<float, int>>::iterator it = spectrum.begin();
404 0 : it != spectrum.end(); ++it) {
405 0 : points.push_back(
406 0 : std::make_pair(it->first, (it->second).first / (it->second).second));
407 : }
408 :
409 0 : return points;
410 : }
411 :
412 : #endif // CRPROPA_HAVE_FFTW3F
413 :
414 :
415 : } // namespace crpropa
|
