FDTD Poynting Vector

Question

I'm attempting to validate my FDTD code against Meep by calculating the Poynting vector field across a simulation consisting of a monochromatic point source within a box (no boundary conditions (Dirichlet) so everything just bounces off the domain boundaries.

In this case, the simulation is 2D TMz (point source in $E{z}$, coupling into $H{x}$ and $H{y}$).

After quite a bit of research and talking with my thesis advisor, I wound up using this for Poynting vector $\hat{P}$ calculation:

$\hat{E} \times \hat{H}^*$ (complex conjugate is the same as the original $H$ vector in this case since everything is in the real domain).

$ \begin{bmatrix} \hat{P}_X\\ \hat{P}_Y\\ \hat{P}_Z \end{bmatrix} = \begin{bmatrix} E_YH_Z - E_zH_Y\\ E_zH_X - E_XH_Z\\ E_XH_Y - E_YH_X\\ \end{bmatrix} $

Given that $\hat{E} = \begin{bmatrix}0 & 0 & E_Z\end{bmatrix}$ and $\hat{H} = \begin{bmatrix}H_X & H_Y & 0\end{bmatrix}$

(since there are no $E_X$, $E_Y$, or $H_Z$ components in a 2D TMz sim)

$ \begin{bmatrix} \hat{P}_X\\ \hat{P}_Y\\ \hat{P}_Z \end{bmatrix} = \begin{bmatrix} - E_ZH_Y\\ E_ZH_X\\ 0 \end{bmatrix} $

When calculating the Poynting vector field, I iterate through the entire $E_Z$ array, and use the above formula for obtaining the $P$ vector at each cell. Then, basic algebra to get the magnitude and direction:

$|P| = \sqrt{P_X^2 + P_Y^2 + P_Z^2} \\ \\ P_{normal} = {P \over{|P|}} $

Now, plotting the magnitude before the source wave hits the boundary of the simulation ($E_Z$ pictured below),

I would expect the Poynting vector magnitude field to resemble a cone, indicating that the magnitude falls off with $1\over{R^2}$ distance $R$ from the source position.

However, I'm seeing this:

Now, my question is: does my math look right? Or, am I missing something here? Any suggestions or links to example code are welcome.

Even a picture of what the Poynting magnitude and direction field should look like for a monochromatic point source would be helpful.

Answer 1

After much deliberation, I tracked down the problem to a few software bugs:

1. An array index is not a field value.
2. Floating point precision is a factor when squaring a VERY small number.
3. NaN really isn't a number, and trying to normalize a zero-length vector requires careful thought.
4. The really big misconception on my part: since the power flows between adjacent E points via connecting H points, and the source is a sinusoid, the magnitude (and Poynting vector) will, of course, oscillate. I've added a pair of RMS calculations that looks at the average and peak values over a period, and I'm getting something more like what I expected.

Ultimately, my math was correct, but my implementation had a few issues. Thankfully it was on the Poynting vector calculation side rather than the FDTD side.