From f5cd590612137e47de788eec28ef3576ec9921f7 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 09:57:26 -0500
Subject: [PATCH 01/16] modify FFI treatment to use dt / tau

---
 core/oscillations.c | 10 +++++++++-
 core/step.c         | 11 ++++-------
 2 files changed, 13 insertions(+), 8 deletions(-)

diff --git a/core/oscillations.c b/core/oscillations.c
index d0a9e40..dde8f06 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -187,11 +187,19 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
 #pragma omp parallel
   {
     struct of_photon *ph = photon_lists[omp_get_thread_num()];
+    double            X[NDIM], Kcov[NDIM], Kcon[NDIM];
     while (ph != NULL) {
       if (ph->type != TYPE_TRACER) {
         int ix1, ix2, icosth[LOCAL_NUM_BASES];
         get_local_angle_bins(ph, &ix1, &ix2, &icosth[0], &icosth[1]);
 
+        // find local FFI time scale nph already computed because of
+        // set_Rmunu which must be called each step
+        int i, j, k;
+        get_X_K_interp(ph, t, P, X, Kcov, Kcon);
+        Xtoijk(X, &i, &j, &k);
+        double tau = 1. / (T_unit * NUFERM * nph[i][j][k] + SMALL);
+
         int    b_osc = local_b_osc[ix1][ix2];
         int    imu   = icosth[b_osc];
         double A     = local_moments[b_osc][MOMENTS_A][ix1][ix2];
@@ -220,7 +228,7 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
           double p_survival =
               in_shallow ? peq : (1 - (1 - peq) * shallow / (deep + SMALL));
 #endif // FORCE_EQUIPARTITION
-          double p_osc = 1. - p_survival;
+          double p_osc = MY_MIN(1.0, dt / (tau + SMALL))*(1. - p_survival);
           if (get_rand() < p_osc) {
             // JMM:
             // Type order is NUE, NUEBAR, NUX, NUXBAR
diff --git a/core/step.c b/core/step.c
index 983c16a..3e32843 100644
--- a/core/step.c
+++ b/core/step.c
@@ -74,15 +74,12 @@ void step() {
   // check_nu_type("after bound"); // DEBUG
 #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS && \
     RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS
-  int oscillations_active = (dt_osc <= dt);
   if (mpi_io_proc()) {
-    printf("\t[Oscillations] Active? %d dt_osc = %.14e\n", oscillations_active,
-        dt_osc);
-  }
-  if (oscillations_active) { // TOOD(JMM): Some safety factor?
-    accumulate_local_angles();
-    oscillate(local_moments, Gnu);
+    printf("\t[Oscillations] tau_ffi = %.14e, dt/tau = %.14e\n",
+           dt_osc, dt / (dt_osc + SMALL));
   }
+  accumulate_local_angles();
+  oscillate(local_moments, Gnu);
   // check_nu_type("after oscillate"); // DEBUG
 #endif // OSCILLATIONS
 #endif

From d0b0eb5424454d4bd26cccbbd12304d669c277bd Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 10:51:36 -0500
Subject: [PATCH 02/16] clarify names for ffi vs CFI

---
 core/decs.h         | 2 +-
 core/oscillations.c | 2 +-
 core/step.c         | 4 ++--
 3 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index ef37661..db0dfd2 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -1041,7 +1041,7 @@ void accumulate_local_angles();
 void compute_local_gnu(grid_local_angles_type local_angles,
     grid_Gnu_type local_Ns, grid_Gnu_type local_wsqr, grid_Gnu_type gnu);
 void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments);
-void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu);
+void oscillate_ffi(grid_local_moment_type local_moments, grid_Gnu_type gnu);
 #endif // RAD_NUM_TYPES >= 4
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 #endif // RADIATION
diff --git a/core/oscillations.c b/core/oscillations.c
index dde8f06..43990d9 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -182,7 +182,7 @@ void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments) {
   }
 }
 
-void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
+void oscillate_ffi(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
   timer_start(TIMER_OSCILLATIONS);
 #pragma omp parallel
   {
diff --git a/core/step.c b/core/step.c
index 3e32843..0aac2e1 100644
--- a/core/step.c
+++ b/core/step.c
@@ -75,11 +75,11 @@ void step() {
 #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS && \
     RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS
   if (mpi_io_proc()) {
-    printf("\t[Oscillations] tau_ffi = %.14e, dt/tau = %.14e\n",
+    printf("\t[Oscillations] tau_ffi = %.14e, dt/tau_ffi = %.14e\n",
            dt_osc, dt / (dt_osc + SMALL));
   }
   accumulate_local_angles();
-  oscillate(local_moments, Gnu);
+  oscillate_ffi(local_moments, Gnu);
   // check_nu_type("after oscillate"); // DEBUG
 #endif // OSCILLATIONS
 #endif

From 93163f53724b0b6a37191cad349b2890cbac1a0a Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 11:00:33 -0500
Subject: [PATCH 03/16] rename more things to ffi to distinguish from cfi

---
 core/decs.h                          | 2 +-
 core/oscillations.c                  | 2 +-
 core/step.c                          | 4 ++--
 prob/{oscillations => ffi}/build.py  | 0
 prob/{oscillations => ffi}/problem.c | 0
 test/{oscillations.py => ffi.py}     | 6 +++---
 6 files changed, 7 insertions(+), 7 deletions(-)
 rename prob/{oscillations => ffi}/build.py (100%)
 rename prob/{oscillations => ffi}/problem.c (100%)
 rename test/{oscillations.py => ffi.py} (97%)

diff --git a/core/decs.h b/core/decs.h
index db0dfd2..c20cb1c 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -1033,7 +1033,7 @@ double alpha_nu_hdf(double nu, int type, const struct of_microphysics *m);
 
 // oscillations.c
 #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS
-double get_dt_oscillations();
+double get_dt_ffi();
 void get_local_angle_bins(
     struct of_photon *ph, int *pi, int *pj, int *pmu1, int *pmu2);
 void accumulate_local_angles();
diff --git a/core/oscillations.c b/core/oscillations.c
index 43990d9..d0ef7d6 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -10,7 +10,7 @@
 #if RADIATION == RADTYPE_NEUTRINOS
 #if LOCAL_ANGULAR_DISTRIBUTIONS
 
-double get_dt_oscillations() {
+double get_dt_ffi() {
   timer_start(TIMER_OSCILLATIONS);
   set_Rmunu(); // So we have Nsph and nph
   double nph_max = 0;
diff --git a/core/step.c b/core/step.c
index 0aac2e1..c317f54 100644
--- a/core/step.c
+++ b/core/step.c
@@ -31,7 +31,7 @@ void step() {
     RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS
   // not used for timestep control. Used to turn oscillations on or
   // off.
-  double dt_osc = get_dt_oscillations();
+  double dt_ffi = get_dt_ffi();
 #endif // oscillations
 #endif // radiation
   dtsave = dt;
@@ -76,7 +76,7 @@ void step() {
     RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS
   if (mpi_io_proc()) {
     printf("\t[Oscillations] tau_ffi = %.14e, dt/tau_ffi = %.14e\n",
-           dt_osc, dt / (dt_osc + SMALL));
+           dt_ffi, dt / (dt_ffi + SMALL));
   }
   accumulate_local_angles();
   oscillate_ffi(local_moments, Gnu);
diff --git a/prob/oscillations/build.py b/prob/ffi/build.py
similarity index 100%
rename from prob/oscillations/build.py
rename to prob/ffi/build.py
diff --git a/prob/oscillations/problem.c b/prob/ffi/problem.c
similarity index 100%
rename from prob/oscillations/problem.c
rename to prob/ffi/problem.c
diff --git a/test/oscillations.py b/test/ffi.py
similarity index 97%
rename from test/oscillations.py
rename to test/ffi.py
index 2da64b3..2df2766 100644
--- a/test/oscillations.py
+++ b/test/ffi.py
@@ -17,7 +17,7 @@
 
 TMP_DIR = 'TMP'
 util.safe_remove(TMP_DIR)
-PROBLEM = 'oscillations'
+PROBLEM = 'ffi'
 AUTO = '-auto' in sys.argv
 MPI = '-mpi' in sys.argv
 gam = 1.4
@@ -175,8 +175,8 @@
 
 plt.tight_layout()
 
-plt.savefig('oscillations_1zone.png', bbox_inches='tight')
-plt.savefig('oscillations_1zone.pdf', bbox_inches='tight')
+plt.savefig('ffi_1zone.png', bbox_inches='tight')
+plt.savefig('ffi_1zone.pdf', bbox_inches='tight')
 
 # clean up
 util.safe_remove(TMP_DIR)

From 06a9a55fcda6197e5842154cd386bf0610651fb7 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 11:00:46 -0500
Subject: [PATCH 04/16] add 1zone ffi test to automated suite

---
 test/test_all.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/test_all.py b/test/test_all.py
index 5784054..aa51e2a 100644
--- a/test/test_all.py
+++ b/test/test_all.py
@@ -45,7 +45,7 @@
 M_NU_TESTS = ['yedecay.py', 'yedecay.py -antinu',
               'yedecay.py -mpi', # 'yedecay.py -mpi -antinu',
               'multiscatt.py',
-              'tracers1d.py']
+              'tracers1d.py', 'ffi.py']
 M_LIGHT_TESTS = ['binning.py', 'brem.py',
                  'thermalization.py', 'thermalization_mpi.py',
                  'comptonization.py']

From ec66cdc9e37b2ce323e363ff172db42853f0cb37 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 22:29:42 -0500
Subject: [PATCH 05/16] in progress

---
 core/constants.h |  5 +++++
 core/decs.h      | 21 +++++++++++++++++++--
 core/defs.h      |  8 ++++++++
 3 files changed, 32 insertions(+), 2 deletions(-)

diff --git a/core/constants.h b/core/constants.h
index d4fe9b0..5d98507 100644
--- a/core/constants.h
+++ b/core/constants.h
@@ -32,6 +32,11 @@
 #define ROOT2 (1.4142135623730951)
 #define NUFERM ((ROOT2*GFERM) / HBAR)
 
+// Cos(2 theta21) where theta_m is the vacuum mixing angle theta21
+#define C2TH21 (0.386)
+// square of mass gap between m1 and m2 for neutrino masses
+#define DM21SQR_EV (7.5e-5) // eV
+
 // Unit conversions
 #define EV (1.60217653e-12)   // Electron-volt
 #define MEV (1.0e6 * EV)      // Mega-Electron-Volt
diff --git a/core/decs.h b/core/decs.h
index c20cb1c..c94bedc 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -171,17 +171,20 @@
 #define ANTINU_HEAVY (3)
 #if MULTISCATT_TEST
 #define RAD_SCATT_TYPES (3)
-#else
+#else // NOT MULTISCATT TEST
 #define RAD_SCATT_TYPES (4) // TODO: Should be 5, including electrons
 #define RSCATT_TYPE_P (0)
 #define RSCATT_TYPE_N (1)
 #define RSCATT_TYPE_A (2)
 #define RSCATT_TYPE_ALPHA (3)
 #define RSCATT_TYPE_E (4) // TOOD: implement me
-#endif
+#endif // NOT MULTISCATT TEST
 #define NRADCOMP (2)
 #define RADG_YE (4)
 #define RADG_YE_EM (5)
+#define CFI_INACTIVE (0)      // For collisional flavor instability
+#define CFI_TYPE_GAPPED (1)   // Omega_- mode
+#define CFI_TYPE_GAPLESS (2)  // Omega_+ mode
 #elif RADIATION == RADTYPE_LIGHT
 #define RAD_SCATT_TYPES (1)
 #define NRADCOMP (0)
@@ -403,6 +406,13 @@ extern grid_local_count_type local_osc_count;
 #endif // #if RAD_NUM_TYPES >= 4
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
+#if DO_CFI
+// Which CFI mode is active, if any
+extern grid_int_type cfi_active_mode;
+// change towards asymptotic state and time scale
+extern grid_double_type cfi_delta_asymp, cfi_tau_asymp;
+#endif // DO_CFI
+
 #endif // RADIATION
 
 // Default initialization is 0, which in this case is
@@ -1033,6 +1043,7 @@ double alpha_nu_hdf(double nu, int type, const struct of_microphysics *m);
 
 // oscillations.c
 #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS
+// FFI
 double get_dt_ffi();
 void get_local_angle_bins(
     struct of_photon *ph, int *pi, int *pj, int *pmu1, int *pmu2);
@@ -1044,6 +1055,12 @@ void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments);
 void oscillate_ffi(grid_local_moment_type local_moments, grid_Gnu_type gnu);
 #endif // RAD_NUM_TYPES >= 4
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
+
+// CFI
+#if DO_CFI
+void compute_cfi_active_mode(grid_int_type cfi_active_mode);
+#endif // DO_CFI
+
 #endif // RADIATION
 
 // passive.c
diff --git a/core/defs.h b/core/defs.h
index 6dd006b..fc0c214 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -67,6 +67,14 @@ grid_local_basis_idx_type local_b_osc;
 grid_local_count_type local_osc_count;
 #endif //  RAD_NUM_TYPES
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
+
+#if DO_CFI
+// Which CFI mode is active, if any
+grid_int_type cfi_active_mode;
+// change towards asymptotic state and time scale
+grid_double_type cfi_delta_asymp, cfi_tau_asymp;
+#endif // DO_CFI
+
 #endif // RADIATION
 
 #if ELECTRONS

From 56db1368be3ed4c1d39d76acb10f656937dc41d6 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 22:36:45 -0500
Subject: [PATCH 06/16] argh

---
 core/decs.h         |  3 ++-
 core/oscillations.c | 19 ++++++++++---------
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index c94bedc..e6fe608 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -1051,7 +1051,8 @@ void accumulate_local_angles();
 #if RAD_NUM_TYPES >= 4
 void compute_local_gnu(grid_local_angles_type local_angles,
     grid_Gnu_type local_Ns, grid_Gnu_type local_wsqr, grid_Gnu_type gnu);
-void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments);
+void compute_local_moments(grid_local_angles_type f,
+                           grid_Gnu_type gnu, grid_local_moment_type moments);
 void oscillate_ffi(grid_local_moment_type local_moments, grid_Gnu_type gnu);
 #endif // RAD_NUM_TYPES >= 4
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
diff --git a/core/oscillations.c b/core/oscillations.c
index d0ef7d6..d22ca44 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -51,7 +51,7 @@ void accumulate_local_angles() {
 #pragma omp atomic
           local_Ns[b][ix1][ix2][icosth[b]] += 1.;
 #pragma omp atomic
-          local_wsqr[b][ix1][ix2][icosth[b]] += (ph->w)*(ph->w);
+          local_wsqr[b][ix1][ix2][icosth[b]] += (ph->w) * (ph->w);
         }
       }
       ph = ph->next;
@@ -106,7 +106,7 @@ void get_local_angle_bins(
   knorm = 1. / (fabs(knorm) + SMALL);
   costh1 *= knorm;
   costh2 *= knorm;
- 
+
   *pi = MY_MAX(
       0, MY_MIN(LOCAL_ANGLES_NX1 - 1, (X[1] - startx_rad[1]) / local_dx1_rad));
   *pj = MY_MAX(
@@ -118,8 +118,8 @@ void get_local_angle_bins(
 }
 
 #if RAD_NUM_TYPES >= 4
-void compute_local_gnu(grid_local_angles_type f, grid_Gnu_type local_Ns,
-                       grid_Gnu_type local_wsqr, grid_Gnu_type gnu) {
+void        compute_local_gnu(grid_local_angles_type f, grid_Gnu_type local_Ns,
+           grid_Gnu_type local_wsqr, grid_Gnu_type gnu) {
 #pragma omp parallel for collapse(4)
   for (int b = 0; b < LOCAL_NUM_BASES; ++b) {
     LOCALXMULOOP {
@@ -151,7 +151,8 @@ void compute_local_gnu(grid_local_angles_type f, grid_Gnu_type local_Ns,
 
 // JMM: We can also compute, e.g., the average bin momentum if we need
 // to, e.g., compute higher moment integrands
-void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments) {
+void compute_local_moments(grid_local_angles_type f, grid_Gnu_type gnu,
+    grid_local_moment_type moments) {
   // We are reducing over mu, but if we just parallel loop over b,i,j,
   // there is no danger of index collisions.
   LOCALMULOOP {
@@ -208,7 +209,7 @@ void oscillate_ffi(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
 
         // gnu == 0 when we activated stddev trigger. Don't oscillate.
         if (((G != 0) || FORCE_EQUIPARTITION) && (A != 0) && (B != 0)) {
-        // if ((A != 0) && (B != 0)) {
+          // if ((A != 0) && (B != 0)) {
           // If A == B then which region we treat as shallow is
           // unimportant. Psurvive = 1/3 for both regions.
           int    A_is_shallow = A < B;
@@ -221,14 +222,14 @@ void oscillate_ffi(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
 
           int in_shallow = (A_is_shallow && g_in_A) || (B_is_shallow && g_in_B);
 
-          double peq = nu_is_heavy(ph->type) ? (2./3.) : (1./3.);
+          double peq = nu_is_heavy(ph->type) ? (2. / 3.) : (1. / 3.);
 #if FORCE_EQUIPARTITION
           double p_survival = peq;
 #else
           double p_survival =
               in_shallow ? peq : (1 - (1 - peq) * shallow / (deep + SMALL));
 #endif // FORCE_EQUIPARTITION
-          double p_osc = MY_MIN(1.0, dt / (tau + SMALL))*(1. - p_survival);
+          double p_osc = MY_MIN(1.0, dt / (tau + SMALL)) * (1. - p_survival);
           if (get_rand() < p_osc) {
             // JMM:
             // Type order is NUE, NUEBAR, NUX, NUXBAR
@@ -237,7 +238,7 @@ void oscillate_ffi(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
             ph->type = (ph->type + (RAD_NUM_TYPES / 2)) % RAD_NUM_TYPES;
             ph->osc_count += 1;
 
-            #pragma omp atomic
+#pragma omp atomic
             local_osc_count[ix1][ix2] += ph->w;
           }
         }

From 0b75c9e840f182cecbd25e1014ceb4b221e458e0 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 22:57:27 -0500
Subject: [PATCH 07/16] update with Gails formula

---
 core/decs.h         |  3 ++-
 core/oscillations.c | 34 +++++++++++++++++++++++++---------
 core/step.c         |  2 +-
 3 files changed, 28 insertions(+), 11 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index ef37661..5863724 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -1041,7 +1041,8 @@ void accumulate_local_angles();
 void compute_local_gnu(grid_local_angles_type local_angles,
     grid_Gnu_type local_Ns, grid_Gnu_type local_wsqr, grid_Gnu_type gnu);
 void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments);
-void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu);
+void oscillate(grid_local_angles_type f,
+               grid_local_moment_type local_moments, grid_Gnu_type gnu);
 #endif // RAD_NUM_TYPES >= 4
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 #endif // RADIATION
diff --git a/core/oscillations.c b/core/oscillations.c
index dde8f06..1595de9 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -182,7 +182,8 @@ void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments) {
   }
 }
 
-void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
+void oscillate(grid_local_angles_type f,
+               grid_local_moment_type local_moments, grid_Gnu_type gnu) {
   timer_start(TIMER_OSCILLATIONS);
 #pragma omp parallel
   {
@@ -193,13 +194,6 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
         int ix1, ix2, icosth[LOCAL_NUM_BASES];
         get_local_angle_bins(ph, &ix1, &ix2, &icosth[0], &icosth[1]);
 
-        // find local FFI time scale nph already computed because of
-        // set_Rmunu which must be called each step
-        int i, j, k;
-        get_X_K_interp(ph, t, P, X, Kcov, Kcon);
-        Xtoijk(X, &i, &j, &k);
-        double tau = 1. / (T_unit * NUFERM * nph[i][j][k] + SMALL);
-
         int    b_osc = local_b_osc[ix1][ix2];
         int    imu   = icosth[b_osc];
         double A     = local_moments[b_osc][MOMENTS_A][ix1][ix2];
@@ -208,7 +202,29 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
 
         // gnu == 0 when we activated stddev trigger. Don't oscillate.
         if (((G != 0) || FORCE_EQUIPARTITION) && (A != 0) && (B != 0)) {
-        // if ((A != 0) && (B != 0)) {
+          // finite volumes cell
+          int i, j, k;
+          get_X_K_interp(ph, t, P, X, Kcov, Kcon);
+          Xtoijk(X, &i, &j, &k);
+
+          // compute crossing depth
+          double n_tot_supercell = 0;
+          TYPELOOP {
+            LOCALMULOOP {
+              // total neutrinos of all flavors
+              n_tot_supercell += fabs(f[0][ix1][ix2][itp][imu]);
+            }
+          }
+          // crossing depth in units of number/cm^3
+          // sqrt(AB)/n_tot_supercell is normalized unitless crossing
+          // depth but units should be in number/cm^3.
+          double crossing_depth = sqrt(fabs(A * B))*nph[i][j][k] / (n_tot_supercell + SMALL);
+
+          // find local FFI time scale nph already computed because of
+          // Rmunu This crossing depth formula is slightly more
+          // sophisticated than just 1 / G_F n
+          double tau = 1. / (T_unit * NUFERM * crossing_depth  + SMALL);
+
           // If A == B then which region we treat as shallow is
           // unimportant. Psurvive = 1/3 for both regions.
           int    A_is_shallow = A < B;
diff --git a/core/step.c b/core/step.c
index 3e32843..5ed6d12 100644
--- a/core/step.c
+++ b/core/step.c
@@ -79,7 +79,7 @@ void step() {
            dt_osc, dt / (dt_osc + SMALL));
   }
   accumulate_local_angles();
-  oscillate(local_moments, Gnu);
+  oscillate(local_angles, local_moments, Gnu);
   // check_nu_type("after oscillate"); // DEBUG
 #endif // OSCILLATIONS
 #endif

From fa0f575613c4cc4de0d703fc859ea0eecf86a7ef Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 23:03:07 -0500
Subject: [PATCH 08/16] oops bad merge fixed now

---
 core/oscillations.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/core/oscillations.c b/core/oscillations.c
index e56b9b6..bd7cac3 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -151,8 +151,8 @@ void        compute_local_gnu(grid_local_angles_type f, grid_Gnu_type local_Ns,
 
 // JMM: We can also compute, e.g., the average bin momentum if we need
 // to, e.g., compute higher moment integrands
-void compute_local_moments(grid_local_angles_type f, grid_Gnu_type gnu,
-    grid_local_moment_type moments) {
+void compute_local_moments(grid_Gnu_type gnu,
+                           grid_local_moment_type moments) {
   // We are reducing over mu, but if we just parallel loop over b,i,j,
   // there is no danger of index collisions.
   LOCALMULOOP {

From cddc6219ed4c4c5ac47acb81d7f67711a1b19261 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 23:06:44 -0500
Subject: [PATCH 09/16] cleanup

---
 core/decs.h | 4 ++--
 core/defs.h | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index c1db1cc..c8706ac 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -409,8 +409,8 @@ extern grid_local_count_type local_osc_count;
 #if DO_CFI
 // Which CFI mode is active, if any
 extern grid_int_type cfi_active_mode;
-// change towards asymptotic state and time scale
-extern grid_double_type cfi_delta_asymp, cfi_tau_asymp;
+// time scale for asymptotic state
+extern grid_double_type cfi_tau_asymp;
 #endif // DO_CFI
 
 #endif // RADIATION
diff --git a/core/defs.h b/core/defs.h
index fc0c214..171e6d2 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -71,8 +71,8 @@ grid_local_count_type local_osc_count;
 #if DO_CFI
 // Which CFI mode is active, if any
 grid_int_type cfi_active_mode;
-// change towards asymptotic state and time scale
-grid_double_type cfi_delta_asymp, cfi_tau_asymp;
+// time scale for asymptotic state
+grid_double_type cfi_tau_asymp;
 #endif // DO_CFI
 
 #endif // RADIATION

From 5e6e2a2d7d67b5e12bbb7ebcaa86e324f78aeeb7 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 23:07:37 -0500
Subject: [PATCH 10/16] cleanup

---
 core/oscillations.c | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/core/oscillations.c b/core/oscillations.c
index bd7cac3..26bddcc 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -118,8 +118,8 @@ void get_local_angle_bins(
 }
 
 #if RAD_NUM_TYPES >= 4
-void        compute_local_gnu(grid_local_angles_type f, grid_Gnu_type local_Ns,
-           grid_Gnu_type local_wsqr, grid_Gnu_type gnu) {
+void compute_local_gnu(grid_local_angles_type f, grid_Gnu_type local_Ns,
+                       grid_Gnu_type local_wsqr, grid_Gnu_type gnu) {
 #pragma omp parallel for collapse(4)
   for (int b = 0; b < LOCAL_NUM_BASES; ++b) {
     LOCALXMULOOP {

From 17c560729ce25896e205e68dca8745b562c92db4 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 23:34:31 -0500
Subject: [PATCH 11/16] rename oscillations -> ffi

---
 core/decs.h             | 59 ++++++++++++++++++++++-------------------
 core/defs.h             | 14 +++++++---
 core/diag.c             | 16 +++++------
 core/io.c               | 10 +++----
 core/step.c             |  4 +--
 core/timing.c           |  2 +-
 prob/ffi/build.py       |  2 +-
 prob/ffi/problem.c      |  6 ++---
 prob/torus_cbc/build.py |  8 +++---
 script/config.py        | 10 ++++---
 10 files changed, 72 insertions(+), 59 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index c8706ac..1982249 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -100,12 +100,6 @@
 // in scattering stability criterion
 #define SCATT_BIAS_SAFETY (0.5) // (0.9/RAD_SCATT_TYPES)
 
-//#define NUMIN     (1.e10)
-//#define NUMAX     (1.e25)
-//#define NU_BINS_SPEC    (200)
-//#define NTH (8)
-//#define NPHI (8)
-
 // Whether to move polar axis slightly off of coordinate singularity
 #define COORDSINGFIX 1
 #define SINGSMALL (1.E-20)
@@ -182,9 +176,9 @@
 #define NRADCOMP (2)
 #define RADG_YE (4)
 #define RADG_YE_EM (5)
-#define CFI_INACTIVE (0)      // For collisional flavor instability
-#define CFI_TYPE_GAPPED (1)   // Omega_- mode
-#define CFI_TYPE_GAPLESS (2)  // Omega_+ mode
+#define CFI_INACTIVE (0)     // For collisional flavor instability
+#define CFI_TYPE_GAPPED (1)  // Omega_- mode
+#define CFI_TYPE_GAPLESS (2) // Omega_+ mode
 #elif RADIATION == RADTYPE_LIGHT
 #define RAD_SCATT_TYPES (1)
 #define NRADCOMP (0)
@@ -345,9 +339,9 @@ extern grid_int_type    Nsph;
 extern grid_double_type nph;
 #if RZ_HISTOGRAMS
 extern rz_hist_type rz_r_orig_hist, rz_z_orig_hist;
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
 extern rz_hist_type osc_rz_r_orig_hist, osc_rz_z_orig_hist;
-#endif // NEUTRINO_OSCILLATIONS
+#endif // NEUTRINO_OSCILLATIONS_FFI
 #endif // RZ_HISTOGRAMS
 
 extern struct of_photon **photon_lists;
@@ -399,18 +393,29 @@ typedef double grid_local_moment_type[LOCAL_NUM_BASES][LOCAL_NUM_MOMENTS]
                                      [LOCAL_ANGLES_NX1][LOCAL_ANGLES_NX2];
 typedef int grid_local_basis_idx_type[LOCAL_ANGLES_NX1][LOCAL_ANGLES_NX2];
 typedef double grid_local_count_type[LOCAL_ANGLES_NX1][LOCAL_ANGLES_NX2];
-extern grid_Gnu_type          Gnu, local_Ns, local_wsqr;
-extern grid_local_moment_type local_moments;
+extern grid_Gnu_type             Gnu, local_Ns, local_wsqr;
+extern grid_local_moment_type    local_moments;
 extern grid_local_basis_idx_type local_b_osc;
-extern grid_local_count_type local_osc_count;
+extern grid_local_count_type     local_osc_count;
 #endif // #if RAD_NUM_TYPES >= 4
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
 #if DO_CFI
+typedef double grid_double_type[N1 + 2 * NG][N2 + 2 * NG];
+typedef int grid_symm_int_type[N1 + 2 * NG][N2 + 2 * NG];
+typedef double grid_symm_radtype_type[N1 + 2 * NG][N2 + 2 * NG][RAD_NUM_TYPES];
+typedef double grid_CFI_Gamma_type[N1 + 2 * NG][N2 + 2 * NG][2];
+
+// number of neutrinos per flavor
+extern grid_symm_radtype_type nph_flavor;
+// individual distribution function-weighted opacity averages
+extern grid_symm_radtype_type kappa_avg;
+extern grid_CFI_Gamma_type cfi_Gamma;
+
 // Which CFI mode is active, if any
-extern grid_int_type cfi_active_mode;
+extern grid_symm_int_type cfi_active_mode;
 // time scale for asymptotic state
-extern grid_double_type cfi_tau_asymp;
+extern grid_symm_double_type cfi_tau_asymp;
 #endif // DO_CFI
 
 #endif // RADIATION
@@ -592,11 +597,11 @@ struct of_photon {
   // radiation type. For neutrinos, flavor. Always active.
   // Not always important.
   // TODO: make sure to always set type when it's needed.
-  int               type;
-  int               nscatt;
-  int               origin[NDIM];
-  double            t0;
-  int               is_tracked;
+  int    type;
+  int    nscatt;
+  int    origin[NDIM];
+  double t0;
+  int    is_tracked;
   // Only relevant for neutrino oscillations
   int               osc_count;
   struct of_photon *next;
@@ -721,8 +726,8 @@ extern int global_stop[NDIM];
 #define JRADLOOP for (int n = 0; n < MAXNSCATT + 2; n++)
 #define NULOOP for (int inu = 0; inu < NU_BINS + 1; inu++)
 
-#define LOCALXLOOP                             \
-  for (int i = 0; i < LOCAL_ANGLES_NX1; ++i)   \
+#define LOCALXLOOP                           \
+  for (int i = 0; i < LOCAL_ANGLES_NX1; ++i) \
     for (int j = 0; j < LOCAL_ANGLES_NX2; ++j)
 #define LOCALMULOOP for (int imu = 0; imu < LOCAL_ANGLES_NMU; ++imu)
 #define LOCALXMULOOP LOCALXLOOP LOCALMULOOP
@@ -1046,11 +1051,11 @@ double alpha_nu_hdf(double nu, int type, const struct of_microphysics *m);
 // FFI
 double get_dt_ffi();
 void get_local_angle_bins(
-    struct of_photon *ph, int *pi, int *pj, int *pmu1, int *pmu2);
+  struct of_photon *ph, int *pi, int *pj, int *pmu1, int *pmu2);
 void accumulate_local_angles();
 #if RAD_NUM_TYPES >= 4
 void compute_local_gnu(grid_local_angles_type local_angles,
-    grid_Gnu_type local_Ns, grid_Gnu_type local_wsqr, grid_Gnu_type gnu);
+  grid_Gnu_type local_Ns, grid_Gnu_type local_wsqr, grid_Gnu_type gnu);
 void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments);
 void oscillate_ffi(grid_local_angles_type f,
                    grid_local_moment_type local_moments, grid_Gnu_type gnu);
@@ -1065,12 +1070,10 @@ void compute_cfi_active_mode(grid_int_type cfi_active_mode);
 #endif // RADIATION
 
 // passive.c
-//#if NVAR_PASSIVE > 0
 void fixup_passive(
     int i, int j, int k, double pv[NVAR], double pv_prefloor[NVAR]);
 void init_passives();
 void name_passives();
-//#endif
 
 // phys.c
 void   primtoflux(double *pr, struct of_state *q, int dir, int magnetic,
@@ -1253,7 +1256,7 @@ double find_median(double *array, int size);
 double interp_1d(double x, const double xmin, const double xmax, const int imin,
     const int imax, const double *restrict tab_x, const double *restrict tab_y);
 int    find_index(double value, const double *array, int size);
-void * safe_malloc(size_t size);
+void  *safe_malloc(size_t size);
 void   safe_system(const char *command);
 void   safe_fscanf(FILE *stream, const char *format, ...);
 int    is_practically_nan(double v);
diff --git a/core/defs.h b/core/defs.h
index 171e6d2..8e68f94 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -29,9 +29,9 @@ grid_int_type      Nsph;
 grid_double_type   nph;
 #if RZ_HISTOGRAMS
 rz_hist_type rz_r_orig_hist, rz_z_orig_hist;
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
 rz_hist_type osc_rz_r_orig_hist, osc_rz_z_orig_hist;
-#endif // NEUTRINO_OSCILLATIONS
+#endif // NEUTRINO_OSCILLATIONS_FFI
 #endif // RZ_HISTOGRAMS
 
 struct of_photon **photon_lists;
@@ -69,10 +69,16 @@ grid_local_count_type local_osc_count;
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
 #if DO_CFI
+// number of neutrinos per flavor
+grid_symm_radtype_type nph_flavor;
+// individual distribution function-weighted opacity averages
+grid_symm_radtype_type kappa_avg;
+grid_CFI_Gamma_type cfi_Gamma;
+
 // Which CFI mode is active, if any
-grid_int_type cfi_active_mode;
+grid_symm_int_type cfi_active_mode;
 // time scale for asymptotic state
-grid_double_type cfi_tau_asymp;
+grid_symm_double_type cfi_tau_asymp;
 #endif // DO_CFI
 
 #endif // RADIATION
diff --git a/core/diag.c b/core/diag.c
index 10c4da6..416ef2e 100644
--- a/core/diag.c
+++ b/core/diag.c
@@ -38,14 +38,14 @@ void reset_dump_variables() {
 #if RZ_HISTOGRAMS
   memset(rz_r_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double));
   memset(rz_z_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double));
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
   memset(osc_rz_r_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double));
   memset(osc_rz_z_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double));
-#endif // NEUTRINO_OSCILLATIONS
+#endif // NEUTRINO_OSCILLATIONS_FFI
 #endif // RZ_HISTOGRAMS
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
   memset(local_osc_count, 0, LOCAL_ANGLES_NX1*LOCAL_ANGLES_NX2*sizeof(double));
-#endif // NEUTRINO_OSCILLATIONS
+#endif // NEUTRINO_OSCILLATIONS_FFI
 #endif // RADIATION
 }
 
@@ -310,7 +310,7 @@ void diag(int call_code) {
 #if ELECTRONS
       fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_ELECTRON));
 #endif
-#if NEUTRINO_OSCILLATIONS || LOCAL_ANGULAR_DISTRIBUTIONS
+#if NEUTRINO_OSCILLATIONS_FFI || LOCAL_ANGULAR_DISTRIBUTIONS
       fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_OSCILLATIONS));
 #endif
       fprintf(ener_file, "\n");
@@ -596,21 +596,21 @@ void generate_rz_histograms() {
 #pragma omp atomic
       rz_z_orig_hist[iz] += ph->w;
 
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
       if (ph->osc_count) {
 #pragma omp atomic
         osc_rz_r_orig_hist[ir] += ph->w;
 #pragma omp atomic
         osc_rz_z_orig_hist[iz] += ph->w;
       }
-#endif // NEUTRINO_OSCILLATIONS
+#endif // NEUTRINO_OSCILLATIONS_FFI
 
       ph = ph->next;
     }
   }
   mpi_dbl_allreduce_array((double *)rz_r_orig_hist, RZ_HISTOGRAMS_N);
   mpi_dbl_allreduce_array((double *)rz_z_orig_hist, RZ_HISTOGRAMS_N);
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
   mpi_dbl_allreduce_array((double *)osc_rz_r_orig_hist, RZ_HISTOGRAMS_N);
   mpi_dbl_allreduce_array((double *)osc_rz_z_orig_hist, RZ_HISTOGRAMS_N);
 #endif
diff --git a/core/io.c b/core/io.c
index b6aa310..9e03cad 100644
--- a/core/io.c
+++ b/core/io.c
@@ -878,8 +878,8 @@ void dump() {
   WRITE_HDR(numin, TYPE_DBL);
   WRITE_HDR(numax, TYPE_DBL);
 
-  int neutrino_oscillations = NEUTRINO_OSCILLATIONS;
-  WRITE_HDR(neutrino_oscillations, TYPE_INT);
+  int neutrino_oscillations_ffi = NEUTRINO_OSCILLATIONS_FFI;
+  WRITE_HDR(neutrino_oscillations_ffi, TYPE_INT);
   int force_equipartition = FORCE_EQUIPARTITION;
   WRITE_HDR(force_equipartition, TYPE_INT);
   int local_angular_distributions = LOCAL_ANGULAR_DISTRIBUTIONS;
@@ -1145,7 +1145,7 @@ void dump() {
 #undef RANK
     }
 
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
     {
       mpi_dbl_allreduce_array(
           (double *)local_osc_count, LOCAL_ANGLES_NX1 * LOCAL_ANGLES_NX2);
@@ -1169,7 +1169,7 @@ void dump() {
           TYPE_DBL);
 #undef RANK
     }
-#endif // NEUTRINO_OSCILLATIONS
+#endif // NEUTRINO_OSCILLATIONS_FFI
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
 #if RZ_HISTOGRAMS
@@ -1186,7 +1186,7 @@ void dump() {
     }
     WRITE_ARRAY(rz_r_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL);
     WRITE_ARRAY(rz_z_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL);
-#if NEUTRINO_OSCILLATIONS
+#if NEUTRINO_OSCILLATIONS_FFI
     WRITE_ARRAY(osc_rz_r_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL);
     WRITE_ARRAY(osc_rz_z_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL);
 #endif
diff --git a/core/step.c b/core/step.c
index aea99aa..25cf4a8 100644
--- a/core/step.c
+++ b/core/step.c
@@ -28,7 +28,7 @@ void step() {
 #if RADIATION
   double dt_cool;
 #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS && \
-    RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS
+    RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS_FFI
   // not used for timestep control. Used to turn oscillations on or
   // off.
   double dt_ffi = get_dt_ffi();
@@ -73,7 +73,7 @@ void step() {
   bound_superphotons(Ph, t, dt);
   // check_nu_type("after bound"); // DEBUG
 #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS && \
-    RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS
+    RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS_FFI
   if (mpi_io_proc()) {
     printf("\t[Oscillations] tau_ffi = %.14e, dt/tau_ffi = %.14e\n",
            dt_ffi, dt / (dt_ffi + SMALL));
diff --git a/core/timing.c b/core/timing.c
index 6a5b3c3..00296a7 100644
--- a/core/timing.c
+++ b/core/timing.c
@@ -91,7 +91,7 @@ void report_performance() {
     fprintf(stdout, "   INTERACT: %8.4g s (%.4g %%)\n",
         times[TIMER_INTERACT] / dnstep,
         100. * times[TIMER_INTERACT] / times[TIMER_ALL]);
-#if NEUTRINO_OSCILLATIONS || LOCAL_ANGULAR_DISTRIBUTIONS
+#if NEUTRINO_OSCILLATIONS_FFI || LOCAL_ANGULAR_DISTRIBUTIONS
     fprintf(stdout, "   OSCILL:   %8.4g s (%.4g %%)\n",
         times[TIMER_OSCILLATIONS] / dnstep,
         100. * times[TIMER_OSCILLATIONS] / times[TIMER_ALL]);
diff --git a/prob/ffi/build.py b/prob/ffi/build.py
index ac6367b..a4c23e2 100644
--- a/prob/ffi/build.py
+++ b/prob/ffi/build.py
@@ -129,7 +129,7 @@
 bhl.config.set_cparm('LOCAL_ANGLES_NX1', NTOT)
 bhl.config.set_cparm('LOCAL_ANGLES_NX2', NTOT)
 bhl.config.set_cparm('RAD_NUM_TYPES', 4)
-bhl.config.set_cparm('NEUTRINO_OSCILLATIONS', True)
+bhl.config.set_cparm('NEUTRINO_OSCILLATIONS_FFI', True)
 bhl.config.set_cparm('FORCE_EQUIPARTITION', False)
 bhl.config.set_cparm("RZ_HISTOGRAMS", True)
 
diff --git a/prob/ffi/problem.c b/prob/ffi/problem.c
index b0d8049..198ded1 100644
--- a/prob/ffi/problem.c
+++ b/prob/ffi/problem.c
@@ -54,14 +54,14 @@ double distmu(int type, double x) {
 // Initialize dynamical variables
 void init_prob() {
   if ((RADIATION != RADTYPE_NEUTRINOS) || !LOCAL_ANGULAR_DISTRIBUTIONS ||
-      !NEUTRINO_OSCILLATIONS || (RAD_NUM_TYPES != 4)) {
+      !NEUTRINO_OSCILLATIONS_FFI || (RAD_NUM_TYPES != 4)) {
     if (mpi_io_proc()) {
       printf("Test needs the following physics active!\n"
              "\tRADIATION == RADTYPE_NEUTRINOS: %d\n"
              "\tLOCAL_ANGULAR_DISTRIBUTIONS: %d\n"
-             "\tNEUTRINO_OSCILLATIONS: %d\n"
+             "\tNEUTRINO_OSCILLATIONS_FFI: %d\n"
              "\tRAD_NUM_TYPES == 4: %d\n",
-          RADIATION, LOCAL_ANGULAR_DISTRIBUTIONS, NEUTRINO_OSCILLATIONS,
+          RADIATION, LOCAL_ANGULAR_DISTRIBUTIONS, NEUTRINO_OSCILLATIONS_FFI,
           RAD_NUM_TYPES);
     }
     exit(1);
diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py
index b39e2d5..1ac3561 100644
--- a/prob/torus_cbc/build.py
+++ b/prob/torus_cbc/build.py
@@ -40,7 +40,7 @@
 TRACERTEST = '-tracertest' in sys.argv
 RESTARTTEST = '-restarttest' in sys.argv
 HDF = '-hdf' in sys.argv
-OSCILLATIONS = "-oscillations" in sys.argv
+FFI = "-ffi" in sys.argv
 N1N2N3CPU_FROM_CLI = '-n1n2n3cpu' in sys.argv
 N1N2N3TOT_FROM_CLI = '-n1n2n3tot' in sys.argv
 
@@ -116,7 +116,7 @@
 if FORTRAN:
     OPACPATH = "opacity.SFHo.nohoro.juo.brem1.bin"
 else:
-    if OSCILLATIONS:
+    if FFI:
         OPACPATH = "NuLib_rho70_temp62_ye50_ng61_ns4_version1.0_20241120_bhlight.h5"
     else:
         OPACPATH = "opacity.SFHo.nohoro.juo.brem1.h5"
@@ -426,13 +426,13 @@
 bhl.config.set_cparm('RZ_HISTOGRAMS', True)
 # bhl.config.set_cparm('RECORD_DT_MIN', True)
 
-if OSCILLATIONS:
+if FFI:
     bhl.config.set_cparm('LOCAL_ANGULAR_DISTRIBUTIONS', True)
     bhl.config.set_cparm('LOCAL_ANGLES_NMU', 32)
     bhl.config.set_cparm('LOCAL_ANGLES_NX1', 32)
     bhl.config.set_cparm('LOCAL_ANGLES_NX2', 32)
     bhl.config.set_cparm('RAD_NUM_TYPES', 4)
-    bhl.config.set_cparm('NEUTRINO_OSCILLATIONS', True)
+    bhl.config.set_cparm('NEUTRINO_OSCILLATIONS_FFI', True)
 
 # Special. Don't turn this on if you don't need to
 if DIAGNOSTIC:
diff --git a/script/config.py b/script/config.py
index 892091d..91ea80f 100644
--- a/script/config.py
+++ b/script/config.py
@@ -345,15 +345,19 @@ def build(PROBLEM, PATHS):
       else:
         set_cparm("RAD_NUM_TYPES", 3)
       print_config("RAD_NUM_TYPES", CPARMS["RAD_NUM_TYPES"])
-    if util.parm_is_active(CPARMS, "NEUTRINO_OSCILLATIONS"):
-      print_config("NEUTRINO_OSCILLATIONS ", CPARMS["NEUTRINO_OSCILLATIONS"])
+    if util.parm_is_active(CPARMS, "NEUTRINO_OSCILLATIONS_FFI"):
+      print_config("NEUTRINO_OSCILLATIONS_FFI ", CPARMS["NEUTRINO_OSCILLATIONS_FFI"])
       if util.parm_is_active(CPARMS, 'FORCE_EQUIPARTITION'):
         print_config("FORCE_EQUIPARTITION ", CPARMS["FORCE_EQUIPARTITION"])
       else:
         set_cparm('FORCE_EQUIPARTITION', 0)
     else:
-      set_cparm('NEUTRINO_OSCILLATIONS', 0)
+      set_cparm('NEUTRINO_OSCILLATIONS_FFI', 0)
       set_cparm('FORCE_EQUIPARTITION', 0)
+    if util.parm_is_active(CPARMS, "NEUTRINO_OSCILLATIONS_CFI"):
+      print_config("NEUTRINO_OSCILLATIONS_CFI ", CPARMS["NEUTRINO_OSCILLATIONS_CFI"])
+    else:
+      set_cparm('NEUTRINO_OSCILLATIONS_CFI', 0)
     if util.parm_is_active(CPARMS, "RZ_HISTOGRAMS"):
       print_config("RZ_HISTOGRAMS", CPARMS["RZ_HISTOGRAMS"])
       if not util.parm_is_active(CPARMS, "RZ_HISTOGRAMS_N"):

From cf3bd2cc47b49f328da096014d3bafd7e5e2c6eb Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Wed, 4 Mar 2026 23:38:00 -0500
Subject: [PATCH 12/16] CFI correct name.

---
 core/decs.h | 10 +++++-----
 core/defs.h |  4 ++--
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index 1982249..2f12807 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -400,8 +400,8 @@ extern grid_local_count_type     local_osc_count;
 #endif // #if RAD_NUM_TYPES >= 4
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
-#if DO_CFI
-typedef double grid_double_type[N1 + 2 * NG][N2 + 2 * NG];
+#if NEUTRINO_OSCILLATIONS_CFI
+typedef double grid_symm_double_type[N1 + 2 * NG][N2 + 2 * NG];
 typedef int grid_symm_int_type[N1 + 2 * NG][N2 + 2 * NG];
 typedef double grid_symm_radtype_type[N1 + 2 * NG][N2 + 2 * NG][RAD_NUM_TYPES];
 typedef double grid_CFI_Gamma_type[N1 + 2 * NG][N2 + 2 * NG][2];
@@ -416,7 +416,7 @@ extern grid_CFI_Gamma_type cfi_Gamma;
 extern grid_symm_int_type cfi_active_mode;
 // time scale for asymptotic state
 extern grid_symm_double_type cfi_tau_asymp;
-#endif // DO_CFI
+#endif // NEUTRINO_OSCILLATIONS_CFI
 
 #endif // RADIATION
 
@@ -1063,9 +1063,9 @@ void oscillate_ffi(grid_local_angles_type f,
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
 // CFI
-#if DO_CFI
+#if NEUTRINO_OSCILLATIONS_CFI
 void compute_cfi_active_mode(grid_int_type cfi_active_mode);
-#endif // DO_CFI
+#endif // NEUTRINO_OSCILLATIONS_CFI
 
 #endif // RADIATION
 
diff --git a/core/defs.h b/core/defs.h
index 8e68f94..eef677a 100644
--- a/core/defs.h
+++ b/core/defs.h
@@ -68,7 +68,7 @@ grid_local_count_type local_osc_count;
 #endif //  RAD_NUM_TYPES
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
-#if DO_CFI
+#if NEUTRINO_OSCILLATIONS_CFI
 // number of neutrinos per flavor
 grid_symm_radtype_type nph_flavor;
 // individual distribution function-weighted opacity averages
@@ -79,7 +79,7 @@ grid_CFI_Gamma_type cfi_Gamma;
 grid_symm_int_type cfi_active_mode;
 // time scale for asymptotic state
 grid_symm_double_type cfi_tau_asymp;
-#endif // DO_CFI
+#endif // NEUTRINO_OSCILLATIONS_CFI
 
 #endif // RADIATION
 

From f2470ad71e13815c9f5e8bdfa2c0a35ddc60c5da Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 5 Mar 2026 00:01:05 -0500
Subject: [PATCH 13/16] add some todos for this project

---
 core/emissivity.c        | 4 ++++
 core/make_superphotons.c | 4 ++++
 core/radiation.c         | 2 ++
 3 files changed, 10 insertions(+)

diff --git a/core/emissivity.c b/core/emissivity.c
index 0ce6e26..8c1166e 100644
--- a/core/emissivity.c
+++ b/core/emissivity.c
@@ -46,6 +46,7 @@ void init_emissivity() {
 #endif
 }
 
+// todo(JMM): Pass in a rescaling option
 double jnu(double nu, int type, const struct of_microphysics *m, double theta) {
   double jnu = 0.;
 
@@ -71,6 +72,7 @@ double jnu(double nu, int type, const struct of_microphysics *m, double theta) {
   return jnu;
 }
 
+// todo(JMM): Pass in a rescaling option
 double Jnu(double nu, int type, const struct of_microphysics *m) {
   double Jnu = 0.;
 
@@ -96,6 +98,7 @@ double Jnu(double nu, int type, const struct of_microphysics *m) {
   return Jnu;
 }
 
+// todo(JMM): Pass in a rescaling option
 double integrandJ(double x, void *params) {
   struct of_J_params *    p    = (struct of_J_params *)params;
   struct of_microphysics *m    = p->microphysics;
@@ -110,6 +113,7 @@ double integrandJ(double x, void *params) {
   return Jsamp;
 }
 
+// todo(JMM): Pass in a rescaling option
 double get_J(struct of_microphysics *m) {
 #if RADIATION == RADTYPE_LIGHT
   {
diff --git a/core/make_superphotons.c b/core/make_superphotons.c
index 4f77cd7..8eaa734 100644
--- a/core/make_superphotons.c
+++ b/core/make_superphotons.c
@@ -153,6 +153,7 @@ void sample_photon(int i, int j, int k, double t, double dt, int type,
   double weight = get_wgt(nu, dtau);
 
   // Sample emissivity in solid angle
+  // todo(JMM): Pass in a rescaling option
   double jmax = jnu(nu, type, m, 0.5 * M_PI);
   do {
     cth[0] = 2. * get_rand() - 1.;
@@ -288,6 +289,7 @@ void get_dndlnu(int i, int j, int k, double dt, double dndlnu[NU_BINS + 1],
 
   double dndlnu_max = -1.e100;
   for (int n = 0; n <= NU_BINS; n++) {
+    // todo(JMM): Pass in a rescaling option
     double Jsamp = Jnu(nusamp[n], type, m);
     Jsamp *= dx[1] * dx[2] * dx[3] * pow(L_unit, 3.) * ggeom[i][j][CENT].g;
 
@@ -330,6 +332,7 @@ void set_weight(grid_prim_type Prad, grid_eosvar_type extra) {
       get_fluid_zone(i, j, k, Prad, extra, &m, Ucon, Ucov, Bcon, Bcov);
       TYPELOOP {
         for (int n = 0; n <= NU_BINS; n++) {
+          // todo(JMM): Pass in a rescaling option
           Jtot += Jnu(nusamp[n], itp, &m) * zoneVol * ggeom[i][j][CENT].g;
         }
       } // TYPELOOP
@@ -353,6 +356,7 @@ double f(double x, void *params) {
   int                     type = p->type;
 
   double nu    = exp(x);
+  // todo(JMM): Pass in a rescaling option
   double Jsamp = Jnu(nu, type, m) * nu;
   double wgt   = get_wgt(nu, get_dtau(nu, type, dt, m));
   if (isinf(Jsamp) || wgt < SMALL) {
diff --git a/core/radiation.c b/core/radiation.c
index aa56f22..7dd393f 100644
--- a/core/radiation.c
+++ b/core/radiation.c
@@ -27,6 +27,7 @@ double Bnu_inv(double nu,
 #endif
 }
 
+// todo(JMM): Pass in a rescaling option
 double jnu_inv(
     double nu, int type, const struct of_microphysics *m, double theta) {
   double j;
@@ -47,6 +48,7 @@ double alpha_inv_scatt(
 }
 
 // Invariant absorption opacity
+// TODO(JMM): Need to pass in a rescaling option
 double alpha_inv_abs(
     double nu, int type, const struct of_microphysics *m, double theta) {
 #if RADIATION == RADTYPE_NEUTRINOS

From cb0f94f1d5441f93b36b4d35fc5dbab0c41cdcab Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 5 Mar 2026 00:01:46 -0500
Subject: [PATCH 14/16] start building the averaging machinery

---
 core/oscillations.c | 53 ++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 50 insertions(+), 3 deletions(-)

diff --git a/core/oscillations.c b/core/oscillations.c
index 26bddcc..3aad444 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -189,10 +189,10 @@ void oscillate_ffi(grid_local_angles_type f,
 #pragma omp parallel
   {
     struct of_photon *ph = photon_lists[omp_get_thread_num()];
-    double            X[NDIM], Kcov[NDIM], Kcon[NDIM];
     while (ph != NULL) {
       if (ph->type != TYPE_TRACER) {
         int ix1, ix2, icosth[LOCAL_NUM_BASES];
+        double X[NDIM], Kcov[NDIM], Kcon[NDIM];
         get_local_angle_bins(ph, &ix1, &ix2, &icosth[0], &icosth[1]);
 
         int    b_osc = local_b_osc[ix1][ix2];
@@ -254,7 +254,7 @@ void oscillate_ffi(grid_local_angles_type f,
             ph->type = (ph->type + (RAD_NUM_TYPES / 2)) % RAD_NUM_TYPES;
             ph->osc_count += 1;
 
-#pragma omp atomic
+            #pragma omp atomic
             local_osc_count[ix1][ix2] += ph->w;
           }
         }
@@ -266,6 +266,53 @@ void oscillate_ffi(grid_local_angles_type f,
 }
 
 #endif // RAD_NUM_TYPES >= 4
-
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
+
+#if NEUTRINO_OSCILLATIONS_CFI
+void compute_cfi_active_mode(grid_int_type cfi_active_mode) {
+  timer_start(TIMER_OSCILLATIONS);
+
+  static const int SYMM_SIZE = (N1 + 2*NG)*(N2 + 2*NG);
+
+  memset(nph_flavor, 0, SYMM_SIZE * RAD_NUM_TYPES * sizeof(double));
+  memset(kappa_avg, 0, SYMM_SIZE * RAD_NUM_TYPES * sizeof(double));
+  memset(cfi_Gamma, 0, 2 * SYMM_SIZE * sizeof(double));
+
+#pragma omp parallel
+  {
+    struct of_photon *ph = photon_lists[omp_get_thread_num()];
+    while (ph != NULL) {
+      if (ph->type != TYPE_TRACER) {
+        double X[NDIM], Kcov[NDIM], Kcon[NDIM];
+        int i, j, k;
+        get_X_K_interp(ph, t, P, X, Kcov, Kcon);
+        Xtoijk(X, &i, &j, &k);
+        
+        #pragma omp atomic
+        nph_flavor[i][j][ph->type] += ph->w;
+
+        double theta = get_bk_angle(
+            X, Kcon, Ucov_grd[i][j][k],
+            Bcov_grd[i][j][k], m_grd[i][j][k].B);
+        double nu = get_fluid_nu(X, Kcov, Ucon_grd[i][j][k]);
+        double kappa = alpha_inv_abs(nu, ph->type,
+                                     &(m_grd[i][j][k]), theta)
+          / (nu + SMALL);
+        #pragma omp atomic
+        kappa_avg[i][j][ph->type] += kappa * ph->w;
+
+        // TODO: DO Gamma
+      }
+      ph = ph->next;
+    }
+  }
+
+  mpi_dbl_allreduce_array((double *)nph_flavor, SYMM_SIZE * RAD_NUM_TYPES);
+  mpi_dbl_allreduce_array((double *)kappa_avg, SYMM_SIZE * RAD_NUM_TYPES);
+  mpi_dbl_allreduce_array((double *)cfi_Gamma, 2 * SYMM_SIZE);
+
+  timer_stop(TIMER_OSCILLATIONS);
+}
+#endif // NEUTRINO_OSCILLATIONS_CFI
+
 #endif // RADIATION

From 1e271204c14f8ba3f54664d6264aeb43fdace95d Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 5 Mar 2026 15:13:38 -0500
Subject: [PATCH 15/16] compute average gammas and kappas for growth rate calc

---
 core/decs.h         |  4 ++++
 core/oscillations.c | 50 +++++++++++++++++++++++++++++++++++++++------
 core/rad_utils.c    | 10 ++++++++-
 3 files changed, 57 insertions(+), 7 deletions(-)

diff --git a/core/decs.h b/core/decs.h
index 2f12807..d149ab9 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -1064,6 +1064,9 @@ void oscillate_ffi(grid_local_angles_type f,
 
 // CFI
 #if NEUTRINO_OSCILLATIONS_CFI
+void compute_cfi_symmetrized_avgs(grid_symm_radtype_type nph_flavor,
+                                  grid_symm_radtype_type kappa_avg,
+                                  grid_CFI_Gamma_type cfi_gamma);
 void compute_cfi_active_mode(grid_int_type cfi_active_mode);
 #endif // NEUTRINO_OSCILLATIONS_CFI
 
@@ -1150,6 +1153,7 @@ void record_lepton_flux(const struct of_photon *ph);
 void check_nu_type(const char *location);
 int  get_lepton_sign(const struct of_photon *ph);
 int  nu_is_heavy(const int radtype);
+int is_antiparticle(const struct of_photon *ph);
 #endif // NEUTRINOS
 #endif // RADIATION
 
diff --git a/core/oscillations.c b/core/oscillations.c
index 3aad444..5b1c23b 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -269,7 +269,9 @@ void oscillate_ffi(grid_local_angles_type f,
 #endif // LOCAL_ANGULAR_DISTRIBUTIONS
 
 #if NEUTRINO_OSCILLATIONS_CFI
-void compute_cfi_active_mode(grid_int_type cfi_active_mode) {
+void compute_cfi_symmetrized_avgs(grid_symm_radtype_type nph_flavor,
+                                  grid_symm_radtype_type kappa_avg,
+                                  grid_CFI_Gamma_type cfi_gamma) {
   timer_start(TIMER_OSCILLATIONS);
 
   static const int SYMM_SIZE = (N1 + 2*NG)*(N2 + 2*NG);
@@ -284,6 +286,7 @@ void compute_cfi_active_mode(grid_int_type cfi_active_mode) {
     while (ph != NULL) {
       if (ph->type != TYPE_TRACER) {
         double X[NDIM], Kcov[NDIM], Kcon[NDIM];
+        double kappa[RAD_NUM_TYPES];
         int i, j, k;
         get_X_K_interp(ph, t, P, X, Kcov, Kcon);
         Xtoijk(X, &i, &j, &k);
@@ -295,13 +298,40 @@ void compute_cfi_active_mode(grid_int_type cfi_active_mode) {
             X, Kcon, Ucov_grd[i][j][k],
             Bcov_grd[i][j][k], m_grd[i][j][k].B);
         double nu = get_fluid_nu(X, Kcov, Ucon_grd[i][j][k]);
-        double kappa = alpha_inv_abs(nu, ph->type,
-                                     &(m_grd[i][j][k]), theta)
-          / (nu + SMALL);
+        TYPELOOP {
+          kappa[itp] = alpha_inv_abs(nu, itp,
+                                     &(m_grd[i][j][k]), theta);
+          kappa[itp] /= (nu + SMALL);
+        }
+
         #pragma omp atomic
-        kappa_avg[i][j][ph->type] += kappa * ph->w;
+        kappa_avg[i][j][ph->type] += kappa[ph->type] * ph->w;
+
+        // electrons and antis positive, heavies and antis negative
+        int gamma_sign = (ph->type < NU_HEAVY) ? 1 : -1;
+
+        // Check for the need to divide neutrino count for heavies by 2
+        double heavy_dupfac = ((RAD_NUM_TYPES < 4) && (ph->type >= NU_HEAVY))? 0.5 : 1.0;
+
+        // check for antiparticles
+        int isanti = is_antiparticle(ph);
+
+        int kappa_e_idx = NU_ELECTRON + isanti; // kappa_e for es, kappa_ebar for ebars
+        // kappa_x for xs, kappa_xbar for xbars. But account for x = xbar
+        int kappa_x_idx = NU_HEAVY + (RAD_NUM_TYPES == 4)*isanti; 
 
-        // TODO: DO Gamma
+        // es contribute + to this sum, xs minus. Need to divide by 2
+        // if heavies and their antis are bundled together.
+        double weight = gamma_sign * heavy_dupfac * ph->w;
+
+        // gamma_0 = (1/2) (kappa_e + kappa_x)
+        // gamma_1 = (1/2) (kappa_ebar + kappa_xbar)
+        //
+        // Method A
+        // avg gamma_0 = \int E^2 dE (kappa_e + kappa_x) (f_e - f_x) / \int E^2 / \int E^2 dE (f_e - f_x)
+        // avg gamma_1 = \int E^2 dE (kappa_ebar + kappa_xbar) (f_ebar - f_xbar) / \int E^2 / \int E^2 dE (f_ebar - f_xbar)
+        #pragma omp atomic
+        cfi_Gamma[i][j][isanti] += 0.5*(kappa[kappa_e_idx] + kappa[kappa_x_idx])*weight;
       }
       ph = ph->next;
     }
@@ -311,6 +341,14 @@ void compute_cfi_active_mode(grid_int_type cfi_active_mode) {
   mpi_dbl_allreduce_array((double *)kappa_avg, SYMM_SIZE * RAD_NUM_TYPES);
   mpi_dbl_allreduce_array((double *)cfi_Gamma, 2 * SYMM_SIZE);
 
+  #pragma omp parallel for collapse(2)
+  ILOOP {
+    JLOOP {
+      cfi_Gamma[i][j][0] /= (nph_flavor[i][j][NU_ELECTRON] - nph_flavor[i][j][NU_HEAVY]);
+      cfi_Gamma[i][j][1] /= (nph_flavor[i][j][ANTINU_ELECTRON] - nph_flavor[i][j][ANTINU_HEAVY]);
+    }
+  }
+
   timer_stop(TIMER_OSCILLATIONS);
 }
 #endif // NEUTRINO_OSCILLATIONS_CFI
diff --git a/core/rad_utils.c b/core/rad_utils.c
index eddcfe5..9bbf6da 100644
--- a/core/rad_utils.c
+++ b/core/rad_utils.c
@@ -749,7 +749,7 @@ double get_min_dt_cool(grid_prim_type P, grid_eosvar_type extra) {
 }
 
 #if RADIATION == RADTYPE_NEUTRINOS
-void        record_lepton_flux(const struct of_photon *ph) {
+void record_lepton_flux(const struct of_photon *ph) {
 #pragma omp atomic
   lepton_lost_local += (ph->w) * get_lepton_sign(ph);
 }
@@ -765,6 +765,14 @@ int nu_is_heavy(const int radtype) {
   return ((radtype == NU_HEAVY) || (radtype == ANTINU_HEAVY));
 }
 
+int is_antiparticle(const struct of_photon *ph) {
+#if RAD_NUM_TYPES == 4
+  return (ph->type == ANTINU_ELECTRON) || (ph->type == ANTINU_HEAVY);
+#else
+  return (ph->type == ANTINU_ELECTRON) || (ph->type == NU_HEAVY);
+#endif
+}
+
 // for debugging
 void        check_nu_type(const char *location) {
 #pragma omp parallel

From eed23c2bfbdb9c6c984bbf6324770d57430ab995 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 5 Mar 2026 23:39:11 -0500
Subject: [PATCH 16/16] Bnu inv now supports Fermi-Dirac

---
 core/radiation.c | 24 +++++++++++++++++-------
 1 file changed, 17 insertions(+), 7 deletions(-)

diff --git a/core/radiation.c b/core/radiation.c
index 7dd393f..d324236 100644
--- a/core/radiation.c
+++ b/core/radiation.c
@@ -9,21 +9,31 @@
 #include "decs.h"
 
 #if RADIATION
-double Bnu_inv(double nu,
-    const struct of_microphysics
-        *m) { // unused for neutrinos. Use tables by Burrows et al.
-#if RADIATION == RADTYPE_LIGHT
+double Bnu_inv(double nu, const struct of_microphysics *m) {
   double x;
+#if RADIATION == RADTYPE_LIGHT
   x = HPL * nu / (ME * CL * CL * m->Thetae);
+#elif RADIATION == RADTYPE_NEUTRINOS
+  x = HPL * nu / (KBOL * m->T);
+#endif
 
+#if RADIATION == RADTYPE_LIGHT
+  // Planck (i.e., Bose-Einstein) distribution
   if (x < 1.e-3) // Taylor-expand small arguments for numerical accuracy
     return (
         (2. * HPL / (CL * CL)) / (x / 24. * (24. + x * (12. + x * (4. + x)))));
   else
     return ((2. * HPL / (CL * CL)) / (exp(x) - 1.));
-#endif
-#if RADIATION == RADTYPE_NEUTRINOS
-  return 0.; // STUB
+#elif RADIATION == RADTYPE_NEUTRINOS
+  // Fermi-Dirac distribution
+  double dist;
+  if (x < 1.e-3) { // Taylor-expand small arguments for numerical accuracy
+    double x2 = x*x;
+    dist = 0.5*(1 - 0.5*x*(1 - (x2/12)*(1 - 0.1*x2)));
+  } else {
+    dist = 1. / (exp(x) + 1.);
+  }
+  return (2. * HPL / (CL * CL)) * dist;
 #endif
 }