move zeta into sneut struct

neutrinos/sneut5.H

@@ -291,9 +291,18 @@ struct sneutf_t {
     Real rmdz;
     Real rmi;
 
+    Real zeta;
+    Real zeta2;
+    Real zeta3;
+
+    Real zetadt;
+    Real zetada;
+    Real zetadz;
+
 };
 
 
+template <int do_derivatives>
 AMREX_GPU_HOST_DEVICE inline
 sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) {
 
@@ -339,6 +348,20 @@ sneutf_t get_sneut_factors(Real den, Real temp, Real abar, Real zbar) {
     sf.rmdz   = den * sf.abari;
     sf.rmi    = 1.0e0_rt / sf.rm;
 
+    Real a0 = sf.rm * 1.0e-9_rt;
+    Real a1 = std::pow(a0, nu_constants::oneth);
+
+    sf.zeta = a1 * sf.xlm1;
+    sf.zeta2 = sf.zeta * sf.zeta;
+    sf.zeta3 = sf.zeta2 * sf.zeta;
+
+    if constexpr (do_derivatives) {
+        Real a2 = nu_constants::oneth * a1*sf.rmi * sf.xlm1;
+        sf.zetadt = -a1 * sf.xlm2 * sf.xldt;
+        sf.zetada = a2 * sf.rmda;
+        sf.zetadz = a2 * sf.rmdz;
+    }
+
     return sf;
 
 }
@@ -806,11 +829,9 @@ void sneut5(const Real temp, const Real den,
 
     // pair production
     Real gl,gldt,
-         zeta,zetadt,zeta2,zeta3,
          xnum,xnumdt,xden,xdendt,
          fpair,fpairdt,qpair,qpairdt,
          fpairda,fpairdz,qpairda,qpairdz,
-         zetada,zetadz,
          xnumda,xnumdz,xdenda,xdendz;
 
     // plasma
@@ -859,21 +880,8 @@ void sneut5(const Real temp, const Real den,
 
     if (temp < 1.0e7_rt) return;
 
-    auto sf = get_sneut_factors(den, temp, abar, zbar);
-
-    a0     = sf.rm * 1.0e-9_rt;
-    a1     = std::pow(a0, nu_constants::oneth);
-    zeta   = a1 * sf.xlm1;
-
-    if constexpr (do_derivatives) {
-        a2     = nu_constants::oneth * a1*sf.rmi * sf.xlm1;
-        zetadt = -a1 * sf.xlm2 * sf.xldt;
-        zetada = a2 * sf.rmda;
-        zetadz = a2 * sf.rmdz;
-    }
+    auto sf = get_sneut_factors<do_derivatives>(den, temp, abar, zbar);
 
-    zeta2 = zeta * zeta;
-    zeta3 = zeta2 * zeta;
 
     // pair neutrino section
     // for reactions like e+ + e- => nu_e + nubar_e
@@ -886,27 +894,27 @@ void sneut5(const Real temp, const Real den,
 
     // equation 2.7
 
-    a1     = 6.002e19_rt + 2.084e20_rt*zeta + 1.872e21_rt*zeta2;
+    a1     = 6.002e19_rt + 2.084e20_rt*sf.zeta + 1.872e21_rt*sf.zeta2;
 
     if (sf.t9 < 10.0_rt) {
-       b1     = std::exp(-5.5924e0_rt*zeta);
+       b1     = std::exp(-5.5924e0_rt*sf.zeta);
        if constexpr (do_derivatives) {
            b2     = -b1*5.5924e0_rt;
        }
     } else {
-       b1     = std::exp(-4.9924e0_rt*zeta);
+       b1     = std::exp(-4.9924e0_rt*sf.zeta);
        if constexpr (do_derivatives) {
            b2     = -b1*4.9924e0_rt;
        }
     }
 
     xnum   = a1 * b1;
     if constexpr (do_derivatives) {
-        a2     = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*zeta;
+        a2     = 2.084e20_rt + 2.0e0_rt*1.872e21_rt*sf.zeta;
         c      = a2*b1 + a1*b2;
-        xnumdt = c*zetadt;
-        xnumda = c*zetada;
-        xnumdz = c*zetadz;
+        xnumdt = c*sf.zetadt;
+        xnumda = c*sf.zetada;
+        xnumdz = c*sf.zetadz;
     }
 
     if (sf.t9 < 10.0_rt) {
@@ -921,12 +929,12 @@ void sneut5(const Real temp, const Real den,
        }
     }
 
-    xden   = zeta3 + a1;
+    xden   = sf.zeta3 + a1;
     if constexpr (do_derivatives) {
-        b1   = 3.0e0_rt*zeta2;
-        xdendt = b1*zetadt + a2*sf.xldt;
-        xdenda = b1*zetada;
-        xdendz = b1*zetadz;
+        b1   = 3.0e0_rt*sf.zeta2;
+        xdendt = b1*sf.zetadt + a2*sf.xldt;
+        xdenda = b1*sf.zetada;
+        xdendz = b1*sf.zetadz;
     }
 
     a1      = 1.0e0_rt/xden;
@@ -1341,30 +1349,30 @@ void sneut5(const Real temp, const Real den,
 
 
     // equation 3.4
-    dum   = a0 + a1*zeta + a2*zeta2;
+    dum   = a0 + a1*sf.zeta + a2*sf.zeta2;
     if constexpr (do_derivatives) {
-        dumdt = f0 + f1*zeta + a1*zetadt + f2*zeta2 + 2.0e0_rt*a2*zeta*zetadt;
-        dumda = a1*zetada + 2.0e0_rt*a2*zeta*zetada;
-        dumdz = a1*zetadz + 2.0e0_rt*a2*zeta*zetadz;
+        dumdt = f0 + f1*sf.zeta + a1*sf.zetadt + f2*sf.zeta2 + 2.0e0_rt*a2*sf.zeta*sf.zetadt;
+        dumda = a1*sf.zetada + 2.0e0_rt*a2*sf.zeta*sf.zetada;
+        dumdz = a1*sf.zetadz + 2.0e0_rt*a2*sf.zeta*sf.zetadz;
     }
 
-    z      = std::exp(-cc*zeta);
+    z      = std::exp(-cc*sf.zeta);
 
     xnum   = dum*z;
     if constexpr (do_derivatives) {
-        xnumdt = dumdt*z - dum*z*cc*zetadt;
-        xnumda = dumda*z - dum*z*cc*zetada;
-        xnumdz = dumdz*z - dum*z*cc*zetadz;
+        xnumdt = dumdt*z - dum*z*cc*sf.zetadt;
+        xnumda = dumda*z - dum*z*cc*sf.zetada;
+        xnumdz = dumdz*z - dum*z*cc*sf.zetadz;
     }
 
-    xden   = zeta3 + 6.290e-3_rt*sf.xlm1 + 7.483e-3_rt*sf.xlm2 + 3.061e-4_rt*sf.xlm3;
+    xden   = sf.zeta3 + 6.290e-3_rt*sf.xlm1 + 7.483e-3_rt*sf.xlm2 + 3.061e-4_rt*sf.xlm3;
 
-    dum    = 3.0e0_rt*zeta2;
+    dum    = 3.0e0_rt*sf.zeta2;
     if constexpr (do_derivatives) {
-        xdendt = dum*zetadt - sf.xldt*(6.290e-3_rt*sf.xlm2
+        xdendt = dum*sf.zetadt - sf.xldt*(6.290e-3_rt*sf.xlm2
                                     + 2.0e0_rt*7.483e-3_rt*sf.xlm3 + 3.0e0_rt*3.061e-4_rt*sf.xlm4);
-        xdenda = dum*zetada;
-        xdendz = dum*zetadz;
+        xdenda = dum*sf.zetada;
+        xdendz = dum*sf.zetadz;
     }
     dum      = 1.0e0_rt/xden;
     fphot   = xnum*dum;