129 lines
No EOL
3.7 KiB
C++
129 lines
No EOL
3.7 KiB
C++
#include <iostream>
|
|
#include <string>
|
|
#include <vector>
|
|
#include <fstream>
|
|
#include <numeric>
|
|
#include <gsl/gsl_errno.h>
|
|
#include <gsl/gsl_math.h>
|
|
#include <gsl/gsl_odeiv2.h>
|
|
#include <gsl/gsl_spline.h>
|
|
|
|
|
|
extern "C" const int gsl_success() { return GSL_SUCCESS; } // It's zero, but just for clarity sake.
|
|
|
|
// Our units are {kiloparsec, solar mass, gigayear}
|
|
constexpr double G = 4.498317481097514e-06;
|
|
|
|
class Interp {
|
|
public:
|
|
Interp(std::vector<double>& x, std::vector<double>& y)
|
|
{
|
|
acc = gsl_interp_accel_alloc();
|
|
spline = gsl_spline_alloc(gsl_interp_cspline, x.size());
|
|
gsl_spline_init(spline, x.data(), y.data(), x.size());
|
|
}
|
|
inline double operator()(double x) const
|
|
{
|
|
return gsl_spline_eval(spline, x, acc);
|
|
}
|
|
|
|
private:
|
|
gsl_interp_accel *acc;
|
|
gsl_spline *spline;
|
|
};
|
|
|
|
class Plummer {
|
|
public:
|
|
Plummer(double M, double b)
|
|
: M(M), b(b) {}
|
|
void calc_acceleration(const double *pos, double *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;
|
|
}
|
|
|
|
private:
|
|
double M, b;
|
|
};
|
|
|
|
class Galaxy {
|
|
public:
|
|
Galaxy(std::string file_name)
|
|
{
|
|
std::vector<double> t_data, M_halo_data, b_halo_data;
|
|
|
|
std::ifstream file(file_name);
|
|
std::string line;
|
|
while (std::getline(file, line)) {
|
|
auto pos = line.find('#');
|
|
if (pos != std::string::npos) line = line.substr(0, pos);
|
|
pos = line.find_first_not_of(" \t");
|
|
if (pos == std::string::npos) continue;
|
|
double data[3];
|
|
sscanf(line.c_str(), "%*s %lf %lf %lf", &data[0], &data[1], &data[2]);
|
|
t_data.push_back(data[0]);
|
|
M_halo_data.push_back(data[1]);
|
|
b_halo_data.push_back(data[2]); // note, this is not half-mass radius
|
|
}
|
|
interp_M_halo = new Interp(t_data, M_halo_data);
|
|
interp_b_halo = new Interp(t_data, b_halo_data);
|
|
}
|
|
int func(double t, const double y[], double f[], void *params)
|
|
{
|
|
double M_halo = (*interp_M_halo)(t);
|
|
double b_halo = (*interp_b_halo)(t);
|
|
/*printf("xxxxxxxxx %e, %e msun\n", t, M_halo);
|
|
printf("xxxxxxxxx %e, %e kpc\n", t, b_halo);
|
|
exit(0);*/
|
|
Plummer plummer(M_halo, b_halo);
|
|
f[0] = y[3]; // vx -> x'
|
|
f[1] = y[4]; // vy -> y'
|
|
f[2] = y[5]; // vz -> z'
|
|
plummer.calc_acceleration(y, f+3); // a -> v'
|
|
return GSL_SUCCESS;
|
|
}
|
|
|
|
private:
|
|
Interp *interp_M_halo;
|
|
Interp *interp_b_halo;
|
|
} galaxy("file.dat");
|
|
// Not very nice to have it as a global variable but GSL will have problem otherwise.
|
|
|
|
int jac(double t, const double y[], double *dfdy, double dfdt[], void *params)
|
|
{
|
|
return GSL_SUCCESS;
|
|
}
|
|
|
|
|
|
int func(double t, const double y[], double f[], void *params)
|
|
{
|
|
return galaxy.func(t, y, f, params);
|
|
}
|
|
|
|
extern "C"
|
|
int integrate(const double y0[], const double t_max, double y[])
|
|
{
|
|
double t = 2.145;
|
|
constexpr double h = 1./4096.;
|
|
constexpr double epsabs = 1e-7;
|
|
constexpr double epsrel = 0;
|
|
const gsl_odeiv2_step_type *T = gsl_odeiv2_step_rk8pd;
|
|
gsl_odeiv2_step *s = gsl_odeiv2_step_alloc(T, 6);
|
|
gsl_odeiv2_evolve *e = gsl_odeiv2_evolve_alloc(6);
|
|
gsl_odeiv2_control *c = gsl_odeiv2_control_y_new(epsabs, 0);
|
|
gsl_odeiv2_system sys = {func, jac, 6, nullptr};
|
|
gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, T, h, epsabs, epsrel);
|
|
|
|
std::copy(y0, y0+6, y);
|
|
int status = gsl_odeiv2_driver_apply(d, &t, t_max, y);
|
|
return status;
|
|
}
|
|
|
|
int main()
|
|
{
|
|
std::cout << "bye" << std::endl;
|
|
} |