Processed subhalo data and modified integrator to read and use new data structure

loadtxt.cpp

@@ -62,7 +62,7 @@ void Loadtxt::line_to_buf(std::vector<int> cols, std::string line, double *buffe
         if (col++ == cols[idx]) buffer[idx++] = std::stod(num_as_string);
     }
     if (col++ == cols[idx]) buffer[idx++] = std::stod(line);
-    if (idx < n_cols) throw;
+    if (idx < n_cols) throw std::runtime_error("Read fewer columns than expected");
 }
 
 // Below is a deomonstration. The file has multiple columns but we are only

main.cpp

@@ -1,6 +1,7 @@
 #include <algorithm>
 #include <array>
 #include <boost/numeric/odeint.hpp>
+#include <Eigen/Geometry> 
 #include <fstream>
 #include <gsl/gsl_spline.h>
 #include <iostream>
@@ -11,7 +12,7 @@
 
 #include "loadtxt.h"
 
-
+using double3 = Eigen::Vector3d;
 
 // Our units are {kiloparsec, solar mass, gigayear}
 constexpr double G = 4.498317481097514e-06;
@@ -41,52 +42,78 @@ private:
     gsl_spline *spline;
 };
 
-class Plummer {
-public:
-    Plummer(double M, double b)
-        : M(M), b(b) {}
-    void calc_acceleration(const double *pos, double *acc)
+void plummer(double M, double b, const double3& pos, double3& acc)
 {
     double r2 = (pos[0]*pos[0] + pos[1]*pos[1] + pos[2]*pos[2] + b*b);
     double r  = sqrt(r2);
     double r3_inv = 1/(r*r2);
-        acc[0] = -G*M*pos[0]*r3_inv;
-        acc[1] = -G*M*pos[1]*r3_inv;
-        acc[2] = -G*M*pos[2]*r3_inv;
+    acc = -G*M*pos*r3_inv;
 }
 
-private:
-    double M, b;
-};
+void miyamoto_nagai(double M, double a, double b, double phi, double theta, const double3& pos, double3& acc)
+{
+    // Construct new z-axis from the angles
+    double cos_theta = cos(theta);
+    double sin_theta = sqrt(1-cos_theta*cos_theta);
+    double cos_phi   = cos(phi);
+    double sin_phi   = sqrt(1-cos_phi*cos_phi);
+    Eigen::Vector3d old_z_axis = {cos_phi*sin_theta, sin_phi*sin_theta, cos_theta};
+
+    // Construct rotation
+    auto rot = Eigen::Quaternion<double>::FromTwoVectors(old_z_axis, Eigen::Vector3d::UnitZ());
+
+    // Rotate the position vector
+    auto new_pos = rot * pos;
+
+    // Calculate acceleration in new frame
+    Eigen::Vector3d acc_in_rotated_frame;
+    auto z_tmp  = sqrt(new_pos[2]*new_pos[2] + b*b);
+    auto r2_tmp = new_pos[0]*new_pos[0] + new_pos[1]*new_pos[1] + (z_tmp + a)*(z_tmp + a);
+    auto r_tmp  = sqrt(r2_tmp);
+    auto tmp = G*M / (r2_tmp*r_tmp);
+    acc_in_rotated_frame[0] = - tmp * new_pos[0];
+    acc_in_rotated_frame[1] = - tmp * new_pos[1];
+    acc_in_rotated_frame[2] = - tmp * new_pos[2] * (z_tmp + a)/z_tmp;
+
+    // Return to original frame
+    acc = rot.inverse() * acc_in_rotated_frame;
+}
 
 class Galaxy {
 public:
     Galaxy(std::string file_name)
     {
-        auto data = Loadtxt("file.dat", {1, 2, 3}).get_cols();
+        auto data = Loadtxt(file_name, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11}).get_cols();
         auto& t_data = data[0];
-        auto& halo_m_data = data[1];
-        auto& halo_b_data = data[2];
         std::transform(begin(t_data), end(t_data), begin(t_data), [t0=t_data[0]](const double& x){ return x-t0; });
-        std::transform(halo_b_data.begin(), halo_b_data.end(), halo_b_data.begin(), [](const double& x){ return x*0.7664209365408798; });
-        interp_halo_m = Interp(t_data, halo_m_data);
-        interp_halo_b = Interp(t_data, halo_b_data);
+        for (int i=0; i<10; i++) interp[i] = Interp(t_data, data[i+1]);
     }
     void func(const std::array<double, 6> &y, std::array<double, 6> &f, const double t)
     {
-        double halo_m = interp_halo_m(t);
-        double halo_b = interp_halo_b(t);
-        Plummer plummer(halo_m, halo_b);
         f[0] = y[3]; // vx -> x'
         f[1] = y[4]; // vy -> y'
         f[2] = y[5]; // vz -> z'
-        plummer.calc_acceleration(y.data(), f.data()+3); // a -> v'
+        double m_disk      = interp[0](t);
+        double phi         = interp[1](t);
+        double theta       = interp[2](t);
+        double a_disk      = interp[3](t);
+        double b_disk      = interp[4](t);
+        double m_halo      = interp[5](t);
+        double b_halo      = interp[6](t);
+        double3 dalembert;
+        // for (int i=0; i<3; i++) dalembert[i] = interp[7+i](t);
+        dalembert = {0,0,0};
+
+        double3 acc_disk, acc_halo;
+        const double3 pos(y.data());
+        miyamoto_nagai(m_disk, a_disk, b_disk, phi, theta, pos, acc_disk);
+        plummer(m_halo, b_halo, pos, acc_halo);
+        for (int i=0; i<3; i++) f[3+i] = acc_disk[i] + acc_halo[i] + dalembert[i];
     }
 
 private:
-    Interp interp_halo_m;
-    Interp interp_halo_b;
-} galaxy("file.dat");
+    Interp interp[10];
+} galaxy("subhalo_411321_parameters_processed.dat");
 
 extern "C"
 int integrate(const double y0[], const double t_max, const double stride_size, double y[])
@@ -113,14 +140,24 @@ int integrate(const double y0[], const double t_max, const double stride_size, d
 }
 
 
+
+
+
 int main()
 {
     std::cout << "Hi" << 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);
+    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];

plot.py

@@ -4,7 +4,7 @@
 import numpy as np
 import ctypes, subprocess
 
-recompile   = True
+recompile = False
 
 if recompile:
     p = subprocess.Popen('g++ -I/home/meiron/boost_1_72_0 -shared -o libmain.so -fPIC loadtxt.cpp main.cpp -lgsl'.split())
@@ -21,28 +21,43 @@ def integrate(y0, t_max, step_size=None):
     y = np.array(y).reshape(size,6)
     return np.array(y), status
 
-def integrate_boost(y0, t_max, step_size=None):
-    y0 = (ctypes.c_double*6)(*y0)
-    if step_size is None: step_size = t_max
-    size = int(t_max // step_size) + 1
-    y = (ctypes.c_double*size*6)()
-    libmain.integrate_boost(y0, ctypes.c_double(t_max), ctypes.c_double(step_size), y)
-    y = np.array(y).reshape(size,6)
-    return np.array(y), 0
-
-#gsl_success = libmain.gsl_success()
-
 from pylab import *
 t_max = 10
+x_max = 85
 ic = [80,0,0,0,80,0]
 res = integrate(ic, t_max, step_size=32/4096)
 x, y, z, vx, vy, vz = res[0].T
-zzz = x[-1]
-plot(x,y)
 
-res = integrate(ic, t_max)
-x, y, z, vx, vy, vz = res[0].T
-plot(x,y,'o')
-print(zzz - x[-1])
+
+figure(figsize=(8,8))
+subplot(221)
+plot(x, y)
+plot(x[0], y[0], 'x')
 gca().set_aspect('equal')
+xlabel('x')
+ylabel('y')
+xlim(-x_max,x_max)
+ylim(-x_max,x_max)
+
+subplot(223)
+plot(x, z)
+plot(x[0], z[0], 'x')
+gca().set_aspect('equal')
+xlabel('x')
+ylabel('z')
+xlim(-x_max,x_max)
+ylim(-x_max,x_max)
+
+subplot(224)
+plot(y, z)
+plot(y[0], z[0], 'x')
+gca().set_aspect('equal')
+xlabel('z')
+ylabel('y')
+xlim(-x_max,x_max)
+ylim(-x_max,x_max)
+
+subplots_adjust(wspace=0.25)
+
+savefig('test_orbit.png')
 show()