/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Alejandro Benitez-Llambay
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see .
*
******************************************************************************/
#ifndef SWIFT_POTENTIAL_NFW_MN_H
#define SWIFT_POTENTIAL_NFW_MN_H
/* Config parameters. */
#include
/* Some standard headers. */
#include
#include
/* Local includes. */
#include "error.h"
#include "gravity.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "space.h"
#include "units.h"
/**
* @brief External Potential Properties - NFW Potential + Miyamoto-Nagai
*
* halo --> rho(r) = rho_0 / ( (r/R_s)*(1+r/R_s)^2 )
* disk --> phi(R,z) = -G * Mdisk / (R^2 + (Rdisk + (z^2+Zdisk^2)^1/2)^2)^(1/2)
*
* We however parameterise this in terms of c and virial_mass, Mdisk, Rdisk
* and Zdisk
*/
struct external_potential {
/*! Position of the centre of potential */
double x[3];
/*! The scale radius of the NFW potential */
double r_s;
/*! The pre-factor \f$ 4 \pi G \rho_0 \r_s^3 \f$ */
double pre_factor;
/*! The critical density of the universe */
double rho_c;
/*! The concentration parameter */
double c_200;
/*! The virial mass */
double M_200;
/*! Disk Size */
double Rdisk;
/*! Disk height */
double Zdisk;
/*! Disk Mass */
double Mdisk;
/*! Time-step condition pre_factor, this factor is used to multiply times the
* orbital time, so in the case of 0.01 we take 1% of the orbital time as
* the time integration steps */
double timestep_mult;
/*! Minimum time step based on the orbital time at the softening times
* the timestep_mult */
double mintime;
/*! Common log term \f$ \ln(1+c_{200}) - \frac{c_{200}}{1 + c_{200}} \f$ */
double log_c200_term;
/*! Softening length */
double eps;
};
/**
* @brief Computes the time-step due to the acceleration from the NFW + MN
* potential as a fraction (timestep_mult) of the circular orbital time of that
* particle.
*
* @param time The current time.
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static float external_gravity_timestep(
double time, const struct external_potential* restrict potential,
const struct phys_const* restrict phys_const,
const struct gpart* restrict g) {
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
const float R2 = dx * dx + dy * dy;
const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps);
const float mr = potential->M_200 *
(logf(1.f + r / potential->r_s) - r / (r + potential->r_s)) /
potential->log_c200_term;
const float Vcirc_NFW = sqrtf((phys_const->const_newton_G * mr) / r);
const float f1 =
potential->Rdisk + sqrtf(potential->Zdisk * potential->Zdisk + dz * dz);
const float Vcirc_MN = sqrtf(phys_const->const_newton_G * potential->Mdisk *
R2 / pow(R2 + f1 * f1, 3.0 / 2.0));
const float Vcirc = sqrtf(Vcirc_NFW * Vcirc_NFW + Vcirc_MN * Vcirc_MN);
const float period = 2 * M_PI * r / Vcirc;
/* Time-step as a fraction of the circular period */
const float time_step = potential->timestep_mult * period;
return max(time_step, potential->mintime);
}
/**
* @brief Computes the gravitational acceleration from an NFW Halo potential +
* MN disk.
*
* Note that the accelerations are multiplied by Newton's G constant
* later on.
*
* a_x = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * x
* a_y = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * y
* a_z = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * z
*
* @param time The current time.
* @param potential The #external_potential used in the run.
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
__attribute__((always_inline)) INLINE static void external_gravity_acceleration(
double time, const struct external_potential* restrict potential,
const struct phys_const* restrict phys_const, struct gpart* restrict g) {
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
/* First for the NFW part */
const float R2 = dx * dx + dy * dy;
const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps);
const float r_inv = 1.f / r;
const float term1 = potential->pre_factor;
const float term2 =
r / (r + potential->r_s) - logf(1.0f + r / potential->r_s);
const float acc_nfw = term1 * term2 * r_inv * r_inv * r_inv;
const float pot_nfw =
-potential->pre_factor * logf(1.0f + r / potential->r_s) * r_inv;
g->a_grav[0] += acc_nfw * dx;
g->a_grav[1] += acc_nfw * dy;
g->a_grav[2] += acc_nfw * dz;
gravity_add_comoving_potential(g, pot_nfw);
/* Now the the MN disk */
const float f1 = sqrtf(potential->Zdisk * potential->Zdisk + dz * dz);
const float f2 = potential->Rdisk + f1;
const float f3 = powf(R2 + f2 * f2, -1.5f);
const float mn_term = potential->Rdisk + sqrtf(potential->Zdisk + dz * dz);
const float pot_mn = -potential->Mdisk / sqrtf(R2 + mn_term * mn_term);
g->a_grav[0] -= potential->Mdisk * f3 * dx;
g->a_grav[1] -= potential->Mdisk * f3 * dy;
g->a_grav[2] -= potential->Mdisk * f3 * (f2 / f1) * dz;
gravity_add_comoving_potential(g, pot_mn);
}
/**
* @brief Computes the gravitational potential energy of a particle in an
* NFW potential + MN potential.
*
* phi = -4 * pi * G * rho_0 * r_s^3 * ln(1+r/r_s) - G * Mdisk / sqrt(R^2 +
* (Rdisk + sqrt(z^2 + Zdisk^2))^2)
*
* @param time The current time (unused here).
* @param potential The #external_potential used in the run.
* @param phys_const Physical constants in internal units.
* @param g Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
external_gravity_get_potential_energy(
double time, const struct external_potential* potential,
const struct phys_const* const phys_const, const struct gpart* g) {
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
/* First for the NFW profile */
const float R2 = dx * dx + dy * dy;
const float r = sqrtf(R2 + dz * dz + potential->eps * potential->eps);
const float term1 = -potential->pre_factor / r;
const float term2 = logf(1.0f + r / potential->r_s);
/* Now for the MN disk */
const float mn_term = potential->Rdisk + sqrtf(potential->Zdisk + dz * dz);
const float mn_pot = -potential->Mdisk / sqrtf(R2 + mn_term * mn_term);
return phys_const->const_newton_G * (term1 * term2 + mn_pot);
}
/**
* @brief Initialises the external potential properties in the internal system
* of units.
*
* @param parameter_file The parsed parameter file
* @param phys_const Physical constants in internal units
* @param us The current internal system of units
* @param potential The external potential properties to initialize
*/
static INLINE void potential_init_backend(
struct swift_params* parameter_file, const struct phys_const* phys_const,
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
/* Read in the position of the centre of potential */
parser_get_param_double_array(parameter_file, "NFW_MNPotential:position", 3,
potential->x);
/* Is the position absolute or relative to the centre of the box? */
const int useabspos =
parser_get_param_int(parameter_file, "NFW_MNPotential:useabspos");
if (!useabspos) {
potential->x[0] += s->dim[0] / 2.;
potential->x[1] += s->dim[1] / 2.;
potential->x[2] += s->dim[2] / 2.;
}
/* Read the other parameters of the model */
potential->timestep_mult =
parser_get_param_double(parameter_file, "NFW_MNPotential:timestep_mult");
potential->c_200 =
parser_get_param_double(parameter_file, "NFW_MNPotential:concentration");
potential->M_200 =
parser_get_param_double(parameter_file, "NFW_MNPotential:M_200");
potential->rho_c = parser_get_param_double(
parameter_file, "NFW_MNPotential:critical_density");
potential->Mdisk =
parser_get_param_double(parameter_file, "NFW_MNPotential:Mdisk");
potential->Rdisk =
parser_get_param_double(parameter_file, "NFW_MNPotential:Rdisk");
potential->Zdisk =
parser_get_param_double(parameter_file, "NFW_MNPotential:Zdisk");
potential->eps = 0.05;
/* Compute R_200 */
const double R_200 =
cbrtf(3.0 * potential->M_200 / (4. * M_PI * 200.0 * potential->rho_c));
/* NFW scale-radius */
potential->r_s = R_200 / potential->c_200;
const double r_s3 = potential->r_s * potential->r_s * potential->r_s;
/* Log(c_200) term appearing in many expressions */
potential->log_c200_term =
log(1. + potential->c_200) - potential->c_200 / (1. + potential->c_200);
const double rho_0 =
potential->M_200 / (4.f * M_PI * r_s3 * potential->log_c200_term);
/* Pre-factor for the accelerations (note G is multiplied in later on) */
potential->pre_factor = 4.0f * M_PI * rho_0 * r_s3;
/* Compute the orbital time at the softening radius */
const double sqrtgm = sqrt(phys_const->const_newton_G * potential->M_200);
const double epslnthing = log(1.f + potential->eps / potential->r_s) -
potential->eps / (potential->eps + potential->r_s);
potential->mintime = 2. * M_PI * potential->eps * sqrtf(potential->eps) *
sqrtf(potential->log_c200_term / epslnthing) / sqrtgm *
potential->timestep_mult;
}
/**
* @brief Prints the properties of the external potential to stdout.
*
* @param potential The external potential properties.
*/
static INLINE void potential_print_backend(
const struct external_potential* potential) {
message(
"External potential is 'NFW + MN disk' with properties are (x,y,z) = "
"(%e, %e, %e), scale radius = %e timestep multiplier = %e, mintime = %e",
potential->x[0], potential->x[1], potential->x[2], potential->r_s,
potential->timestep_mult, potential->mintime);
}
#endif /* SWIFT_POTENTIAL_NFW_MN_H */