Improved initial conditions

scattering_experiment.cpp

@@ -13,6 +13,9 @@
 
 using double3 = Eigen::Vector3d;
 
+// 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[]);
 
@@ -81,6 +84,53 @@ private:
     double target_d, target_z, target_vd, target_vt, target_vz;
 };
 
+class Uniform_sphere_escape_velocity_cutoff {
+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<double>(0, 1);
+    }
+    template<typename Generator>
+    std::array<double,6> operator()(Generator& generator)
+    {
+        std::array<double,6> result;
+
+        double r         = pow(uniform_distribution(generator), 1./3.)*r_max;
+        double cos_theta = uniform_distribution(generator)*2-1;
+        double sin_theta = sqrt(1 - cos_theta*cos_theta);
+        double phi       = uniform_distribution(generator)*2*M_PI;
+        result[0] = r*cos(phi)*sin_theta;
+        result[1] = r*sin(phi)*sin_theta;
+        result[2] = r*cos_theta;
+
+        double potential = nfw_potential(r);
+        if (potential > 0) throw std::runtime_error("Cannot deal with a positive potential");
+        double v_max = sqrt(-2*potential);
+        double v_magnitude;
+        do {
+            for (int i=3; i<6; i++) result[i] = uniform_distribution(generator);
+            v_magnitude = sqrt(result[3]*result[3]+result[4]*result[4]+result[5]*result[5]);
+        } while (v_magnitude >= v_max);
+
+        return result;
+    }
+    friend std::ostream& operator<<(std::ostream& stream, const Uniform_sphere_escape_velocity_cutoff& obj)
+    {
+        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';
+    }
+private:
+    inline double nfw_potential(double r)
+    {
+        return -4*M_PI*G*rho_0*b*b*b*log(1+r/b)/r;
+    }
+    double r_max, rho_0, b;
+    std::uniform_real_distribution<double> uniform_distribution;
+};
+
 class Box_initial_conditions {
 public:
     Box_initial_conditions(double x_max, double v_max)
@@ -162,8 +212,6 @@ public:
 
                 if (accept(y+6)) {
                     std::lock_guard<std::mutex> lock(mtx);
-                    char buffer[512];
-                    file << std::setw(3) << tid;
                     for (int j=0; j<12; j++) file << std::setw(24) << y[j];
                     file << std::endl;
                 }
@@ -215,10 +263,10 @@ int main()
     target_pos *= 0.001; // Now in kpc
     target_vel *= 1.0226911647958985; // Now in kpc/Gyr
 
-    In_box_cyl accept(2, 20, target_pos, target_vel);
+    In_box_cyl accept(5, 50, target_pos, target_vel);
-    Box_initial_conditions box_initial_conditions(75, 250);
+    Uniform_sphere_escape_velocity_cutoff uniform_sphere_escape_velocity_cutoff(75, 1.82777387E+07, 9.86242602E+00);
 
-    Scattering_experiment scattering_experiment("results.dat", t_max, 1000000000L, box_initial_conditions, accept);
+    Scattering_experiment scattering_experiment("results.dat", t_max, 1000000000L, uniform_sphere_escape_velocity_cutoff, accept);
     scattering_experiment.launch();
 
     std::cout << "Bye\n";