tutorial for non-axisymmetric sources in cylindrical coordinates

[stevengj]

It would be nice to have a tutorial showing how to use non-axisymmetric sources in cylindrical coordinates by expanding the source in a exp(imφ) Fourier series. For example:

• An off-axis planewave can be handled with the Jacobi–Anger expansion, as mentioned in tutorial with cylindrical coordinates and planewaves #902. Another example can be found in this paper eqs (5–10) and many other papers.

• LDOS for off-axis points. This is a point source, so when expanded in a Fourier series it is just a constant as a function of m. Power orthogonality means that you can then just sum the LDOS's for different m's.

[stevengj]

In particular, you can do two calculations

1. One with point dipole at the origin, with a single LDOS calculation
2. One with the point dipole at r > 0, computing the LDOS at m=0,1,2,3,4,… .

If you sum the LDOS(m) from (2), multiplying all of the m ≠ 0 terms by 2 (to account for the -m solutions), then it should converge to the answer from (1). (There might be a constant factor of 2 or π or something.)

The farther off-axis you go, the more m's it should require to converge. Would be good to have a plot of LDOS(m) as a function of m from (2). I'm guessing you'll need m's up to about radius×ω.

[oskooi]

Two source polarizations/configurations to check:

1. $E_z$ at $r=0$ with $m=0$
2. $E_r$ at $r=0$ with $m=1$

[oskooi]

Here is an initial attempt (not working) to demonstrate that the radiated flux in air from an $E_r$ point source in a dielectric layer above a lossless ground plane (same setup as the LED extraction efficiency tutorial) is the same for two different configurations: (1) $r=0$ and $m=-1$ (one simulation) and (2) $r\neq 0$ summed over $m=0,1,2,3...M-1$ ($M$ simulations). For (2), a plot of flux vs. $m$ shows that for this particular source location the flux is negligible (<1e-6) for $m&gt;8$. I also verified that the $+m$ and $-m$ simulations produce the same flux in (2) as expected.

However, there is a rather large discrepancy in the radiated flux for (1) and (2): 0.33 vs. 1014.46. This could mean there are scaling coefficients of the individual flux values in (2) that need to be included when computing their sum.

"""Verifies that the radiated flux from a point dipole in a                                                       
   dielectric layer above a lossless ground plane computed in                                                     
   cylindrical coordinates is the same for a dipole placed at                                                     
   either r=0 or r≠0.                                                                                             
"""

import numpy as np
import matplotlib
matplotlib.use("agg")
import matplotlib.pyplot as plt
import meep as mp

resolution = 50  # pixels/μm                                                                                      
dpml = 0.5       # thickness of PML                                                                               
dair = 1.0       # thickness of air padding                                                                       
L = 6.0          # length of non-PML region                                                                       
n = 2.4          # refractive index of surrounding medium                                                         
wvl = 1.0        # wavelength (in vacuum)                                                                         

fcen = 1 / wvl   # center frequency of source/monitor                                                             

# source properties (cylindrical)                                                                                 
df = 0.05 * fcen
cutoff = 10.0
src = mp.GaussianSource(fcen, fwidth=df, cutoff=cutoff)

# field decay termination criteria                                                                                
tol = 1e-8

def radiated_flux(dmat: float, h: float, rpos: float, m: int) -> float:
    """Computes the radiated flux of a point dipole embedded                                              
       within a dielectric layer above a lossless ground plane in                                                 
       cylindrical coordinates.                                                                                   
                                                                                                                  
    Args:                                                                                                         
      dmat: thickness of dielectric layer.                                                                        
      h: height of dipole above ground plane as fraction of dmat.                                                 
      rpos: position of source in r direction.                                                                    
      m: angular φ dependence of the fields exp(imφ).                                                             
    """
    sr = L + dpml
    sz = dmat + dair + dpml
    cell_size = mp.Vector3(sr, 0, sz)

    boundary_layers = [
        mp.PML(dpml, direction=mp.R),
        mp.PML(dpml, direction=mp.Z, side=mp.High),
    ]

    src_cmpt = mp.Er
    src_pt = mp.Vector3(rpos, 0, -0.5 * sz + h * dmat)
    sources = [mp.Source(src=src, component=src_cmpt, center=src_pt)]

    geometry = [
        mp.Block(
            material=mp.Medium(index=n),
            center=mp.Vector3(0, 0, -0.5 * sz + 0.5 * dmat),
            size=mp.Vector3(mp.inf, mp.inf, dmat),
        )
    ]

    sim = mp.Simulation(
        resolution=resolution,
        cell_size=cell_size,
        dimensions=mp.CYLINDRICAL,
        m=m,
        boundary_layers=boundary_layers,
        sources=sources,
        geometry=geometry,
    )

    flux_air = sim.add_flux(
        fcen,
        0,
        1,
        mp.FluxRegion(
            center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
            size=mp.Vector3(L, 0, 0),
        ),
        mp.FluxRegion(
            center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * dair),
            size=mp.Vector3(0, 0, dair),
        ),
    )

    sim.run(
        until_after_sources=mp.stop_when_fields_decayed(
            20,
            src_cmpt,
            src_pt,
            tol,
        ),
    )                                    

    out_flux = mp.get_fluxes(flux_air)[0]

    return out_flux


if __name__ == "__main__":
    layer_thickness = 0.7 * wvl / n
    dipole_height = 0.5

    # (1) r = 0                                                                                                   
    rpos = 0
    m = -1
    ref_out_flux = radiated_flux(
        layer_thickness,
        dipole_height,
        rpos,
        m,
    )
    print(f"flux-ref:, {ref_out_flux}")

    # (2) r ≠ 0                                                                                                   
    rpos = 0.215
    num_m = 20
    out_flux = np.zeros(num_m)
    for m in range(num_m):
        out_flux[m] = radiated_flux(
            layer_thickness,
            dipole_height,
            rpos,
            m,
        )
        print(f"flux:, {m}, {out_flux[m]}")

    sum_out_flux = out_flux[0] + 2*np.sum(out_flux[1:])

    print(f"flux-compare:, {ref_out_flux:.4f} (r=0), {sum_out_flux:.4f} (r≠0)")

[kalosu]

Hi Meep community,

I was wondering if in this case(using cylindrical coordinates) is it possible to have a point like source object which is meant to be polarized along x. I know that for a plane wave excitation that is possible (based on the Tutorial on cylindrical coordinates) but for a point-like (dipole emitter) I am not fully sure if this works as well?

I am interested in simulating such a dipole located at r=0 with a polarization given a single cartesian component.

[oskooi]

I was wondering if in this case(using cylindrical coordinates) is it possible to have a point like source object which is meant to be polarized along x.

In cylindrical coordinates, the $E_r$ polarization corresponds to the "in-plane" direction for a rotationally symmetric structure.

I am interested in simulating such a dipole located at r=0 with a polarization given a single cartesian component.

We already have a tutorial for point dipole emission at $r=0$: https://meep.readthedocs.io/en/latest/Python_Tutorials/Local_Density_of_States/#extraction-efficiency-of-a-light-emitting-diode-led.

[stevengj]

Can you try two different nonzero radii (e.g. double the radius), and see if the answer just change proportional to r or r^2?

[stevengj]

You could also try Ep vs Er — they should be the same?

[oskooi]

Can you try two different nonzero radii (e.g. double the radius), and see if the answer just change proportional to r or r^2?

Looks like the radiated flux varies as $r^2$ based on the data for an $E_r$ source at $r$, $2r$, and $4r$ where $r=0.215$.

>>> flux = [1077.0635, 4479.4474, 19244.9381]  # flux @ R = r, 2r, 4r; resolution = 100
>>> flux[2] / flux[1]  # flux(4r) / flux(2r)
4.296275049462574
>>> flux[1] / flux[0]  # flux(2r) / flux(r)
4.15894457476277
>>> flux[2] / flux[0]  # flux(4r) / flux(r)
17.867969808651022

This means there is likely a $Cr^2$ factor that we need to divide the $r\neq 0$ flux values when comparing to the $r=0$ result.

You could also try Ep vs Er — they should be the same?

The radiated flux is the (roughly the) same for $E_\phi$ and $E_r$ at $r$, $2r$, and $4r$:

>>> flux_Ep = [1107.9242, 4369.4396] # flux @ R = r, 2r, 4r; resolution = 50
>>> flux_Er = [1014.4596, 4359.5257] # flux @ R = r, 2r, 4r; resolution = 50

[oskooi]

It appears that the further the point source is positioned from $r=0$, the more Fourier orders of $\exp(im\phi)$ are necessary for the summed radiated flux to converge. As a demonstration, for a point source positioned at $r$, $2r$, and $3r$ where $r=0.215$, the required number of Fourier orders is 12, 22, and 33.

It would be useful to know apriori how many terms to compute in order to minimize the number of simulations necessary to obtain an accurate result using this approach.

[stevengj]

@kalosu, please don't post unrelated discussions in this issue. Start a new discussion thread if needed.

[stevengj]

I think there is a factor of 1/2πr that we are missing in the source normalization, which is why you are seeing (2πr)^2 scaling in the power (~ source^2). For r=0, I think this corresponds to r=0.5/resolution, so there will be a resolution-dependent factor for the r=0 source as we already saw in the LDOS calculations.

(The fact that the Fourier series takes longer to converge with m as the radius increases is expected.)

[oskooi]

I think there is a factor of 1/2πr that we are missing in the source normalization, which is why you are seeing (2πr)^2 scaling in the power (~ source^2).

In order to empirically determine the constant factor which relates the radiated flux from (a) a point source at $r=0$ to (b) a point source at $r\neq 0$, I made a plot of the ratio of the flux from (a) to the flux from (b) divided by $(2\pi r)^2$ as a function of $r$. This flux ratio seems to be converging to a value ~0.0005 as $r \rightarrow \infty$. A value of 0.0005 is roughly equivalent to $2(\pi/res)^3$ where $res$ is the resolution (50 pixels/μm in this example). For reference, we used $\Delta V = \pi/(res)^3=0.000025$ in #2148 (comment) when comparing the radiated power from a point source at $r=0$ in cylindrical coordinates to a point source in 3d.

[stevengj]

I think the correct factor for LDOS is probably ΔV = Δr * Δz * 2πr, where for r=0 you treat r as Δr/2. This is where you get the Δr * Δz * 2πΔr/2 = π/resolution^3 factor in the r=0 LDOS.

Here, however, you aren't computing LDOS, you are computing total power, so I think maybe you're getting a missing 2πr factor from the current source normalization (since we normalize the currents as a delta function but IIRC we don't include the 2πr from a δ function in cylindrical coordinates?).

Try plotting flux(r=0)/flux(r≠0) * (2r/Δr)^2?

[oskooi]

From the convergence study, it looks like the radiated flux does converge with (1) s (non-PML thickness) and (2) PML thickness in the $r$ direction but not (3) PML thickness in the $z$ direction.

[stevengj]

If you are still having a problem where the (m=+1, fcen=+1) simulation is not exactly a complex conjugate of (m=-1, fcen=-1), with the same flux, it would be good to dig into this. Maybe the (+1,+1) and (-1,-1) simulations for just a few timesteps — before the fields reach the PML — and see if the complex fields are exactly conjugates. If they aren't, then there is a complex factor we are missing in the timestepping? If they are, run it longer (until the fields reach the boundaries) to see if the PML is the source of the problem.

[stevengj]

Note that in order to give yourself the freedom to change the time profile as you change the radius, you might want to normalize the power by the Fourier transform of the current source. One way to get that is from the dft_ldos, which provides this in the Jdata field (#2148).

[oskooi]

edit (1/16): improved the formatting of the results.

see if the complex fields are exactly conjugates. If they aren't, then there is a complex factor we are missing in the timestepping? If they are, run it longer (until the fields reach the boundaries) to see if the PML is the source of the problem.

As previously reported, computing the radiated flux from an $E_r$ point source at $r=z=0$ in free space for a cell surrounded by PML on all sides for various configurations of m (angular dependence of the fields) and fcen (center frequency of the pulse) produces an unexpected result:

m	fcen	flux
-1	+1	A
-1	-1	A
+1	+1	B
+1	-1	B

Note: the third column (the radiated flux) consists of two values represented by constants A and B.

As discussed, we would expect row 1 (m=-1, fcen=+1) and row 4 (m=+1, fcen=-1) to produce the same flux because the fields $e^{-i\omega t}e^{im\phi}$ are complex conjugates. (The same would be expected for row 2 and row 3.) However, the results in the table seem to indicate that row 1 and row 2 are complex conjugate (as are row 3 and row 4) because their flux values are equivalent.

For reference, the correct set of results would be:

m	fcen	flux
-1	+1	A
-1	-1	B
+1	+1	B
+1	-1	A

As a second experiment, we compare the complex fields (rather than the flux) obtained from sim.get_array from two different runs ( $F_1$ and $F_2$ ) by computing their norm $|| F_1 - F_2^* ||$ where $^*$ is the complex conjugate as well as $|| F_1 - F_2 ||$.

$F_1$	$F_2$	$||F_1-F_2^*||$	$||F_1-F_2||$
-1,+1	-1,-1	0	C
-1,+1	+1,+1	C	0
-1,+1	+1,-1	0	C
-1,-1	+1,+1	0	C
-1,-1	+1,-1	C	0
+1,+1	+1,-1	0	C

Note: the third and fourth columns consists of two values: 0 and C≠0.

Row 1 indicates that (m=-1, fcen=+1) and (m=-1, fcen=-1) are indeed complex conjugates (i.e., the norm is zero). (Row 6 indicates that (m=+1, fcen=+1) and (m=+1, fcen=-1) are also indeed conjugates.) This result is the same with and without the PMLs. These findings are consistent with the first experiment.

However, in the second experiment, row 3 indicates that (m=-1, fcen=+1) and (m=+1, fcen=-1) are complex conjugates. (Row 4 indicates that (m=-1, fcen=-1) and (m=+1, fcen=+1) are complex conjugates.) But these results are inconsistent with the first experiment since row 1 (m=-1, fcen=+1) and row 4 (m=+1, fcen=-1) have flux values which are not the same.

--> Is it possible that the Poynting flux is somehow not being computed correctly in cylindrical coordinates?

link to experiment 1 script and results

link to experiment 2 script and results

[smartalecH]

Is it possible that the Poynting flux is somehow not being computed correctly in cylindrical coordinates?

So it looks like you verified that the time-domain fields are complex conjugates. Maybe check that the raw dft fields are as expected. Then you can narrow it down to the Poynting flux routine.

[stevengj]

The time domain fields should be conjugates … 

The weird thing is that the conclusion from the table is that:

1. Flipping the sign of m does nothing
2. Conjugating the source therefore conjugates the solution, regardless of the sign of m.

It's weird to me that flipping m has no effect since we clearly use the sign of m in the update equations:

meep/src/step_db.cpp

Lines 266 to 271 in 19b8abe

	else // no PML - fu - conductivity
	for (int ir = ir0 == 0; ir <= gv.nr(); ++ir) {
	realnum rinv = the_m / (ir + ir0);
	for (int iz = 0; iz <= gv.nz(); ++iz) {
	ptrdiff_t idx = ir * sr + iz;
	the_f[idx] += rinv * g[idx];

Can you double-check that m is indeed negative? i.e. check sim.fields.m

[oskooi]

Can you double-check that m is indeed negative? i.e. check sim.fields.m

I verified that m is indeed passed in correctly to sim.fields.m and that internally in the C++ code the_m changes sign with m.

Following a request from @smartalecH, I verified that that DFT fields have the same properties as the time-domain fields shown in the second table of results above (as expected). ref: link to experiment 3 script and results.

These results indicate that we likely need to look elsewhere for the cause of the discrepancy in the results among experiments 1 and 2/3.

[stevengj]

Maybe try a custom source that is purely real, so we can eliminate the variable of the source term, e.g. a source sin(ωt). Then the only imaginary parts of the fields should come from the m terms in the curl equations.

[oskooi]

Maybe try a custom source that is purely real, so we can eliminate the variable of the source term, e.g. a source sin(ωt).

Unfortunately, this approach does not work. This is because a purely real source implemented as:

    sources = [
        mp.Source(
            src=mp.CustomSource(src_func=lambda t: np.sin(2*np.pi*fcen*t)),
            center=mp.Vector3(),
            component=mp.Er,
	),
    ]

produces purely real fields for a m=±1 simulation everywhere in the cell.

If instead we use a purely imaginary source (1j*np.sin(2*np.pi*fcen*t))), the results for the time-domain and DFT fields are identical to the second table shown above with six data rows.

[oskooi]

Reviewing Maxwell's equations in cylindrical coordinates (reproduced below), a real current for $E_r$ (which enters on the RHS of eqn. 4) would produce (purely) real $E_r$ fields which would in turn produce purely imaginary $H_z$ fields for $m=\pm 1$ based on eqn. 3. In my findings $H_z$ is indeed purely imaginary whereas $E_r$ and $E_z$ are purely real. This is consistent with these equations.

[stevengj]

From the equations, if you have a purely real Er source in equation (4), you will get purely real {Er, Ez, Hφ} and purely imaginary {Eφ, Hr, Hz}. (This is something we could exploit to get a factor-of-2 speedup for any m≠0, by the way?)

If you flip the sign of m, then it should complex-conjugate the solution, which flips the signs of {Eφ, Hr, Hz}.

It would be worth checking this for both m=1 and m=2 (it should be true for both), just in case there is a bug in the special handling at the origin for m=1.

[oskooi]

I also discovered from rerunning experiment #1 involving the radiated flux that |m|=2, 3, ... always produces the same flux for all four test cases:

m	fcen	flux
-2	+1	A
-2	-1	A
+2	+1	A
+2	-1	A

This could be an indication that there is a bug only for the m=±1 case.

[stevengj]

In particular, this might be a bug in the m=1 code at r=0 here which looks like should include a * m factor

realnum f_m_mult = ft == D_stuff ? 2 : (1 - 2 * cmp) * m;

i.e. it was previously only correct for m = +1, not m = -1. However, since the error was only at r=0, it would still converge to the correct answer for m = -1, just with a first-order error.

(Note that this affects both the PML and non-PML regions, regardless of the location of the source. It may have contributed to the degradation of the PML that we previously observed in #2148?)

[oskooi]

In particular, this might be a bug in the m=1 code at r=0 here which looks like should include a * m factor

realnum f_m_mult = ft == D_stuff ? 2 : (1 - 2 * cmp) * m;

Looks like this was indeed the bug. With this fix, the results for the m=±1 simulations are now consistent.

Here is a summary of the updated results using this bug fix.

From the equations, if you have a purely real Er source in equation (4), you will get purely real {Er, Ez, Hφ} and purely imaginary {Eφ, Hr, Hz}.

Confirmed.

If you flip the sign of m, then it should complex-conjugate the solution, which flips the signs of {Eφ, Hr, Hz}.

Confirmed.

Also, the first experiment now produces the same flux values for all four test cases which is similar to m = ±2, ±3, ...:

m	fcen	flux
-1	+1	A
-1	-1	A
+1	+1	A
+1	-1	A

These results can be the basis for a unit test for this feature which currently does not exist.

Finally, the bug fix now resolves the slowly decaying fields at $r=0$ reported in #2148 which means it is no longer necessary to increase the cutoff value of the GaussianSource. We will therefore also need to update Tutorial/Extraction Efficiency of a Light-Emitting Diode (LED).

[stevengj]

Note that if you vary the radial length L, the (time-averaged/frequency-domain) flux in the R and Z directions should go to nonzero constants since you are exciting a guided mode trapped in the waveguide.

(Since you have PEC boundary conditions at the bottom, this corresponds to the odd-symmetry modes of the waveguide. For the in-plane polarization, that includes the fundamental mode, so I would expect waveguide modes to exist regardless of the thickness — there is no cutoff.)

As to whether the R or Z flux is bigger, that depends. On the one hand, the waveguide modes have a large density of states than vacuum so point sources tend to preferentially radiate into them. On the other hand, you would couple better into the guided mode using a z-directed point source than an r-directed point source.

[stevengj]

Note also that you should only expect Poynting's theorem to hold up to discretization error, because we are using the Poynting flux ExH derived for the exact Maxwell equations, which is not quite the same as the discrete Poynting flux.

[oskooi]

Note also that you should only expect Poynting's theorem to hold up to discretization error

Unfortunately, this discretization error does not seem to be the cause of the discrepancy in the results for the different $L$ calculations. This is because doubling the grid resolution to 100 pixels/μm and then again to 200 pixels/μm still produces flux values which depend on $L$. See the two plots in the figures below.

I also tried using accurate_fields_near_cylorigin=True as an alternative to the default heuristic approach of manually zeroing the field values near $r=0$. This also did not change the results.

In summary: even though we have fixed a couple of bugs in the $z$-PML at $r=0$ for $m=0, \pm 1$ (#2382, #2383), we are still not able to demonstrate that the total radiated flux for a point source at $r &gt; 0$ normalized by $2\pi r^2$ is a constant independent of $r$ and $L$.

[oskooi]

On 12/8/2022, @stevengj wrote:

I think the correct factor for LDOS is probably ΔV = Δr * Δz * 2πr, where for r=0 you treat r as Δr/2. This is where you get the Δr * Δz * 2πΔr/2 = π/resolution^3 factor in the r=0 LDOS.

Here, however, you aren't computing LDOS, you are computing total power, so I think maybe you're getting a missing 2πr factor from the current source normalization (since we normalize the currents as a delta function but IIRC we don't include the 2πr from a δ function in cylindrical coordinates?).

As specified in fields::add_volume_source, the normalization for a δ-function source in cylindrical coordinates is Δr * Δz:

meep/src/sources.cpp

Lines 476 to 479 in a6a1192

	LOOP_OVER_DIRECTIONS(gv.dim, d) {
	if (where.in_direction(d) == 0.0 && !nosize_direction(d)) // delta-fun
	data.amp *= gv.a; // correct units for J delta-function amplitude
	}

(The if statement evaluates to true for directions R and Z which each add a factor of Δr or Δz, respectively, where Δr = Δz.)

There is no 2πr term, as you suggested. The normalization of the δ-function source is identical in 2d and cylindrical. Could this be a bug? It seems the δ-function source in cylindrical coordinates should be handled as a special case which it is currently not.

[stevengj]

My suggestion is to remove the waveguide (it slows down convergence with resolution) for now, fix the source position (at r ≠ 0), and plot the power P(L) - P(L+1) vs resolution. What we'd like to see is that this difference goes to zero. (You want a bunch of points here because as you know sometimes the convergence can be a bit oscillatory, so two data points might be deceptive.)

If you want you can fix m as well, since Poynting's theorem should hold for every m individually.

(If the difference stops decreasing with resolution at some point, check the runtime and the PML thickness.)

[oskooi]

My suggestion is to remove the waveguide (it slows down convergence with resolution) for now, fix the source position (at r ≠ 0), and plot the power P(L) - P(L+1) vs resolution. What we'd like to see is that this difference goes to zero.

I performed a series of tests for a single m as well as a sum over all ms and it seems that the difference P(2L) - P(L) for a source at $r=0.5$ does indeed converge to zero quadratically with resolution for different L (as expected).

Additionally, here are two plots showing the radiated flux for m = 1 (top) and sum over all ms (bottom) vs. the source position at a resolution of 200 (the largest value that I have simulated). The radiated flux has not been normalized by $(2\pi r)^2$ as in previous plots. The right inset in each plot shows the relative error defined by $(\max(P(L)) - \min(P(L)) / mean(P(L))$ for $L = 5, 10, 20$ at each source position. Note that the error changes with source position but is always less than ~3%.

[oskooi]

Putting all the recent results and bug fixes together, it seems we could be getting closer to wrapping this up.

If we plot $P(r)$ vs. $r$ for $r &gt; 0$ on a log-log scale (where $P(r)$ is obtained by summing over ms), we can clearly see the $r^2$ dependence:

Taking the above data and normalizing $P(r)$ by $r^2$, we don't quite see a constant value as expected but the $r$-dependence of the data is small (less than ~5%) plotted on a linear scale:

Based on this data, let's specify the normalized flux at all $r &gt; 0$ to be $P(r) / r^2 \approx 7600$. We need to determine a scale factor $C$ such that $C \times P(r=0) = P(r) / r^2$. From the table below, $P(r=0) \approx 0.1$. Putting this all together yields $C \approx \pi \times (\pi / \Delta r)^2$ where $\Delta r = \Delta z = 1 / resolution$. In this example, the resolution is 50.

$L$	flux($r=0$), m=+1
5	0.09769383
10	0.09739900
20	0.09805311

[stevengj]

Looks like the L dependence is conclusively a resolution effect. Is it converging with resolution or resolution^2?

Similarly, would be good to plot the flux(2r) - flux(r) (normalized), for a fixed L, vs resolution, to see if that converges as well.

[oskooi]

Looks like the L dependence is conclusively a resolution effect. Is it converging with resolution or resolution^2?

P(2L) - P(L) is converging to zero quadratically with resolution (as expected). I have updated the log-log plot in the comment above to include the quadratic reference $1/(resolution)^2$ in order to demonstrate this scaling clearly.

Also, based on an empirical analysis of the results involving the radiated flux $P$ for a dipole at $r=0$ and $r &gt; 0$, the scale factor $C$ in $C \times P(r=0) = P(r) / r^2$ is $C \approx \pi \times (\pi / \Delta r)^2$ where $\Delta r = \Delta z = 1 / resolution$. This analytic expression is what we need to complete the tutorial although it would be good to explain where it comes from.

[oskooi]

The constant factor $C$ necessary to convert the radiated flux $P_{cyl}$ at $r=0$ in cylindrical coordinates to the radiated flux $P_{3d}$ at $(x,y) = (0,0)$ in 3d is $C = 2 / (\pi \Delta r)^2$ where $\Delta r = \Delta z = 1/resolution$. That is, $C \times P_{cyl} = P_{3d}$. To verify this relationship empirically, here are the results for dipole emission in vacuum at four different resolutions.

resolution	$P_{cyl}$	$P_{3d}$	$F = P_{3d} / P_{cyl}$	$C = 2 / (\pi \Delta r)^2$	$error = (F - C) / C$
25	0.40757553	52.345147	128.4305439	126.6514795	0.01404693
50	0.10293927	52.355973	508.6103000	506.6059182	0.00395649
100	0.02581597	52.358913	2028.1598174	2026.4236728	0.00085675
200	0.00646006	52.359600	8105.1316045	8105.6946914	0.00006947

Based on the results in this table, a plot of the relative error in $C$ (last/sixth column) vs. resolution (first column) shows a ~quadratic scaling, as expected.

Schematic of cell layout in cylindrical coordinates:

scripts used to generate results:
freespace_3d_flux.py
freespace_cyl_flux_r0.py

[oskooi]

Similarly, would be good to plot the flux(2r) - flux(r) (normalized), for a fixed L, vs resolution, to see if that converges as well.

The quantity $P(r)/r^2 - P(2r)/(2r)^2$ [where $P(r)$ is the radiated flux (summed over ms) for a point source at $r$] does seem to be decreasing monotonically to zero with increasing resolution. The plot shows this trend for three sets of data.

However, based on this data, the convergence rate seems to be much smaller than $1/resolution^2$. The slow convergence is troubling for the purpose of these types of calculations where the normalized radiated flux is expected to be independent of $r$ at modest grid resolutions.

Here is a plot for just the $r = 0.05$ case.

[oskooi]

Looks like $P(r)/r^2$ approaches a constant only for $r \rightarrow \infty$. I simply needed to increase the maximum $r$ in my plots from 4 to 8 to see this trend clearly (see figures below). Using the aymptotic value for $P(r)/r^2$ in the empirical analysis (rather than a best guess), the scale factor $C$ in $C \times P(0) = P(r) / r^2$ is exactly $C = \pi \times (\pi / \Delta r)^2$. It is not clear why $P(r)/r^2$ is not a constant for $r \lesssim 5$.

Putting all the results together, to convert the Poynting flux $P(r)$ in cylindrical coordinates to 3d requires: $P_{3d} = P(r) / r^2 \times 2 / \pi^5$.

[stevengj]

Would be nice to do a tutorial for the other case mentioned at the beginning, an off-axis planewave.