diff --git a/core/decs.h b/core/decs.h
index 8eb3680..5914fd1 100644
--- a/core/decs.h
+++ b/core/decs.h
@@ -292,9 +292,10 @@
 #define TIMER_MAKE (7)
 #define TIMER_PUSH (8)
 #define TIMER_INTERACT (9)
-#define TIMER_MICRO (10)
-#define TIMER_ALL (11)
-#define NUM_TIMERS (12)
+#define TIMER_OSCILLATIONS (10)
+#define TIMER_MICRO (11)
+#define TIMER_ALL (12)
+#define NUM_TIMERS (13)
 
 // Units
 #define NEED_UNITS (RADIATION || EOS == EOS_TYPE_TABLE)
diff --git a/core/diag.c b/core/diag.c
index e326bdd..949d021 100644
--- a/core/diag.c
+++ b/core/diag.c
@@ -298,6 +298,9 @@ void diag(int call_code) {
           get_time_per_step(TIMER_ALL));
 #if ELECTRONS
       fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_ELECTRON));
+#endif
+#if NEUTRINO_OSCILLATIONS || LOCAL_ANGULAR_DISTRIBUTIONS
+      fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_OSCILLATIONS));
 #endif
       fprintf(ener_file, "\n");
       fflush(ener_file);
diff --git a/core/io.c b/core/io.c
index 0bea6f6..43d3907 100644
--- a/core/io.c
+++ b/core/io.c
@@ -873,6 +873,13 @@ void dump() {
   WRITE_HDR(nubins_spec, TYPE_INT);
   WRITE_HDR(numin, TYPE_DBL);
   WRITE_HDR(numax, TYPE_DBL);
+
+  int neutrino_oscillations = NEUTRINO_OSCILLATIONS;
+  WRITE_HDR(neutrino_oscillations, TYPE_INT);
+  int force_equipartition = FORCE_EQUIPARTITION;
+  WRITE_HDR(force_equipartition, TYPE_INT);
+  int local_angular_distributions = LOCAL_ANGULAR_DISTRIBUTIONS;
+  WRITE_HDR(local_angular_distributions, TYPE_INT);
 #endif
 
 #if NVAR_PASSIVE > 0
diff --git a/core/oscillations.c b/core/oscillations.c
index 240907d..4574b14 100644
--- a/core/oscillations.c
+++ b/core/oscillations.c
@@ -11,6 +11,7 @@
 #if LOCAL_ANGULAR_DISTRIBUTIONS
 
 double get_dt_oscillations() {
+  timer_start(TIMER_OSCILLATIONS);
   set_Rmunu(); // So we have Nsph and nph
   double nph_max = 0;
 #pragma omp parallel for reduction(max : nph_max) collapse(3)
@@ -21,10 +22,12 @@ double get_dt_oscillations() {
   double dt_osc = 1. / (NUFERM * nph_max + SMALL);
   dt_osc /= T_unit; // code units
   dt_osc = mpi_min(dt_osc);
+  timer_stop(TIMER_OSCILLATIONS);
   return dt_osc;
 }
 
 void accumulate_local_angles() {
+  timer_start(TIMER_OSCILLATIONS);
   static const int LOCAL_ANGLES_SIZE = LOCAL_NUM_BASES * LOCAL_ANGLES_NX1 *
                                        LOCAL_ANGLES_NX2 * LOCAL_ANGLES_NMU *
                                        RAD_NUM_TYPES;
@@ -64,6 +67,7 @@ void accumulate_local_angles() {
   compute_local_gnu(local_angles, local_Ns, local_wsqr, Gnu);
   compute_local_moments(Gnu, local_moments);
 #endif // RAD_NUM_TYPES >= 4
+  timer_stop(TIMER_OSCILLATIONS);
 }
 
 void get_local_angle_bins(
@@ -182,7 +186,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_moment_type local_moments, grid_Gnu_type gnu) {
+  timer_start(TIMER_OSCILLATIONS);
 #pragma omp parallel
   {
     struct of_photon *ph = photon_lists[omp_get_thread_num()];
@@ -233,6 +238,7 @@ void        oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) {
       ph = ph->next;
     }
   } // omp parallel
+  timer_stop(TIMER_OSCILLATIONS);
 }
 
 #endif // RAD_NUM_TYPES >= 4
diff --git a/core/timing.c b/core/timing.c
index 605fce9..6a5b3c3 100644
--- a/core/timing.c
+++ b/core/timing.c
@@ -91,6 +91,11 @@ 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
+    fprintf(stdout, "   OSCILL:   %8.4g s (%.4g %%)\n",
+        times[TIMER_OSCILLATIONS] / dnstep,
+        100. * times[TIMER_OSCILLATIONS] / times[TIMER_ALL]);
+#endif
     fprintf(stdout, "   MICRPHYS: %8.4g s (%.4g %%)\n",
         times[TIMER_MICRO] / dnstep,
         100. * times[TIMER_MICRO] / times[TIMER_ALL]);
diff --git a/script/analysis/hdf5_to_dict.py b/script/analysis/hdf5_to_dict.py
index 6ff11e1..dc2a856 100644
--- a/script/analysis/hdf5_to_dict.py
+++ b/script/analysis/hdf5_to_dict.py
@@ -115,6 +115,9 @@ def read_hdr_var(hdrname,varname,default = 0):
   read_hdr_var('nubins_spec','nubins_spec',200)
   read_hdr_var('diagnostics_use_radtypes','diagnostics_use_radtypes',0)
   read_hdr_var('FULL_DUMP','FULL_DUMP',1)
+  read_hdr-var("LOCAL_ANGULAR_DISTRIBUTIONS", "local_angular_distributions", 0)
+  read_hdr_var("NEUTRINO_OSCILLATIONS", "neutrino_oscillations", 0)
+  read_hdr_var("FORCE_EQUIPARTITION", "force_equipartition", 0)
 
   hdr['vnams'] = [h5_to_string(i) for i in dfile['P'].attrs['vnams']]
 
@@ -383,6 +386,16 @@ def load_dump(fname, geom=None, nulegacy=False):
         dump['Nabs_x'] = dump['Nabs_phys'][:,:,:,2]
       else:
         dump['Nabs_ph'] = dump['Nabs_phys'][:,:,:,0]
+    if hdr['NEUTRINO_OSCILLATIONS'] or hdr['LOCAL_ANGULAR_DISTRIBUTIONS']:
+      dump['local_angles'] = dfile['local_angles'][:]
+      dump['Gnu'] = dfile['Gnu'][:]
+      dump['local_Ns'] = dfile['local_Ns'][:]
+      dump['local_wsqr'] = dfile['local_wsqr'][:]
+      wmean = dump['local_angles'].sum(axis=3) / dump['local_Ns']
+      Ns = dump['local_Ns'][:]
+      wb2N = Ns*wmean*wmean
+      w2 = dump['local_wsqr'][:]
+      dump['local_stddev'] = np.sqrt(wb2N + N2*(w2 - wb2N)/(Ns - 1))
 
   ucon, ucov, bcon, bcov = get_state(dump, geom)
 
@@ -561,6 +574,9 @@ def load_diag(path, hdr = None, nulegacy = False, timers = True):
     if hdr['ELECTRONS']:
       diag['TIMER_ELECTRON'] = dfile[nbase + 1]
       nbase += 1
+    if hdr['NEUTRINO_OSCILLATIONS'] or hdr['LOCAL_ANGULAR_DISTRIBUTIONS']:
+      diag['TIMER_OSCILLATIONS'] = dfile[nbase + 1]
+      nbase += 1
 
   # Purge vanishing data due to restarts
   restart_mask = np.ones(len(diag['t']),dtype=bool)