From 1dd87f26a5c722e756db763560924dc1f84e5318 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Fri, 28 Feb 2020 17:14:08 -0700
Subject: [PATCH 01/14] changes to meep_mpb_eps so that s, o, and p are
 rotated. Still need to figure out how to rotate the fields object f

---
 src/mpb.cpp | 33 ++++++++++++++++++++++++++++++---
 1 file changed, 30 insertions(+), 3 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index ea1637801..08db40acf 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -45,15 +45,42 @@ typedef struct {
 
 static void meep_mpb_eps(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
                          void *eps_data_) {
+  int coordcycle = 0;
+  return meep_mpb_eps_coordcycle(eps, eps_inv, r, eps_data_, coordcycle);
+}
+
+static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
+                         void *eps_data_, int coordcycle = 0) {
   meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_;
-  const double *s = eps_data->s;
-  const double *o = eps_data->o;
+  double *s_temp = eps_data->s;
+  double *o_temp = eps_data->o;
   double omega = eps_data->omega;
   vec p(eps_data->dim == D3 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1], o[2] + r[2] * s[2])
                             : (eps_data->dim == D2 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1]) :
                                                    /* D1 */ vec(o[2] + r[2] * s[2])));
-  const fields *f = eps_data->f;
+  fields *f_temp = eps_data->f;
+
+  if (eps_date->dim == D3 && coordcycle) {
+    int i;
+    for (i = 0; i < coordcycle; i++) {
+      double s_copy[3] = *s_temp;
+      double o_copy[3] = *o_temp;
+      s_temp[0] = s_copy[1];
+      s_temp[1] = s_copy[2];
+      s_temp[2] = s_copy[0];
+      o_temp[0] = o_copy[1];
+      o_temp[1] = o_copy[2];
+      o_temp[2] = o_copy[0];
+      p = vec(p.y(), p.z(), p.x());
+      //operation on f ?
+    }
+  }
+
+  const double *s = s_temp;
+  const double *o = o_temp;
+  const fields *f = f_temp;
 
+  // need to find out what eps_inv is doing exactly, but I believe all the necessary changes are made above this
   eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
   eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));
   eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega));

From 724efda16e145c0258a44e409759c1865addf72e Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Fri, 28 Feb 2020 22:47:49 -0700
Subject: [PATCH 02/14] Changes to get_eigenmode for coordinate system cycling

---
 src/mpb.cpp | 56 ++++++++++++++++++++++++++++++++++++++++++++---------
 1 file changed, 47 insertions(+), 9 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 08db40acf..94233df39 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -51,6 +51,9 @@ static void meep_mpb_eps(symmetric_matrix *eps, symmetric_matrix *eps_inv, const
 
 static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
                          void *eps_data_, int coordcycle = 0) {
+  if (coordcycle != 0 && coordcycle != 1 && coordcycle != 2) {
+    abort("unsupported coordcycle value");
+  }
   meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_;
   double *s_temp = eps_data->s;
   double *o_temp = eps_data->o;
@@ -282,6 +285,14 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
                             int band_num, const vec &_kpoint, bool match_frequency, int parity,
                             double resolution, double eigensolver_tol, double *kdom,
                             void **user_mdata) {
+  return get_eigenmode_coordcycle(omega_src, d, where, eig_vol, band_num, &_kpoint, match_frequency, parity, resolution,
+                                  eigensolver_tol, kdom, user_mdata, 0);
+}
+
+void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volume where, const volume eig_vol,
+                            int band_num, const vec &_kpoint, bool match_frequency, int parity,
+                            double resolution, double eigensolver_tol, double *kdom,
+                            void **user_mdata, int coordcycle) {
   /*--------------------------------------------------------------*/
   /*- part 1: preliminary setup for calling MPB  -----------------*/
   /*--------------------------------------------------------------*/
@@ -325,15 +336,42 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
 
   switch (gv.dim) {
     case D3:
-      o[0] = eig_vol.in_direction_min(X);
-      o[1] = eig_vol.in_direction_min(Y);
-      o[2] = eig_vol.in_direction_min(Z);
-      s[0] = eig_vol.in_direction(X);
-      s[1] = eig_vol.in_direction(Y);
-      s[2] = eig_vol.in_direction(Z);
-      kcart[0] = kpoint.in_direction(X);
-      kcart[1] = kpoint.in_direction(Y);
-      kcart[2] = kpoint.in_direction(Z);
+      switch (coordcycle) {
+        case 0:
+          o[0] = eig_vol.in_direction_min(X);
+          o[1] = eig_vol.in_direction_min(Y);
+          o[2] = eig_vol.in_direction_min(Z);
+          s[0] = eig_vol.in_direction(X);
+          s[1] = eig_vol.in_direction(Y);
+          s[2] = eig_vol.in_direction(Z);
+          kcart[0] = kpoint.in_direction(X);
+          kcart[1] = kpoint.in_direction(Y);
+          kcart[2] = kpoint.in_direction(Z);
+          break;
+        case 1:
+          o[0] = eig_vol.in_direction_min(Y);
+          o[1] = eig_vol.in_direction_min(Z);
+          o[2] = eig_vol.in_direction_min(X);
+          s[0] = eig_vol.in_direction(Y);
+          s[1] = eig_vol.in_direction(Z);
+          s[2] = eig_vol.in_direction(X);
+          kcart[0] = kpoint.in_direction(Y);
+          kcart[1] = kpoint.in_direction(Z);
+          kcart[2] = kpoint.in_direction(X);
+          break;
+        case 2:
+          o[0] = eig_vol.in_direction_min(Z);
+          o[1] = eig_vol.in_direction_min(X);
+          o[2] = eig_vol.in_direction_min(Y);
+          s[0] = eig_vol.in_direction(Z);
+          s[1] = eig_vol.in_direction(X);
+          s[2] = eig_vol.in_direction(Y);
+          kcart[0] = kpoint.in_direction(Z);
+          kcart[1] = kpoint.in_direction(X);
+          kcart[2] = kpoint.in_direction(Y);
+          break;
+        default: abort("unsupported coordcycle value");
+      }
       break;
     case D2:
       o[0] = eig_vol.in_direction_min(X);

From a1099806b1116d2be61f8eae76d892c89f9c96d2 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Fri, 28 Feb 2020 23:44:54 -0700
Subject: [PATCH 03/14] Change to meep_mpb_eps_coordcycle so that s, o, and f
 are simply no longer const

---
 src/mpb.cpp | 26 +++++++++++---------------
 1 file changed, 11 insertions(+), 15 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 94233df39..f8f3e0630 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -55,34 +55,30 @@ static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps
     abort("unsupported coordcycle value");
   }
   meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_;
-  double *s_temp = eps_data->s;
-  double *o_temp = eps_data->o;
+  double *s = eps_data->s;
+  double *o = eps_data->o;
   double omega = eps_data->omega;
   vec p(eps_data->dim == D3 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1], o[2] + r[2] * s[2])
                             : (eps_data->dim == D2 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1]) :
                                                    /* D1 */ vec(o[2] + r[2] * s[2])));
-  fields *f_temp = eps_data->f;
+  fields *f = eps_data->f;
 
   if (eps_date->dim == D3 && coordcycle) {
     int i;
     for (i = 0; i < coordcycle; i++) {
-      double s_copy[3] = *s_temp;
-      double o_copy[3] = *o_temp;
-      s_temp[0] = s_copy[1];
-      s_temp[1] = s_copy[2];
-      s_temp[2] = s_copy[0];
-      o_temp[0] = o_copy[1];
-      o_temp[1] = o_copy[2];
-      o_temp[2] = o_copy[0];
+      double s_copy[3] = *s;
+      double o_copy[3] = *s;
+      s[0] = s_copy[1];
+      s[1] = s_copy[2];
+      s[2] = s_copy[0];
+      o[0] = o_copy[1];
+      o[1] = o_copy[2];
+      o[2] = o_copy[0];
       p = vec(p.y(), p.z(), p.x());
       //operation on f ?
     }
   }
 
-  const double *s = s_temp;
-  const double *o = o_temp;
-  const fields *f = f_temp;
-
   // need to find out what eps_inv is doing exactly, but I believe all the necessary changes are made above this
   eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
   eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));

From d9ec3aa58650631078b361a0d2e3ba1d2de58988 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Fri, 28 Feb 2020 23:54:50 -0700
Subject: [PATCH 04/14] Add TODOs, fix typos

---
 src/mpb.cpp | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index f8f3e0630..7c4c46cba 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -75,11 +75,10 @@ static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps
       o[1] = o_copy[2];
       o[2] = o_copy[0];
       p = vec(p.y(), p.z(), p.x());
-      //operation on f ?
+      //TODO: operation on f ?
     }
   }
 
-  // need to find out what eps_inv is doing exactly, but I believe all the necessary changes are made above this
   eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
   eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));
   eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega));
@@ -281,7 +280,7 @@ void *fields::get_eigenmode(double omega_src, direction d, const volume where, c
                             int band_num, const vec &_kpoint, bool match_frequency, int parity,
                             double resolution, double eigensolver_tol, double *kdom,
                             void **user_mdata) {
-  return get_eigenmode_coordcycle(omega_src, d, where, eig_vol, band_num, &_kpoint, match_frequency, parity, resolution,
+  return get_eigenmode_coordcycle(omega_src, d, where, eig_vol, band_num, _kpoint, match_frequency, parity, resolution,
                                   eigensolver_tol, kdom, user_mdata, 0);
 }
 
@@ -453,6 +452,7 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
       kmatch = kcart_len;
   }
   else {
+    // TODO: rotate
     kmatch = G[d - X][d - X] * k[d - X]; // k[d] in cartesian
     kdir[d - X] = 1;                     // kdir = unit vector in d direction
   }
@@ -468,14 +468,17 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
       if (gv.dim == D2) k[2] = beta;
     }
     else {
+      // TODO: rotate
       k[d - X] = kmatch * R[d - X][d - X]; // convert to reciprocal basis
       if (eig_vol.in_direction(d) > 0 &&
           fabs(k[d - X]) > 0.4) // ensure k is well inside the Brillouin zone
         k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;
+      // end rotate to do
     }
     if (verbosity > 1) master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
   }
 
+  // TODO: rotate
   set_maxwell_data_parity(mdata, parity);
   update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);
 
@@ -562,6 +565,7 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
           if (gv.dim == D2) k[2] = beta;
         }
         else {
+          // TODO: rotate
           k[d - X] = kmatch * R[d - X][d - X];
         }
         update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);
@@ -703,6 +707,7 @@ double get_group_velocity(void *vedata) {
 
 vec get_k(void *vedata) {
   eigenmode_data *edata = (eigenmode_data *)vedata;
+  // TODO: rotate
   return vec(edata->Gk[0], edata->Gk[1], edata->Gk[2]);
 }
 

From 85bdf99f44ebd539cef6f0fb0bdd840d5e1a179e Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Tue, 10 Mar 2020 16:23:54 -0600
Subject: [PATCH 05/14] Fixed meep_mpb_eps_coordcycle

---
 src/mpb.cpp | 57 ++++++++++++++++++++++++++++++-----------------------
 1 file changed, 32 insertions(+), 25 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 7c4c46cba..953bdfced 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -55,37 +55,43 @@ static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps
     abort("unsupported coordcycle value");
   }
   meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_;
-  double *s = eps_data->s;
-  double *o = eps_data->o;
+  const double *s = eps_data->s;
+  const double *o = eps_data->o;
   double omega = eps_data->omega;
   vec p(eps_data->dim == D3 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1], o[2] + r[2] * s[2])
                             : (eps_data->dim == D2 ? vec(o[0] + r[0] * s[0], o[1] + r[1] * s[1]) :
                                                    /* D1 */ vec(o[2] + r[2] * s[2])));
-  fields *f = eps_data->f;
-
-  if (eps_date->dim == D3 && coordcycle) {
-    int i;
-    for (i = 0; i < coordcycle; i++) {
-      double s_copy[3] = *s;
-      double o_copy[3] = *s;
-      s[0] = s_copy[1];
-      s[1] = s_copy[2];
-      s[2] = s_copy[0];
-      o[0] = o_copy[1];
-      o[1] = o_copy[2];
-      o[2] = o_copy[0];
-      p = vec(p.y(), p.z(), p.x());
-      //TODO: operation on f ?
-    }
+  const fields *f = eps_data->f;
+
+  // operations only need to happen in chi1inv
+  switch(coordcycle) {
+    case 0:
+      eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
+      eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));
+      eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega));
+      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Z, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, omega)), 0);
+      break;
+    case 1:
+      eps_inv->m00 = real(f->get_chi1inv(Ex, Y, p, omega));
+      eps_inv->m11 = real(f->get_chi1inv(Ey, Z, p, omega));
+      eps_inv->m22 = real(f->get_chi1inv(Ez, X, p, omega));
+      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Z, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, X, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, X, p, omega)), 0);
+      break;
+    case 2:
+      eps_inv->m00 = real(f->get_chi1inv(Ex, Z, p, omega));
+      eps_inv->m11 = real(f->get_chi1inv(Ey, X, p, omega));
+      eps_inv->m22 = real(f->get_chi1inv(Ez, Y, p, omega));
+      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, X, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Y, p, omega)), 0);
+      break;
+    default: abort("unsupported coordcycle value");
   }
 
-  eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
-  eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));
-  eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega));
-
-  ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
-  ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Z, p, omega)), 0);
-  ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, omega)), 0);
   /*
   master_printf("m11(%g,%g) = %g\n", p.x(), p.y(), eps_inv->m00);
   master_printf("m22(%g,%g) = %g\n", p.x(), p.y(), eps_inv->m11);
@@ -329,6 +335,7 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
   if (!eig_vol.contains(where))
     abort("invalid grid_volume in add_eigenmode_source (WHERE must be in EIG_VOL)");
 
+  // TODO: change here
   switch (gv.dim) {
     case D3:
       switch (coordcycle) {

From f5c825f6c037a07526a7f716ecde8a130cef1e83 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Tue, 10 Mar 2020 16:25:46 -0600
Subject: [PATCH 06/14] Removed extra coordcycle check

---
 src/mpb.cpp | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 953bdfced..d78ab6560 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -51,9 +51,6 @@ static void meep_mpb_eps(symmetric_matrix *eps, symmetric_matrix *eps_inv, const
 
 static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
                          void *eps_data_, int coordcycle = 0) {
-  if (coordcycle != 0 && coordcycle != 1 && coordcycle != 2) {
-    abort("unsupported coordcycle value");
-  }
   meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_;
   const double *s = eps_data->s;
   const double *o = eps_data->o;

From 3bd5579f11419795d298fb49616d432e075ca4b2 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Tue, 10 Mar 2020 20:18:49 -0600
Subject: [PATCH 07/14] Rotations in place, added to header

---
 src/meep.hpp |   6 ++
 src/mpb.cpp  | 154 +++++++++++++++++++++++++++++++++++++--------------
 2 files changed, 117 insertions(+), 43 deletions(-)

diff --git a/src/meep.hpp b/src/meep.hpp
index 4d95e8289..4a3da6857 100644
--- a/src/meep.hpp
+++ b/src/meep.hpp
@@ -1662,6 +1662,11 @@ class fields {
   // that can be passed to eigenmode_amplitude() to get
   // values of field components at arbitrary points in space.
   // call destroy_eigenmode_data() to deallocate it when finished.
+  void *get_eigenmode_coordcycle(double omega_src, direction d, const volume where, const volume eig_vol,
+                      int band_num, const vec &kpoint, bool match_frequency, int parity,
+                      double resolution, double eigensolver_tol, double *kdom = 0,
+                      void **user_mdata = 0, int coordcycle = 0);
+
   void *get_eigenmode(double omega_src, direction d, const volume where, const volume eig_vol,
                       int band_num, const vec &kpoint, bool match_frequency, int parity,
                       double resolution, double eigensolver_tol, double *kdom = 0,
@@ -1982,6 +1987,7 @@ void green3d(std::complex<double> *EH, const vec &x, double freq, double eps, do
 void destroy_eigenmode_data(void *vedata, bool destroy_mdata = true);
 std::complex<double> eigenmode_amplitude(void *vedata, const vec &p, component c);
 double get_group_velocity(void *vedata);
+vec get_k_coordcycle(void *vedata, int coordcycle = 0);
 vec get_k(void *vedata);
 
 realnum linear_interpolate(realnum rx, realnum ry, realnum rz, realnum *data, int nx, int ny,
diff --git a/src/mpb.cpp b/src/mpb.cpp
index d78ab6560..090c7c1a3 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -43,12 +43,6 @@ typedef struct {
   const fields *f;
 } meep_mpb_eps_data;
 
-static void meep_mpb_eps(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
-                         void *eps_data_) {
-  int coordcycle = 0;
-  return meep_mpb_eps_coordcycle(eps, eps_inv, r, eps_data_, coordcycle);
-}
-
 static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
                          void *eps_data_, int coordcycle = 0) {
   meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_;
@@ -60,7 +54,7 @@ static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps
                                                    /* D1 */ vec(o[2] + r[2] * s[2])));
   const fields *f = eps_data->f;
 
-  // operations only need to happen in chi1inv
+  // TODO: rotate meep -> mpb
   switch(coordcycle) {
     case 0:
       eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
@@ -71,20 +65,20 @@ static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps
       ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, omega)), 0);
       break;
     case 1:
-      eps_inv->m00 = real(f->get_chi1inv(Ex, Y, p, omega));
-      eps_inv->m11 = real(f->get_chi1inv(Ey, Z, p, omega));
-      eps_inv->m22 = real(f->get_chi1inv(Ez, X, p, omega));
-      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Z, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, X, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, X, p, omega)), 0);
+      eps_inv->m00 = real(f->get_chi1inv(Ey, Y, p, omega));
+      eps_inv->m11 = real(f->get_chi1inv(Ez, Z, p, omega));
+      eps_inv->m22 = real(f->get_chi1inv(Ex, X, p, omega));
+      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ey, Z, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ey, X, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ez, X, p, omega)), 0);
       break;
     case 2:
-      eps_inv->m00 = real(f->get_chi1inv(Ex, Z, p, omega));
-      eps_inv->m11 = real(f->get_chi1inv(Ey, X, p, omega));
-      eps_inv->m22 = real(f->get_chi1inv(Ez, Y, p, omega));
-      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, X, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Y, p, omega)), 0);
+      eps_inv->m00 = real(f->get_chi1inv(Ez, Z, p, omega));
+      eps_inv->m11 = real(f->get_chi1inv(Ex, X, p, omega));
+      eps_inv->m22 = real(f->get_chi1inv(Ey, Y, p, omega));
+      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ez, X, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ez, Y, p, omega)), 0);
+      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
       break;
     default: abort("unsupported coordcycle value");
   }
@@ -100,6 +94,12 @@ static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps
   maxwell_sym_matrix_invert(eps, eps_inv);
 }
 
+static void meep_mpb_eps(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
+                         void *eps_data_) {
+  int coordcycle = 0;
+  return meep_mpb_eps_coordcycle(eps, eps_inv, r, eps_data_, coordcycle);
+}
+
 /**************************************************************/
 /* prototype for position-dependent amplitude function passed */
 /* to add_volume_source                                       */
@@ -279,14 +279,7 @@ static double dot_product(const mpb_real a[3], const mpb_real b[3]) {
 /* field components at arbitrary points in space.               */
 /* call destroy_eigenmode_data() to deallocate when finished.   */
 /****************************************************************/
-void *fields::get_eigenmode(double omega_src, direction d, const volume where, const volume eig_vol,
-                            int band_num, const vec &_kpoint, bool match_frequency, int parity,
-                            double resolution, double eigensolver_tol, double *kdom,
-                            void **user_mdata) {
-  return get_eigenmode_coordcycle(omega_src, d, where, eig_vol, band_num, _kpoint, match_frequency, parity, resolution,
-                                  eigensolver_tol, kdom, user_mdata, 0);
-}
-
+// TODO: create get_eigenmode_coordcycle and wrap get_eigenmode around it
 void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volume where, const volume eig_vol,
                             int band_num, const vec &_kpoint, bool match_frequency, int parity,
                             double resolution, double eigensolver_tol, double *kdom,
@@ -332,9 +325,9 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
   if (!eig_vol.contains(where))
     abort("invalid grid_volume in add_eigenmode_source (WHERE must be in EIG_VOL)");
 
-  // TODO: change here
   switch (gv.dim) {
     case D3:
+      // TODO: rotate  meep -> mpb
       switch (coordcycle) {
         case 0:
           o[0] = eig_vol.in_direction_min(X);
@@ -456,9 +449,22 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
       kmatch = kcart_len;
   }
   else {
-    // TODO: rotate
-    kmatch = G[d - X][d - X] * k[d - X]; // k[d] in cartesian
-    kdir[d - X] = 1;                     // kdir = unit vector in d direction
+    // TODO: rotate meep -> mpb
+    switch(coordcycle) {
+      case 0:
+        kmatch = G[d - X][d - X] * k[d - X]; // k[d] in cartesian
+        kdir[d - X] = 1;                     // kdir = unit vector in d direction
+        break;
+      case 1:
+        kmatch = G[((d + 1) % 3) - X][((d + 1) % 3) - X] * k[((d + 1) % 3) - X];
+        kdir[((d + 1) % 3) - X] = 1;
+        break;
+      case 2:
+        kmatch = G[((d + 2) % 3) - X][((d + 2) % 3) - X] * k[((d + 2) % 3) - X];
+        kdir[((d + 2) % 3) - X] = 1;
+        break;
+      default: abort("unsupported coordcycle value");
+    }
   }
 
   // if match_frequency is true, we need at least a crude guess for kmatch;
@@ -472,18 +478,46 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
       if (gv.dim == D2) k[2] = beta;
     }
     else {
-      // TODO: rotate
-      k[d - X] = kmatch * R[d - X][d - X]; // convert to reciprocal basis
-      if (eig_vol.in_direction(d) > 0 &&
-          fabs(k[d - X]) > 0.4) // ensure k is well inside the Brillouin zone
-        k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;
-      // end rotate to do
+      // TODO: rotate meep -> mpb
+      switch(coordcycle) {
+        case 0:
+          k[d - X] = kmatch * R[d - X][d - X]; // convert to reciprocal basis
+          if (eig_vol.in_direction(d) > 0 &&
+            fabs(k[d - X]) > 0.4) // ensure k is well inside the Brillouin zone
+          k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;
+          break;
+        case 1:
+          k[((d + 1) % 3) - X] = kmatch * R[((d + 1) % 3) - X][((d + 1) % 3) - X];
+          if (eig_vol.in_direction(((d + 1) % 3)) > 0 &&
+            fabs(k[((d + 1) % 3) - X]) > 0.4)
+          k[((d + 1) % 3) - X] = k[((d + 1) % 3) - X] > 0 ? 0.4 : -0.4;
+          break;
+        case 2:
+          k[((d + 2) % 3) - X] = kmatch * R[((d + 2) % 3) - X][((d + 2) % 3) - X];
+          if (eig_vol.in_direction(((d + 2) % 3)) > 0 &&
+            fabs(k[((d + 2) % 3) - X]) > 0.4)
+          k[((d + 2) % 3) - X] = k[((d + 2) % 3) - X] > 0 ? 0.4 : -0.4;
+          break;
+        default: abort("unsupported coordcycle value");
+      }
+
     }
     if (verbosity > 1) master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
   }
 
   // TODO: rotate
-  set_maxwell_data_parity(mdata, parity);
+  switch(coordcycle) {
+    case 0:
+      set_maxwell_data_parity(mdata, parity);
+      break;
+    case 1:
+      set_maxwell_data_parity(mdata, (parity + 1) % 3);
+      break;
+    case 2:
+      set_maxwell_data_parity(mdata, (parity + 2) % 3);
+      break;
+    default: abort("unsupported coordcycle value");
+  }
   update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);
 
   if (k[0] == 0 && k[1] == 0 && k[2] == 0) {
@@ -569,8 +603,19 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
           if (gv.dim == D2) k[2] = beta;
         }
         else {
-          // TODO: rotate
-          k[d - X] = kmatch * R[d - X][d - X];
+          // TODO: rotate unsure which direction, currently meep -> mpb
+          switch(coordcycle) {
+            case 0:
+              k[d - X] = kmatch * R[d - X][d - X];
+              break;
+            case 1:
+              k[((d + 1) % 3) - X] = kmatch * R[((d + 1) % 3) - X][((d + 1) % 3) - X];
+              break;
+            case 2:
+              k[((d + 2) % 3) - X] = kmatch * R[((d + 2) % 3) - X][((d + 2) % 3) - X];
+              break;
+            default: abort("unsupported coordcycle value");
+          }
         }
         update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);
       }
@@ -696,6 +741,14 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
   return (void *)edata;
 }
 
+void *fields::get_eigenmode(double omega_src, direction d, const volume where, const volume eig_vol,
+                            int band_num, const vec &_kpoint, bool match_frequency, int parity,
+                            double resolution, double eigensolver_tol, double *kdom,
+                            void **user_mdata) {
+  return get_eigenmode_coordcycle(omega_src, d, where, eig_vol, band_num, _kpoint, match_frequency, parity, resolution,
+                                  eigensolver_tol, kdom, user_mdata, 0);
+}
+
 void destroy_eigenmode_data(void *vedata, bool destroy_mdata) {
   eigenmode_data *edata = (eigenmode_data *)vedata;
   destroy_evectmatrix(edata->H);
@@ -709,10 +762,25 @@ double get_group_velocity(void *vedata) {
   return edata->group_velocity;
 }
 
-vec get_k(void *vedata) {
+vec get_k_coordcycle(void *vedata, int coordcycle = 0) {
   eigenmode_data *edata = (eigenmode_data *)vedata;
-  // TODO: rotate
-  return vec(edata->Gk[0], edata->Gk[1], edata->Gk[2]);
+  // TODO: rotate mpb -> meep
+  switch(coordcycle) {
+    case 0:
+      return vec(edata->Gk[0], edata->Gk[1], edata->Gk[2]);
+      break;
+    case 1:
+      return vec(edata->Gk[2], edata->Gk[0], edata->Gk[1]);
+      break;
+    case 2:
+      return vec(edata->Gk[1], edata->Gk[2], edata->Gk[0]);
+      break;
+    default: abort("unsupported coordcycle value");
+  }
+}
+
+vec get_k(void *vedata) {
+  return get_k_coordcycle(vedata, 0);
 }
 
 /***************************************************************/

From 18d4e8416576d6fa12aa48f7bb90816784b0a3cd Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Wed, 11 Mar 2020 19:11:15 -0600
Subject: [PATCH 08/14] attempt to cast ints as direction

---
 src/mpb.cpp | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 090c7c1a3..500b62529 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -456,12 +456,12 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
         kdir[d - X] = 1;                     // kdir = unit vector in d direction
         break;
       case 1:
-        kmatch = G[((d + 1) % 3) - X][((d + 1) % 3) - X] * k[((d + 1) % 3) - X];
-        kdir[((d + 1) % 3) - X] = 1;
+        kmatch = G[((d + (direction)1) % (direction)3) - X][((d + (direction)1) % (direction)3) - X] * k[((d + (direction)1) % (direction)3) - X];
+        kdir[((d + (direction)1) % (direction)3) - X] = 1;
         break;
       case 2:
-        kmatch = G[((d + 2) % 3) - X][((d + 2) % 3) - X] * k[((d + 2) % 3) - X];
-        kdir[((d + 2) % 3) - X] = 1;
+        kmatch = G[((d + (direction)2) % (direction)3) - X][((d + (direction)2) % (direction)3) - X] * k[((d + (direction)2) % (direction)3) - X];
+        kdir[((d + (direction)2) % (direction)3) - X] = 1;
         break;
       default: abort("unsupported coordcycle value");
     }
@@ -487,16 +487,16 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
           k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;
           break;
         case 1:
-          k[((d + 1) % 3) - X] = kmatch * R[((d + 1) % 3) - X][((d + 1) % 3) - X];
-          if (eig_vol.in_direction(((d + 1) % 3)) > 0 &&
-            fabs(k[((d + 1) % 3) - X]) > 0.4)
-          k[((d + 1) % 3) - X] = k[((d + 1) % 3) - X] > 0 ? 0.4 : -0.4;
+          k[((d + (direction)1) % (direction)3) - X] = kmatch * R[((d + (direction)1) % (direction)3) - X][((d + (direction)1) % (direction)3) - X];
+          if (eig_vol.in_direction(((d + (direction)1) % (direction)3)) > 0 &&
+            fabs(k[((d + (direction)1) % (direction)3) - X]) > 0.4)
+          k[((d + (direction)1) % (direction)3) - X] = k[((d + (direction)1) % (direction)3) - X] > 0 ? 0.4 : -0.4;
           break;
         case 2:
-          k[((d + 2) % 3) - X] = kmatch * R[((d + 2) % 3) - X][((d + 2) % 3) - X];
-          if (eig_vol.in_direction(((d + 2) % 3)) > 0 &&
-            fabs(k[((d + 2) % 3) - X]) > 0.4)
-          k[((d + 2) % 3) - X] = k[((d + 2) % 3) - X] > 0 ? 0.4 : -0.4;
+          k[((d + (direction)2) % (direction)3) - X] = kmatch * R[((d + (direction)2) % (direction)3) - X][((d + (direction)2) % (direction)3) - X];
+          if (eig_vol.in_direction(((d + (direction)2) % (direction)3)) > 0 &&
+            fabs(k[((d + (direction)2) % (direction)3) - X]) > 0.4)
+          k[((d + (direction)2) % (direction)3) - X] = k[((d + (direction)2) % (direction)3) - X] > 0 ? 0.4 : -0.4;
           break;
         default: abort("unsupported coordcycle value");
       }
@@ -609,10 +609,10 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
               k[d - X] = kmatch * R[d - X][d - X];
               break;
             case 1:
-              k[((d + 1) % 3) - X] = kmatch * R[((d + 1) % 3) - X][((d + 1) % 3) - X];
+              k[((d + (direction)1) % (direction)3) - X] = kmatch * R[((d + (direction)1) % (direction)3) - X][((d + (direction)1) % (direction)3) - X];
               break;
             case 2:
-              k[((d + 2) % 3) - X] = kmatch * R[((d + 2) % 3) - X][((d + 2) % 3) - X];
+              k[((d + (direction)2) % (direction)3) - X] = kmatch * R[((d + (direction)2) % (direction)3) - X][((d + (direction)2) % (direction)3) - X];
               break;
             default: abort("unsupported coordcycle value");
           }

From 77773fc9c28c8ba600a2ecee1e77e98a2bd59dc5 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Wed, 11 Mar 2020 19:19:37 -0600
Subject: [PATCH 09/14] finally found definition of direction, it's an enum, so
 changed d-X values to by d-Y or d-Z based on coordcycle value

---
 src/mpb.cpp | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 500b62529..c228bcf76 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -456,12 +456,12 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
         kdir[d - X] = 1;                     // kdir = unit vector in d direction
         break;
       case 1:
-        kmatch = G[((d + (direction)1) % (direction)3) - X][((d + (direction)1) % (direction)3) - X] * k[((d + (direction)1) % (direction)3) - X];
-        kdir[((d + (direction)1) % (direction)3) - X] = 1;
+        kmatch = G[d - Y][d - Y] * k[d - Y];
+        kdir[d - Y] = 1;
         break;
       case 2:
-        kmatch = G[((d + (direction)2) % (direction)3) - X][((d + (direction)2) % (direction)3) - X] * k[((d + (direction)2) % (direction)3) - X];
-        kdir[((d + (direction)2) % (direction)3) - X] = 1;
+        kmatch = G[d - Z][d - Z] * k[d - Z];
+        kdir[d - Z] = 1;
         break;
       default: abort("unsupported coordcycle value");
     }
@@ -487,16 +487,16 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
           k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;
           break;
         case 1:
-          k[((d + (direction)1) % (direction)3) - X] = kmatch * R[((d + (direction)1) % (direction)3) - X][((d + (direction)1) % (direction)3) - X];
-          if (eig_vol.in_direction(((d + (direction)1) % (direction)3)) > 0 &&
-            fabs(k[((d + (direction)1) % (direction)3) - X]) > 0.4)
-          k[((d + (direction)1) % (direction)3) - X] = k[((d + (direction)1) % (direction)3) - X] > 0 ? 0.4 : -0.4;
+          k[d - Y] = kmatch * R[d - Y][d - Y];
+          if (eig_vol.in_direction(d) > 0 &&
+            fabs(k[d - Y]) > 0.4)
+          k[d - Y] = k[d - Y] > 0 ? 0.4 : -0.4;
           break;
         case 2:
-          k[((d + (direction)2) % (direction)3) - X] = kmatch * R[((d + (direction)2) % (direction)3) - X][((d + (direction)2) % (direction)3) - X];
-          if (eig_vol.in_direction(((d + (direction)2) % (direction)3)) > 0 &&
-            fabs(k[((d + (direction)2) % (direction)3) - X]) > 0.4)
-          k[((d + (direction)2) % (direction)3) - X] = k[((d + (direction)2) % (direction)3) - X] > 0 ? 0.4 : -0.4;
+          k[d - Z] = kmatch * R[d - Z][d - Z];
+          if (eig_vol.in_direction(d) > 0 &&
+            fabs(k[d - Z]) > 0.4)
+          k[d - Z] = k[d - Z] > 0 ? 0.4 : -0.4;
           break;
         default: abort("unsupported coordcycle value");
       }
@@ -609,10 +609,10 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
               k[d - X] = kmatch * R[d - X][d - X];
               break;
             case 1:
-              k[((d + (direction)1) % (direction)3) - X] = kmatch * R[((d + (direction)1) % (direction)3) - X][((d + (direction)1) % (direction)3) - X];
+              k[d - Y] = kmatch * R[d - Y][d - Y];
               break;
             case 2:
-              k[((d + (direction)2) % (direction)3) - X] = kmatch * R[((d + (direction)2) % (direction)3) - X][((d + (direction)2) % (direction)3) - X];
+              k[d - Z] = kmatch * R[d - Z][d - Z];
               break;
             default: abort("unsupported coordcycle value");
           }

From a245065be13699efcdd71641c029c6e6d9b085d9 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Wed, 11 Mar 2020 19:28:39 -0600
Subject: [PATCH 10/14] removed default parameter from mpb.cpp

---
 src/mpb.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index c228bcf76..720941306 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -762,7 +762,7 @@ double get_group_velocity(void *vedata) {
   return edata->group_velocity;
 }
 
-vec get_k_coordcycle(void *vedata, int coordcycle = 0) {
+vec get_k_coordcycle(void *vedata, int coordcycle) {
   eigenmode_data *edata = (eigenmode_data *)vedata;
   // TODO: rotate mpb -> meep
   switch(coordcycle) {

From c9c36281efcad0fb933c5e7accf62a2fa66c6093 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Wed, 11 Mar 2020 20:02:31 -0600
Subject: [PATCH 11/14] cleanup parity

---
 src/mpb.cpp | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 720941306..0320840e5 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -505,16 +505,16 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
     if (verbosity > 1) master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
   }
 
-  // TODO: rotate
+  // TODO: rotate, not sure how to do this on parity
   switch(coordcycle) {
     case 0:
       set_maxwell_data_parity(mdata, parity);
       break;
     case 1:
-      set_maxwell_data_parity(mdata, (parity + 1) % 3);
+      set_maxwell_data_parity(mdata, parity);
       break;
     case 2:
-      set_maxwell_data_parity(mdata, (parity + 2) % 3);
+      set_maxwell_data_parity(mdata, parity);
       break;
     default: abort("unsupported coordcycle value");
   }

From e49c9a44c148ca165c6053b0c3c572744f8449d4 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Thu, 12 Mar 2020 01:28:54 -0600
Subject: [PATCH 12/14] Restructure coordcycle rotations; direction and
 component casting and recasting should be in order. Need to rotate parity
 bitfield

---
 src/mpb.cpp | 142 ++++++++++++----------------------------------------
 1 file changed, 31 insertions(+), 111 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 0320840e5..169e1eed9 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -45,6 +45,9 @@ typedef struct {
 
 static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps_inv, const mpb_real r[3],
                          void *eps_data_, int coordcycle = 0) {
+  if (coordcycle != 0 && coordcycle != 1 && coordcycle != 2) {
+    abort("unsupported coordcycle value");
+  }
   meep_mpb_eps_data *eps_data = (meep_mpb_eps_data *)eps_data_;
   const double *s = eps_data->s;
   const double *o = eps_data->o;
@@ -55,33 +58,12 @@ static void meep_mpb_eps_coordcycle(symmetric_matrix *eps, symmetric_matrix *eps
   const fields *f = eps_data->f;
 
   // TODO: rotate meep -> mpb
-  switch(coordcycle) {
-    case 0:
-      eps_inv->m00 = real(f->get_chi1inv(Ex, X, p, omega));
-      eps_inv->m11 = real(f->get_chi1inv(Ey, Y, p, omega));
-      eps_inv->m22 = real(f->get_chi1inv(Ez, Z, p, omega));
-      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ex, Z, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ey, Z, p, omega)), 0);
-      break;
-    case 1:
-      eps_inv->m00 = real(f->get_chi1inv(Ey, Y, p, omega));
-      eps_inv->m11 = real(f->get_chi1inv(Ez, Z, p, omega));
-      eps_inv->m22 = real(f->get_chi1inv(Ex, X, p, omega));
-      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ey, Z, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ey, X, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ez, X, p, omega)), 0);
-      break;
-    case 2:
-      eps_inv->m00 = real(f->get_chi1inv(Ez, Z, p, omega));
-      eps_inv->m11 = real(f->get_chi1inv(Ex, X, p, omega));
-      eps_inv->m22 = real(f->get_chi1inv(Ey, Y, p, omega));
-      ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv(Ez, X, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv(Ez, Y, p, omega)), 0);
-      ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv(Ex, Y, p, omega)), 0);
-      break;
-    default: abort("unsupported coordcycle value");
-  }
+  eps_inv->m00 = real(f->get_chi1inv((component)((Ex - Ex + coordcycle) % 3), (direction)((X - X + coordcycle) % 3), p, omega));
+  eps_inv->m11 = real(f->get_chi1inv((component)((Ey - Ex + coordcycle) % 3), (direction)((Y - X + coordcycle) % 3), p, omega));
+  eps_inv->m22 = real(f->get_chi1inv((component)((Ez - Ex + coordcycle) % 3), (direction)((Z - X + coordcycle) % 3), p, omega));
+  ASSIGN_ESCALAR(eps_inv->m01, real(f->get_chi1inv((component)((Ex - Ex + coordcycle) % 3), (direction)((Y - X + coordcycle) % 3), p, omega)), 0);
+  ASSIGN_ESCALAR(eps_inv->m02, real(f->get_chi1inv((component)((Ex - Ex + coordcycle) % 3), (direction)((Z - X + coordcycle) % 3), p, omega)), 0);
+  ASSIGN_ESCALAR(eps_inv->m12, real(f->get_chi1inv((component)((Ey - Ex + coordcycle) % 3), (direction)((Z - X + coordcycle) % 3), p, omega)), 0);
 
   /*
   master_printf("m11(%g,%g) = %g\n", p.x(), p.y(), eps_inv->m00);
@@ -279,6 +261,9 @@ static double dot_product(const mpb_real a[3], const mpb_real b[3]) {
 /* field components at arbitrary points in space.               */
 /* call destroy_eigenmode_data() to deallocate when finished.   */
 /****************************************************************/
+//#define EVEN_X_PARITY (1<<4)
+//#define ODD_X_PARITY (1<<5)
+
 // TODO: create get_eigenmode_coordcycle and wrap get_eigenmode around it
 void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volume where, const volume eig_vol,
                             int band_num, const vec &_kpoint, bool match_frequency, int parity,
@@ -287,6 +272,9 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
   /*--------------------------------------------------------------*/
   /*- part 1: preliminary setup for calling MPB  -----------------*/
   /*--------------------------------------------------------------*/
+  if (coordcycle != 0 && coordcycle != 1 && coordcycle != 2) {
+    abort("unsupported coordcycle value");
+  }
 
   // if the mode region extends over the full computational grid and we are bloch-periodic
   // in any direction, set the corresponding component of the eigenmode initial-guess
@@ -328,42 +316,15 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
   switch (gv.dim) {
     case D3:
       // TODO: rotate  meep -> mpb
-      switch (coordcycle) {
-        case 0:
-          o[0] = eig_vol.in_direction_min(X);
-          o[1] = eig_vol.in_direction_min(Y);
-          o[2] = eig_vol.in_direction_min(Z);
-          s[0] = eig_vol.in_direction(X);
-          s[1] = eig_vol.in_direction(Y);
-          s[2] = eig_vol.in_direction(Z);
-          kcart[0] = kpoint.in_direction(X);
-          kcart[1] = kpoint.in_direction(Y);
-          kcart[2] = kpoint.in_direction(Z);
-          break;
-        case 1:
-          o[0] = eig_vol.in_direction_min(Y);
-          o[1] = eig_vol.in_direction_min(Z);
-          o[2] = eig_vol.in_direction_min(X);
-          s[0] = eig_vol.in_direction(Y);
-          s[1] = eig_vol.in_direction(Z);
-          s[2] = eig_vol.in_direction(X);
-          kcart[0] = kpoint.in_direction(Y);
-          kcart[1] = kpoint.in_direction(Z);
-          kcart[2] = kpoint.in_direction(X);
-          break;
-        case 2:
-          o[0] = eig_vol.in_direction_min(Z);
-          o[1] = eig_vol.in_direction_min(X);
-          o[2] = eig_vol.in_direction_min(Y);
-          s[0] = eig_vol.in_direction(Z);
-          s[1] = eig_vol.in_direction(X);
-          s[2] = eig_vol.in_direction(Y);
-          kcart[0] = kpoint.in_direction(Z);
-          kcart[1] = kpoint.in_direction(X);
-          kcart[2] = kpoint.in_direction(Y);
-          break;
-        default: abort("unsupported coordcycle value");
-      }
+      o[0] = eig_vol.in_direction_min((direction)((X - X + coordcycle) % 3));
+      o[1] = eig_vol.in_direction_min((direction)((Y - X + coordcycle) % 3));
+      o[2] = eig_vol.in_direction_min((direction)((Z - X + coordcycle) % 3));
+      s[0] = eig_vol.in_direction((direction)((X - X + coordcycle) % 3));
+      s[1] = eig_vol.in_direction((direction)((Y - X + coordcycle) % 3));
+      s[2] = eig_vol.in_direction((direction)((Z - X + coordcycle) % 3));
+      kcart[0] = kpoint.in_direction((direction)((X - X + coordcycle) % 3));
+      kcart[1] = kpoint.in_direction((direction)((Y - X + coordcycle) % 3));
+      kcart[2] = kpoint.in_direction((direction)((Z - X + coordcycle) % 3));
       break;
     case D2:
       o[0] = eig_vol.in_direction_min(X);
@@ -450,21 +411,8 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
   }
   else {
     // TODO: rotate meep -> mpb
-    switch(coordcycle) {
-      case 0:
-        kmatch = G[d - X][d - X] * k[d - X]; // k[d] in cartesian
-        kdir[d - X] = 1;                     // kdir = unit vector in d direction
-        break;
-      case 1:
-        kmatch = G[d - Y][d - Y] * k[d - Y];
-        kdir[d - Y] = 1;
-        break;
-      case 2:
-        kmatch = G[d - Z][d - Z] * k[d - Z];
-        kdir[d - Z] = 1;
-        break;
-      default: abort("unsupported coordcycle value");
-    }
+    kmatch = G[(d - X + coordcycle) % 3][(d - X + coordcycle) % 3] * k[(d - X + coordcycle) % 3]; // k[d] in cartesian
+    kdir[(d - X + coordcycle) % 3] = 1;                                                           // kdir = unit vector in d direction
   }
 
   // if match_frequency is true, we need at least a crude guess for kmatch;
@@ -479,27 +427,10 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
     }
     else {
       // TODO: rotate meep -> mpb
-      switch(coordcycle) {
-        case 0:
-          k[d - X] = kmatch * R[d - X][d - X]; // convert to reciprocal basis
-          if (eig_vol.in_direction(d) > 0 &&
-            fabs(k[d - X]) > 0.4) // ensure k is well inside the Brillouin zone
-          k[d - X] = k[d - X] > 0 ? 0.4 : -0.4;
-          break;
-        case 1:
-          k[d - Y] = kmatch * R[d - Y][d - Y];
-          if (eig_vol.in_direction(d) > 0 &&
-            fabs(k[d - Y]) > 0.4)
-          k[d - Y] = k[d - Y] > 0 ? 0.4 : -0.4;
-          break;
-        case 2:
-          k[d - Z] = kmatch * R[d - Z][d - Z];
-          if (eig_vol.in_direction(d) > 0 &&
-            fabs(k[d - Z]) > 0.4)
-          k[d - Z] = k[d - Z] > 0 ? 0.4 : -0.4;
-          break;
-        default: abort("unsupported coordcycle value");
-      }
+      k[(d - X + coordcycle) % 3] = kmatch * R[(d - X + coordcycle) % 3][(d - X + coordcycle) % 3]; // convert to reciprocal basis
+          if (eig_vol.in_direction((direction)((d - X + coordcycle) % 3)) > 0 &&
+            fabs(k[(d - X + coordcycle) % 3]) > 0.4) // ensure k is well inside the Brillouin zone
+          k[(d - X + coordcycle) % 3] = k[(d - X + coordcycle) % 3] > 0 ? 0.4 : -0.4;
 
     }
     if (verbosity > 1) master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
@@ -604,18 +535,7 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
         }
         else {
           // TODO: rotate unsure which direction, currently meep -> mpb
-          switch(coordcycle) {
-            case 0:
-              k[d - X] = kmatch * R[d - X][d - X];
-              break;
-            case 1:
-              k[d - Y] = kmatch * R[d - Y][d - Y];
-              break;
-            case 2:
-              k[d - Z] = kmatch * R[d - Z][d - Z];
-              break;
-            default: abort("unsupported coordcycle value");
-          }
+          k[(d - X + coordcycle) % 3] = kmatch * R[(d - X + coordcycle) % 3][(d - X + coordcycle) % 3];
         }
         update_maxwell_data_k(mdata, k, G[0], G[1], G[2]);
       }

From 93d88cd4c3f99ada00c33ed3423618f7883c8ee2 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Mon, 16 Mar 2020 16:09:04 -0600
Subject: [PATCH 13/14] added parity rotation

---
 src/mpb.cpp | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 169e1eed9..20a265d31 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -261,8 +261,8 @@ static double dot_product(const mpb_real a[3], const mpb_real b[3]) {
 /* field components at arbitrary points in space.               */
 /* call destroy_eigenmode_data() to deallocate when finished.   */
 /****************************************************************/
-//#define EVEN_X_PARITY (1<<4)
-//#define ODD_X_PARITY (1<<5)
+#define EVEN_X_PARITY (1<<4)
+#define ODD_X_PARITY (1<<5)
 
 // TODO: create get_eigenmode_coordcycle and wrap get_eigenmode around it
 void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volume where, const volume eig_vol,
@@ -442,9 +442,17 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
       set_maxwell_data_parity(mdata, parity);
       break;
     case 1:
+      // Rotating x,y,z parity to y,z,x involves rotating the bits 2 bits to the left
+      unsigned int num_sig_bits = 6;
+      unsigned int rotation = 2;
+      parity = ((parity >> rotation)|(parity << (num_sig_bits - rotation))) % 64;
       set_maxwell_data_parity(mdata, parity);
       break;
     case 2:
+      // Rotating x,y,z parity to z,x,y involves rotating the bits 4 bits to the left
+      unsigned int num_sig_bits = 6;
+      unsigned int rotation = 4;
+      parity = ((parity >> rotation)|(parity << (num_sig_bits - rotation))) % 64;
       set_maxwell_data_parity(mdata, parity);
       break;
     default: abort("unsupported coordcycle value");

From ed1b446667b33cbf398ac249d4fd4c77e651f4e5 Mon Sep 17 00:00:00 2001
From: danielwboyce <danielwboyce@gmail.com>
Date: Mon, 16 Mar 2020 16:48:05 -0600
Subject: [PATCH 14/14] couple typos

---
 src/mpb.cpp | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/mpb.cpp b/src/mpb.cpp
index 20a265d31..d9c4c4173 100644
--- a/src/mpb.cpp
+++ b/src/mpb.cpp
@@ -436,22 +436,22 @@ void *fields::get_eigenmode_coordcycle(double omega_src, direction d, const volu
     if (verbosity > 1) master_printf("NEW KPOINT: %g, %g, %g\n", k[0], k[1], k[2]);
   }
 
-  // TODO: rotate, not sure how to do this on parity
+  // TODO: rotate meep -> mpb
+  unsigned int num_sig_bits = 6;
+  unsigned int rotation;
   switch(coordcycle) {
     case 0:
       set_maxwell_data_parity(mdata, parity);
       break;
     case 1:
       // Rotating x,y,z parity to y,z,x involves rotating the bits 2 bits to the left
-      unsigned int num_sig_bits = 6;
-      unsigned int rotation = 2;
+      rotation = 2;
       parity = ((parity >> rotation)|(parity << (num_sig_bits - rotation))) % 64;
       set_maxwell_data_parity(mdata, parity);
       break;
     case 2:
       // Rotating x,y,z parity to z,x,y involves rotating the bits 4 bits to the left
-      unsigned int num_sig_bits = 6;
-      unsigned int rotation = 4;
+      rotation = 4;
       parity = ((parity >> rotation)|(parity << (num_sig_bits - rotation))) % 64;
       set_maxwell_data_parity(mdata, parity);
       break;