diff --git a/src/const.h b/src/const.h
index fbd46ac488a6e63d8dcc49640c25344ebd360f10..7e08e6b266c0d3ee0c51052daf1d69298abb4518 100644
--- a/src/const.h
+++ b/src/const.h
@@ -20,6 +20,8 @@
#ifndef SWIFT_CONST_H
#define SWIFT_CONST_H
+#include "physical_constants.h"
+
/* Hydrodynamical constants. */
#define const_hydro_gamma (5.0f / 3.0f)
@@ -56,7 +58,6 @@
0.57735f /* Opening criteria, which is the ratio of the \
cell distance over the cell width. */
-#define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */
#define const_iepsilon 714.285714286f /* Inverse gravity blending distance. */
#define const_iepsilon2 (const_iepsilon* const_iepsilon)
@@ -79,8 +80,18 @@
/* System of units */
-#define const_unit_length_in_cgs 1 /* 3.08567810e16 /\* 1Mpc *\/ */
-#define const_unit_mass_in_cgs 1 /* 1.9891e33 /\* 1 M_sun *\/ */
-#define const_unit_velocity_in_cgs 1 /* 1e5 /\* km s^-1 *\/ */
+#define const_unit_length_in_cgs (1000 * PARSEC_IN_CGS ) /* kpc */
+#define const_unit_mass_in_cgs (SOLAR_MASS_IN_CGS) /* solar mass */
+#define const_unit_velocity_in_cgs (1e5) /* km/s */
+
+/* Derived constants */
+#define const_unit_time_in_cgs (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
+#define const_G (NEWTON_GRAVITY_CGS * const_unit_mass_in_cgs*const_unit_time_in_cgs*const_unit_time_in_cgs/ (const_unit_length_in_cgs*const_unit_length_in_cgs*const_unit_length_in_cgs))
+
+/* External Potential Constants */
+#define External_Potential_X (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs)
+#define External_Potential_Y (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs)
+#define External_Potential_Z (50000 * PARSEC_IN_CGS / const_unit_length_in_cgs)
+#define External_Potential_Mass (1e10 * SOLAR_MASS_IN_CGS / const_unit_mass_in_cgs)
#endif /* SWIFT_CONST_H */
diff --git a/src/engine.c b/src/engine.c
index b269f1b94fc4550e278a38a7d402541b27d92d71..d84d74c4391fc56e016e8f77655a4032e66aa50d 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -119,6 +119,14 @@ void engine_mkghosts(struct engine *e, struct cell *c, struct cell *super) {
/* Add the kick task. */
c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0, c,
NULL, 0);
+
+ /* /\* Add the gravity tasks *\/ */
+ /* c->grav_external = scheduler_addtask(s, task_type_grav_external, task_subtype_none, 0, 0, */
+ /* c, NULL, 0); */
+
+ /* /\* Enforce gravity calculated before kick *\/ */
+ /* scheduler_addunlock(s, c->grav_external, c->kick); */
+
}
}
@@ -775,6 +783,7 @@ void engine_maketasks(struct engine *e) {
}
}
+#ifdef GRAVITY
#ifdef DEFAULT_GRAVITY
/* Add the gravity mm tasks. */
for (i = 0; i < nr_cells; i++)
@@ -786,6 +795,20 @@ void engine_maketasks(struct engine *e) {
scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
&cells[i], &cells[j], 0);
}
+#endif
+/* #ifdef EXTERNAL_POTENTIAL */
+/* /\* Add the external potential gravity *\/ */
+/* for (i = 0; i < nr_cells; i++) */
+/* if (cells[i].gcount > 0) { */
+/* scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0, */
+/* &cells[i], NULL, 0); */
+/* for (j = i + 1; j < nr_cells; j++) */
+/* if (cells[j].gcount > 0) */
+/* scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0, */
+/* &cells[i], &cells[j], 0); */
+/* } */
+
+/* #endif */
#endif
scheduler_print_tasks(sched, "sched.txt");
diff --git a/src/gravity/Default/runner_iact_grav.h b/src/gravity/Default/runner_iact_grav.h
index e62be446e8263bf02e3fd73f902b28cb1c3b16cf..11730e1148652d5b4a512b380294e9176cec3f3a 100644
--- a/src/gravity/Default/runner_iact_grav.h
+++ b/src/gravity/Default/runner_iact_grav.h
@@ -21,9 +21,10 @@
#define SWIFT_RUNNER_IACT_GRAV_H
/* Includes. */
-#include "const.h"
-#include "kernel.h"
-#include "vector.h"
+#include "../../const.h"
+#include "../../kernel.h"
+#include "../../vector.h"
+#include "../../timers.h"
/**
* @file runner_iact_grav.h
@@ -31,6 +32,8 @@
*
*/
+
+
/**
* @brief Gravity potential
*/
diff --git a/src/gravity/ExternalPotential/runner_iact_grav.h b/src/gravity/ExternalPotential/runner_iact_grav.h
index e62be446e8263bf02e3fd73f902b28cb1c3b16cf..16ea59c599b0352cd8add7a5efa1a3dbdbfb91c1 100644
--- a/src/gravity/ExternalPotential/runner_iact_grav.h
+++ b/src/gravity/ExternalPotential/runner_iact_grav.h
@@ -26,7 +26,7 @@
#include "vector.h"
/**
- * @file runner_iact_grav.h
+ * @file runner_iact_grav.c
* @brief Gravity interaction functions.
*
*/
diff --git a/src/runner.c b/src/runner.c
index 1cd52e1dd89f08172c3a3b679f15a096a0007b1b..a33bbd18b8cadb7432b52ba743205855b8802ff7 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -68,6 +68,11 @@ const float runner_shift[13 * 3] = {
const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+#define OUT if(CELL_ID == MY_CELL) {message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time);}
+//#define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time);
+//#define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d %f\n",c->loc[0],c->loc[1],c->loc[2], CELL_ID, c->count, r->e->time);
+
+
/* Import the density loop functions. */
#define FUNCTION density
#include "runner_doiact.h"
@@ -80,6 +85,69 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
/* Import the gravity loop functions. */
#include "runner_doiact_grav.h"
+
+/**
+ * @brief Calculate gravity acceleration from external potential
+ *
+ * @param runner task
+ * @param cell
+ */
+void runner_dograv_external(struct runner *r, struct cell *c) {
+
+ struct part *p, *parts = c->parts;
+ float rinv;
+ float L[3], E;
+ int i, k, count = c->count;
+ float t_end = r->e->time;
+ //struct space *s = r->e->s;
+ //double CentreOfPotential[3];
+ TIMER_TIC
+
+ /* /\* location of external gravity point mass - should pass in as paraneter *\/ */
+ /* CentreOfPotential[0] = 0.5 * s->dim[0]; */
+ /* CentreOfPotential[1] = 0.5 * s->dim[1]; */
+ /* CentreOfPotential[2] = 0.5 * s->dim[2]; */
+
+ /* Recurse? */
+ if (c->split) {
+ for (k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) runner_dograv_external(r, c->progeny[k]);
+ return;
+ }
+
+#ifdef TASK_VERBOSE
+ OUT
+#endif
+ /* /\* Loop over the parts in this cell. *\/ */
+ /* for (i = 0; i < count; i++) { */
+
+ /* /\* Get a direct pointer on the part. *\/ */
+ /* p = &parts[i]; */
+
+ /* /\* Is this part within the time step? *\/ */
+ /* if (p->t_end <= t_end) { */
+ /* rinv = 1 / sqrtf((p->x[0]-External_Potential_X)*(p->x[0]-External_Potential_X) + (p->x[1]-External_Potential_Y)*(p->x[1]-External_Potential_Y) + (p->x[2]-External_Potential_Z)*(p->x[2]-External_Potential_Z)); */
+
+ /* /\* check for energy and angular momentum conservation *\/ */
+ /* E = 0.5 * ((p->v[0]*p->v[0]) + (p->v[1]*p->v[1]) + (p->v[2]*p->v[2])) - const_G * External_Potential_Mass * rinv; */
+ /* L[0] = (p->x[1] - External_Potential_X)*p->v[2] - (p->x[2] - External_Potential_X)*p->v[1]; */
+ /* L[1] = (p->x[2] - External_Potential_Y)*p->v[0] - (p->x[0] - External_Potential_Y)*p->v[2]; */
+ /* L[2] = (p->x[0] - External_Potential_Z)*p->v[1] - (p->x[1] - External_Potential_Z)*p->v[0]; */
+ /* if(p->id == 0) { */
+ /* message("update %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f\n", r->e->time, 1./rinv, p->x[0], p->x[1], p->x[2], E, L[0], L[1], L[2]); */
+ /* message(" ... %f %f %f\n", p->v[0], p->v[1], p->v[2]); */
+ /* message(" ... %e %e\n", const_G, External_Potential_Mass); */
+ /* } */
+ /* p->grav_accel[0] = - const_G * External_Potential_Mass * (p->x[0] - External_Potential_X) * rinv * rinv * rinv; */
+ /* p->grav_accel[1] = - const_G * External_Potential_Mass * (p->x[1] - External_Potential_Y) * rinv * rinv * rinv; */
+ /* p->grav_accel[2] = - const_G * External_Potential_Mass * (p->x[2] - External_Potential_Z) * rinv * rinv * rinv; */
+
+ /* } */
+ /* } */
+ /* TIMER_TOC(timer_dograv_external); */
+}
+
+
/**
* @brief Sort the entries in ascending order using QuickSort.
*
diff --git a/src/task.c b/src/task.c
index 69109f9e6d4fe8730a317db46ea3862e65ab90b2..e93ce131a8823c567fae7cb9f8024e23154e9e4b 100644
--- a/src/task.c
+++ b/src/task.c
@@ -45,7 +45,7 @@
const char *taskID_names[task_type_count] = {
"none", "sort", "self", "pair", "sub", "init",
"ghost", "drift", "kick", "send", "recv", "grav_pp",
- "grav_mm", "grav_up", "grav_down", "psort", "split_cell", "rewait"};
+ "grav_mm", "grav_up", "grav_down", "grav_external", "psort", "split_cell", "rewait"};
const char *subtaskID_names[task_type_count] = {"none", "density",
"force", "grav"};
diff --git a/src/task.h b/src/task.h
index b86631cc49bfad102302e3bab380bfb5eb8ed1e0..27b1b1fe9b0f899a368597c937e3074b411b24a0 100644
--- a/src/task.h
+++ b/src/task.h
@@ -25,6 +25,7 @@
#include "cycle.h"
/* Some constants. */
+//#define TASK_VERBOSE
#define task_maxwait 3
#define task_maxunlock 15
@@ -45,6 +46,7 @@ enum task_types {
task_type_grav_mm,
task_type_grav_up,
task_type_grav_down,
+ task_type_grav_external,
task_type_psort,
task_type_split_cell,
task_type_rewait,