Added Dehnen potential and moved external potentials to a new compilation unit

Makefile

@@ -25,7 +25,7 @@ CPPFLAGS     += -DHAS_HDF5
 LIB          += -lhdf5 -lz -ldl
 
 default:
-	$(MPICXX) $(CPPFLAGS) $(CXXFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) io.cpp config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB)
+	$(MPICXX) $(CPPFLAGS) $(CXXFLAGS) -DETICS_DTSCF=$(ETICS_DTSCF) $(INC) external.cpp io.cpp config.cpp phigrape.cpp -o $(EXECUTABLE) $(LIB)
 
 yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS))
 yebisu: default

config.cpp

@@ -3,9 +3,12 @@
 #include <cstdio>
 #include <fstream>
 #include <algorithm>
+#include <limits>
 
 using Dictionary = std::unordered_map<std::string,std::string>;
 
+static constexpr double nix = -std::numeric_limits<double>::max(); // avoid nans
+
 std::string Config::strip(const std::string str)
 {
     std::string str_new = str;
@@ -138,6 +141,8 @@ void Config::error_checking()
         throw std::runtime_error("To use external disk gravity, please specify positive ext_m_disk, ext_a_disk and ext_b_disk");
     if (((ext_log_halo_r > 0) && (ext_log_halo_v <= 0)) || ((ext_log_halo_r <= 0) && (ext_log_halo_v > 0)))
         throw std::runtime_error("To use external logarithmic halo gravity, please specify positive ext_log_halo_r and ext_log_halo_v");
+    if ((ext_dehnen_m > 0) && ((ext_dehnen_r <= 0) || (ext_dehnen_gamma <= 0)))
+        throw std::runtime_error("To use external Dehnen model, please specify positive ext_dehnen_r and ext_dehnen_gamma");
 }
 
 Config::Config(std::string file_name)
@@ -186,7 +191,9 @@ Config::Config(std::string file_name)
     ext_b_halo_plummer = get_parameter<double>(dictionary, "ext_b_halo_plummer", -1);
     ext_log_halo_v = get_parameter<double>(dictionary, "ext_log_halo_v", 0);
     ext_log_halo_r = get_parameter<double>(dictionary, "ext_log_halo_r", 0);
-
+    ext_dehnen_m = get_parameter<double>(dictionary, "ext_dehnen_m", 0);
+    ext_dehnen_r = get_parameter<double>(dictionary, "ext_dehnen_r", -1);
+    ext_dehnen_gamma = get_parameter<double>(dictionary, "ext_dehnen_gamma", -1);
 
     error_checking();
 }

config.h

@@ -41,6 +41,9 @@ public:
     double           ext_b_halo_plummer;
     double           ext_log_halo_r;
     double           ext_log_halo_v;
+    double           ext_dehnen_m;
+    double           ext_dehnen_r;
+    double           ext_dehnen_gamma;
 
 private:
     using Dictionary = std::unordered_map<std::string,std::string>;

external.cpp

@@ -0,0 +1,86 @@
+#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);
+}

external.h

@@ -0,0 +1,59 @@
+#pragma once
+#include "double3.h"
+
+class External_gravity {
+public:
+    void apply(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[])
+    {
+        for (int i=0; i<n_act; i++) {
+            this->set_coordinates(x[i], v[i]);
+            this->calc_gravity();
+            pot[i]  += potential;
+            a[i]    += acceleration;
+            adot[i] += jerk;
+        }
+    }
+    virtual void calc_gravity() = 0;
+    bool is_active;
+protected:
+    double potential;
+    double3 acceleration, jerk;
+    double3 x, v;
+    void set_coordinates(double3 x, double3 v)
+    {
+        this->x = x;
+        this->v = v;
+    }
+};
+
+class Plummer : public External_gravity {
+public:
+    Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);}
+    void calc_gravity();
+private:
+    double m, b;
+};
+
+class Miyamoto_Nagai : public External_gravity {
+public:
+    Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0);}
+    void calc_gravity();
+private:
+    double m, a, b;
+};
+
+class Logarithmic_halo : public External_gravity {
+public:
+    Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);}
+    void calc_gravity();
+private:
+    double v2_halo, r2_halo;
+};
+
+class Dehnen : public External_gravity {
+public:
+    Dehnen(double m, double r, double gamma) : m(m), r(r), gamma(gamma) {is_active=(m>0);}
+    void calc_gravity();
+private:
+    double m, r, gamma;
+};

phigrape.conf

@@ -49,6 +49,7 @@ input_file_name =                       data.con
 # Whether to output a warning on the screen when the minimum time step is encountered. [default: false]
 #dt_min_warning =                        false
 
+
 ####################
 # EXTERNAL GRAVITY #
 ####################
@@ -84,6 +85,10 @@ input_file_name =                       data.con
 #ext_log_halo_v =                        240 # km/s
 #ext_log_halo_r =                        1   # kpc
 
+# This is a spherical Dehnen model.
+#ext_dehnen_m =                          1E11 # MSun
+#ext_dehnen_r =                          2    # kpc
+#ext_dehnen_gamma =                      0.5
 
 
 ###################################
@@ -107,6 +112,7 @@ input_file_name =                       data.con
 # Number of nearest neighbours to the SMBH (or SMBHs) to include in output. [default: 10]
 #live_smbh_neighbor_number =             10
 
+
 ##################################
 # BINARY SUPERMASSIVE BLACK HOLE #
 ##################################
@@ -129,6 +135,7 @@ input_file_name =                       data.con
 # The speed of light in N-body units [default: 500]
 #pn_c =                                  477.12
 
+
 ####################################
 # Negative powers of two           #
 ####################################

phigrape.cpp

@@ -56,16 +56,10 @@ Last redaction : 2019.04.16 12:55
 #include "double3.h"
 #include "config.h"
 #include "io.h"
+#include "external.h"
 
 Config *config;
 
-#define NORM            // Physical normalization
-
-//#define EXTPOT_BH        // BH - usually NB units
-
-
-//#define EXTPOT_GAL_DEH      // Dehnen Galactic - usually physical units
-
 #ifdef ETICS
 #include "grapite.h"
 // why do we need CEP as a compilaion flag... just have it always on when ETICS is on. IF there is no CEP, there should be a graceful skipping of those operations.
@@ -131,35 +125,6 @@ Config *config;
 #define TWOPi       6.283185307179
 #define sqrt_TWOPi  2.506628274631
 
-#ifdef NORM
-//http://pdg.lbl.gov/2015/reviews/rpp2015-rev-astrophysical-constants.pdf
-
-#define G      6.67388E-11        // (m/s^2) * (m^2/kg)
-#define Msol   1.988489E+30        // kg
-#define Rsol   6.957E+08        // m
-#define AU     149597870700.0        // m
-#define pc     3.08567758149E+16    // m
-#define Year   31556925.2        // s
-#define c_feny 299792458.0        // m/s
-
-#define kpc    (1.0E+03*pc)        // m
-#define km     1.0E+03            // km   -> m
-#define cm3    1.0E-06            // cm^3 -> m^3
-#define Myr    (1.0E+06*Year)        // s
-#define Gyr    (1.0E+09*Year)        // s
-#define R_gas  8.31447215        // J/(K*mol)
-#define k_gas  1.380650424E-23        // J/K
-#define N_A    6.022141510E+23      // 1/mol
-#define mu     1.6605388628E-27        // kg
-#define mp     1.67262163783E-27     // kg
-#define me     9.1093821545E-31     // kg
-
-#define pc2    (pc*pc)
-#define pc3    (pc*pc*pc)
-#define kpc2   (kpc*kpc)
-#define kpc3   (kpc*kpc*kpc)
-#endif
-
 /*
     1KB =    1024
     2KB =    2048
@@ -230,24 +195,8 @@ double eps_BH=0.0;
 
 /* external potential... */
 
-
-#ifdef EXTPOT_GAL_DEH
-double m_ext, r_ext, g_ext,
-    tmp_r2, tmp_r3, dum, dum2, dum3, dum_g, tmp_g;
-#endif
-
-#ifdef EXTPOT_BH
-double  r2, rv_ij,
-        m_bh, b_bh, eps_bh; // NOTE different from eps_BH
-        // NOTE there are other m_bh and b_bh defined outside this ifdef block
-        // NOTE this is a mess. Looks like eps_bh is never used, and the m_bh from outside the block is never used.
-#endif
-
-// double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3];
 double3 x_bh1, x_bh2, v_bh1, v_bh2;
 
-// double pot_bh1, a_bh1[3], adot_bh1[3],
-//     pot_bh2, a_bh2[3], adot_bh2[3];
 double pot_bh1, pot_bh2;
 double3 a_bh1, adot_bh1, a_bh2, adot_bh2;
 
@@ -401,106 +350,6 @@ void write_bh_nb_data(double time_cur, int N, double m[], double3 x[], double3 v
     fclose(out);
 }
 
-class External_gravity {
-public:
-    void apply(const int n_act, const double3 x[], const double3 v[], double pot[], double3 a[], double3 adot[])
-    {
-        for (int i=0; i<n_act; i++) {
-            this->set_coordinates(x[i], v[i]);
-            this->calc_gravity();
-            pot[i]  += potential;
-            a[i]    += acceleration;
-            adot[i] += jerk;
-        }
-    }
-    virtual void calc_gravity() = 0;
-    bool is_active;
-protected:
-    double potential;
-    double3 acceleration, jerk;
-    double3 x, v;
-    void set_coordinates(double3 x, double3 v)
-    {
-        this->x = x;
-        this->v = v;
-    }
-};
-
-class Plummer : public External_gravity {
-public:
-    Plummer(double m, double b) : m(m), b(b) {is_active=(m>0);}
-private:
-    void calc_gravity()
-    {
-        double r2 = SQR(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);
-    }
-    double m, b;
-};
-
-class Miyamoto_Nagai : public External_gravity {
-public:
-    Miyamoto_Nagai(double m, double a, double b) : m(m), a(a), b(b) {is_active=(m>0);}
-private:
-    void 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 + SQR(b);
-        auto z_tmp  = sqrt(z2_tmp);
-        auto r2_tmp = x_ij*x_ij + y_ij*y_ij + SQR(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*SQR(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*SQR(z_tmp + a));
-        jerk[2] = tmp * (- vz_ij*(z_tmp + a)*(x_ij*x_ij*(z2_tmp*z_tmp + a*SQR(b)) +
-                y_ij*y_ij*(z2_tmp*z_tmp + a*SQR(b)) -
-                (2.0*a*(SQR(z_ij) - SQR(b))*z_tmp +
-                2.0*SQR(z_ij*z_ij) +
-                SQR(b)*SQR(z_ij) -
-                SQR(b)*(SQR(a) + SQR(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;
-    }
-    double m, a, b;
-};
-
-class Logarithmic_halo : public External_gravity {
-public:
-    Logarithmic_halo(double v_halo, double r_halo) : v2_halo(v_halo*v_halo), r2_halo(r_halo*r_halo) {is_active=(r_halo>0);}
-private:
-    void 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);
-    }
-    double v2_halo, r2_halo;
-};
-
 void calc_self_grav(double t, double eps2, double &g6_calls, int n_loc,
                     int n_act, int ind_act[],
                     double3 x_act_new[], double3 v_act_new[],
@@ -555,90 +404,10 @@ int ni = n_act; // TODO redundant?
 
 #ifdef TIMING
 get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0);
-#endif
-
-#ifdef EXTPOT_GAL_DEH
-
-for (i=0; i<ni; i++)
-    {
-
-    x_ij = x_act_new[i][0];
-    y_ij = x_act_new[i][1];
-    z_ij = x_act_new[i][2];
-
-    vx_ij = v_act_new[i][0];
-    vy_ij = v_act_new[i][1];
-    vz_ij = v_act_new[i][2];
-
-    tmp_r = sqrt( SQR(x_ij) + SQR(y_ij) + SQR(z_ij) );
-    tmp_r2 = SQR(tmp_r);
-    tmp_r3 = POW3(tmp_r);
-
-    dum   = tmp_r/(tmp_r + r_ext);
-    dum2  = SQR(dum);
-    dum3  = POW3(dum);
-    dum_g = pow( dum , -g_ext );
-
-    pot_act_ext[i] -= ((m_ext/r_ext) / (2.0-g_ext) ) * ( 1.0 - dum2 * dum_g );
-
-    tmp = (m_ext/tmp_r3) * dum3 * dum_g;
-
-    a_act_new[i][0] -= tmp * x_ij;
-    a_act_new[i][1] -= tmp * y_ij;
-    a_act_new[i][2] -= tmp * z_ij;
-
-    rv_ij = ( vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij );
-
-    tmp = ( m_ext/POW3(tmp_r+r_ext) ) * dum_g;
-    tmp_g = ((r_ext*g_ext + 3.0*tmp_r)/(tmp_r2*(tmp_r+r_ext)) ) * rv_ij;
-
-    adot_act_new[i][0] += tmp * ( tmp_g * x_ij - vx_ij );
-    adot_act_new[i][1] += tmp * ( tmp_g * y_ij - vy_ij );
-    adot_act_new[i][2] += tmp * ( tmp_g * z_ij - vz_ij );
-
-    } /* i */
-
 #endif
 
     /* For simple Plummer potential... */
 
-#ifdef EXTPOT_BH
-
-for (i=0; i<ni; i++)
-    {
-
-    x_ij = x_act_new[i][0];
-    y_ij = x_act_new[i][1];
-    z_ij = x_act_new[i][2];
-
-    vx_ij = v_act_new[i][0];
-    vy_ij = v_act_new[i][1];
-    vz_ij = v_act_new[i][2];
-
-    r2 = SQR(b_bh);
-
-    r2 += x_ij*x_ij + y_ij*y_ij + z_ij*z_ij;   r = sqrt(r2);
-
-    rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
-
-    tmp = m_bh / r;
-
-    pot_act_ext[i] -= tmp;
-
-    tmp /= r2;
-
-    a_act_new[i][0] -= tmp * x_ij;
-    a_act_new[i][1] -= tmp * y_ij;
-    a_act_new[i][2] -= tmp * z_ij;
-
-    adot_act_new[i][0] -= tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
-    adot_act_new[i][1] -= tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
-    adot_act_new[i][2] -= tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
-
-    } /* i */
-
-#endif
-
 #ifdef TIMING
 get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst);
 DT_EXT_GRAV += (CPU_tmp_user - CPU_tmp_user0);
@@ -895,16 +664,6 @@ int main(int argc, char *argv[])
         inp_data : name of the input file (data.inp)
         */
 
-        /* Read the data for EXTernal POTential */
-
-#ifdef EXTPOT_BH
-        fscanf(inp,"%lE %lE", &m_bh, &b_bh);
-#endif
-
-#ifdef EXTPOT_GAL_DEH
-        fscanf(inp,"%lE %lE %lE", &m_ext, &r_ext, &g_ext);
-#endif
-
         printf("\n");
         printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc);
         printf("\n");
@@ -915,25 +674,9 @@ int main(int argc, char *argv[])
         printf("eta     = %.6E \n", eta);
         printf("\n");
 
-#ifdef NORM
-        printf("Normalization: \n");
-        printf("\n");
-        printf("m_norm = %.4E [Msol]  r_norm = %.4E [pc] \n", m_norm/Msol, r_norm/pc);
-        printf("v_morm = %.4E [km/s]  t_morm = %.4E [Myr] \n", v_norm/km, t_norm/Myr);
-        printf("\n");
-#endif
-
         printf("External Potential: \n");
         printf("\n");
 
-#ifdef EXTPOT_BH
-        printf("m_bh = %.4E  b_bh = %.4E \n", m_bh, b_bh);
-#endif
-
-#ifdef EXTPOT_GAL_DEH
-        printf("m_ext  = %.6E r_ext  = %.6E \t g_ext = %.3E \n", m_ext, r_ext, g_ext);
-#endif
-
         if ((diskstep == 0) && (time_cur == 0.0)) {
             out = fopen("contr.dat", "w");
             fclose(out);
@@ -977,24 +720,6 @@ int main(int argc, char *argv[])
     MPI_Bcast(&dt_bh,    1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
     MPI_Bcast(&time_cur, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
 
-#ifdef NORM
-    MPI_Bcast(&m_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Bcast(&r_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Bcast(&v_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Bcast(&t_norm, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-#endif        // NORM
-
-#ifdef EXTPOT_BH
-    MPI_Bcast(&m_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Bcast(&b_bh, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-#endif
-
-#ifdef EXTPOT_GAL_DEH
-    MPI_Bcast(&m_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Bcast(&r_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-    MPI_Bcast(&g_ext, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD);
-#endif
-
     /* Wait to all processors to finish his works... */
     MPI_Barrier(MPI_COMM_WORLD);
 
@@ -1008,17 +733,20 @@ int main(int argc, char *argv[])
     Miyamoto_Nagai ext_disk(config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length);
     Plummer ext_halo_plummer(config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
     Logarithmic_halo ext_log_halo(config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);
+    Dehnen ext_dehnen(config->ext_dehnen_m*normalization_mass, config->ext_dehnen_r*normalization_length, config->ext_dehnen_gamma);
 
     std::vector<External_gravity*> external_gravity_components;
     external_gravity_components.push_back(&ext_bulge);
     external_gravity_components.push_back(&ext_disk);
     external_gravity_components.push_back(&ext_halo_plummer);
     external_gravity_components.push_back(&ext_log_halo);
+    external_gravity_components.push_back(&ext_dehnen);
 
     if (ext_bulge.is_active)        printf("m_bulge = %.4E                        b_bulge = %.4E\n", config->ext_m_bulge*normalization_mass, config->ext_b_bulge*normalization_length);
     if (ext_disk.is_active)         printf("m_disk  = %.4E  a_disk  = %.4E  b_disk  = %.4E\n", config->ext_m_disk*normalization_mass, config->ext_a_disk*normalization_length, config->ext_b_disk*normalization_length);
     if (ext_halo_plummer.is_active) printf("m_halo  = %.4E                        b_halo = %.4E\n", config->ext_m_halo_plummer*normalization_mass, config->ext_b_halo_plummer*normalization_length);
     if (ext_log_halo.is_active)     printf("v_halo = %.6E r_halo = %.6E \n", config->ext_log_halo_v*normalization_velocity, config->ext_log_halo_r*normalization_length);
+    if (ext_dehnen.is_active)       printf("m_ext  = %.6E r_ext  = %.6E \t g_ext = %.3E \n", config->ext_dehnen_m*normalization_mass, config->ext_dehnen_r*normalization_length, config->ext_dehnen_gamma);
     printf("\n");
     fflush(stdout);