From 3d30e1cb600b3fb2f60159a5a3477ecb34a6b9ae Mon Sep 17 00:00:00 2001
From: Michael Zingale <michael.zingale@stonybrook.edu>
Date: Sun, 10 Sep 2023 12:53:09 -0400
Subject: [PATCH] add some comments on the Fermi integrals

---
 neutrinos/sneut5.H | 63 ++++++++++++++++++++++++++++++----------------
 1 file changed, 41 insertions(+), 22 deletions(-)

diff --git a/neutrinos/sneut5.H b/neutrinos/sneut5.H
index 6b066a2895..44b5d93cff 100644
--- a/neutrinos/sneut5.H
+++ b/neutrinos/sneut5.H
@@ -60,36 +60,55 @@ Real ifermi12(const Real f)
 
     if (f < 4.0e0_rt) {
 
-       rn  = f + a1(m1);
+        // build sum_{i=0, m1} a_i x**i.  This is the numerator
+        // in Eq. 4 of Antia.
+        //
+        // with our 1-based indexing, this expression is
+        // a[1] + a[2] * f + a[3] * f**2 + ... a[m1+1] * f**m1
+        //
+        // we do the sum starting with the largest term and working
+        // on a single power of f each iteration.
+        //
+        // in the starting rn here, the leading f is actually
+        // a1(m1+1) * f, but that element is 1
+        rn  = f + a1(m1);
+
+        for (int i = m1 - 1; i >= 1; --i) {
+            rn  = rn*f + a1(i);
+        }
+
+        // now we do the denominator in Eq. 4.  None of the coefficients
+        // are 1, so we loop over all
+        den = b1(k1+1);
+
+        for (int i = k1; i >= 1; --i) {
+            den = den*f + b1(i);
+        }
+
+        // Eq. 6 of Antia
+
+        ifermi12r = std::log(f * rn/den);
 
-       for (int i = m1 - 1; i >= 1; --i) {
-          rn  = rn*f + a1(i);
-       }
-
-       den = b1(k1+1);
-
-       for (int i = k1; i >= 1; --i) {
-          den = den*f + b1(i);
-       }
+    } else {
 
-       ifermi12r = std::log(f * rn/den);
+        ff = 1.0e0_rt/std::pow(f, 1.0e0_rt/(1.0e0_rt + an));
 
-    } else {
+        // this construction is the same as above, but using the
+        // second set of coefficients
 
-       ff = 1.0e0_rt/std::pow(f, 1.0e0_rt/(1.0e0_rt + an));
-       rn = ff + a2(m2);
+        rn = ff + a2(m2);
 
-       for (int i = m2 - 1; i >= 1; --i) {
-          rn = rn*ff + a2(i);
-       }
+        for (int i = m2 - 1; i >= 1; --i) {
+            rn = rn*ff + a2(i);
+        }
 
-       den = b2(k2+1);
+        den = b2(k2+1);
 
-       for (int i = k2; i >= 1; --i) {
-          den = den*ff + b2(i);
-       }
+        for (int i = k2; i >= 1; --i) {
+            den = den*ff + b2(i);
+        }
 
-       ifermi12r = rn/(den*ff);
+        ifermi12r = rn/(den*ff);
 
     }