diff --git a/camb/tests/camb_test.py b/camb/tests/camb_test.py
index 2c58fb12..35446f9c 100644
--- a/camb/tests/camb_test.py
+++ b/camb/tests/camb_test.py
@@ -706,22 +706,32 @@ def testTimeTransfers(self):
             cls += [results.get_total_cls(CMB_unit="muK")]
         self.assertTrue(np.allclose((cls[1] - cls[0])[2:300, 2] * 2, (cls[2] - cls[0])[2:300, 2], rtol=1e-3))
 
-    def testDarkEnergy(self):
-        pars = camb.CAMBparams()
-        pars.set_cosmology(H0=71)
-        pars.InitPower.set_params(ns=0.965, r=0)
-        for m in ["fluid", "ppf"]:
+    def testDarkEnergy(self):
+        pars = camb.CAMBparams()
+        pars.set_cosmology(H0=71)
+        pars.InitPower.set_params(ns=0.965, r=0)
+        for m in ["fluid", "ppf"]:
             pars.set_dark_energy(w=-0.7, wa=0.2, dark_energy_model=m)
             C1 = camb.get_results(pars).get_cmb_power_spectra()
             a = np.logspace(-5, 0, 1000)
             w = -0.7 + 0.2 * (1 - a)
             pars2 = pars.copy()
             pars2.set_dark_energy_w_a(a, w, dark_energy_model=m)
-            C2 = camb.get_results(pars2).get_cmb_power_spectra()
-            for f in ["lens_potential", "lensed_scalar"]:
-                self.assertTrue(np.allclose(C1[f][2:, 0], C2[f][2:, 0]))
-            pars3 = pars2.copy()
-            self.assertAlmostEqual(-0.7, pars3.DarkEnergy.w)
+            C2 = camb.get_results(pars2).get_cmb_power_spectra()
+            for f in ["lens_potential", "lensed_scalar"]:
+                self.assertTrue(np.allclose(C1[f][2:, 0], C2[f][2:, 0]))
+            pars3 = pars2.copy()
+            self.assertAlmostEqual(-0.7, pars3.DarkEnergy.w)
+
+        w0 = -0.7
+        wa = 0.2
+        pars.set_dark_energy(w=w0, wa=wa, dark_energy_model="ppf")
+        rho, w = camb.get_background(pars, no_thermo=True).get_dark_energy_rho_w(np.array([0.3, 0.5, 0.8, 1.0]))
+        a = np.array([0.3, 0.5, 0.8, 1.0])
+        expected_rho = a ** (-3 * (1 + w0 + wa)) * np.exp(-3 * wa * (1 - a))
+        expected_w = w0 + wa * (1 - a)
+        np.testing.assert_allclose(rho, expected_rho, rtol=1e-10)
+        np.testing.assert_allclose(w, expected_w, rtol=1e-12)
 
     def testInitialPower(self):
         pars = camb.CAMBparams()
diff --git a/fortran/DarkEnergyInterface.f90 b/fortran/DarkEnergyInterface.f90
index 0912cb77..97bf5893 100644
--- a/fortran/DarkEnergyInterface.f90
+++ b/fortran/DarkEnergyInterface.f90
@@ -220,6 +220,8 @@ function TDarkEnergyEqnOfState_grho_de(this, a) result(grho_de) !relative densit
     real(dl), intent(IN) :: a
 
     if(.not. this%use_tabulated_w) then
+        ! grho_de is proportional to a^4 rho_de/rho_de0, so for CPL w(a)=w0+wa(1-a)
+        ! it differs from rho_de/rho_de0 by an extra factor of a^4.
         grho_de = a ** (1._dl - 3. * this%w_lam - 3. * this%wa)
         if (this%wa/=0) grho_de=grho_de*exp(-3. * this%wa * (1._dl - a))
     else
diff --git a/fortran/halofit.f90 b/fortran/halofit.f90
index 1dc8fe81..df0737d3 100644
--- a/fortran/halofit.f90
+++ b/fortran/halofit.f90
@@ -516,6 +516,7 @@ end subroutine wint
 
     function omega_m(aa,om_m0,om_v0,wval,waval)
     real(dl) omega_m,omega_t,om_m0,om_v0,aa,wval,waval,Qa2
+    ! Qa2 is a^2*rho_de/rho_de0, matching the a^2 H^2 form used below.
     Qa2= aa**(-1.0-3.0*(wval+waval))*dexp(-3.0*(1-aa)*waval)
     omega_t=1.0+(om_m0+om_v0-1.0)/(1-om_m0-om_v0+om_v0*Qa2+om_m0/aa)
     omega_m=omega_t*om_m0/(om_m0+om_v0*aa*Qa2)
@@ -527,6 +528,7 @@ end function omega_m
 
     function omega_v(aa,om_m0,om_v0,wval,waval)
     real(dl) aa,omega_v,om_m0,om_v0,omega_t,wval,waval,Qa2
+    ! Qa2 is a^2*rho_de/rho_de0, matching the a^2 H^2 form used below.
     Qa2= aa**(-1.0-3.0*(wval+waval))*dexp(-3.0*(1-aa)*waval)
     omega_t=1.0+(om_m0+om_v0-1.0)/(1-om_m0-om_v0+om_v0*Qa2+om_m0/aa)
     omega_v=omega_t*om_v0*Qa2/(om_v0*Qa2+om_m0/aa)
@@ -2464,6 +2466,7 @@ FUNCTION X_de(z,cosm)
     REAL(dl) :: a
 
     a=1./(1.+z)
+    ! X_de is the physical rho_de/rho_de0 for CPL w(a)=w0+wa(1-a).
     X_de=(a**(-3*(1+cosm%w+cosm%wa)))*exp(-3*cosm%wa*(1-a))
 
     END FUNCTION X_de