Inserted ccosmological factors in the Gadget2 force interactions.

src/adiabatic_index.h

@@ -437,7 +437,6 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) {
  * Computes \f$x^{3\gamma - 2}\f$.
  *
  * @param x Argument
- * @return Argument to the power one over the adiabatic index
  */
 __attribute__((always_inline)) INLINE static float pow_three_gamma_minus_two(
     float x) {
@@ -466,4 +465,39 @@ __attribute__((always_inline)) INLINE static float pow_three_gamma_minus_two(
 #endif
 }
 
+/**
+ * @brief Return the argument to the power three adiabatic index minus five over
+ * two.
+ *
+ * Computes \f$x^{(3\gamma - 5)/2}\f$.
+ *
+ * @param x Argument
+ */
+__attribute__((always_inline)) INLINE static float
+pow_three_gamma_minus_five_over_two(float x) {
+
+#if defined(HYDRO_GAMMA_5_3)
+
+  return 1.f; /* x^(0) */
+
+#elif defined(HYDRO_GAMMA_7_5)
+
+  return powf(x, -0.4f); /* x^(-2/5) */
+
+#elif defined(HYDRO_GAMMA_4_3)
+
+  return 1.f / sqrtf(x); /* x^(-1/2) */
+
+#elif defined(HYDRO_GAMMA_2_1)
+
+  return sqrtf(x); /* x^(1/2) */
+
+#else
+
+  error("The adiabatic index is not defined !");
+  return 0.f;
+
+#endif
+}
+
 #endif /* SWIFT_ADIABATIC_INDEX_H */

src/hydro/Gadget2/hydro_iact.h

@@ -436,9 +436,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
 
   float wi, wj, wi_dx, wj_dx;
 
-  /* Will change with cosmological integration */
-  const float fac_mu = 1.f;
-  const float a2_Hubble = 0.f;
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
 
   const float r = sqrtf(r2);
   const float r_inv = 1.0f / r;
@@ -555,9 +555,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
 
   float wi, wj, wi_dx, wj_dx;
 
-  /* Will change with cosmological integration */
-  const float fac_mu = 1.f;
-  const float a2_Hubble = 0.f;
+  /* Cosmological factors entering the EoMs */
+  const float fac_mu = pow_three_gamma_minus_five_over_two(a);
+  const float a2_Hubble = a * a * H;
 
   const float r = sqrtf(r2);
   const float r_inv = 1.0f / r;