diff --git a/src/disc/disc_derivatives.c b/src/disc/disc_derivatives.c
index 01026016a457704fdfe13e6454f977338bdcf825..b2997b3c09244eb1920285e4b20559ccea838eb3 100644
--- a/src/disc/disc_derivatives.c
+++ b/src/disc/disc_derivatives.c
@@ -713,8 +713,8 @@ void disc_derivatives(struct disc_t * const disc,
      * Appendix B3: Angular momentum flux correction because
      * of mass loss
      */
-    if(1)disc_edge_loss_angular_momentum_flux(disc,
-                                              binary);
+    disc_edge_loss_angular_momentum_flux(disc,
+                                         binary);
 }
 
 #endif // DISCS
diff --git a/src/disc/disc_parameters.h b/src/disc/disc_parameters.h
index bcb19915b2ad69c9d816b87a873c3a3d1f282bf3..c71420ace2e72deb9e13a53ca363446574a3cbb9 100644
--- a/src/disc/disc_parameters.h
+++ b/src/disc/disc_parameters.h
@@ -33,7 +33,7 @@
 /*
  * Disc equation checks
  */
-//#define DISC_EQUATION_CHECKS
+#define DISC_EQUATION_CHECKS
 
 /*
  * Log equation checks
diff --git a/src/disc/disc_structure_functions.c b/src/disc/disc_structure_functions.c
index 7dad8579ab1f5df66130b5a59e95f7b83ebc9250..31e2dc1e52a20c5868e7bc7c6bebe3aa901eccb1 100644
--- a/src/disc/disc_structure_functions.c
+++ b/src/disc/disc_structure_functions.c
@@ -800,7 +800,7 @@ int disc_n_valid_zones(const struct disc_t * const disc)
 void disc_equation_checker(const double r,
                            struct disc_t * const disc,
                            const struct binary_system_t * const binary,
-                           double * const model,
+                           double * const model MAYBE_UNUSED,
                            double * const star)
 {
     /*
@@ -823,7 +823,7 @@ void disc_equation_checker(const double r,
     /*
      * Disc temperature as a function of radius
      */
-    double T_LHS = disc_temperature(r,disc);
+    const double T_LHS = disc_temperature(r,disc);
 
     /*
      * Evaluate T^4 RHS from disc model
@@ -837,27 +837,32 @@ void disc_equation_checker(const double r,
     const double Hlow = disc_scale_height(rlow,disc,binary);
     const double dHdR = (H-Hlow) / (r - rlow);
     const double dlnHdlnR = (r/H) * dHdR;
-    const double A = 27.0 * disc->kappa * nu * Pow2(Sigma * omega) / (64.0 * STEFAN_BOLTZMANN_CONSTANT);
+
     double mdot = 0.0;
-    if(1)
+    /*
+     * Use non-edge mass loss
+     */
+    for(int i=0;i<DISC_LOSS_N;i++)
     {
-        /*
-         * Use non-edge mass loss
-         */
-        int i;
-        for(i=0;i<DISC_LOSS_N;i++)
+        if(!Disc_loss_strip_from_edge(i))
         {
-            if(!Disc_loss_strip_from_edge(i))
-            {
-                mdot += disc->loss[i].mdot;
-            }
+            mdot += disc->loss[i].mdot;
         }
     }
-    const double B = GRAVITATIONAL_CONSTANT * binary->mtot * mdot /
+
+    /*
+     * A, B, C, D are Eqs. 4, 6, 9 and 10
+     */
+    const double A = 27.0 * disc->kappa * nu * Pow2(Sigma * omega) / (64.0 * STEFAN_BOLTZMANN_CONSTANT);
+    const double B = 2.0/(3.0*PI) * binary->flux * (1.0 - disc->albedo) * Pow3(binary->Rstar / r) / STEFAN_BOLTZMANN_CONSTANT;
+    const double C = 3.0 * PI * H / (4.0 * binary->Rstar) * (dlnHdlnR - 1.0);
+    const double D = GRAVITATIONAL_CONSTANT * binary->mtot * mdot /
         (2.0 * PI * binary->separation * Pow2(r) * STEFAN_BOLTZMANN_CONSTANT);
-    const double C = 2.0/(3.0*PI) * binary->flux * (1.0 - disc->albedo) * Pow3(binary->Rstar / r) / STEFAN_BOLTZMANN_CONSTANT;
-    const double D = 3.0 * PI * H / (4.0 * binary->Rstar) * (dlnHdlnR - 1.0);
-    const double T4_RHS = A + 0.0*B + C * (1.0 + D);
+
+    /*
+     * The usual RHS : note that the D term does not contribute
+     */
+    const double T4_RHS = A + 0.0*D + B * (1.0 + C);
 #ifdef DISC_EQUATION_CHECKS_LOG
     const double T4_LHS = Pow4(T_LHS);
 #endif // DISC_EQUATION_CHECKS_LOG
@@ -868,10 +873,10 @@ void disc_equation_checker(const double r,
     disc->equations_T_error_pc = Max(disc->equations_T_error_pc,
                                      fabs(dT/T_LHS)*100.0);
     disc->equations_mass_loss_term_pc = Max(disc->equations_mass_loss_term_pc,
-                                            fabs(B/(A+C*(1+D)))*100.0);
+                                            fabs(D/(A+B*(1+C)))*100.0);
 
 #ifdef DISC_EQUATION_CHECKS_LOG
-    printf("EQN2 r=%12g T_LHS=%12g T^4(LHS) = %12g == (A=%12g) + (B=%12g) + (C=%12g) * (1 + (D=%12g)) = %12g :: T4 ratio %12g : dT = %12g : dT/T = %5.2f %% : B/(A+C*(1+D)) = %5.2f %%\n",
+    printf("EQN2 r=%12g T_LHS=%12g T^4(LHS) = %12g == (A=%12g) + (B=%12g) * (1 + (C=%12g)) + (D=%12g) = %12g :: T4 ratio %12g : dT = %12g : dT/T = %5.2f %% : D/(A+B*(1+C)) = %5.2f %%\n",
            r,
            T_LHS,
            T4_LHS,
diff --git a/src/setup/version.c b/src/setup/version.c
index 138dca4e669f153f413bce1d8dd3f1b13821b206..acd05b5aed7a5d61059300026ef391ecb1c9d971 100644
--- a/src/setup/version.c
+++ b/src/setup/version.c
@@ -1435,6 +1435,7 @@ void version(struct stardata_t * RESTRICT const stardata)
     Show_int_macro(DISCS_NUMBER_OF_POWER_LAWS);
     Show_int_macro(NDISCS);
     Show_int_macro(DISC_DEBUG);
+    Macrotest(DISC_EQUATION_CHECKS);
     Show_float_macro(DISC_MAX_LIFETIME);
     Show_float_macro(DISC_MIN_TIMESTEP);
     Show_float_macro(DISC_MAX_TIMESTEP);