85 lines
2.5 KiB
C++
85 lines
2.5 KiB
C++
#include "external.h"
|
|
|
|
inline double square(double x) {return x*x;}
|
|
|
|
void Plummer::calc_gravity()
|
|
{
|
|
double r2 = square(b);
|
|
r2 += x.norm2();
|
|
double r = sqrt(r2);
|
|
double rv_ij = v*x;
|
|
double tmp = m / r;
|
|
potential = -tmp;
|
|
tmp /= r2;
|
|
acceleration = - tmp * x;
|
|
jerk = - tmp * (v - 3*rv_ij * x/r2);
|
|
}
|
|
|
|
void Miyamoto_Nagai::calc_gravity()
|
|
{
|
|
double x_ij=x[0], y_ij=x[1], z_ij=x[2];
|
|
double vx_ij=v[0], vy_ij=v[1], vz_ij=v[2];
|
|
auto z2_tmp = z_ij*z_ij + square(b);
|
|
auto z_tmp = sqrt(z2_tmp);
|
|
auto r2_tmp = x_ij*x_ij + y_ij*y_ij + square(z_tmp + a);
|
|
auto r_tmp = sqrt(r2_tmp);
|
|
potential = - m / r_tmp;
|
|
auto tmp = m / (r2_tmp*r_tmp);
|
|
acceleration[0] = - tmp * x_ij;
|
|
acceleration[1] = - tmp * y_ij;
|
|
acceleration[2] = - tmp * z_ij * (z_tmp + a)/z_tmp;
|
|
tmp = m / (z_tmp*r2_tmp*r2_tmp*r_tmp);
|
|
jerk[0] = tmp * (- vx_ij*z_tmp*r2_tmp
|
|
+ 3.0*x_ij*vx_ij*x_ij*z_tmp
|
|
+ 3.0*x_ij*vy_ij*y_ij*z_tmp
|
|
+ 3.0*x_ij*vz_ij*z_ij*square(z_tmp + a));
|
|
jerk[1] = tmp * (- vy_ij*z_tmp*r2_tmp
|
|
+ 3.0*y_ij*vx_ij*x_ij*z_tmp
|
|
+ 3.0*y_ij*vy_ij*y_ij*z_tmp
|
|
+ 3.0*y_ij*vz_ij*z_ij*square(z_tmp + a));
|
|
jerk[2] = tmp * (- vz_ij*(z_tmp + a)*(x_ij*x_ij*(z2_tmp*z_tmp + a*square(b)) +
|
|
y_ij*y_ij*(z2_tmp*z_tmp + a*square(b)) -
|
|
(2.0*a*(square(z_ij) - square(b))*z_tmp +
|
|
2.0*square(z_ij*z_ij) +
|
|
square(b)*square(z_ij) -
|
|
square(b)*(square(a) + square(b))))
|
|
+ 3.0*vx_ij*x_ij*z_ij*z2_tmp*(z_tmp + a)
|
|
+ 3.0*vy_ij*y_ij*z_ij*z2_tmp*(z_tmp + a)) / z2_tmp;
|
|
}
|
|
|
|
void Logarithmic_halo::calc_gravity()
|
|
{
|
|
auto r2 = x.norm2();
|
|
auto rv_ij = 2.0*(v*x);
|
|
auto r2_r2_halo = (r2 + r2_halo);
|
|
potential = - 0.5*v2_halo * log(1.0 + r2/r2_halo);
|
|
auto tmp = v2_halo/r2_r2_halo;
|
|
acceleration = - tmp * x;
|
|
tmp /= (r2 + r2_halo);
|
|
jerk = tmp * (rv_ij * x - r2_r2_halo * v);
|
|
}
|
|
|
|
void Dehnen::calc_gravity()
|
|
{
|
|
|
|
auto tmp_r = x.norm();
|
|
auto tmp_r2 = tmp_r*tmp_r;
|
|
auto tmp_r3 = tmp_r2*tmp_r;
|
|
|
|
auto dum = tmp_r/(tmp_r + r);
|
|
auto dum2 = dum*dum;
|
|
auto dum3 = dum2*dum;
|
|
auto dum_g = pow(dum, -gamma);
|
|
|
|
potential = - ( (m/r) / (2.0-gamma) ) * ( 1.0 - dum2 * dum_g );
|
|
|
|
auto tmp = (m/tmp_r3) * dum3 * dum_g;
|
|
|
|
acceleration = - tmp * x;
|
|
|
|
auto rv_ij = v*x;
|
|
tmp = ( m/((tmp_r+r)*(tmp_r+r)*(tmp_r+r)) ) * dum_g;
|
|
auto tmp_g = ( (r*gamma + 3.0*tmp_r)/(tmp_r2*(tmp_r+r)) ) * rv_ij;
|
|
|
|
jerk = tmp * (tmp_g * x - v);
|
|
}
|