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

python/examples/dipole_in_vacuum_1D.py

@@ -0,0 +1,357 @@
+"""Radiation pattern of a dipole in vacuum obtained using a 1D calculation."""
+
+from enum import Enum
+from typing import Tuple
+import math
+import warnings
+
+import meep as mp
+import numpy as np
+
+
+Current = Enum("Current", "Jx Jy Kx Ky")
+
+RESOLUTION_UM = 25
+WAVELENGTH_UM = 1.0
+NUM_K = 40
+NUM_POLAR = 50
+NUM_AZIMUTHAL = 50
+PML_UM = 1.0
+AIR_UM = 10.0
+DEBUG_OUTPUT = False
+
+frequency = 1 / WAVELENGTH_UM
+
+
+def dipole_in_vacuum(
+    dipole_component: Current, kx: float, ky: float
+) -> Tuple[complex, complex, complex, complex]:
+    """                                                                                               
+    Returns the near fields of a dipole in vacuum.                                                    
+
+    Args:                                                                                             
+        dipole_component: the component of the dipole.                                                
+        kx, ky: the wavevector components.                                                            
+
+    Returns:                                                                                          
+        The surface-tangential electric and magnetic near fields as a 4-tuple.                        
+    """
+    size_z_um = PML_UM + AIR_UM + PML_UM
+    cell_size = mp.Vector3(0, 0, size_z_um)
+    pml_layers = [mp.PML(thickness=PML_UM)]
+
+    if dipole_component.name == "Jx":
+        src_cmpt = mp.Ex
+    elif dipole_component.name == "Jy":
+        src_cmpt = mp.Ey
+    elif dipole_component.name == "Kx":
+        src_cmpt = mp.Hx
+    elif dipole_component.name == "Ky":
+        src_cmpt = mp.Hy
+
+    sources = [
+        mp.Source(
+            src=mp.GaussianSource(frequency, fwidth=0.05 * frequency),
+            component=src_cmpt,
+            center=mp.Vector3(0, 0, -0.5 * AIR_UM),
+        )
+    ]
+
+    sim = mp.Simulation(
+        resolution=RESOLUTION_UM,
+        force_complex_fields=True,
+        cell_size=cell_size,
+        sources=sources,
+        boundary_layers=pml_layers,
+        k_point=mp.Vector3(kx, ky, 0),
+    )
+
+    dft_fields = sim.add_dft_fields(
+        [mp.Ex, mp.Ey, mp.Hx, mp.Hy],
+        frequency,
+        0,
+        1,
+        center=mp.Vector3(0, 0, 0.5 * AIR_UM),
+        size=mp.Vector3()
+    )
+
+    sim.run(
+        until_after_sources=mp.stop_when_fields_decayed(
+            10, src_cmpt, mp.Vector3(0, 0, 0.5 * AIR_UM), 1e-6
+        )
+    )
+
+    ex_dft = sim.get_dft_array(dft_fields, mp.Ex, 0)
+    ey_dft = sim.get_dft_array(dft_fields, mp.Ey, 0)
+    hx_dft = sim.get_dft_array(dft_fields, mp.Hx, 0)
+    hy_dft = sim.get_dft_array(dft_fields, mp.Hy, 0)
+
+    return ex_dft, ey_dft, hx_dft, hy_dft
+
+
+def equivalent_currents(
+    ex_dft: complex, ey_dft: complex, hx_dft: complex, hy_dft: complex
+) -> Tuple[Tuple[complex, complex], Tuple[complex, complex]]:
+    """Computes the equivalent electric and magnetic currents on a surface.                           
+
+    Args:                                                                                             
+        ex_dft, ey_dft, hx_dft, hy_dft: the surface tangential DFT fields.                            
+
+    Returns:                                                                                          
+        A 2-tuple of the electric and magnetic sheet currents in x and y.                             
+    """
+
+    electric_current = (-hy_dft, hx_dft)
+    magnetic_current = (ey_dft, -ex_dft)
+
+    return (electric_current, magnetic_current)
+
+
+def farfield_amplitudes(
+    kx: float,
+    kz: float,
+    current_amplitude: complex,
+    current_component: Current,
+) -> Tuple[complex, complex, complex]:
+    """Computes the S- or P-polarized far-field amplitudes from a sheet current.                      
+
+    Assumes ky=0 such that wavevector is in xz plane.                                                 
+
+    Args:                                                                                             
+        kx, kz: wavevector of the outgoing planewave in the x,z direction. Units of 2π.               
+        current_amplitude: amplitude of the sheet current.                                            
+        current_component: component of the sheet current.                                            
+
+    Returns:                                                                                          
+        A 3-tuple of the electric and magnetic far fields.                                            
+    """
+    if current_component.name == "Jx":
+        # Jx --> (Ex, Hy, Ez) [P pol.]                                                                
+        ex0 = -kz * current_amplitude / (2 * frequency)
+        hy0 = frequency * ex0 / kz
+    elif current_component.name == "Jy":
+        # Jy --> (Hx, Ey, Hz) [S pol.]                                                                
+        ey0 = -frequency * current_amplitude / (2 * kz)
+        hx0 = -kz * ey0 / frequency
+    elif current_component.name == "Kx":
+        # Kx --> (Hx, Ey, Hz) [S pol.]                                                                
+        hx0 = kz * current_amplitude / (2 * frequency)
+        ey0 = -frequency * hx0 / kz
+    elif current_component.name == "Ky":
+        # Ky --> (Ex, Hy, Ez) [P pol.]                                                                
+        hy0 = frequency * current_amplitude / (2 * kz)
+        ex0 = kz * hy0 / frequency
+
+    if current_component.name == "Jx" or current_component.name == "Ky":
+        ez0 = -kx * hy0 / frequency
+        return (ex0, hy0, ez0)
+    elif current_component.name == "Jy" or current_component.name == "Kx":
+        hz0 = kx * ey0 / frequency
+        return (hx0, ey0, hz0)
+
+
+def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, float]:
+    """Converts a far point in spherical to Cartesian coordinates.                                    
+
+    Args:                                                                                             
+        polar_rad: polar angle of the point.                                                          
+        azimuthal_rad: azimuthal angle of the point.                                                  
+
+    Returns:                                                                                          
+        The x,y,z coordinates of the far point as a 3-tuple.                                          
+    """
+    x = np.sin(polar_rad) * np.cos(azimuthal_rad)
+    y = np.sin(polar_rad) * np.sin(azimuthal_rad)
+    z = np.cos(polar_rad)
+
+    return (x, y, z)
+
+
+if __name__ == "__main__":
+    dipole_component = Current.Jx
+
+    kxs = np.linspace(-frequency, frequency, NUM_K)
+    kys = np.linspace(-frequency, frequency, NUM_K)
+
+    # Far fields are defined on the surface of a hemisphere.                                          
+    polar_rad = np.linspace(
+        0,
+        0.5 * math.pi,
+        NUM_POLAR
+    )
+    azimuthal_rad = np.linspace(0, 2 * np.pi, NUM_AZIMUTHAL)
+    delta_azimuthal_rad = 2 * np.pi / (NUM_AZIMUTHAL - 1)
+
+    current_components = [Current.Jx, Current.Jy, Current.Kx, Current.Ky]
+    farfield_components = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"]
+
+    poynting_flux = np.zeros((NUM_POLAR, NUM_AZIMUTHAL))
+
+    # Brillouin-zone integration over 2D grid of (kx, ky).                                            
+    for kx in kxs:
+        for ky in kys:
+
+            # Skip wavevectors which are outside the light cone.                                      
+            if np.sqrt(kx**2 + ky**2) >= (0.945 * frequency):
+                continue
+
+            kz = (frequency**2 - kx**2 - ky**2) ** 0.5
+
+            if DEBUG_OUTPUT:
+                print(f"(kx, ky, kz) = ({kx:.6f}, {ky:.6f}, {kz:.6f})")
+
+            ex_dft, ey_dft, hx_dft, hy_dft = dipole_in_vacuum(dipole_component, kx, ky)
+
+            if DEBUG_OUTPUT:
+                print(f"ex_dft:, {ex_dft}")
+                print(f"ey_dft:, {ey_dft}")
+                print(f"hx_dft:, {hx_dft}")
+                print(f"hy_dft:, {hy_dft}")
+
+            # Rotation angle around z axis to force ky=0. 0 radians is +x.                            
+            rotation_rad = -math.atan2(ky, kx)
+
+            if DEBUG_OUTPUT:
+                print(f"rotation angle:, {math.degrees(rotation_rad):.2f}°")
+
+            rotation_matrix = np.array(
+                [
+                    [np.cos(rotation_rad), -np.sin(rotation_rad)],
+                    [np.sin(rotation_rad), np.cos(rotation_rad)],
+                ]
+            )
+
+            k_rotated = rotation_matrix @ np.transpose(np.array([kx, ky]))
+            if k_rotated[1] > 1e-10:
+                raise ValueError(f"rotated ky is nonzero: {k_rotated[1]}")
+
+            electric_current, magnetic_current = equivalent_currents(
+                ex_dft, ey_dft, hx_dft, hy_dft
+            )
+
+            electric_current_rotated = rotation_matrix @ np.transpose(
+                np.array([electric_current[0], electric_current[1]])
+            )
+
+            magnetic_current_rotated = rotation_matrix @ np.transpose(
+                np.array([magnetic_current[0], magnetic_current[1]])
+            )
+
+            if DEBUG_OUTPUT:
+                print(f"electric_current:, {electric_current}")
+                print(f"magnetic_current:, {magnetic_current}")
+                print(f"electric_current_rotated:, {electric_current_rotated}")
+                print(f"magnetic_current_rotated:, {magnetic_current_rotated}")
+
+            current_amplitudes = [
+                electric_current_rotated[0],
+                electric_current_rotated[1],
+                magnetic_current_rotated[0],
+                magnetic_current_rotated[1],
+            ]
+
+            # Obtain the far fields over a hemisphere given a sheet current                           
+            # with linear in-plane polarization.                                                      
+
+            farfields = {}
+            for component in farfield_components:
+                farfields[component] = np.zeros(
+                    (NUM_POLAR, NUM_AZIMUTHAL),
+                    dtype=np.complex128
+                )
+
+            for i in range(NUM_POLAR):
+                for j in range(NUM_AZIMUTHAL):
+
+                    # Obtain the farfields for each type of sheet current                             
+                    # and keep a running sum.                                                         
+                    for current_component, current_amplitude in zip(
+                        current_components, current_amplitudes
+                    ):
+                        if abs(current_amplitude) == 0:
+                            continue
+
+                        farfield_pol = farfield_amplitudes(
+                            k_rotated[0],
+                            kz,
+                            current_amplitude,
+                            current_component,
+                        )
+
+                        if (
+                            current_component.name == "Jx" or
+                            current_component.name == "Ky"
+                        ):
+                            # P polarization                                                          
+                            farfields["Ex"][i, j] += farfield_pol[0]
+                            farfields["Hy"][i, j] += farfield_pol[1]
+                            farfields["Ez"][i, j] += farfield_pol[2]
+                        elif (
+                            current_component.name == "Jy" or
+                            current_component.name == "Kx"
+                        ):
+                            # S polarization                                                          
+                            farfields["Hx"][i, j] += farfield_pol[0]
+                            farfields["Ey"][i, j] += farfield_pol[1]
+                            farfields["Hz"][i, j] += farfield_pol[2]
+
+
+            inverse_rotation_matrix = np.transpose(rotation_matrix)
+
+            for i in range(NUM_POLAR):
+                for j in range(NUM_AZIMUTHAL):
+                    ex_rotated = (
+                        inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] +
+                        inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j]
+                    )
+                    ey_rotated = (
+                        inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] +
+                        inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j]
+                    )
+
+                    hx_rotated = (
+                        inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] +
+                        inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j]
+                    )
+                    hy_rotated = (
+                        inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] +
+                        inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j]
+                    )
+
+                    flux_x = np.real(
+                        np.conj(ey_rotated) * farfields["Hz"][i, j] -
+                        np.conj(farfields["Ez"][i, j]) * hy_rotated
+                    )
+
+                    flux_y = np.real(
+                        np.conj(farfields["Ez"][i, j]) * hx_rotated -
+                        np.conj(ex_rotated) * farfields["Hz"][i, j]
+                    )
+
+                    flux_z = np.real(
+                        np.conj(ex_rotated) * hy_rotated -
+                        np.conj(ey_rotated) * hx_rotated
+                    )
+
+                    rx, ry, rz = spherical_to_cartesian(polar_rad[i], azimuthal_rad[j])
+                    poynting_flux[i, j] += rx * flux_x + ry * flux_y + rz * flux_z
+
+
+    if mp.am_master():
+        np.savez(
+            "dipole_farfields.npz",
+            RESOLUTION_UM=RESOLUTION_UM,
+            WAVELENGTH_UM=WAVELENGTH_UM,
+            NUM_K=NUM_K,
+            NUM_POLAR=NUM_POLAR,
+            NUM_AZIMUTHAL=NUM_AZIMUTHAL,
+            PML_UM=PML_UM,
+            AIR_UM=AIR_UM,
+            polar_rad=polar_rad,
+            azimuthal_rad=azimuthal_rad,
+            kxs=kxs,
+            kys=kys,
+            dipole_component=dipole_component,
+            poynting_flux=poynting_flux
+        )
+