diff --git a/examples/IsothermalPotential/GravityOnly/test.pro b/examples/IsothermalPotential/GravityOnly/test.pro
index ac1f5e915c38c107408ffed9579c52944295f079..edfa50121d2e5adb7e039f3c38d6d4c0b4d5e34f 100644
--- a/examples/IsothermalPotential/GravityOnly/test.pro
+++ b/examples/IsothermalPotential/GravityOnly/test.pro
@@ -127,7 +127,7 @@ print,' relative Lz change: (per cent) ',minmax(dl) * 100.
if(iplot eq 1) then begin
; plot results on energy conservation for some particles
nplot = min(10, nfollow)
- wset,0
+ win,0
xr = [min(tout), max(tout)]
yr = [-2,2]*1d-2 ; in percent
plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='time',ytitle='dE/E, dL/L (%)'
@@ -137,7 +137,7 @@ if(iplot eq 1) then begin
screen_to_png,'e-time.png'
; plot orbits of those particles
- wset,2
+ win,2
xr = [-100,100]
yr = xr
plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/iso,/nodata,xtitle='x',ytitle='y'
@@ -146,7 +146,7 @@ if(iplot eq 1) then begin
screen_to_png,'orbit.png'
; plot radial position of these particles
- wset,4
+ win,4
xr = [min(tout), max(tout)]
yr = [0,80]
plot,[0],[0],xr=xr,yr=yr,/xs,/ys,/nodata,xtitle='t',ytitle='r'
@@ -155,7 +155,7 @@ for i=0,nplot-1 do begin dr = sqrt(reform(xout[i,*])^2 + reform(yout[i,*])^2) &
screen_to_png,'r-time.png'
; make histogram of energy changes at end
- wset,6
+ win,6
ohist,de,x,y,-0.05,0.05,0.001
plot,x,y,psym=10,xtitle='de (%)'
screen_to_png,'de-hist.png'
diff --git a/src/const.h b/src/const.h
index 498b60b711e88e8054b5e0e73a18475d98277448..34db85df64e023ca929b70e193747a3d472b98ff 100644
--- a/src/const.h
+++ b/src/const.h
@@ -65,8 +65,9 @@
#define const_gravity_eta 0.025f
/* External gravity properties */
-#define EXTERNAL_POTENTIAL_POINTMASS
+//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+#define EXTERNAL_POTENTIAL_DISK_PATCH
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 5903ce1384f616a123c74b579539fe53747d17e4..74a68fdc335b95e393c1cfad6cd8d07c32ce5458 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -46,7 +46,9 @@ gravity_compute_timestep_external(const struct external_potential* potential,
dt = fminf(dt, external_gravity_isothermalpotential_timestep(potential,
phys_const, gp));
#endif
-
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ dt = fmin(dt, external_gravity_disk_patch_timestep(potential, phys_const, gp));
+#endif
return dt;
}
@@ -110,7 +112,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
* @param const_G Newton's constant
*/
__attribute__((always_inline)) INLINE static void gravity_end_force(
- struct gpart* gp, double const_G) {
+ struct gpart* gp, const double const_G) {
/* Let's get physical... */
gp->a_grav[0] *= const_G;
@@ -137,6 +139,9 @@ __attribute__((always_inline)) INLINE static void external_gravity(
#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
external_gravity_isothermalpotential(potential, phys_const, gp);
#endif
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ external_gravity_disk_patch_potential(potential, phys_const, gp);
+#endif
}
/**
diff --git a/src/potentials.c b/src/potentials.c
index 5cbe05e008b7de17310b1e738e032af90684c25e..9a4480adadf8186e34174fa1c4d0b766a4a03d4c 100644
--- a/src/potentials.c
+++ b/src/potentials.c
@@ -65,6 +65,17 @@ void potential_init(const struct swift_params* parameter_file,
parameter_file, "IsothermalPotential:timestep_mult");
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ potential -> disk_patch_potential.surface_density =
+ parser_get_param_double(parameter_file, "Disk-PatchPotential:surface_density");
+ potential -> disk_patch_potential.scale_height =
+ parser_get_param_double(parameter_file, "Disk-PatchPotential:scale_height");
+ potential -> disk_patch_potential.z_disk =
+ parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
+ potential -> disk_patch_potential.timestep_mult =
+ parser_get_param_double(parameter_file, "Disk-PatchPotential:timestep_mult");
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
+
}
/**
@@ -94,4 +105,11 @@ void potential_print(const struct external_potential* potential) {
potential->isothermal_potential.timestep_mult);
#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ message("Disk-patch potential properties are surface_density = %e disk height= %e scale height= %e timestep multiplier= %e",
+ potential->disk_patch_potential.surface_density,
+ potential->disk_patch_potential.z_disk,
+ potential->disk_patch_potential.scale_height,
+ potential->disk_patch_potential.timestep_mult);
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
}
diff --git a/src/potentials.h b/src/potentials.h
index 5373cc1f3f55b0c6d9876cbe6348fdcefbe242aa..83bdaf6fc817344f97612a8e8ddfc42c2d7ffcfa 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -34,6 +34,7 @@
#include "part.h"
#include "physical_constants.h"
#include "units.h"
+#include "math.h"
/* External Potential Properties */
struct external_potential {
@@ -53,9 +54,66 @@ struct external_potential {
float timestep_mult;
} isothermal_potential;
#endif
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+ struct {
+ double surface_density;
+ double scale_height;
+ double z_disk;
+ double timestep_mult;
+ } disk_patch_potential;
+#endif
};
-/* Include external isothermal potential */
+/* external potential: disk_patch */
+#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
+/**
+ * @brief Computes the time-step due to the acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
+ *
+ * @param phys_cont The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
+ */
+__attribute__((always_inline))
+ INLINE static float external_gravity_disk_patch_timestep(
+ const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const g) {
+
+ const float dz = fabs(g->x[2] - potential->disk_patch_potential.z_disk);
+ const float z_accel = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density
+ * tanh(dz / potential->disk_patch_potential.scale_height);
+ const float dz_accel_over_dt = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density /
+ potential->disk_patch_potential.scale_height / cosh(dz / potential->disk_patch_potential.scale_height) / cosh(dz / potential->disk_patch_potential.scale_height) * fabs(g->v_full[2]);
+
+ /* demand that time step * derivative of acceleration < fraction of acceleration */
+ const float dt1 = potential->disk_patch_potential.timestep_mult * z_accel / dz_accel_over_dt;
+
+ /* demand that dt^2 * acceleration < fraction of scale height of disk */
+ const float dt2 = potential->disk_patch_potential.timestep_mult * potential->disk_patch_potential.scale_height / fabs(z_accel);
+
+ float dt = dt1;
+ if(dt2 < dt) dt = dt2;
+ return dt2;
+}
+/**
+ * @brief Computes the gravitational acceleration from a hydrostatic disk (Creasy, Theuns & Bower, 2013)
+ *
+ * @param phys_cont The physical constants in internal units.
+ * @param gp Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_disk_patch_potential(const struct external_potential* potential, const struct phys_const* const phys_const, struct gpart* g) {
+
+ const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
+ g->a_grav[0] = 0;
+ g->a_grav[1] = 0;
+ const float z_accel = 2 * M_PI * phys_const->const_newton_G * potential->disk_patch_potential.surface_density
+ * tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
+ if (dz > 0)
+ g->a_grav[2] -= z_accel;
+ if (dz < 0)
+ g->a_grav[2] += z_accel;
+}
+#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
+
#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
/**