Tutorial for nonaxisymmetric point-dipole sources in cylindrical coordinates

doc/docs/Python_Tutorials/Cylindrical_Coordinates.md

@@ -585,3 +585,202 @@ Note that the volume specified in `get_farfields` via `center` and `size` is in
 Shown below is the far-field energy-density profile around the focal length for both the *r* and *z* coordinate directions for three lens designs with $N$ of 25, 50, and 100. The focus becomes sharper with increasing $N$ due to the enhanced constructive interference of the diffracted beam. As the number of zones $N$ increases, the size of the focal spot (full width at half maximum) at $z = 200$ μm decreases as $1/\sqrt{N}$ (see eq. 17 of the [reference](http://zoneplate.lbl.gov/theory)). This means that doubling the resolution (halving the spot width) requires quadrupling the number of zones.
 
 ![](../images/zone_plate_farfield.png#center)
+
+Nonaxisymmetric Dipole Sources
+-------------------------------
+
+In [Tutorial/Local Density of States/Extraction Efficiency of a Light-Emitting Diode (LED)](Local_Density_of_States.md#extraction-efficiency-of-a-light-emitting-diode-led), the extraction efficiency of an LED was computed using an axisymmetric point-dipole source at $r = 0$. This involved a single simulation with $m = \pm 1$. Simulating a point-dipole source at $r > 0$ (as shown in the schematic below) is more challenging because it is nonaxisymmetric: any point source at $r > 0$ is equivalent to a ring source in cylindrical coordinates. In order to simulate a point-dipole source at $r > 0$, it is necessary to represent the $\delta(r)$ (Dirac delta function) current source as a Fourier-series expansion of $\phi$ (the field angular dependence $\exp(im\phi)$). This computational approach involves two parts: (1) performing a series of simulations for $m = 0, 1, 2, ..., M$ for some cutoff $M$ of the Fourier series (described below), and (2) because of power orthogonality, summing the results from each $m$-simulation in post processing where the $m \neq 0$ terms are multiplied by two to account for the $-m$ solutions.
+
+![](../images/cyl_nonaxisymmetric_source_layout.png#center)
+
+Two features of this method may provide a significant speedup compared to an identical 3d simulation:
+
+1. Convergence of the Fourier series may require only a small number ($M + 1$) of simulations. For a given source position $r$, $M \approx r \omega$ where $\omega$ is the angular frequency of the source within the source medium. For $m > M$, the field oscillations tend to be too rapid and the current source therefore cannot radiate any power into the far field. As an example, a point-dipole source with wavelength of $1.0 \mu m$ at a radial position of $r = 1.0 \mu m$ within a medium of $n = 2.4$ would require roughly $M=16$ simulations. (In practice, however, we usually truncate the Fourier-series expansion earlier: whenever the radiated flux at some $m$ has dropped to some small fraction of its maximum value in the summation.) The plot below shows the radiated flux vs. $m$ for three different source positions. Generally, the farther the point source is from $r = 0$, the more simulations are required for the Fourier-series summation to converge.
+
+2. Each $m$-simulation in the Fourier-series expansion is independent of the others. The simulations can therefore be executed simultaneously using an embarassingly parallel approach.
+
+![](../images/cyl_nonaxisymmetric_source_flux_vs_m.png#center)
+
+The Poynting flux computed in cylindrical coordinates can be converted to 3d Cartesian coordinates with the use of various scaling factors:
+
+- To convert the radiated flux from a point dipole at $r=0$ in cylindrical coordinates to 3d, the result must be multiplied by $2 / (\pi \Delta r)^2$ where $\Delta r$ is the grid size ($1/resolution$).
+
+- The relationship between the radiated flux $P(r)$ from a point dipole at $r = 0$ and $r > 0$ is $(\pi / (\Delta r)^2) P(0) = P(r) / (\pi r)^2$.
+
+- Combining these two sets of results, to convert the radiated flux at $r > 0$ in cylindrical coordinates to 3d, the result must be multiplied by $2 / (\pi^5 * r^2)$.
+
+As a demonstration, we compute the radiated flux from a point dipole at three different locations with $r > 0$ within the LED. These results are compared to the radiated flux from a point dipole at $r = 0$. The results are shown in the table below. At this resolution, the relative error is at most ~4%. The error decreases with increasing resolution.
+
+| `rpos` |  **flux** | **error** |
+|:------:|:---------:|:---------:|
+|    0   | 48.390787 |    N/A    |
+|   3.3  | 49.038307 |   0.013   |
+|   7.5  | 49.868279 |   0.031   |
+|  12.1  | 50.124093 |   0.036   |
+
+The simulation script is in [examples/noaxisym_point_dipole_cyl.py](https://github.com/NanoComp/meep/blob/master/python/examples/nonaxisym_point_dipole_cyl.py).
+
+```py
+import meep as mp
+import numpy as np
+
+
+resolution = 30  # pixels/μm
+n = 2.4  # refractive index of dielectric layer
+wvl = 1.0  # wavelength (in vacuum)
+fcen = 1 / wvl  # center frequency of source/monitor
+
+
+def radiated_flux_cyl(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 a fraction of dmat.
+       rpos: position of source in radial direction.
+       m: angular φ dependence of the fields exp(imφ).
+
+    Returns:
+       The radiated flux in cylindrical coordinates.
+    """
+    L = 20  # length of non-PML region in radial direction
+    dair = 1.0  # thickness of air padding
+    dpml = 1.0  # PML thickness
+    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=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
+            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_z = 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),
+        ),
+    )
+
+    flux_air_r = sim.add_flux(
+        fcen,
+        0,
+        1,
+        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(
+            50.0,
+            src_cmpt,
+            src_pt,
+            1e-8,
+        ),
+    )
+
+    flux_z = mp.get_fluxes(flux_air_z)[0]
+    flux_r = mp.get_fluxes(flux_air_r)[0]
+    flux_tot = flux_r + flux_z
+
+    return flux_tot
+
+
+def radiated_flux_3d(rpos: float) -> float:
+    """Computes the radiated flux in 3d Cartesian coordinates using a
+       point-dipole source in cylindrical coordinates.
+
+    Args:
+       rpos: position of source in radial direction.
+
+    Returns:
+       The radiated flux in 3d Cartesian coordinates.
+    """
+    layer_thickness = 0.7 * wvl / n
+    dipole_height = 0.5
+
+    if rpos == 0:
+        # r = 0 source requires a single simulation with m = ±1
+        m = 1
+        flux_cyl = radiated_flux_cyl(
+            layer_thickness,
+            dipole_height,
+            rpos,
+            m,
+        )
+
+        flux_3d = flux_cyl * 2 * (resolution / np.pi) ** 2
+
+    else:
+        # r > 0 source requires Fourier-series expansion of φ
+        flux_tol = 1e-6  # relative tolerance for flux for truncating expansion
+        cutoff_M = int(2 * rpos * 2 * np.pi * fcen * n)
+        ms = range(cutoff_M + 1)
+        flux_cyl_tot = 0
+        flux_cyl_max = 0
+        for m in ms:
+            flux_cyl = radiated_flux_cyl(
+                layer_thickness,
+                dipole_height,
+                rpos,
+                m,
+            )
+            print(f"flux-m:, {rpos}, {m}, {flux_cyl:.6f}")
+            flux_cyl_tot += flux_cyl if m == 0 else 2 * flux_cyl
+            if flux_cyl > flux_cyl_max:
+                flux_cyl_max = flux_cyl
+            if m > 0 and (flux_cyl / flux_cyl_max) < flux_tol:
+                break
+
+        flux_3d = (flux_cyl_tot / rpos**2) * (2 / np.pi**5)
+
+    print(f"flux-3d:, {rpos:.2f}, {flux_3d:.6f}")
+
+    return flux_3d
+
+
+if __name__ == "__main__":
+    # reference result for dipole at r = 0
+    Pcyl_ref = radiated_flux_3d(0)
+
+    rs = [3.3, 7.5, 12.1]
+    for r in rs:
+        Pcyl = radiated_flux_3d(r)
+        err = abs(Pcyl - Pcyl_ref) / Pcyl_ref
+        print(f"err:, {r}, {Pcyl:.6f}, {err:.6f}")
+```

doc/docs/Python_Tutorials/Local_Density_of_States.md

@@ -327,7 +327,7 @@ A more efficient approach to computing the total emitted power from the dipole i
 
 To demonstrate this feature of the LDOS, we will compute the extraction efficiency of an LED-like structure consisting of a point dipole embedded within a dielectric layer ($n$=2.4 and thickness 0.7 $\lambda/n$) surrounded by vacuum and positioned above an infinite ground plane of a lossless metal. We will compute the extraction efficiency (at a wavelength $\lambda$ = 1.0 μm for a dipole with polarization parallel to the ground plane) as a function of the height of the dipole source above the ground plane using two different coordinate systems — 3D Cartesian and cylindrical — and demonstrate that the results are equivalent (apart from discretization error). The key difference is that simulation using cylindrical (2D) coordinates is significantly faster and more memory efficient than 3D.
 
-The simulation setup is shown in the figures below for 3D Cartesian (cross section in $xz$) and cylindrical coordinates. (Note that this single-dipole calculation differs from the somewhat related flux calculation in [Tutorials/Custom Source/Stochastic Dipole Emission in Light Emitting Diodes](Custom_Source.md#stochastic-dipole-emission-in-light-emitting-diodes) which involved modeling a *collection* of dipoles.)
+The simulation setup is shown in the figures below for 3D Cartesian (cross section in $xz$) and cylindrical coordinates. (Note that this single-dipole calculation differs from the somewhat related flux calculation in [Tutorials/Custom Source/Stochastic Dipole Emission in Light Emitting Diodes](Custom_Source.md#stochastic-dipole-emission-in-light-emitting-diodes) which involved modeling a *collection* of dipoles.) In this example, the point-dipole source is positioned at $r=0$ which is relatively straightforward to set up. Non-axisymmetric dipoles positioned at $r>0$ are more challenging because they involve combining the results of *multiple* simulations. For a demonstration, see [Cylindrical Coordinates/Non-axisymmetric Dipole Sources](Cylindrical_Cordinates.md#non-axisymmetric-dipole-sources).
 
 ![](../images/dipole_extraction_eff_3D.png#center)
 

python/examples/nonaxisym_point_dipole_cyl.py

@@ -0,0 +1,172 @@
+"""Tutorial example for non-axisymmetric point sources in cylindrical coords.
+
+This example demonstrates 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 anywhere along the radial
+direction.
+
+Reference: https://meep.readthedocs.io/en/latest/Python_Tutorials/Cylindrical_Coordinates/#nonaxisymmetric-dipole-sources
+"""
+
+import meep as mp
+import numpy as np
+
+
+resolution = 30  # pixels/μm
+n = 2.4  # refractive index of dielectric layer
+wvl = 1.0  # wavelength (in vacuum)
+fcen = 1 / wvl  # center frequency of source/monitor
+
+
+def radiated_flux_cyl(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 a fraction of dmat.
+       rpos: position of source in radial direction.
+       m: angular φ dependence of the fields exp(imφ).
+
+    Returns:
+       The radiated flux in cylindrical coordinates.
+    """
+    L = 20  # length of non-PML region in radial direction
+    dair = 1.0  # thickness of air padding
+    dpml = 1.0  # PML thickness
+    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=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
+            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_z = 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),
+        ),
+    )
+
+    flux_air_r = sim.add_flux(
+        fcen,
+        0,
+        1,
+        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(
+            50.0,
+            src_cmpt,
+            src_pt,
+            1e-8,
+        ),
+    )
+
+    flux_z = mp.get_fluxes(flux_air_z)[0]
+    flux_r = mp.get_fluxes(flux_air_r)[0]
+    flux_tot = flux_r + flux_z
+
+    return flux_tot
+
+
+def radiated_flux_3d(rpos: float) -> float:
+    """Computes the radiated flux in 3d Cartesian coordinates using a
+       point-dipole source in cylindrical coordinates.
+
+    Args:
+       rpos: position of source in radial direction.
+
+    Returns:
+       The radiated flux in 3d Cartesian coordinates.
+    """
+    layer_thickness = 0.7 * wvl / n
+    dipole_height = 0.5
+
+    if rpos == 0:
+        # r = 0 source requires a single simulation with m = ±1
+        m = 1
+        flux_cyl = radiated_flux_cyl(
+            layer_thickness,
+            dipole_height,
+            rpos,
+            m,
+        )
+
+        flux_3d = flux_cyl * 2 * (resolution / np.pi) ** 2
+
+    else:
+        # r > 0 source requires Fourier-series expansion of φ
+        flux_tol = 1e-6  # relative tolerance for flux for truncating expansion
+        cutoff_M = int(2 * rpos * 2 * np.pi * fcen * n)
+        ms = range(cutoff_M + 1)
+        flux_cyl_tot = 0
+        flux_cyl_max = 0
+        for m in ms:
+            flux_cyl = radiated_flux_cyl(
+                layer_thickness,
+                dipole_height,
+                rpos,
+                m,
+            )
+            print(f"flux-m:, {rpos}, {m}, {flux_cyl:.6f}")
+            flux_cyl_tot += flux_cyl if m == 0 else 2 * flux_cyl
+            if flux_cyl > flux_cyl_max:
+                flux_cyl_max = flux_cyl
+            if m > 0 and (flux_cyl / flux_cyl_max) < flux_tol:
+                break
+
+        flux_3d = (flux_cyl_tot / rpos**2) * (2 / np.pi**5)
+
+    print(f"flux-3d:, {rpos:.2f}, {flux_3d:.6f}")
+
+    return flux_3d
+
+
+if __name__ == "__main__":
+    # reference result for dipole at r = 0
+    Pcyl_ref = radiated_flux_3d(0)
+
+    rs = [3.3, 7.5, 12.1]
+    for r in rs:
+        Pcyl = radiated_flux_3d(r)
+        err = abs(Pcyl - Pcyl_ref) / Pcyl_ref
+        print(f"err:, {r}, {Pcyl:.6f}, {err:.6f}")