diff --git a/docs/source/changelog.rst b/docs/source/changelog.rst
index 34fc942a7..74fc7d8a3 100644
--- a/docs/source/changelog.rst
+++ b/docs/source/changelog.rst
@@ -35,6 +35,12 @@ MESA no longer stops when reactions for which special rates are set are not in t
 Bug Fixes
 ---------
 
+The plasma neutrino cooling rate was using an hardcoded prefactor calculated with a Weinberg angle of 0.2319, while all other neutrino cooling processes used calculated prefactors taking the Weinberg angle as input, with default value 0.22290.
+Thus, modifying the value of the Weinberg angle resulted in changes to neutrino cooling processes except for the plasma neutrinos.
+This was the case in all previous MESA versions, and was found and fixed by user Garv Chauhan, see `gh-998 https://github.com/MESAHub/mesa/pull/998`_. Plasma neutrinos now use the same Weinberg angle as all other processes and changing its value will affect the corresponding cooling rate.
+The change in default Weinberg angle results in small numerical differences for stars where plasma neutrino cooling is significant.
+
+
 The parameter ``report_max_infall_inside_fe_core`` was ignored in versions r25.12.1 and r26.4.1 and always had it's default value. See `gh-981 https://github.com/MESAHub/mesa/pull/981`_.
 
 ``fe_core_infall_limit`` now obeys ``when_to_stop_rtol`` and ``when_to_stop_atol`` again (broken since r11532).
diff --git a/neu/private/mod_neu.f90 b/neu/private/mod_neu.f90
index 5108b3a6b..1a992e35e 100644
--- a/neu/private/mod_neu.f90
+++ b/neu/private/mod_neu.f90
@@ -35,6 +35,7 @@ module mod_neu
       real(dp), parameter :: cvp   = 1.0d0 - cv
       real(dp), parameter :: ca    = 0.5d0
       real(dp), parameter :: cap   = 1.0d0 - ca
+      real(dp), parameter :: tfac0 = cv*cv + (num_neu_fam-1.0d0) * (cvp*cvp)
       real(dp), parameter :: tfac1 = cv*cv + ca*ca + (num_neu_fam-1.0d0) * (cvp*cvp+cap*cap)
       real(dp), parameter :: tfac2 = cv*cv - ca*ca + (num_neu_fam-1.0d0) * (cvp*cvp - cap*cap)
       real(dp), parameter :: tfac3 = tfac2/tfac1
@@ -237,13 +238,13 @@ subroutine neutrinos(T, logT, Rho, logRho, abar, zbar, log10_Tlim,  &
          temp = T
       end if
 
-      if (T <= 0) then
+      if (temp <= 0) then
          info = -1
          return
       end if
 
       if (logT == arg_not_provided) then
-         logtemp = log10(T)
+         logtemp = log10(temp)
       else
          logtemp = logT
       end if
@@ -263,13 +264,13 @@ subroutine neutrinos(T, logT, Rho, logRho, abar, zbar, log10_Tlim,  &
          den = Rho
       end if
 
-      if (Rho <= 0) then
+      if (den <= 0) then
          info = -1
          return
       end if
 
       if (logRho == arg_not_provided) then
-         logden = log10(Rho)
+         logden = log10(den)
       else
          logden = logRho
       end if
@@ -1579,8 +1580,8 @@ subroutine plas_neu(splas,splasdt,splasdd,splasda,splasdz, input)
          splasdz = a2*splasdz + a3*gl2dz*a1
 
 
-         a2      = 0.93153d0 * 3.0d21 * input% xl9
-         a3      = 0.93153d0 * 3.0d21 * 9.0d0*input% xl8*input% xldt
+         a2      = tfac0 * 3.0d21 * input% xl9
+         a3      = tfac0 * 3.0d21 * 9.0d0*input% xl8*input% xldt
 
          a1      = splas
          splas   = a2*a1