yohai/cosmo-replay
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Changed the integrate interface to use std::array instead of C-style

Browse source
This commit is contained in:
Yohai Meiron 2020-05-09 23:29:42 -04:00
parent 1eeab159cc
commit b86625ea13
2 changed files with 12 additions and 87 deletions
Download patch file Download diff file Expand all files Collapse all files

94
replay.cpp
View file

@ -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<double,6> 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<double, 6>;
@ -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<stride_count; i++) {
@ -144,81 +141,10 @@ void integrate(const double y0[], const double t_max, const double stride_size,
}
}
#include <random>
#include <chrono>
#include <unistd.h>
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<double> position_distribution(-x_max, x_max);
auto random_position_component = [&]() { return position_distribution(generator); };
std::uniform_real_distribution<double> 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<double> w = {8, 0.12, 0.04, -0.08, 200, -0.001};
std::cout << "Bye" << std::endl;
return 0;
std::array<double,6> y0_array;
std::copy(y0, y0+6, begin(y0_array));
integrate(y0_array, t_max, stride_size, y);
}

5
scattering_experiment.cpp
View file

@ -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<double,6> 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<std::mutex> lock(mtx);