diff --git a/main.cpp b/main.cpp
index f10ef51095708ce8e187e68e310c5aa721de5ca7..c66983d48ee620d879e46f7b06fa462798404f10 100644
--- a/main.cpp
+++ b/main.cpp
@@ -3,6 +3,7 @@
#include <cstring>
#include "sound_column.h"
#include <execution>
+#include <algorithm>
#define DEFAULT_DISTANCE 20'000.0f
#define DEFAULT_DEPTH 5'000.0f
@@ -67,15 +68,15 @@ int run_simulation(float r_max, float z_max, float f_min, float f_max){
auto nf = size_t(f_max / f_min);
float delta_r = r_max/10'000;
+ std::vector<sound_column*> column(nf, nullptr);
- std::vector<sound_column*> column;
- column.reserve(nf);
+ //std::ranges::iota_view indexes((size_t) 0, nf);
- std::ranges::iota_view indexes((size_t) 0, nf);
- std::for_each(std::execution::par_unseq, indexes.begin(), indexes.end(),
- [&](size_t i) {column[i] = new sound_column(delta_r, (size_t) 1.10 * z_max*(4*(i+1)*f_min)/WATER_SPEED, 7, r_max, z_max, (i+1)*f_min, WATER_SPEED);
- column[i]->init_source(gaussian_init);});
- std::for_each(std::execution::par_unseq, indexes.begin(), indexes.end(),
- [&](size_t i) {while (column[i]->step_r(media));});
+ std::for_each(std::execution::par_unseq, column.begin(), column.end(),
+ [&](sound_column*& c) {int i = &c - &column[0];
+ c = new sound_column(delta_r, (size_t) 1.10 * z_max*(4*(i+1)*f_min)/WATER_SPEED, 7, r_max, z_max, (i+1)*f_min, WATER_SPEED);
+ c->init_source(gaussian_init);});
+ std::for_each(std::execution::par_unseq, column.begin(), column.end(),
+ [&](sound_column*& c) {while (c->step_r(media));});
return 0;
-}
\ No newline at end of file
+}
diff --git a/sound_column.cpp b/sound_column.cpp
index 85d3e72e688744e55023a084d5a4046d1cc112f3..03724ad53080fb73db9de52c5dedbaabe84f91ec 100644
--- a/sound_column.cpp
+++ b/sound_column.cpp
@@ -14,13 +14,13 @@ sound_column::sound_column(float dr, size_t nz, size_t np, float dm, float zm, f
n_p=np;
freq=f;
c0=speed;
- pressures.reserve(nz);
- _u_.reserve(nz-1);
- _l_.reserve(nz-1);
- _d_.reserve(nz);
- n_next.reserve(nz);
- fapbj.reserve(np);
- fbmaj.reserve(np);
+ pressures.resize(nz);
+ _u_.resize(nz-1);
+ _l_.resize(nz-1);
+ _d_.resize(nz);
+ n_next.resize(nz);
+ fapbj.resize(np);
+ fbmaj.resize(np);
for(int i=0; i < np; i++) {
faj = std::complex<float>(0, 4* M_1_PIf32 * freq / c0) / float((2.f * n_p + 1.f) *
powf(sinf(float(i+1) * M_1_PIf32 / (2.f * n_p + 1.f)), 2) * delta_r);
@@ -36,34 +36,35 @@ sound_column::~sound_column() = default;
bool sound_column::step_r(float (*fun)(float, float)) {
if (r >= d_max)
return true;
- std::ranges::iota_view z_ind((size_t) 0, n_z);
+ //std::ranges::iota_view z_ind((size_t) 0, n_z);
r += delta_r;
- std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
- [&](size_t i) {n_next[i] = (*fun)(r, z_max/n_z*i);});
+ std::for_each(std::execution::par_unseq, n_next.begin(), n_next.end(),
+ [&](float& n) {int i = &n - &n_next[0]; n = (*fun)(r, z_max/n_z*i);});
for(int i=0; i < n_p;i++){
_step_r_(i);
}
// Hank
- std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
- [&](size_t i) {pressures[i] *= 2.f / (M_1_PIf32 * M_1_PIf32 * freq / c0 * r) *
+ std::for_each(std::execution::par_unseq, pressures.begin(), pressures.end(),
+ [&](std::complex<float>& p) {p *= static_cast<std::complex<float>>(2.f / (M_1_PIf32 * M_1_PIf32 * freq / c0 * r)) *
std::exp(std::complex<float>(0, 2.f* M_1_PIf32 * freq / c0 * r - M_1_PIf32 / 4));});
return false;
}
void sound_column::init_source(std::complex<float> (*fun)(float, float)) {
- std::ranges::iota_view z_ind((size_t) 0, n_z);
- std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
- [&](size_t i) {pressures[i] = (*fun)(z_max/n_z*i, 2*M_1_PIf32 * freq / c0);});
+// std::ranges::iota_view z_ind((size_t) 0, n_z);
+ std::for_each(std::execution::par_unseq, pressures.begin(), pressures.end(),
+ [&](std::complex<float>& p) {int i = &p - &pressures[0];
+ p = (*fun)(z_max/n_z*i, 2*M_1_PIf32 * freq / c0);});
}
void sound_column::_step_r_(int ip) {
- std::ranges::iota_view z_ind((size_t) 0, n_z);
+// std::ranges::iota_view z_ind((size_t) 0, n_z);
//tri diag init
std::fill(std::execution::par_unseq, _u_.begin(), _u_.end(), fbmaj[ip]/h2);
std::fill(std::execution::par_unseq, _l_.begin(), _l_.end(), fbmaj[ip]/h2);
- std::for_each(std::execution::par_unseq, z_ind.begin(), z_ind.end(),
- [&](size_t i) {_d_[i] = 4* M_1_PIf32 * freq / c0 + fbmaj[ip] /
+ std::for_each(std::execution::par_unseq, _d_.begin(), _d_.end(),
+ [&](std::complex<float>& d) {int i = &d - &_d_[0];d = static_cast<std::complex<float>>(4* M_1_PIf32 * freq / c0) + fbmaj[ip] /
static_cast<std::complex<float>>(-2.f/h2 + powf(4* M_1_PIf32 * freq / c0, 2) * (powf(n_next[i], 2) - 1.f));});
// first matmul
diff --git a/sound_column.h b/sound_column.h
index aa44e654cfbc3d3f42fa2baec0e59cf72c5215d3..23049d87a0b584eb9849c1306543cb1274c6a78a 100644
--- a/sound_column.h
+++ b/sound_column.h
@@ -6,8 +6,9 @@
#define PARABOLIC___SOUND_COLUMN_H
#include <vector>
#include <complex>
-#include <ranges>
+//#include <ranges>
#include <execution>
+#include <algorithm>
# define M_1_PIf32 (0.318309886183790671537767526745028724)
@@ -30,6 +31,7 @@ private:
std::vector<std::complex<float>> _l_;
std::vector<std::complex<float>> _d_;
std::vector<float> n_next;
+ std::complex<float> (*fun)(float, float) source_fun;
void _step_r_(int ip);
@@ -38,6 +40,7 @@ public:
~sound_column();
bool step_r(float (*fun)(float, float));
void init_source(std::complex<float> (*fun)(float, float));
+ void init_source();
};