From b86625ea13fc173e4d5466930a606f957527363d Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Sat, 9 May 2020 23:29:42 -0400 Subject: [PATCH] Changed the integrate interface to use std::array instead of C-style --- replay.cpp | 94 +++++---------------------------------- scattering_experiment.cpp | 5 +-- 2 files changed, 12 insertions(+), 87 deletions(-) diff --git a/replay.cpp b/replay.cpp index 5d4e40b..00eda41 100644 --- a/replay.cpp +++ b/replay.cpp @@ -36,7 +36,6 @@ public: { return gsl_spline_eval(spline, x, acc); } - private: gsl_interp_accel *acc; gsl_spline *spline; @@ -54,7 +53,7 @@ double3 nfw(const double rho_0, const double b, const double3& pos) { double r2 = pos.squaredNorm(); double r = sqrt(r2); - double tmp = -4*M_PI*G*rho_0*b*b*b*(log(1+r/b) - r/(b+r))/(r2*r); + double tmp = -4*M_PI*G*rho_0*b*b*b*(log1p(r/b) - r/(b+r))/(r2*r); return pos*tmp; } @@ -114,14 +113,12 @@ public: double3 acc_halo = nfw(rho_0_nfw, b_nfw, pos); for (int i=0; i<3; i++) f[3+i] = acc_disk[i] + acc_halo[i]; } - - private: Interp interp[10]; } galaxy("subhalo_411321_parameters_processed.dat"); -extern "C" -void integrate(const double y0[], const double t_max, const double stride_size, double y[]) +void integrate(const std::array y0, const double t_max, const double stride_size, double y[]) +// NOTE we're not responsible to the memory of the output parameter y. Caller should allocate { using namespace boost::numeric::odeint; using Coordinates = std::array; @@ -129,8 +126,8 @@ void integrate(const double y0[], const double t_max, const double stride_size, auto function_wrapper = [](const Coordinates &x, Coordinates &dxdt, const double t) { return galaxy.func(x, dxdt, t); }; const int stride_count = t_max / stride_size; Coordinates y_current; - std::copy(y0, y0+6, begin(y_current)); - std::copy(y0, y0+6, y); + std::copy(begin(y0), end(y0), begin(y_current)); + std::copy(begin(y0), end(y0), y); double t = 0; const double h = std::min(1./4096., stride_size); for (int i=0; i -#include -#include -int mainXXX() +extern "C" +void integrate(const double y0[], const double t_max, const double stride_size, double y[]) { - std::cout << "Hi" << std::endl; - - double x_max = 100; - double v_max = 200; - double t_max = 10.; - - double3 target_pos = {-9.1817914423447837E+03, -2.3161972143758221E+03, 4.2321813890923086E+03}; // pc - double3 target_vel = { 1.6995805137819784E+02, 1.3476658021349482E+02, -1.3342192079702110E+02}; // km/s - - // Unit adjustment - target_pos *= 0.001; // Now in kpc - target_vel *= 1.0226911647958985; // Now in kpc/Gyr - - double tolrance = 2; - double max_dx = target_pos.norm()*tolrance; - double max_dv = target_vel.norm()*tolrance; - - const auto current_time = std::chrono::system_clock::now().time_since_epoch().count(); - const unsigned long big_prime = 840580612873; - const auto pid = getpid(); - const auto seed = current_time + big_prime*pid; // TODO add thread number - - std::default_random_engine generator(seed); - - std::uniform_real_distribution position_distribution(-x_max, x_max); - auto random_position_component = [&]() { return position_distribution(generator); }; - - std::uniform_real_distribution velocity_distribution(-v_max, v_max); - auto random_velocity_component = [&]() { return velocity_distribution(generator); }; - - for (int i=0; i<100000; i++) { - if (i%1000==0) std::cout << i << "\n"; - double y[12], y0[6]; - for (int i=0; i<3; i++) y0[i] = random_position_component(); - for (int i=3; i<6; i++) y0[i] = random_velocity_component(); - - integrate(y0, t_max, t_max, y); - double3 new_pos(y+6); - double3 new_vel(y+9); - - - if (((new_pos - target_pos).norm() < max_dx) && ((new_vel - target_vel).norm() < max_dv)) - std::cout << "*"; - } - - - - - // double y[12]; - // double y0[] = {80,0,0,0,80,0}; - // integrate(y0, 10, 10, y); - // for (int i=0; i<12; i++) std::cout << y[i] << std::endl; - - - - // double y[12]; - // double y0[] = {80,0,0,0,80,0}; - // for (int i=0; i<8000; i++) - // integrate(y0, 10, 10, y); - -/*double y0[] = {80,0,0,0,80,0}; -double y[12]; -integrate(y0, 10, 10, y); -for (int i=0; i<12; i++) std::cout << y[i] << std::endl;*/ - - //std::vector w = {8, 0.12, 0.04, -0.08, 200, -0.001}; - - - - std::cout << "Bye" << std::endl; - return 0; + std::array y0_array; + std::copy(y0, y0+6, begin(y0_array)); + integrate(y0_array, t_max, stride_size, y); } \ No newline at end of file diff --git a/scattering_experiment.cpp b/scattering_experiment.cpp index 84417ac..383efdd 100644 --- a/scattering_experiment.cpp +++ b/scattering_experiment.cpp @@ -17,8 +17,7 @@ using Generator = std::default_random_engine; // Our units are {kiloparsec, solar mass, gigayear} constexpr double G = 4.498317481097514e-06; -extern "C" -void integrate(const double y0[], const double t_max, const double stride_size, double y[]); +void integrate(const std::array y0, const double t_max, const double stride_size, double y[]); std::string git_version() { @@ -236,7 +235,7 @@ public: auto ic = initial_conditions(generator); double y[12]; - integrate(ic.data(), t_max, t_max, y); + integrate(ic, t_max, t_max, y); if (accept(y+6)) { std::lock_guard lock(mtx);