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 2 : void loadGridFromTxt(ref_ptr<Grid3f> grid, std::string filename, double c) {
239 2 : std::ifstream fin(filename.c_str());
240 2 : 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 3 : while (fin.peek() == '#')
247 1 : fin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
248 :
249 8 : for (int ix = 0; ix < grid->getNx(); ix++) {
250 21 : for (int iy = 0; iy < grid->getNy(); iy++) {
251 66 : for (int iz = 0; iz < grid->getNz(); iz++) {
252 : Vector3f &b = grid->get(ix, iy, iz);
253 51 : fin >> b.x >> b.y >> b.z;
254 51 : b *= c;
255 51 : if (fin.eof())
256 0 : throw std::runtime_error("load Grid3f: file too short");
257 : }
258 : }
259 : }
260 2 : fin.close();
261 2 : }
262 :
263 1 : ref_ptr<Grid3f> loadGrid3fFromTxt(std::string filename, double c) {
264 1 : std::ifstream fin(filename.c_str());
265 1 : if (!fin) {
266 0 : std::stringstream ss;
267 0 : ss << "load Grid3f: " << filename << " not found";
268 0 : throw std::runtime_error(ss.str());
269 0 : }
270 :
271 : // search in header lines for GridProperties
272 1 : while (fin.peek() == '#') {
273 : std::string line;
274 1 : std::getline(fin, line);
275 :
276 : // find gridproperties in the header
277 1 : if (line.rfind("GridProperties:") == 2) {
278 : GridProperties gp(Vector3d(0.), 1, 1, 1, 1.); // simple grid properties for default
279 1 : std::stringstream ss(line);
280 :
281 : // skip first names and check type
282 : std::string name, type;
283 1 : ss >> name >> name >> name >> type;
284 1 : if (type != "Grid3f")
285 0 : throw std::runtime_error("Tried to load Grid3f, but Gridproperties assume grid type " + type);
286 :
287 : // grid origin
288 : double x, y, z;
289 1 : ss >> name >> x >> y >> z ;
290 : gp.origin = Vector3d(x, y, z);
291 :
292 : // grid size
293 1 : ss >> name >> gp.Nx >> gp.Ny >> gp.Nz;
294 :
295 : // spacing
296 : double dX, dY, dZ;
297 1 : ss >> name >> dX >> dY >> dZ;
298 : gp.spacing = Vector3d(dX, dY, dZ);
299 :
300 : // reflective
301 1 : ss >> name >> gp.reflective;
302 :
303 : // clip volume
304 : bool clip;
305 1 : ss >> name >> clip;
306 1 : gp.setClipVolume(clip);
307 :
308 : // interpolation type
309 1 : ss >> name >> type;
310 1 : if (type == "TRICUBIC")
311 : gp.setInterpolationType(TRICUBIC);
312 1 : else if (type == "NEAREST_NEIGHBOUR")
313 : gp.setInterpolationType(NEAREST_NEIGHBOUR);
314 : else
315 : gp.setInterpolationType(TRILINEAR);
316 :
317 : // create new grid
318 1 : ref_ptr<Grid3f> grid = new Grid3f(gp);
319 1 : fin.close();
320 :
321 : // load data for grid
322 2 : loadGridFromTxt(grid, filename, c);
323 :
324 1 : return grid;
325 1 : }
326 : }
327 0 : throw std::runtime_error("could not find GridProperties in file " + filename);
328 1 : }
329 :
330 :
331 1 : void loadGridFromTxt(ref_ptr<Grid1f> grid, std::string filename, double c) {
332 1 : std::ifstream fin(filename.c_str());
333 1 : if (!fin) {
334 0 : std::stringstream ss;
335 0 : ss << "load Grid1f: " << filename << " not found";
336 0 : throw std::runtime_error(ss.str());
337 0 : }
338 :
339 : // skip header lines
340 2 : while (fin.peek() == '#')
341 1 : fin.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
342 :
343 4 : for (int ix = 0; ix < grid->getNx(); ix++) {
344 9 : for (int iy = 0; iy < grid->getNy(); iy++) {
345 30 : for (int iz = 0; iz < grid->getNz(); iz++) {
346 : float &b = grid->get(ix, iy, iz);
347 : fin >> b;
348 24 : b *= c;
349 24 : if (fin.eof())
350 0 : throw std::runtime_error("load Grid1f: file too short");
351 : }
352 : }
353 : }
354 1 : fin.close();
355 1 : }
356 :
357 1 : ref_ptr<Grid1f> loadGrid1fFromTxt(std::string filename, double c) {
358 1 : std::ifstream fin(filename.c_str());
359 1 : if (!fin) {
360 0 : std::stringstream ss;
361 0 : ss << "load Grid1f: " << filename << " not found";
362 0 : throw std::runtime_error(ss.str());
363 0 : }
364 :
365 : // search in header lines for GridProperties
366 1 : while (fin.peek() == '#') {
367 : std::string line;
368 1 : std::getline(fin, line);
369 :
370 : // find gridproperties in the header
371 1 : if (line.rfind("GridProperties:") == 2) {
372 : GridProperties gp(Vector3d(0.), 1, 1, 1, 1.); // simple grid properties for default
373 1 : std::stringstream ss(line);
374 :
375 : // skip first names and check type
376 : std::string name, type;
377 1 : ss >> name >> name >> name >> type;
378 1 : if (type != "Grid1f")
379 0 : throw std::runtime_error("Tried to load Grid1f, but Gridproperties assume grid type " + type);
380 :
381 : // grid origin
382 : double x, y, z;
383 1 : ss >> name >> x >> y >> z ;
384 : gp.origin = Vector3d(x, y, z);
385 :
386 : // grid size
387 1 : ss >> name >> gp.Nx >> gp.Ny >> gp.Nz;
388 :
389 : // spacing
390 : double dX, dY, dZ;
391 1 : ss >> name >> dX >> dY >> dZ;
392 : gp.spacing = Vector3d(dX, dY, dZ);
393 :
394 : // reflective
395 1 : ss >> name >> gp.reflective;
396 :
397 : // clip volume
398 : bool clip;
399 1 : ss >> name >> clip;
400 1 : gp.setClipVolume(clip);
401 :
402 : // interpolation type
403 1 : ss >> name >> type;
404 1 : if (type == "TRICUBIC")
405 : gp.setInterpolationType(TRICUBIC);
406 0 : else if (type == "NEAREST_NEIGHBOUR")
407 : gp.setInterpolationType(NEAREST_NEIGHBOUR);
408 : else
409 : gp.setInterpolationType(TRILINEAR);
410 :
411 : // create new grid
412 1 : ref_ptr<Grid1f> grid = new Grid1f(gp);
413 1 : fin.close();
414 :
415 : // load data for grid
416 2 : loadGridFromTxt(grid, filename, c);
417 :
418 1 : return grid;
419 1 : }
420 : }
421 0 : throw std::runtime_error("could not find GridProperties in file " + filename);
422 1 : }
423 :
424 2 : void dumpGridToTxt(ref_ptr<Grid3f> grid, std::string filename, double c, bool saveProp) {
425 2 : std::ofstream fout(filename.c_str());
426 2 : if (!fout) {
427 0 : std::stringstream ss;
428 0 : ss << "dump Grid3f: " << filename << " not found";
429 0 : throw std::runtime_error(ss.str());
430 0 : }
431 :
432 : // store the properties in the file as header information
433 2 : if (saveProp) {
434 : fout << "# GridProperties: Type Grid3f"
435 1 : << "\t" << "origin: " << grid -> getOrigin()
436 : << "\t" << "gridsize: " << grid -> getNx() << " " << grid -> getNy() << " " << grid -> getNz()
437 2 : << "\t" << "spacing: " << grid -> getSpacing ()
438 : << "\t" << "reflective: " << grid -> isReflective()
439 : << "\t" << "clipVolume: " << grid -> getClipVolume()
440 3 : << "\t" << "interpolation: " << grid -> getInterpolationTypeName() << "\n";
441 : }
442 :
443 8 : for (int ix = 0; ix < grid->getNx(); ix++) {
444 21 : for (int iy = 0; iy < grid->getNy(); iy++) {
445 66 : for (int iz = 0; iz < grid->getNz(); iz++) {
446 51 : Vector3f b = grid->get(ix, iy, iz) * c;
447 51 : fout << b << "\n";
448 : }
449 : }
450 : }
451 2 : fout.close();
452 2 : }
453 :
454 1 : void dumpGridToTxt(ref_ptr<Grid1f> grid, std::string filename, double c, bool saveProp) {
455 1 : std::ofstream fout(filename.c_str());
456 1 : if (!fout) {
457 0 : std::stringstream ss;
458 0 : ss << "dump Grid1f: " << filename << " not found";
459 0 : throw std::runtime_error(ss.str());
460 0 : }
461 :
462 : // save properties as header information
463 1 : if (saveProp) {
464 : fout << "# GridProperties: Type Grid1f"
465 1 : << "\t" << "origin: " << grid -> getOrigin()
466 : << "\t" << "gridsize: " << grid -> getNx() << " " << grid -> getNy() << " " << grid -> getNz()
467 2 : << "\t" << "spacing: " << grid -> getSpacing ()
468 : << "\t" << "reflective: " << grid -> isReflective()
469 : << "\t" << "clipVolume: " << grid -> getClipVolume()
470 3 : << "\t" << "interpolation: " << grid -> getInterpolationTypeName() << "\n";
471 : }
472 :
473 4 : for (int ix = 0; ix < grid->getNx(); ix++) {
474 9 : for (int iy = 0; iy < grid->getNy(); iy++) {
475 30 : for (int iz = 0; iz < grid->getNz(); iz++) {
476 24 : float b = grid->get(ix, iy, iz) * c;
477 24 : fout << b << "\n";
478 : }
479 : }
480 : }
481 1 : fout.close();
482 1 : }
483 :
484 : #ifdef CRPROPA_HAVE_FFTW3F
485 :
486 0 : std::vector<std::pair<int, float>> gridPowerSpectrum(ref_ptr<Grid3f> grid) {
487 :
488 0 : double rms = rmsFieldStrength(grid);
489 : size_t n = grid->getNx(); // size of array
490 :
491 : // arrays to hold the complex vector components of the B(k)-field
492 : fftwf_complex *Bkx, *Bky, *Bkz;
493 0 : Bkx = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n);
494 0 : Bky = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n);
495 0 : Bkz = (fftwf_complex *)fftwf_malloc(sizeof(fftwf_complex) * n * n * n);
496 :
497 : fftwf_complex *Bx = (fftwf_complex *)Bkx;
498 : fftwf_complex *By = (fftwf_complex *)Bky;
499 : fftwf_complex *Bz = (fftwf_complex *)Bkz;
500 :
501 : // save to temp
502 : int i;
503 0 : for (size_t ix = 0; ix < n; ix++) {
504 0 : for (size_t iy = 0; iy < n; iy++) {
505 0 : for (size_t iz = 0; iz < n; iz++) {
506 0 : i = ix * n * n + iy * n + iz;
507 : Vector3<float> &b = grid->get(ix, iy, iz);
508 0 : Bx[i][0] = b.x / rms;
509 0 : By[i][0] = b.y / rms;
510 0 : Bz[i][0] = b.z / rms;
511 : }
512 : }
513 : }
514 :
515 : // in-place, real to complex, inverse Fourier transformation on each component
516 : // note that the last elements of B(x) are unused now
517 : fftwf_plan plan_x =
518 0 : fftwf_plan_dft_3d(n, n, n, Bx, Bkx, FFTW_FORWARD, FFTW_ESTIMATE);
519 0 : fftwf_execute(plan_x);
520 0 : fftwf_destroy_plan(plan_x);
521 :
522 : fftwf_plan plan_y =
523 0 : fftwf_plan_dft_3d(n, n, n, By, Bky, FFTW_FORWARD, FFTW_ESTIMATE);
524 0 : fftwf_execute(plan_y);
525 0 : fftwf_destroy_plan(plan_y);
526 :
527 : fftwf_plan plan_z =
528 0 : fftwf_plan_dft_3d(n, n, n, Bz, Bkz, FFTW_FORWARD, FFTW_ESTIMATE);
529 0 : fftwf_execute(plan_z);
530 0 : fftwf_destroy_plan(plan_z);
531 :
532 : float power;
533 : std::map<size_t, std::pair<float, int>> spectrum;
534 : int k;
535 :
536 0 : for (size_t ix = 0; ix < n; ix++) {
537 0 : for (size_t iy = 0; iy < n; iy++) {
538 0 : for (size_t iz = 0; iz < n; iz++) {
539 0 : i = ix * n * n + iy * n + iz;
540 0 : k = static_cast<int>(
541 0 : std::floor(std::sqrt(ix * ix + iy * iy + iz * iz)));
542 0 : if (k > n / 2. || k == 0)
543 0 : continue;
544 0 : power = ((Bkx[i][0] * Bkx[i][0] + Bkx[i][1] * Bkx[i][1]) +
545 0 : (Bky[i][0] * Bky[i][0] + Bky[i][1] * Bky[i][1]) +
546 0 : (Bkz[i][0] * Bkz[i][0] + Bkz[i][1] * Bkz[i][1]));
547 0 : if (spectrum.find(k) == spectrum.end()) {
548 0 : spectrum[k].first = power;
549 0 : spectrum[k].second = 1;
550 : } else {
551 0 : spectrum[k].first += power;
552 0 : spectrum[k].second += 1;
553 : }
554 : }
555 : }
556 : }
557 :
558 0 : fftwf_free(Bkx);
559 0 : fftwf_free(Bky);
560 0 : fftwf_free(Bkz);
561 :
562 : std::vector<std::pair<int, float>> points;
563 0 : for (std::map<size_t, std::pair<float, int>>::iterator it = spectrum.begin();
564 0 : it != spectrum.end(); ++it) {
565 0 : points.push_back(
566 0 : std::make_pair(it->first, (it->second).first / (it->second).second));
567 : }
568 :
569 0 : return points;
570 0 : }
571 :
572 : #endif // CRPROPA_HAVE_FFTW3F
573 :
574 :
575 : } // namespace crpropa
|
