From 1eeab159cc52405923360cd4660b04792d6b9b29 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Sat, 9 May 2020 22:48:29 -0400 Subject: [PATCH] Changed the main class from a template to one that uses polymorphism --- scattering_experiment.cpp | 114 ++++++++++++++++++++++++-------------- 1 file changed, 71 insertions(+), 43 deletions(-) diff --git a/scattering_experiment.cpp b/scattering_experiment.cpp index de6cb77..84417ac 100644 --- a/scattering_experiment.cpp +++ b/scattering_experiment.cpp @@ -12,6 +12,7 @@ #include using double3 = Eigen::Vector3d; +using Generator = std::default_random_engine; // Our units are {kiloparsec, solar mass, gigayear} constexpr double G = 4.498317481097514e-06; @@ -19,7 +20,8 @@ constexpr double G = 4.498317481097514e-06; extern "C" void integrate(const double y0[], const double t_max, const double stride_size, double y[]); -std::string git_version() { +std::string git_version() +{ std::string result; std::ifstream git_log(".git/logs/HEAD"); if (git_log) { @@ -43,7 +45,14 @@ std::string iso_time(std::chrono::system_clock::time_point time) return std::string(time_string); } -class In_box_cyl { +class Acceptor { +// This is a base class for functors that accept or decline the final state of the orbit +public: + virtual bool operator()(const double w[6]) const = 0; + virtual void output_information(std::ostream& stream) const = 0; +}; + +class In_box_cyl : public Acceptor { public: In_box_cyl(const double pos_tol, const double vel_tol, const double3 target_pos, const double3 target_vel) : pos_tol(pos_tol), vel_tol(vel_tol) @@ -54,45 +63,62 @@ public: target_vt = (target_pos[0]*target_vel[1] - target_pos[1]*target_vel[0])/target_d; target_vz = target_vel[2]; } - bool operator()(const double w[6]) { - const double d = sqrt(w[0]*w[0] + w[1]*w[1]); - const double& z = w[2]; - const double vd = (w[0]*w[3] + w[1]*w[4])/d; - const double vt = (w[0]*w[4] - w[1]*w[3])/d; - const double& vz = w[5]; + bool operator()(const double w[6]) const override + { + double d = sqrt(w[0]*w[0] + w[1]*w[1]); + double z = w[2]; + double vd = (w[0]*w[3] + w[1]*w[4])/d; + double vt = (w[0]*w[4] - w[1]*w[3])/d; + double vz = w[5]; + + if (((z<0)&&(target_z>0)) || ((z>0)&&(target_z<0))) { // Flip up-down + z = -z; + vz = -vz; + } + + if (((vt<0)&&(target_vt>0)) || ((vt>0)&&(target_vt<0))) { // Flip left-right + vt = -vt; + } + return (abs(d - target_d) < pos_tol) && (abs(z - target_z) < pos_tol) && (abs(vd - target_vd) < vel_tol) && (abs(vt - target_vt) < vel_tol) && (abs(vz - target_vz) < vel_tol); } - friend std::ostream& operator<<(std::ostream& stream, const In_box_cyl& obj) + void output_information(std::ostream& stream) const override { stream.setf(std::ios::scientific); stream << "# Acceptor type: box in cylindrical coordinates\n"; stream << "# Acceptor parameters:\n"; - stream << "# pos_tol = " << std::setw(13) << obj.pos_tol << '\n'; - stream << "# vel_tol = " << std::setw(13) << obj.vel_tol << '\n'; - stream << "# target_d = " << std::setw(13) << obj.target_d << '\n'; - stream << "# target_z = " << std::setw(13) << obj.target_z << '\n'; - stream << "# target_vd = " << std::setw(13) << obj.target_vd << '\n'; - stream << "# target_vt = " << std::setw(13) << obj.target_vt << '\n'; - stream << "# target_vz = " << std::setw(13) << obj.target_vz << '\n'; + stream << "# pos_tol = " << std::setw(13) << pos_tol << '\n'; + stream << "# vel_tol = " << std::setw(13) << vel_tol << '\n'; + stream << "# target_d = " << std::setw(13) << target_d << '\n'; + stream << "# target_z = " << std::setw(13) << target_z << '\n'; + stream << "# target_vd = " << std::setw(13) << target_vd << '\n'; + stream << "# target_vt = " << std::setw(13) << target_vt << '\n'; + stream << "# target_vz = " << std::setw(13) << target_vz << '\n'; } private: double pos_tol, vel_tol; double target_d, target_z, target_vd, target_vt, target_vz; }; -class Uniform_sphere_escape_velocity_cutoff { +class Initial_conditions { +public: + virtual std::array operator()(Generator& generator) const = 0; + virtual void output_information(std::ostream& stream) const = 0; +}; + +class Uniform_sphere_escape_velocity_cutoff : public Initial_conditions { public: Uniform_sphere_escape_velocity_cutoff(const double r_max, const double rho_0, const double b) - : r_max(r_max), rho_0(rho_0), b(b) - { - uniform_distribution = std::uniform_real_distribution(0, 1); - } - template - std::array operator()(Generator& generator) + : r_max(r_max), + rho_0(rho_0), + b(b), + uniform_distribution(std::uniform_real_distribution(0, 1)) + {} + std::array operator()(Generator& generator) const override { std::array result; @@ -115,23 +141,23 @@ public: return result; } - friend std::ostream& operator<<(std::ostream& stream, const Uniform_sphere_escape_velocity_cutoff& obj) + void output_information(std::ostream& stream) const override { stream.setf(std::ios::scientific); stream << "# Initial conditions type: uniform sphere with escape velocity cutoff\n"; stream << "# Initial conditions parameters:\n"; - stream << "# r_max = " << std::setw(13) << obj.r_max << '\n'; + stream << "# r_max = " << std::setw(13) << r_max << '\n'; } private: - inline double nfw_potential(double r) + inline double nfw_potential(double r) const { - return -4*M_PI*G*rho_0*b*b*b*log(1+r/b)/r; + return -4*M_PI*G*rho_0*b*b*b*log1p(r/b)/r; } double r_max, rho_0, b; - std::uniform_real_distribution uniform_distribution; + mutable std::uniform_real_distribution uniform_distribution; }; -class Box_initial_conditions { +class Box_initial_conditions : public Initial_conditions { public: Box_initial_conditions(double x_max, double v_max) : x_max(x_max), v_max(v_max) @@ -139,33 +165,35 @@ public: position_distribution = std::uniform_real_distribution(-x_max, x_max); velocity_distribution = std::uniform_real_distribution(-v_max, v_max); } - template - std::array operator()(Generator& generator) + std::array operator()(Generator& generator) const override { std::array result; for (int i=0; i<3; i++) result[i] = position_distribution(generator); for (int i=3; i<6; i++) result[i] = velocity_distribution(generator); return result; } - friend std::ostream& operator<<(std::ostream& stream, const Box_initial_conditions& obj) + void output_information(std::ostream& stream) const override { stream.setf(std::ios::scientific); stream << "# Initial conditions type: uniform box\n"; stream << "# Initial conditions parameters:\n"; - stream << "# x_max = " << std::setw(13) << obj.x_max << '\n'; - stream << "# v_max = " << std::setw(13) << obj.v_max << '\n'; + stream << "# x_max = " << std::setw(13) << x_max << '\n'; + stream << "# v_max = " << std::setw(13) << v_max << '\n'; } private: double x_max, v_max; - std::uniform_real_distribution<> position_distribution, velocity_distribution; + mutable std::uniform_real_distribution<> position_distribution, velocity_distribution; }; - -template class Scattering_experiment { public: - Scattering_experiment(const std::string file_name, const double t_max, const unsigned long long n_experiments, Initial_conditions& initial_conditions, Acceptor& accept) - : bunch_size(256), n_experiments(n_experiments), report_interval_seconds(300), t_max(t_max), accept(accept), initial_conditions(initial_conditions) + Scattering_experiment(const std::string& file_name, const double t_max, const unsigned long long n_experiments, const Initial_conditions& initial_conditions, const Acceptor& accept) + : n_experiments(n_experiments), + bunch_size(256), + report_interval_seconds(300), + t_max(t_max), + accept(accept), + initial_conditions(initial_conditions) { auto time = std::chrono::system_clock::now(); file = std::ofstream(file_name/*, std::ofstream::app*/); @@ -185,8 +213,8 @@ public: file << "# Unique identifier: " << unique_id << '\n'; n_threads = std::thread::hardware_concurrency(); file << "# Concurrency: " << n_threads << '\n'; - file << initial_conditions; - file << accept; + initial_conditions.output_information(file); + accept.output_information(file); file.flush(); file.setf(std::ios::scientific | std::ios::right); file.precision(16); @@ -248,8 +276,8 @@ public: int bunch_size; int report_interval_seconds; double t_max; - Acceptor accept; - Initial_conditions initial_conditions; + const Acceptor& accept; + const Initial_conditions& initial_conditions; }; int main()