diff --git a/src/black_holes/EAGLE/black_holes.h b/src/black_holes/EAGLE/black_holes.h
index 2720058ea48e053d2cb5b2ee5f115373fd8e1b24..67333955e8216401573b0841414e617e5adc1079 100644
--- a/src/black_holes/EAGLE/black_holes.h
+++ b/src/black_holes/EAGLE/black_holes.h
@@ -20,9 +20,11 @@
 #define SWIFT_EAGLE_BLACK_HOLES_H
 
 #include <float.h>
+#include "cosmology.h"
 #include "dimension.h"
 #include "kernel_hydro.h"
 #include "minmax.h"
+#include "physical_constants.h"
 
 /**
  * @brief Computes the gravity time-step of a given black hole particle.
@@ -47,6 +49,8 @@ __attribute__((always_inline)) INLINE static void black_holes_first_init_bpart(
     struct bpart* bp) {
 
   bp->time_bin = 0;
+  bp->subgrid_mass = bp->mass;
+  bp->energy_reservoir = 0.;
 }
 
 /**
@@ -65,6 +69,12 @@ __attribute__((always_inline)) INLINE static void black_holes_init_bpart(
 
   bp->density.wcount = 0.f;
   bp->density.wcount_dh = 0.f;
+  bp->rho_gas = 0.f;
+  bp->sound_speed_gas = 0.f;
+  bp->velocity_gas[0] = 0.f;
+  bp->velocity_gas[1] = 0.f;
+  bp->velocity_gas[2] = 0.f;
+  bp->ngb_mass = 0.f;
 }
 
 /**
@@ -112,6 +122,11 @@ __attribute__((always_inline)) INLINE static void black_holes_end_density(
   /* Finish the calculation by inserting the missing h-factors */
   bp->density.wcount *= h_inv_dim;
   bp->density.wcount_dh *= h_inv_dim_plus_one;
+  bp->rho_gas *= h_inv_dim;
+  bp->sound_speed_gas *= h_inv_dim;
+  bp->velocity_gas[0] *= h_inv_dim;
+  bp->velocity_gas[1] *= h_inv_dim;
+  bp->velocity_gas[2] *= h_inv_dim;
 }
 
 /**
@@ -135,6 +150,58 @@ black_holes_bpart_has_no_neighbours(struct bpart* restrict bp,
   bp->density.wcount_dh = 0.f;
 }
 
+/**
+ *
+ */
+__attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
+    struct bpart* restrict bp, const struct phys_const* constants,
+    const struct cosmology* cosmo, const double dt) {
+
+  const double G = constants->const_newton_G;
+  const double c = constants->const_speed_light_c;
+
+  /* (Subgrid) mass of the BH (internal units) */
+  const double BH_mass = bp->subgrid_mass;
+
+  /* Convert the quantities we gathered to physical frame (all internal units)
+   */
+  const double gas_rho_phys = bp->rho_gas * cosmo->a3_inv;
+  const double gas_c_phys = bp->sound_speed_gas * cosmo->a_factor_sound_speed;
+  const double gas_v_peculiar[3] = {bp->velocity_gas[0] * cosmo->a_inv,
+                                    bp->velocity_gas[1] * cosmo->a_inv,
+                                    bp->velocity_gas[2] * cosmo->a_inv};
+
+  const double bh_v_peculiar[3] = {bp->v[0] * cosmo->a_inv,
+                                   bp->v[1] * cosmo->a_inv,
+                                   bp->v[2] * cosmo->a_inv};
+
+  /* Difference in velocity between the gas and the BH */
+  const double v_diff_peculiar[3] = {gas_v_peculiar[0] - bh_v_peculiar[0],
+                                     gas_v_peculiar[1] - bh_v_peculiar[1],
+                                     gas_v_peculiar[2] - bh_v_peculiar[2]};
+  const double v_diff_norm2 = v_diff_peculiar[0] * v_diff_peculiar[0] +
+                              v_diff_peculiar[1] * v_diff_peculiar[1] +
+                              v_diff_peculiar[2] * v_diff_peculiar[2];
+
+  /* We can now compute the Bondi accretion rate (internal units) */
+  const double gas_c_phys2 = gas_c_phys * gas_c_phys;
+  const double denominator2 = v_diff_norm2 + gas_c_phys2;
+  const double denominator_inv = 1. / sqrt(denominator2);
+  const double Bondi_rate = 4. * M_PI * G * G * BH_mass * BH_mass *
+                            gas_rho_phys * denominator_inv * denominator_inv *
+                            denominator_inv;
+
+  /* Factor in the radiative efficiency */
+  const double epsilon_r = 0.1;
+  const double mass_rate = (1. - epsilon_r) * Bondi_rate;
+  const double luminosity = epsilon_r * Bondi_rate * c * c;
+
+  /* Integrate forward in time */
+  const double epsilon_f = 0.15;
+  bp->subgrid_mass += mass_rate * dt;
+  bp->energy_reservoir += luminosity * epsilon_f * dt;
+}
+
 /**
  * @brief Reset acceleration fields of a particle
  *
diff --git a/src/black_holes/EAGLE/black_holes_iact.h b/src/black_holes/EAGLE/black_holes_iact.h
index a85580a0eca5804f71fa0c30d27f7fba3b969db0..c93fdc6de21141eebbe3f298f645d433d7cd2dbf 100644
--- a/src/black_holes/EAGLE/black_holes_iact.h
+++ b/src/black_holes/EAGLE/black_holes_iact.h
@@ -19,6 +19,8 @@
 #ifndef SWIFT_EAGLE_BH_IACT_H
 #define SWIFT_EAGLE_BH_IACT_H
 
+#include "hydro.h"
+
 /**
  * @brief Density interaction between two particles (non-symmetric).
  *
@@ -51,6 +53,29 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_bh_density(
   bi->density.wcount += wi;
   bi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
 
+  /* Neighbour gas mass */
+  const float mj = hydro_get_mass(pj);
+
+  /* Contribution to the BH gas density */
+  bi->rho_gas += mj * wi;
+
+  /* Contribution to the total neighbour mass */
+  bi->ngb_mass += mj;
+
+  /* Neighbour sounds speed */
+  const float cj = hydro_get_comoving_soundspeed(pj);
+
+  /* Contribution to the smoothed sound speed */
+  bi->sound_speed_gas += cj * wi;
+
+  /* Neighbour peculiar drifted velocity */
+  const float v[3] = {pj->v[0], pj->v[1], pj->v[2]};
+
+  /* Contribution to the smoothed velocity */
+  bi->velocity_gas[0] += v[0] * wi;
+  bi->velocity_gas[1] += v[1] * wi;
+  bi->velocity_gas[2] += v[2] * wi;
+
 #ifdef DEBUG_INTERACTIONS_BH
   /* Update ngb counters */
   if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_BH)
diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h
index 360eb45d9e2ccb14e8ff55b6d286e7bb91ef89cc..0aa311ff8084b87e607f6ccbf3d1cce7b15b53ce 100644
--- a/src/black_holes/EAGLE/black_holes_part.h
+++ b/src/black_holes/EAGLE/black_holes_part.h
@@ -16,11 +16,8 @@
  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  *
  ******************************************************************************/
-#ifndef SWIFT_DEFAULT_BLACK_HOLE_PART_H
-#define SWIFT_DEFAULT_BLACK_HOLE_PART_H
-
-/* Some standard headers. */
-#include <stdlib.h>
+#ifndef SWIFT_EAGLE_BLACK_HOLE_PART_H
+#define SWIFT_EAGLE_BLACK_HOLE_PART_H
 
 /**
  * @brief Particle fields for the black hole particles.
@@ -63,6 +60,24 @@ struct bpart {
 
   } density;
 
+  /*! Subgrid mass of the black hole */
+  float subgrid_mass;
+
+  /*! Energy reservoir for feedback */
+  float energy_reservoir;
+
+  /*! Density of the gas surrounding the black hole. */
+  float rho_gas;
+
+  /*! Smoothed sound speed of the gas surrounding the black hole. */
+  float sound_speed_gas;
+
+  /*! Smoothed velocity (peculiar) of the gas surrounding the black hole */
+  float velocity_gas[3];
+
+  /*! Total mass of the gas neighbours. */
+  float ngb_mass;
+
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Time of the last drift */
@@ -89,4 +104,4 @@ struct bpart {
 
 } SWIFT_STRUCT_ALIGN;
 
-#endif /* SWIFT_DEFAULT_BLACK_HOLE_PART_H */
+#endif /* SWIFT_EAGLE_BLACK_HOLE_PART_H */
diff --git a/src/runner.c b/src/runner.c
index e6be52e3dc5a9e04d36434ef8a43ddc02d7ba976..7385f0c99815b12cf00aef3eab691044c8444dd4 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -524,6 +524,8 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
   /* Running value of the maximal smoothing length */
   double h_max = c->black_holes.h_max;
 
+  double dt = 0.001;
+
   TIMER_TIC;
 
   /* Anything to do here? */
@@ -630,6 +632,11 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
           if (((bp->h >= black_holes_h_max) && (f < 0.f)) ||
               ((bp->h <= black_holes_h_min) && (f > 0.f))) {
 
+            /* Compute variables required for the feedback loop */
+            black_holes_prepare_feedback(bp, e->physical_constants,
+                                         e->cosmology, dt);
+
+            /* Reset quantities computed by the feedback loop */
             black_holes_reset_feedback(bp);
 
             /* Ok, we are done with this particle */
@@ -723,6 +730,11 @@ void runner_do_black_holes_ghost(struct runner *r, struct cell *c, int timer) {
         /* Check if h_max has increased */
         h_max = max(h_max, bp->h);
 
+        /* Compute variables required for the feedback loop */
+        black_holes_prepare_feedback(bp, e->physical_constants, e->cosmology,
+                                     dt);
+
+        /* Reset quantities computed by the feedback loop */
         black_holes_reset_feedback(bp);
       }