diff --git a/libpympb/pympb.cpp b/libpympb/pympb.cpp
index b758b4102..e3d03b3de 100644
--- a/libpympb/pympb.cpp
+++ b/libpympb/pympb.cpp
@@ -29,6 +29,18 @@
           int xyz_index = ((i2_ * n1 + i1) * n3 + i3);
 #  endif /* HAVE_MPI */
 
+#ifdef CASSIGN_MULT
+  #undef CASSIGN_MULT
+#endif
+
+/* a = b * c */
+#define CASSIGN_MULT(a, b, c) { \
+  mpb_real bbbb_re = (b).re, bbbb_im = (b).im; \
+  mpb_real cccc_re = (c).re, cccc_im = (c).im; \
+  CASSIGN_SCALAR(a, bbbb_re * cccc_re - bbbb_im * cccc_im, \
+                 bbbb_re * cccc_im + bbbb_im * cccc_re); \
+}
+
 // TODO: Support MPI
 #define mpi_allreduce(sb, rb, n, ctype, t, op, comm) { \
      CHECK((sb) != (rb), "MPI_Allreduce doesn't work for sendbuf == recvbuf");\
@@ -87,18 +99,20 @@ static int mean_epsilon_func(symmetric_matrix* meps, symmetric_matrix *meps_inv,
 /* a couple of utilities to convert libctl data types to the data
    types of the eigensolver & maxwell routines: */
 
-void vector3_to_arr(mpb_real arr[3], vector3 v)
-{
-     arr[0] = v.x;
-     arr[1] = v.y;
-     arr[2] = v.z;
+void vector3_to_arr(mpb_real arr[3], vector3 v) {
+  arr[0] = v.x;
+  arr[1] = v.y;
+  arr[2] = v.z;
 }
 
-void matrix3x3_to_arr(mpb_real arr[3][3], matrix3x3 m)
-{
-     vector3_to_arr(arr[0], m.c0);
-     vector3_to_arr(arr[1], m.c1);
-     vector3_to_arr(arr[2], m.c2);
+void matrix3x3_to_arr(mpb_real arr[3][3], matrix3x3 m) {
+  vector3_to_arr(arr[0], m.c0);
+  vector3_to_arr(arr[1], m.c1);
+  vector3_to_arr(arr[2], m.c2);
+}
+
+cnumber cscalar2cnumber(scalar_complex cs) {
+  return make_cnumber(CSCALAR_RE(cs), CSCALAR_IM(cs));
 }
 
 // Return a string describing the current parity, used for frequency and filename
@@ -1202,6 +1216,69 @@ void mode_solver::get_epsilon() {
                       eps_high == eps_low ? 100.0 : 100.0 * (eps_mean-eps_low) / (eps_high-eps_low));
 }
 
+/* get the mu function, and compute some statistics */
+void mode_solver::get_mu() {
+  mpb_real eps_mean = 0;
+  mpb_real mu_inv_mean = 0;
+  mpb_real eps_high = -1e20;
+  mpb_real eps_low = 1e20;
+  int fill_count = 0;
+
+  if (!mdata) {
+    meep::master_fprintf(stderr, "mode_solver.init must be called before get-mu!\n");
+    return;
+  }
+
+  curfield = (scalar_complex *) mdata->fft_data;
+  mpb_real *mu = (mpb_real *) curfield;
+  curfield_band = 0;
+  curfield_type = mu_CURFIELD_TYPE;
+
+  /* get mu.  Recall that we actually have an inverse
+     dielectric tensor at each point; define an average index by
+     the inverse of the average eigenvalue of the 1/eps tensor.
+     i.e. 3/(trace 1/eps). */
+
+  int N = mdata->fft_output_size;
+
+  for (int i = 0; i < N; ++i) {
+    if (mdata->mu_inv == NULL) {
+      mu[i] = 1.0;
+    }
+    else {
+      mu[i] = mean_medium_from_matrix(mdata->mu_inv + i);
+    }
+
+    if (mu[i] < eps_low) {
+      eps_low = mu[i];
+    }
+    if (mu[i] > eps_high) {
+      eps_high = mu[i];
+    }
+
+    eps_mean += mu[i];
+    mu_inv_mean += 1/mu[i];
+    if (mu[i] > 1.0001) {
+      ++fill_count;
+    }
+  }
+
+  mpi_allreduce_1(&eps_mean, mpb_real, SCALAR_MPI_TYPE, MPI_SUM, mpb_comm);
+  mpi_allreduce_1(&mu_inv_mean, mpb_real, SCALAR_MPI_TYPE, MPI_SUM, mpb_comm);
+  mpi_allreduce_1(&eps_low, mpb_real, SCALAR_MPI_TYPE, MPI_MIN, mpb_comm);
+  mpi_allreduce_1(&eps_high, mpb_real, SCALAR_MPI_TYPE, MPI_MAX, mpb_comm);
+  mpi_allreduce_1(&fill_count, int, MPI_INT, MPI_SUM, mpb_comm);
+
+  N = mdata->nx * mdata->ny * mdata->nz;
+  eps_mean /= N;
+  mu_inv_mean = N/mu_inv_mean;
+
+  meep::master_printf("mu: %g-%g, mean %g, harm. mean %g, %g%% > 1, %g%% \"fill\"\n",
+                      eps_low, eps_high, eps_mean, mu_inv_mean, (100.0 * fill_count) / N,
+                      eps_high == eps_low ? 100.0 : 100.0 *
+                      (eps_mean-eps_low) / (eps_high-eps_low));
+}
+
 void mode_solver::curfield_reset() {
   curfield = NULL;
   curfield_type = '-';
@@ -1595,6 +1672,168 @@ std::vector<mpb_real> mode_solver::get_output_k() {
   return output_k;
 }
 
+mpb_real mode_solver::get_val(int ix, int iy, int iz, int nx, int ny, int nz,
+                              int last_dim_size, mpb_real *data, int stride, int conjugate) {
+// #ifdef HAVE_MPI
+//      CHECK(0, "get-*-point not yet implemented for MPI!");
+// #else
+  (void)nx;
+  (void)last_dim_size;
+  (void)conjugate;
+  return data[(((ix * ny) + iy) * nz + iz) * stride];
+// #endif
+}
+
+mpb_real mode_solver::interp_val(vector3 p, int nx, int ny, int nz, int last_dim_size,
+                                 mpb_real *data, int stride, int conjugate) {
+  double ipart;
+  mpb_real rx, ry, rz, dx, dy, dz;
+  int x, y, z, x2, y2, z2;
+
+  mpb_real latx = geometry_lattice.size.x == 0 ? 1e-20 : geometry_lattice.size.x;
+  mpb_real laty = geometry_lattice.size.y == 0 ? 1e-20 : geometry_lattice.size.y;
+  mpb_real latz = geometry_lattice.size.z == 0 ? 1e-20 : geometry_lattice.size.z;
+
+  rx = modf(p.x/latx + 0.5, &ipart); if (rx < 0) rx += 1;
+  ry = modf(p.y/laty + 0.5, &ipart); if (ry < 0) ry += 1;
+  rz = modf(p.z/latz + 0.5, &ipart); if (rz < 0) rz += 1;
+
+  /* get the point corresponding to r in the grid: */
+  x = rx * nx;
+  y = ry * ny;
+  z = rz * nz;
+
+  /* get the difference between (x,y,z) and the actual point */
+  dx = rx * nx - x;
+  dy = ry * ny - y;
+  dz = rz * nz - z;
+
+  /* get the other closest point in the grid, with periodic boundaries: */
+  x2 = (nx + (dx >= 0.0 ? x + 1 : x - 1)) % nx;
+  y2 = (ny + (dy >= 0.0 ? y + 1 : y - 1)) % ny;
+  z2 = (nz + (dz >= 0.0 ? z + 1 : z - 1)) % nz;
+
+  /* take abs(d{xyz}) to get weights for {xyz} and {xyz}2: */
+  dx = fabs(dx);
+  dy = fabs(dy);
+  dz = fabs(dz);
+
+#define D(x,y,z) (get_val(x,y,z,nx,ny,nz,last_dim_size, data,stride,conjugate))
+
+  return(((D(x,y,z)*(1.0-dx) + D(x2,y,z)*dx) * (1.0-dy) +
+          (D(x,y2,z)*(1.0-dx) + D(x2,y2,z)*dx) * dy) * (1.0-dz) +
+         ((D(x,y,z2)*(1.0-dx) + D(x2,y,z2)*dx) * (1.0-dy) +
+          (D(x,y2,z2)*(1.0-dx) + D(x2,y2,z2)*dx) * dy) * dz);
+
+#undef D
+}
+
+scalar_complex mode_solver::interp_cval(vector3 p, int nx, int ny, int nz,
+                                        int last_dim_size, mpb_real *data, int stride) {
+  scalar_complex cval;
+  cval.re = interp_val(p, nx,ny,nz,last_dim_size, data, stride, 0);
+  cval.im = interp_val(p, nx,ny,nz,last_dim_size,data + 1, stride, 1);
+  return cval;
+}
+
+#define f_interp_val(p,f,data,stride,conj) interp_val(p,f->nx,f->ny,f->nz,f->last_dim_size,data,stride,conj)
+#define f_interp_cval(p,f,data,stride) interp_cval(p,f->nx,f->ny,f->nz,f->last_dim_size,data,stride)
+
+symmetric_matrix mode_solver::interp_eps_inv(vector3 p) {
+  int stride = sizeof(symmetric_matrix) / sizeof(mpb_real);
+  symmetric_matrix eps_inv;
+
+  eps_inv.m00 = f_interp_val(p, mdata, &mdata->eps_inv->m00, stride, 0);
+  eps_inv.m11 = f_interp_val(p, mdata, &mdata->eps_inv->m11, stride, 0);
+  eps_inv.m22 = f_interp_val(p, mdata, &mdata->eps_inv->m22, stride, 0);
+#ifdef WITH_HERMITIAN_EPSILON
+  eps_inv.m01 = f_interp_cval(p, mdata, &mdata->eps_inv->m01.re, stride);
+  eps_inv.m02 = f_interp_cval(p, mdata, &mdata->eps_inv->m02.re, stride);
+  eps_inv.m12 = f_interp_cval(p, mdata, &mdata->eps_inv->m12.re, stride);
+#else
+  eps_inv.m01 = f_interp_val(p, mdata, &mdata->eps_inv->m01, stride, 0);
+  eps_inv.m02 = f_interp_val(p, mdata, &mdata->eps_inv->m02, stride, 0);
+  eps_inv.m12 = f_interp_val(p, mdata, &mdata->eps_inv->m12, stride, 0);
+#endif
+  return eps_inv;
+}
+
+mpb_real mode_solver::get_epsilon_point(vector3 p) {
+  symmetric_matrix eps_inv;
+  eps_inv = interp_eps_inv(p);
+  return mean_medium_from_matrix(&eps_inv);
+}
+
+cmatrix3x3 mode_solver::get_epsilon_inverse_tensor_point(vector3 p) {
+  symmetric_matrix eps_inv;
+  eps_inv = interp_eps_inv(p);
+
+#ifdef WITH_HERMITIAN_EPSILON
+  return make_hermitian_cmatrix3x3(eps_inv.m00,eps_inv.m11,eps_inv.m22,
+                                   cscalar2cnumber(eps_inv.m01),
+                                   cscalar2cnumber(eps_inv.m02),
+                                   cscalar2cnumber(eps_inv.m12));
+#else
+  return make_hermitian_cmatrix3x3(eps_inv.m00,eps_inv.m11,eps_inv.m22,
+                                   make_cnumber(eps_inv.m01,0),
+                                   make_cnumber(eps_inv.m02,0),
+                                   make_cnumber(eps_inv.m12,0));
+#endif
+}
+
+mpb_real mode_solver::get_energy_point(vector3 p) {
+  CHECK(curfield && strchr("DHBR", curfield_type),
+        "compute-field-energy must be called before get-energy-point");
+  return f_interp_val(p, mdata, (mpb_real *) curfield, 1, 0);
+}
+
+cvector3 mode_solver::get_bloch_field_point(vector3 p) {
+  scalar_complex field[3];
+  cvector3 F;
+
+  CHECK(curfield && strchr("dhbecv", curfield_type),
+        "field must be must be loaded before get-*field*-point");
+  field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6);
+  field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6);
+  field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6);
+  F.x = cscalar2cnumber(field[0]);
+  F.y = cscalar2cnumber(field[1]);
+  F.z = cscalar2cnumber(field[2]);
+  return F;
+}
+
+cvector3 mode_solver::get_field_point(vector3 p) {
+  scalar_complex field[3], phase;
+  cvector3 F;
+
+  CHECK(curfield && strchr("dhbecv", curfield_type),
+       "field must be must be loaded before get-*field*-point");
+  field[0] = f_interp_cval(p, mdata, &curfield[0].re, 6);
+  field[1] = f_interp_cval(p, mdata, &curfield[1].re, 6);
+  field[2] = f_interp_cval(p, mdata, &curfield[2].re, 6);
+
+  if (curfield_type != 'v') {
+    mpb_real latx = geometry_lattice.size.x == 0 ? 1e-20 : geometry_lattice.size.x;
+    mpb_real laty = geometry_lattice.size.y == 0 ? 1e-20 : geometry_lattice.size.y;
+    mpb_real latz = geometry_lattice.size.z == 0 ? 1e-20 : geometry_lattice.size.z;
+
+    double phase_phi = TWOPI * (cur_kvector.x * (p.x / latx) +
+                                cur_kvector.y * (p.y / laty) +
+                                cur_kvector.z * (p.z / latz));
+
+    CASSIGN_SCALAR(phase, cos(phase_phi), sin(phase_phi));
+    CASSIGN_MULT(field[0], field[0], phase);
+    CASSIGN_MULT(field[1], field[1], phase);
+    CASSIGN_MULT(field[2], field[2], phase);
+  }
+
+  F.x = cscalar2cnumber(field[0]);
+  F.y = cscalar2cnumber(field[1]);
+  F.z = cscalar2cnumber(field[2]);
+
+  return F;
+}
+
 void mode_solver::multiply_bloch_phase() {
 
   std::vector<mpb_real> kvector = get_output_k();
@@ -1936,6 +2175,108 @@ std::vector<mpb_real> mode_solver::compute_group_velocity_component(vector3 d) {
   return group_v;
 }
 
+/* as above, but only computes for given band */
+mpb_real mode_solver::compute_1_group_velocity_component(vector3 d, int b) {
+  mpb_real u[3];
+  int ib = b - 1;
+  mpb_real group_v;
+  mpb_real scratch;
+
+  curfield_reset();
+
+  if (!mdata) {
+    meep::master_fprintf(stderr, "mode_solver.init must be called first!\n");
+    return group_v;
+  }
+
+  if (!kpoint_index) {
+    meep::master_fprintf(stderr, "mode_solver.solve_kpoint must be called first!\n");
+    return group_v;
+  }
+
+  /* convert d to unit vector in Cartesian coords: */
+  d = unit_vector3(matrix3x3_vector3_mult(Gm, d));
+  u[0] = d.x;
+  u[1] = d.y;
+  u[2] = d.z;
+
+  evectmatrix_resize(&W[0], 1, 0);
+  CHECK(nwork_alloc > 1, "eigensolver-nwork is too small");
+  evectmatrix_resize(&W[1], 1, 0);
+
+  maxwell_compute_H_from_B(mdata, H, W[1], (scalar_complex *) mdata->fft_data, ib, 0, 1);
+  maxwell_ucross_op(W[1], W[0], mdata, u);
+  evectmatrix_XtY_diag_real(W[1], W[0], &group_v, &scratch);
+
+  /* Reset scratch matrix sizes: */
+  evectmatrix_resize(&W[1], W[1].alloc_p, 0);
+  evectmatrix_resize(&W[0], W[0].alloc_p, 0);
+
+  if (freqs[ib] == 0) {  /* v is undefined in this case */
+    group_v = 0.0;  /* just set to zero */
+  }
+  else {
+    group_v /= negative_epsilon_ok ? sqrt(fabs(freqs[ib])) : freqs[ib];
+  }
+  return group_v;
+}
+
+/* returns group velocity for band b, in Cartesian coordinates */
+vector3 mode_solver::compute_1_group_velocity(int b) {
+  vector3 v;
+  vector3 d;
+  matrix3x3 RmT = matrix3x3_transpose(Rm);
+  d.x = 1;
+  d.y = 0;
+  d.z = 0;
+  v.x = compute_1_group_velocity_component(matrix3x3_vector3_mult(RmT,d),b);
+  d.y = 1;
+  d.x = 0;
+  d.z = 0;
+  v.y = compute_1_group_velocity_component(matrix3x3_vector3_mult(RmT,d),b);
+  d.z = 1;
+  d.y = 0;
+  d.x = 0;
+  v.z = compute_1_group_velocity_component(matrix3x3_vector3_mult(RmT,d),b);
+
+  return v;
+}
+
+/* as above, but returns "group velocity" given by gradient of
+   frequency with respect to k in reciprocal coords ... this is useful
+   for band optimization. */
+vector3 mode_solver::compute_1_group_velocity_reciprocal(int b) {
+  return matrix3x3_vector3_mult(matrix3x3_transpose(Gm),
+                                compute_1_group_velocity(b));
+}
+
+/* compute the fraction of the field energy that is located in the
+   given range of dielectric constants: */
+mpb_real mode_solver::compute_energy_in_dielectric(mpb_real eps_low, mpb_real eps_high) {
+  mpb_real *energy = (mpb_real *) curfield;
+  mpb_real epsilon = 0.0;
+  mpb_real energy_sum = 0.0;
+
+  if (!curfield || !strchr("DHBR", curfield_type)) {
+    meep::master_fprintf(stderr, "The D or H energy density must be loaded first.\n");
+    return 0.0;
+  }
+
+  int N = mdata->fft_output_size;
+
+  for (int i = 0; i < N; ++i) {
+    epsilon = mean_medium_from_matrix(mdata->eps_inv +i);
+    if (epsilon >= eps_low && epsilon <= eps_high) {
+      energy_sum += energy[i];
+    }
+  }
+  mpi_allreduce_1(&energy_sum, mpb_real, SCALAR_MPI_TYPE, MPI_SUM, mpb_comm);
+  energy_sum *= vol / H.N;
+  return energy_sum;
+}
+
+
+
 bool mode_solver::with_hermitian_epsilon() {
 #ifdef WITH_HERMITIAN_EPSILON
   return true;
diff --git a/libpympb/pympb.hpp b/libpympb/pympb.hpp
index 97d966a39..bc93c3fe1 100644
--- a/libpympb/pympb.hpp
+++ b/libpympb/pympb.hpp
@@ -15,6 +15,7 @@ namespace py_mpb {
 struct mode_solver {
   static const int MAX_NWORK = 10;
   static const char epsilon_CURFIELD_TYPE = 'n';
+  static const char mu_CURFIELD_TYPE = 'm';
   static const int NUM_FFT_BANDS = 20;
 
   int num_bands;
@@ -90,6 +91,7 @@ struct mode_solver {
   int get_kpoint_index();
   void set_kpoint_index(int i);
   void get_epsilon();
+  void get_mu();
   void get_epsilon_tensor(int c1, int c2, int imag, int inv);
   void get_material_pt(meep_geom::material_type &material, vector3 p);
   void material_epsmu(meep_geom::material_type material, symmetric_matrix *epsmu,
@@ -132,12 +134,30 @@ struct mode_solver {
   std::vector<int> get_dims();
   std::vector<mpb_real> get_output_k();
 
+  mpb_real get_val(int ix, int iy, int iz, int nx, int ny, int nz, int last_dim_size,
+                   mpb_real *data, int stride, int conjugate);
+  mpb_real interp_val(vector3 p, int nx, int ny, int nz, int last_dim_size,
+                      mpb_real *data, int stride, int conjugate);
+  scalar_complex interp_cval(vector3 p, int nx, int ny, int nz, int last_dim_size,
+                             mpb_real *data, int stride);
+  symmetric_matrix interp_eps_inv(vector3 p);
+
+  mpb_real get_epsilon_point(vector3 p);
+  cmatrix3x3 get_epsilon_inverse_tensor_point(vector3 p);
+  mpb_real get_energy_point(vector3 p);
+  cvector3 get_field_point(vector3 p);
+  cvector3 get_bloch_field_point(vector3 p);
+
   void multiply_bloch_phase();
   void fix_field_phase();
   void compute_field_divergence();
   std::vector<mpb_real> compute_zparities();
   std::vector<mpb_real> compute_yparities();
   std::vector<mpb_real> compute_group_velocity_component(vector3 d);
+  mpb_real compute_1_group_velocity_component(vector3 d, int b);
+  vector3 compute_1_group_velocity(int b);
+  vector3 compute_1_group_velocity_reciprocal(int b);
+  mpb_real compute_energy_in_dielectric(mpb_real eps_low, mpb_real eps_high);
   bool with_hermitian_epsilon();
 
 private:
diff --git a/python/mpb.i b/python/mpb.i
index 598626d35..2dd0a6a98 100644
--- a/python/mpb.i
+++ b/python/mpb.i
@@ -104,6 +104,65 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) {
 
     return 1;
 }
+
+static PyObject* v3_to_pyv3(vector3 *v) {
+    PyObject *geom_mod = PyImport_ImportModule("meep.geom");
+    PyObject *v3_class = PyObject_GetAttrString(geom_mod, "Vector3");
+    PyObject *args = Py_BuildValue("(ddd)", v->x, v->y, v->z);
+    PyObject *py_v = PyObject_Call(v3_class, args, NULL);
+
+    Py_DECREF(geom_mod);
+    Py_DECREF(args);
+    Py_DECREF(v3_class);
+
+    return py_v;
+}
+
+static PyObject* cv3_to_pyv3(cvector3 *cv) {
+    PyObject *geom_mod = PyImport_ImportModule("meep.geom");
+    PyObject *v3_class = PyObject_GetAttrString(geom_mod, "Vector3");
+
+    vector3 r = cvector3_re(*cv);
+    vector3 i = cvector3_im(*cv);
+
+    Py_complex x, y, z;
+    x.real = r.x;
+    x.imag = i.x;
+    y.real = r.y;
+    y.imag = i.y;
+    z.real = r.z;
+    z.imag = i.z;
+
+    PyObject *args = Py_BuildValue("(DDD)", &x, &y, &z);
+    PyObject *py_v = PyObject_Call(v3_class, args, NULL);
+
+    Py_DECREF(geom_mod);
+    Py_DECREF(args);
+    Py_DECREF(v3_class);
+
+    return py_v;
+}
+
+static PyObject* cmatrix3x3_to_pymatrix(cmatrix3x3 *m) {
+    PyObject *c1 = cv3_to_pyv3(&m->c0);
+    PyObject *c2 = cv3_to_pyv3(&m->c1);
+    PyObject *c3 = cv3_to_pyv3(&m->c2);
+
+    PyObject *geom_mod = PyImport_ImportModule("meep.geom");
+    PyObject *matrix_class = PyObject_GetAttrString(geom_mod, "Matrix");
+
+    PyObject *args = Py_BuildValue("(OOO)", c1, c2, c3);
+    PyObject *res = PyObject_Call(matrix_class, args, NULL);
+
+    Py_DECREF(c1);
+    Py_DECREF(c2);
+    Py_DECREF(c3);
+    Py_DECREF(geom_mod);
+    Py_DECREF(matrix_class);
+    Py_DECREF(args);
+
+    return res;
+}
 %}
 
 %include "std_string.i"
@@ -132,6 +191,8 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) {
     meep_geom::material_data*
 };
 
+%apply double { mpb_real };
+
 %typemap(in) lattice {
     if (!pylattice_to_lattice($input, &$1)) {
         PyErr_PrintEx(0);
@@ -169,6 +230,30 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) {
     }
 }
 
+%typemap(out) vector3 {
+    $result = v3_to_pyv3(&$1);
+
+    if (!$result) {
+        SWIG_fail;
+    }
+}
+
+%typemap(out) cvector3 {
+    $result = cv3_to_pyv3(&$1);
+
+    if (!$result) {
+        SWIG_fail;
+    }
+}
+
+%typemap(out) cmatrix3x3 {
+    $result = cmatrix3x3_to_pymatrix(&$1);
+
+    if (!$result) {
+        SWIG_fail;
+    }
+}
+
 %include "pympb.hpp"
 
 %pythoncode %{
@@ -193,6 +278,8 @@ static int pylattice_to_lattice(PyObject *py_lat, lattice *l) {
         output_charge_density,
         output_bpwr,
         output_dpwr,
+        output_tot_pwr,
+        output_dpwr_in_objects,
         output_poynting,
         output_poynting_x,
         output_poynting_y,
diff --git a/python/solver.py b/python/solver.py
index 588796ae4..7575db258 100644
--- a/python/solver.py
+++ b/python/solver.py
@@ -1,5 +1,6 @@
 from __future__ import division, print_function
 
+import functools
 import os
 import numbers
 import re
@@ -146,6 +147,9 @@ def ExH(e, h):
 
         return flat_res
 
+    def get_epsilon(self):
+        self.mode_solver.get_epsilon()
+
     def get_bfield(self, which_band):
         return self._get_field('b', which_band)
 
@@ -180,6 +184,44 @@ def _get_field(self, f, band):
 
         return field
 
+    def get_epsilon_point(self, p):
+        return self.mode_solver.get_epsilon_point(p)
+
+    def get_epsilon_inverse_tensor_point(self, p):
+        return self.mode_solver.get_epsilon_inverse_tensor_point(p)
+
+    def get_energy_point(self, p):
+        return self.mode_solver.get_energy_point(p)
+
+    def get_field_point(self, p):
+        return self.mode_solver.get_field_point(p)
+
+    def get_bloch_field_point(self, p):
+        return self.mode_solver.get_bloch_field_point(p)
+
+    def get_tot_pwr(self, which_band):
+        self.get_dfield(which_band)
+        self.compute_field_energy()
+
+        dims = self.mode_solver.get_dims()
+        epwr = np.zeros(np.prod(dims))
+        self.mode_solver.get_curfield(epwr)
+        epwr = np.reshape(epwr, dims)
+
+        self.get_bfield(which_band)
+        self.compute_field_energy()
+
+        hpwr = np.zeros(np.prod(dims))
+        self.mode_solver.get_curfield(hpwr)
+        hpwr = np.reshape(hpwr, dims)
+
+        tot_pwr = epwr + hpwr
+
+        self.mode_solver.set_curfield(tot_pwr.ravel())
+        self.mode_solver.set_curfield_type('R')
+
+        return tot_pwr
+
     # The band-range-data is a list of tuples, each consisting of a (min, k-point)
     # tuple and a (max, k-point) tuple, with each min/max pair describing the
     # frequency range of a band and the k-points where it achieves its minimum/maximum.
@@ -288,10 +330,8 @@ def output_epsilon(self):
         self.output_field_to_file(mp.ALL, self.get_filename_prefix())
 
     def output_mu(self):
-        print("output_mu: Not yet supported")
-        # TODO
-        # self.mode_solver.get_mu()
-        # self.mode_solver.output_field_to_file(-1, self.get_filename_prefix)
+        self.mode_solver.get_mu()
+        self.output_field_to_file(mp.ALL, self.get_filename_prefix())
 
     def output_field_to_file(self, component, fname_prefix):
         curfield_type = self.mode_solver.get_curfield_type()
@@ -445,9 +485,15 @@ def _create_fname(self, fname, prefix, parity_suffix):
     def compute_field_energy(self):
         return self.mode_solver.compute_field_energy()
 
+    def compute_field_divergence(self):
+        mode_solver.compute_field_divergence()
+
     def compute_energy_in_objects(self, objs):
         return self.mode_solver.compute_energy_in_objects(objs)
 
+    def compute_energy_in_dielectric(self, eps_low, eps_high):
+        return self.mode_solver.compute_energy_in_dielectric(eps_low, eps_high)
+
     def compute_group_velocities(self):
         xarg = mp.cartesian_to_reciprocal(mp.Vector3(1), self.geometry_lattice)
         vx = self.mode_solver.compute_group_velocity_component(xarg)
@@ -458,6 +504,16 @@ def compute_group_velocities(self):
 
         return [mp.Vector3(x, y, z) for x, y, z in zip(vx, vy, vz)]
 
+    def compute_group_velocity_component(self, direction):
+        return self.mode_solver.compute_group_velocity_component(direction)
+
+    def compute_one_group_velocity(self, which_band):
+        return self.mode_solver.compute_1_group_velocity(which_band)
+
+    def compute_one_group_velocity_component(self, direction, which_band):
+        return self.mode_solver.compute_1_group_velocity_component(direction,
+                                                                   which_band)
+
     def randomize_fields(self):
         self.mode_solver.randomize_fields()
 
@@ -715,6 +771,42 @@ def bfunc(ms, b_prime):
 
         return ks
 
+    def first_brillouin_zone(self, k):
+        """
+        Function to convert a k-point k into an equivalent point in the
+        first Brillouin zone (not necessarily the irreducible Brillouin zone)
+        """
+        def n(k):
+            return mp.reciprocal_to_cartesian(k, self.geometry_lattice).norm()
+
+        def try_plus(k, v):
+            if n(k + v) < n(k):
+                return try_plus(k + v, v)
+            else:
+                return k
+
+        def _try(k, v):
+            return try_plus(try_plus(k, v), mp.Vector3() - v)
+
+        try_list = [
+            mp.Vector3(1, 0, 0), mp.Vector3(0, 1, 0), mp.Vector3(0, 0, 1),
+            mp.Vector3(0, 1, 1), mp.Vector3(1, 0, 1), mp.Vector3(1, 1, 0),
+            mp.Vector3(0, 1, -1), mp.Vector3(1, 0, -1), mp.Vector3(1, -1, 0),
+            mp.Vector3(1, 1, 1), mp.Vector3(-1, 1, 1), mp.Vector3(1, -1, 1),
+            mp.Vector3(1, 1, -1),
+        ]
+
+        def try_all(k):
+            return functools.reduce(_try, try_list, k)
+
+        def try_all_and_repeat(k):
+            knew = try_all(k)
+            return try_all_and_repeat(knew) if n(knew) < n(k) else k
+
+        k0 = k - mp.Vector3(*[round(x) for x in k])
+
+        return try_all_and_repeat(k0) if n(k0) < n(k) else try_all_and_repeat(k)
+
 
 # Predefined output functions (functions of the band index), for passing to `run`
 
@@ -810,6 +902,32 @@ def output_dpwr(ms, which_band):
     ms.output_field()
 
 
+def output_tot_pwr(ms, which_band):
+    ms.get_tot_pwr(which_band)
+    ms.output_field_to_file(-1, ms.get_filename_prefix() + 'tot.')
+
+
+def output_dpwr_in_objects(output_func, min_energy, objects=[]):
+    """
+    The following function returns an output function that calls output_func for
+    bands with D energy in objects > min-energy. For example,
+    output_dpwr_in_objects(output_dfield, 0.20, some_object) would return an
+    output function that would spit out the D field for bands with at least %20
+    of their D energy in some-object.
+    """
+
+    def _output(ms, which_band):
+        ms.get_dfield(which_band)
+        ms.compute_field_energy()
+        energy = ms.compute_energy_in_objects(objects)
+        fmt = "dpwr:, {}, {}, {} "
+        print(fmt.format(which_band, ms.freqs[which_band - 1], energy))
+        if energy >= min_energy:
+            apply_band_func(ms, output_func, which_band)
+
+    return _output
+
+
 def output_charge_density(ms, which_band):
     ms.get_charge_density(which_band)
     ms.output_field_to_file(-1, ms.get_filename_prefix())
diff --git a/python/tests/data/tutorial-d.k01.b02.te.h5 b/python/tests/data/tutorial-d.k01.b02.te.h5
new file mode 100644
index 000000000..82daee09f
Binary files /dev/null and b/python/tests/data/tutorial-d.k01.b02.te.h5 differ
diff --git a/python/tests/data/tutorial-d.k16.b02.te.h5 b/python/tests/data/tutorial-d.k16.b02.te.h5
new file mode 100644
index 000000000..82daee09f
Binary files /dev/null and b/python/tests/data/tutorial-d.k16.b02.te.h5 differ
diff --git a/python/tests/data/tutorial-mu.h5 b/python/tests/data/tutorial-mu.h5
new file mode 100644
index 000000000..bbba29cdf
Binary files /dev/null and b/python/tests/data/tutorial-mu.h5 differ
diff --git a/python/tests/data/tutorial-tot.rpwr.k16.b08.te.h5 b/python/tests/data/tutorial-tot.rpwr.k16.b08.te.h5
new file mode 100644
index 000000000..fa0d674f2
Binary files /dev/null and b/python/tests/data/tutorial-tot.rpwr.k16.b08.te.h5 differ
diff --git a/python/tests/mpb.py b/python/tests/mpb.py
index 22cc06da5..f6c061e60 100644
--- a/python/tests/mpb.py
+++ b/python/tests/mpb.py
@@ -5,6 +5,7 @@
 import os
 import re
 import sys
+import time
 import unittest
 
 import h5py
@@ -25,6 +26,8 @@ def setUp(self):
         """Store the test name and register a function to clean up all the
         generated h5 files."""
 
+        self.start = time.time()
+
         self.filename_prefix = self.id().split('.')[-1]
         print()
         print(self.filename_prefix)
@@ -38,6 +41,10 @@ def rm_h5():
 
         self.addCleanup(rm_h5)
 
+    def tearDown(self):
+        end = time.time() - self.start
+        print("{}: {:.2f}s".format(self.filename_prefix, end))
+
     def init_solver(self, geom=True):
         num_bands = 8
         k_points = [
@@ -109,6 +116,35 @@ def test_list_split(self):
         for res, exp in zip(k_split[1], expected_list):
             self.assertTrue(res.close(exp))
 
+    def test_first_brillouin_zone(self):
+        ms = self.init_solver()
+
+        res = []
+        for k in ms.k_points:
+            res.append(ms.first_brillouin_zone(k))
+
+        expected = [
+            mp.Vector3(0.0, 0.0, 0.0),
+            mp.Vector3(0.10000000000000003, 0.0, 0.0),
+            mp.Vector3(0.20000000000000004, 0.0, 0.0),
+            mp.Vector3(0.30000000000000004, 0.0, 0.0),
+            mp.Vector3(0.4, 0.0, 0.0),
+            mp.Vector3(0.5, 0.0, 0.0),
+            mp.Vector3(0.5, 0.10000000000000003, 0.0),
+            mp.Vector3(0.5, 0.20000000000000004, 0.0),
+            mp.Vector3(0.5, 0.30000000000000004, 0.0),
+            mp.Vector3(0.5, 0.4, 0.0),
+            mp.Vector3(0.5, 0.5, 0.0),
+            mp.Vector3(0.4, 0.4, 0.0),
+            mp.Vector3(0.30000000000000004, 0.30000000000000004, 0.0),
+            mp.Vector3(0.2, 0.2, 0.0),
+            mp.Vector3(0.1, 0.1, 0.0),
+            mp.Vector3(0.0, 0.0, 0.0),
+        ]
+
+        for e, r in zip(expected, res):
+            self.assertTrue(e.close(r))
+
     def check_band_range_data(self, expected_brd, result, places=3, tol=1e-3):
         for exp, res in zip(expected_brd, result):
             # Compare min freqs
@@ -233,8 +269,6 @@ def test_run_te(self):
                           0.9229965195855876, 1.0001711176882833, 1.0001720032257042,
                           1.092820931747826]
 
-        # Obtained by passing NULL to set_maxwell_dielectric for epsilon_mean_func
-        # in mpb/mpb/medium.c.
         expected_brd = [
             ((0.0, mp.Vector3(0.0, 0.0, 0.0)),
              (0.49683586474489877, mp.Vector3(0.5, 0.5, 0.0))),
@@ -269,6 +303,9 @@ def test_run_te(self):
 
         self.check_gap_list(expected_gap_list, ms.gap_list)
 
+        pt = ms.get_epsilon_point(mp.Vector3(0.5, 0.5))
+        self.assertEqual(pt, 1.0)
+
     def test_run_tm(self):
         expected_brd = [
             ((0.0, mp.Vector3(0.0, 0.0, 0.0)),
@@ -333,12 +370,70 @@ def test_compute_field_energy(self):
         ms = self.init_solver()
         ms.run_te()
         ms.get_dfield(8)
-        res = ms.compute_field_energy()
+        field_pt = ms.get_field_point(mp.Vector3(0.5, 0.5))
+        bloch_field_pt = ms.get_bloch_field_point(mp.Vector3(0.5, 0.5))
+        eps_inv_tensor = ms.get_epsilon_inverse_tensor_point(mp.Vector3(0.5, 0.5))
+
+        energy = ms.compute_field_energy()
+        pt = ms.get_energy_point(mp.Vector3(0.5, 0.5))
+
+        self.assertAlmostEqual(pt, 1.330368347216153e-9)
+
+        expected_fp = mp.Vector3(
+            2.5823356723958247e-5 + 6.713243287584132e-12j,
+            -2.575955745071957e-5 - 6.696552990958943e-12j,
+            0.0 - 0.0j
+        )
+
+        expected_bloch_fp = mp.Vector3(
+            2.5823356723958247e-5 + 6.713243287584132e-12j,
+            -2.575955745071957e-5 - 6.696552990958943e-12j,
+            -0.0 - 0.0j
+        )
+
+        expected_eps_inv_tensor = mp.Matrix(
+            mp.Vector3(1.0 + 0.0j, 0.0 - 0.0j, 0.0 - 0.0j),
+            mp.Vector3(0.0 + 0.0j, 1.0 + 0.0j, 0.0 - 0.0j),
+            mp.Vector3(0.0 + 0.0j, 0.0 + 0.0j, 1.0 + 0.0j)
+        )
+
+        self.assertTrue(expected_fp.close(field_pt))
+        self.assertTrue(expected_bloch_fp.close(bloch_field_pt))
+        self.assertEqual(expected_eps_inv_tensor.c1, eps_inv_tensor.c1)
+        self.assertEqual(expected_eps_inv_tensor.c2, eps_inv_tensor.c2)
+        self.assertEqual(expected_eps_inv_tensor.c3, eps_inv_tensor.c3)
+
+        energy_in_dielectric = ms.compute_energy_in_dielectric(0, 1)
+
+        expected_energy = [1.0000000000000002, 1.726755206037815e-5, 0.4999827324479414,
+                           1.7267552060375955e-5, 0.4999827324479377, 0.0, 0.0]
 
-        expected = [1.0000000000000002, 1.726755206037815e-5, 0.4999827324479414,
-                    1.7267552060375955e-5, 0.4999827324479377, 0.0, 0.0]
+        expected_energy_in_dielectric = 0.6990769686037558
 
-        self.compare_arrays(np.array(expected), np.array(res))
+        self.compare_arrays(np.array(expected_energy), np.array(energy))
+        self.assertAlmostEqual(expected_energy_in_dielectric, energy_in_dielectric, places=3)
+
+    def test_compute_group_velocity(self):
+        ms = self.init_solver()
+        ms.run_te()
+        res1 = ms.compute_group_velocity_component(mp.Vector3(0.5, 0.5))
+        res2 = ms.compute_one_group_velocity(8)
+        res3 = ms.compute_one_group_velocity_component(mp.Vector3(0.5, 0.5), 8)
+
+        expected1 = [
+            0.0, 1.470762578355642e-4, -1.4378185933055663e-4, 1.1897503996483383e-4,
+            -4.892687048681629e-4, 1.1240346140784176e-4, 1.5842474585356007e-4,
+            4.496945573323881e-5,
+        ]
+
+        expected2 = mp.Vector3(3.180010979062989e-5, 3.179611968757397e-5)
+        expected3 = 4.496932512216992e-5
+
+        for e, r in zip(expected1, res1):
+            self.assertAlmostEqual(e, r, places=4)
+
+        self.assertTrue(expected2.close(res2, tol=3))
+        self.assertAlmostEqual(expected3, res3, places=3)
 
     def test_output_efield_z(self):
         ms = self.init_solver()
@@ -351,6 +446,22 @@ def test_output_efield_z(self):
         res_path = re.sub('tutorial', ms.filename_prefix, ref_fname)
         self.compare_h5_files(ref_path, res_path)
 
+    def test_output_dpwr_in_objects(self):
+        ms = self.init_solver()
+        ms.run_te(mpb.output_dpwr_in_objects(mpb.output_dfield, 0.85, ms.geometry))
+
+        ref_fname1 = 'tutorial-d.k01.b02.te.h5'
+        ref_fname2 = 'tutorial-d.k16.b02.te.h5'
+
+        ref_path1 = os.path.join(self.data_dir, ref_fname1)
+        ref_path2 = os.path.join(self.data_dir, ref_fname2)
+
+        res_path1 = re.sub('tutorial', ms.filename_prefix, ref_fname1)
+        res_path2 = re.sub('tutorial', ms.filename_prefix, ref_fname2)
+
+        self.compare_h5_files(ref_path1, res_path1)
+        self.compare_h5_files(ref_path2, res_path2)
+
     def test_triangular_lattice(self):
 
         expected_brd = [
@@ -904,5 +1015,16 @@ def test_run_te_with_mu_material(self):
         ms.run_te()
         self.check_band_range_data(expected_brd, ms.band_range_data)
 
+    def test_output_tot_pwr(self):
+        ms = self.init_solver()
+        ms.run_te()
+        mpb.output_tot_pwr(ms, 8)
+
+        ref_fname = 'tutorial-tot.rpwr.k16.b08.te.h5'
+        ref_path = os.path.join(self.data_dir, ref_fname)
+        res_path = re.sub('tutorial', self.filename_prefix, ref_fname)
+
+        self.compare_h5_files(ref_path, res_path)
+
 if __name__ == '__main__':
     unittest.main()