Tutorial example for radiation pattern of a dipole using 1D calculation and Brillouin-zone integration #2917

oskooi · 2024-10-03T06:28:10Z

Closes #2866.

This is still a work in progress. This implementation is based on the procedure described in #2866 (comment) for computing the radiation pattern of a dipole in vacuum using a 1D calculation and Brillouin-zone integration.

The radiation pattern of a dipole should be similar to a "donut" with rotational symmetry around the axis defined by the dipole orientation. However, as shown in the image below, the radiation pattern generated by this script is rotationally symmetric around the $z$ axis. This seems to be incorrect because the dipole is oriented in the $x$ direction.

For reference, the script used to generate the radiation pattern from the farfield data produced as the output from the script is provided below.

from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np


def spherical_to_cartesian(
	polar_rad,
	azimuthal_rad
) -> Tuple[float, float, float]:
    """Converts a point in spherical to Cartesian coordinates."""

    rx = np.sin(polar_rad) * np.cos(azimuthal_rad)
    ry = np.sin(polar_rad) * np.sin(azimuthal_rad)
    rz = np.cos(polar_rad)

    return rx, ry, rz

if __name__ == "__main__":
    data = np.load("dipole_farfields.npz")

    farfield_radius_um = data['FARFIELD_RADIUS_UM']
    polar_rad = data['polar_rad']
    azimuthal_rad = data['azimuthal_rad']
    num_polar = data['NUM_POLAR']
    num_azimuthal = data['NUM_AZIMUTHAL']

    flux_z = np.zeros((num_polar, num_azimuthal))
    xs = np.zeros((num_polar, num_azimuthal))
    ys = np.zeros((num_polar, num_azimuthal))
    zs = np.zeros((num_polar, num_azimuthal))

    for i in range(num_polar):
	for j in range(num_azimuthal):

            flux_z[i, j] = np.real(
		np.conj(data['Ex'][i, j]) * data['Hy'][i, j] -
		np.conj(data['Ey'][i, j]) * data['Hx'][i, j]
            )

            xs[i, j], ys[i, j], zs[i, j] = spherical_to_cartesian(
		polar_rad[i], azimuthal_rad[j]
            )

            xs[i, j] *= flux_z[i, j]
            ys[i, j] *= flux_z[i, j]
            zs[i, j] *= flux_z[i, j]

    fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(6, 6))
    ax.plot_surface(xs, ys, zs, cmap="inferno")
    ax.set_title("radiation pattern in 3d")
    ax.set_box_aspect((np.amax(xs), np.amax(ys), np.amax(zs)))
    ax.set_zlabel("z flux (a.u.)", labelpad=30)
    ax.tick_params(axis="z", which="major", pad=15)
    ax.set(xticklabels=[], yticklabels=[])

    fig.savefig(
	"radiation_pattern_3d.png",
	dpi=150,
	bbox_inches="tight",
    )

codecov-commenter · 2024-10-03T06:52:10Z

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 73.71%. Comparing base (f29a8c7) to head (af17883).
Report is 35 commits behind head on master.

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #2917      +/-   ##
==========================================
- Coverage   73.81%   73.71%   -0.10%     
==========================================
  Files          18       18              
  Lines        5423     5449      +26     
==========================================
+ Hits         4003     4017      +14     
- Misses       1420     1432      +12

stevengj · 2024-10-03T19:47:49Z

You need to take into account $k_x, k_y$ in the near-to-far calculation, since the current sheets are oscillating in $(x,y)$. (Note that this will also give you nonzero $z$ components in the far field.)

stevengj · 2024-10-03T19:52:06Z

(Note that you might run into an issue with the 1d PML, because when you are right on the light line the fields won't decay.)

oskooi · 2024-10-08T20:35:06Z

You need to take into account $k_x, k_y$ in the near-to-far calculation, since the current sheets are oscillating in $(x,y)$. (Note that this will also give you nonzero $z$ components in the far field.)

Here are the equations for the S- and P-polarized farfields given a current sheet oscillating in $(x,y)$.

oskooi · 2024-10-10T21:00:54Z

After updating the script to include the corrected formulas for obtaining the farfields given a current sheet, the farfields given an $x$-directed dipole is not rotationally symmetric but is still not quite "donut" like.

oskooi · 2024-10-10T21:25:14Z

from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np


def spherical_to_cartesian(
       polar_rad,
       azimuthal_rad
) -> Tuple[float, float, float]:
   """Converts a point in spherical to Cartesian coordinates."""

   rx = np.sin(polar_rad) * np.cos(azimuthal_rad)
   ry = np.sin(polar_rad) * np.sin(azimuthal_rad)
   rz = np.cos(polar_rad)

   return rx, ry, rz


if __name__ == "__main__":
   data = np.load("dipole_farfields.npz")

   farfield_radius_um = data['FARFIELD_RADIUS_UM']
   polar_rad = data['polar_rad']
   azimuthal_rad = data['azimuthal_rad']
   num_polar = data['NUM_POLAR']
   num_azimuthal = data['NUM_AZIMUTHAL']

   flux_x = np.zeros((num_polar, num_azimuthal))
   flux_y = np.zeros((num_polar, num_azimuthal))
   flux_z = np.zeros((num_polar, num_azimuthal))
   xs = np.zeros((num_polar, num_azimuthal))
   ys = np.zeros((num_polar, num_azimuthal))
   zs = np.zeros((num_polar, num_azimuthal))

   for i in range(num_polar):
       for j in range(num_azimuthal):

           flux_x[i, j] = np.real(
               np.conj(data['Ey'][i, j]) * data['Hz'][i, j] -
               np.conj(data['Ez'][i, j]) * data['Hy'][i, j]
           )

           flux_y[i, j] = np.real(
               np.conj(data['Ez'][i, j]) * data['Hx'][i, j] -
               np.conj(data['Ex'][i, j]) * data['Hz'][i, j]
           )

           flux_z[i, j] = np.real(
               np.conj(data['Ex'][i, j]) * data['Hy'][i, j] -
               np.conj(data['Ey'][i, j]) * data['Hx'][i, j]
           )

           flux_r = np.sqrt(flux_x[i, j]**2 + flux_y[i, j]**2 + flux_z[i, j]**2)

           rx, ry, rz = spherical_to_cartesian(polar_rad[i], azimuthal_rad[j])

           xs[i, j] = rx * flux_x[i, j]
           ys[i, j] = ry * flux_y[i, j]
           zs[i, j] = rz * flux_z[i, j]

   fig, ax = plt.subplots(subplot_kw={"projection": "3d"}, figsize=(6, 6))
   ax.plot_surface(xs, ys, zs, cmap="inferno")
   ax.set_title("radiation pattern in 3d")
   ax.set_box_aspect((np.amax(xs), np.amax(ys), np.amax(zs)))
   ax.set_zlabel("z flux (a.u.)", labelpad=30)
   ax.tick_params(axis="z", which="major", pad=15)
   ax.set(xticklabels=[], yticklabels=[])

   fig.savefig(
       "radiation_pattern_3d.png",
       dpi=150,
       bbox_inches="tight",
   )

stevengj · 2024-10-10T21:36:45Z

I think you want use

xs[i, j] = rx
ys[i, j] = ry
zs[i, j] = rz
flux_r[i,j] = rx * flux_x[i, j] + ry * flux_y[i, j] + rz * flux_z[i, j]

with the surface given by (xs,ys,zs) but the color given by flux_r. Here is an example in matplotlib — you have to use the facecolors argument to plot_surface. Something like:

c = matplotlib.cm.magma(flux_r)
plot_surface(xs,ys,zs, facecolors=c)

or just

pcolor(xs,ys, flux_r, cmap="magma")

…nd add missing rotation of farfield points

oskooi · 2024-10-15T02:53:45Z

I fixed a couple of bugs in the calculation of the far fields. However, even with these fixes the radiation pattern for $\phi = 0$ (i.e., a slice of $\phi$ in the range $[0, 2\pi]$) is quite different than the expected $\cos^2(\theta)$. (A fixed $\phi$ is intended to simplify the analysis from a 3D contour which can be difficult to visualize correctly to a 2D polar plot.) The script used to generate the 2D radiation pattern (polar plot of the radial flux) using the simulation data is provided below. For reference, the radiation pattern for a dipole in 3D resembling a "donut" is shown as a separate figure below.

At this point, it might be worth investigating whether setting up the 3D simulation using a 1D grid itself is bug free (e.g., see #2916).

import math
from typing import Tuple

import matplotlib
import matplotlib.pyplot as plt
import numpy as np


if __name__ == "__main__":
    data = np.load("dipole_farfields.npz")

    polar_rad = data['polar_rad']
    num_polar = data['NUM_POLAR']
    flux_x = np.zeros(num_polar)
    flux_z = np.zeros(num_polar)

    for i in range(num_polar):
        flux_x[i] = np.real(
            np.conj(data['Ey'][i, 0]) * data['Hz'][i, 0] -
            np.conj(data['Ez'][i, 0]) * data['Hy'][i, 0]
	)

	flux_z[i] = np.real(
            np.conj(data['Ex'][i, 0]) * data['Hy'][i, 0] -
            np.conj(data['Ey'][i, 0]) * data['Hx'][i, 0]
        )

    flux_r = np.sqrt(np.square(flux_x) + np.square(flux_z))

    np.savez(
        "dipole_radiation_pattern.npz",
        flux_r=flux_r,
        polar_rad=polar_rad,
    )

    fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6))
    ax.plot(polar_rad, flux_r / np.max(flux_r), 'b-', label="Meep")
    ax.plot(polar_rad, np.square(np.cos(polar_rad)), 'r--', label="cos$^2$(θ)")
    ax.set_theta_direction(-1)
    ax.set_theta_offset(0.5 * math.pi)
    ax.set_thetalim(0, 0.5 * math.pi)
    ax.set_rmax(1)
    ax.set_rticks([0,0.5,1])
    ax.grid(True)
    ax.set_rlabel_position(22)
    ax.legend()

    fig.savefig(
        "radiation_pattern_jx.png",
        dpi=150,
        bbox_inches="tight",
    )

python/examples/dipole_in_vacuum_1D.py

…nsure fields decay away

stevengj · 2024-11-05T18:54:40Z

python/examples/dipole_in_vacuum_1D.py

+    """
+    if current_component.name == "Jx":
+        # Jx --> (Ex, Hy, Ez) [P pol.]                                                                       
+        ex0 = frequency * current_amplitude / (2 * kz)


Suggested change

ex0 = frequency * current_amplitude / (2 * kz)

ex0 = kz * current_amplitude / (2 * frequency)

stevengj · 2024-11-05T18:59:37Z

python/examples/dipole_in_vacuum_1D.py

+            for i in range(NUM_POLAR):
+                for j in range(NUM_AZIMUTHAL):
+                    # Map the rotated point to the closest azimuthal angle.                                  
+                    azimuthal_index = (j + azimuthal_index_shift) % NUM_AZIMUTHAL


shifting the index seems wrong. the point index hasn't changed in this process, you only rotated the fields and r.

…FIELD_RADIUS_UM

oskooi · 2024-11-15T07:42:50Z

The results for the radiation pattern of a dipole in vacuum are improved after making these two modifications to the script:

Correcting the analytic formulas for the far fields from a sheet current. These errors were discovered rederiving the formulas.
Specifying a radius of 1 rather than 1e6 for the spherical surface for the far fields. It's not obvious why the radius of the far field surface should affet the results but it does so in a significant way.

The radiation pattern of the dipole computed using the Brillouin-zone integration is now similar to the analytic result. However, the results are not yet identical. The discrepancy should go away when using a finer and finer grid for the Brillouin-zone integration. This is something I am currently investigating.

stevengj · 2024-11-15T21:45:49Z

Probably for large R the far fields E(kx,ky) become very oscillatory functions of (kx,ky), so that your integral is not being computed accurately. Would be nice to plot E(kx,ky) for a few different R to confirm this.
For small R, the angle is not being computed correctly — the angle should be from the location of your dipole source to the far-field point, not from the near-field monitor to the far-field point.

add Python script for tutorial example

774bfa1

oskooi added the documentation label Oct 3, 2024

update script to include integration over kx and ky

73090bd

fix bug in calculation of phase in farfields due to missing ky term a…

27b91fd

…nd add missing rotation of farfield points

oskooi commented Oct 17, 2024

View reviewed changes

python/examples/dipole_in_vacuum_1D.py Outdated Show resolved Hide resolved

oskooi commented Oct 17, 2024

View reviewed changes

python/examples/dipole_in_vacuum_1D.py Outdated Show resolved Hide resolved

Ardavan Oskooi and others added 2 commits October 18, 2024 12:07

perform BZ integration over wavegrid grid just inside light cone to e…

4fb1ca9

…nsure fields decay away

correctly rotate the far fields back into the original coordinate frame

8ce770b

stevengj reviewed Nov 5, 2024

View reviewed changes

use corrected formulas for far fields of sheet current and remove FAR…

af17883

…FIELD_RADIUS_UM

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Tutorial example for radiation pattern of a dipole using 1D calculation and Brillouin-zone integration #2917

Tutorial example for radiation pattern of a dipole using 1D calculation and Brillouin-zone integration #2917

oskooi commented Oct 3, 2024

codecov-commenter commented Oct 3, 2024 •

edited

Loading

stevengj commented Oct 3, 2024

stevengj commented Oct 3, 2024

oskooi commented Oct 8, 2024

oskooi commented Oct 10, 2024

oskooi commented Oct 10, 2024

stevengj commented Oct 10, 2024 •

edited

Loading

oskooi commented Oct 15, 2024

stevengj Nov 5, 2024

stevengj Nov 5, 2024

oskooi commented Nov 15, 2024

stevengj commented Nov 15, 2024

	ex0 = frequency * current_amplitude / (2 * kz)
	ex0 = kz * current_amplitude / (2 * frequency)

Tutorial example for radiation pattern of a dipole using 1D calculation and Brillouin-zone integration #2917

Are you sure you want to change the base?

Tutorial example for radiation pattern of a dipole using 1D calculation and Brillouin-zone integration #2917

Conversation

oskooi commented Oct 3, 2024

codecov-commenter commented Oct 3, 2024 • edited Loading

Codecov Report

stevengj commented Oct 3, 2024

stevengj commented Oct 3, 2024

oskooi commented Oct 8, 2024

oskooi commented Oct 10, 2024

oskooi commented Oct 10, 2024

stevengj commented Oct 10, 2024 • edited Loading

oskooi commented Oct 15, 2024

stevengj Nov 5, 2024

Choose a reason for hiding this comment

stevengj Nov 5, 2024

Choose a reason for hiding this comment

oskooi commented Nov 15, 2024

stevengj commented Nov 15, 2024

codecov-commenter commented Oct 3, 2024 •

edited

Loading

stevengj commented Oct 10, 2024 •

edited

Loading