#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);
}