provide dominant planewave component

[stevengj]

For various calculations, it would be convenient to have a function to return the dominant planewave component. That is, find the G vector corresponding to the biggest amplitude planewave (in the planewave expansion of an eigenmode), and return k+G (the corresponding planewave wavevector).

Has to be implemented at the C level, as a new function in src/maxwell, but should be straightforward since we already have all the requisite information.

(Would help with NanoComp/meep#476, for example.)

[stevengj]

In particular, I'm thinking of a function like:

void maxwell_dominant_planewave(maxwell_data *d, evectmatrix H, int band, double kdom[3])
{
  /* check d != NULL, 1 <= band <= H.p */

  /* find biggest planewave component in band b of H on current process: */
  double max_amp = 0, max_i = 0;
  for (i = 0; i < H.localN; ++i) {
    double amp = SCALAR_NORMSQR(H.data[(i*2+0)*H.p + band-1]) + SCALAR_NORMSQR(H.data[(i*2+1)*H.p + band-1]);
    if (amp > max_amp) {
      max_amp = amp;
      max_i = i;
    }
  }

  k_data cur_k = d->k_plus_G[max_i];
  /* set kdom to cur_k.kmag * (cross product of cur_k.m and cur_k.n) */ 
  ... vector math ...

  /* parallel only: find kdom of biggest max_amp across all processes */
  ... some MPI stuff ...
}

Upon return, kdom is the wavevector (k+G), in Cartesian coordinates, of the dominant planewave component for band (as determined by the amplitudes in H, which should normally be the eigenfunctions).

cc @ChristopherHogan

[stevengj]

(Alternatively we could have a function to return the average planewave vector, weighted by the mode amplitudes. For eigenfunctions in uniform media, these should the same as the maximum. The code is very similar except for summing instead of taking the maximum.)