From 774bfa13b6b824ba9437f227c3293d6ba4372742 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Wed, 2 Oct 2024 22:49:23 -0700 Subject: [PATCH 1/7] add Python script for tutorial example --- python/examples/dipole_in_vacuum_1D.py | 328 +++++++++++++++++++++++++ 1 file changed, 328 insertions(+) create mode 100644 python/examples/dipole_in_vacuum_1D.py diff --git a/python/examples/dipole_in_vacuum_1D.py b/python/examples/dipole_in_vacuum_1D.py new file mode 100644 index 000000000..1c76a1c5f --- /dev/null +++ b/python/examples/dipole_in_vacuum_1D.py @@ -0,0 +1,328 @@ +"""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 = 35 +WAVELENGTH_UM = 1.0 +FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM +NUM_K = 11 +NUM_POLAR = 20 +NUM_AZIMUTHAL = 20 +DEBUG_OUTPUT = True + +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. + """ + pml_um = 1.0 + air_um = 10.0 + 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 + else: + src_cmpt = mp.Hy + + sources = [ + mp.Source( + src=mp.GaussianSource(frequency, fwidth=0.1 * frequency), + component=src_cmpt, + center=mp.Vector3(), + ) + ] + + 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 * size_z_um - pml_um), + size=mp.Vector3(), + ) + + sim.run( + until_after_sources=mp.stop_when_fields_decayed( + 10, src_cmpt, mp.Vector3(0, 0, 0.5423), 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 as 2-tuples. + """ + + electric_current = (-hy_dft, hx_dft) + magnetic_current = (ey_dft, -ex_dft) + + return electric_current, magnetic_current + + +def far_fields( + kz: float, + rz: float, + current_amplitude: complex, + current_component: Current, +) -> Tuple[complex, complex]: + """Computes the S- or P-polarized far fields from a sheet current. + + Args: + kz: wavevector of the outgoing planewave in the z direction. + rz: height of the far-field point on the surface of a hemisphere. + current_amplitude: amplitude of the sheet current. + current_component: component of the sheet current. + + Returns: + A 2-tuple of the electric and magnetic far fields. + """ + + if current_component.name == "Jx": + # Jx --> (Ex, Hy) [P pol.] + e_far = -0.5 * current_amplitude + h_far = -0.5 * current_amplitude + elif current_component.name == "Jy": + # Jy --> (Hx, Ey) [S pol.] + e_far = -0.5 * current_amplitude + h_far = 0.5 * current_amplitude + elif current_component.name == "Kx": + # Kx --> (Hx, Ey) [S pol.] + e_far = 0.5 * current_amplitude + h_far = -0.5 * current_amplitude + else: + # Ky --> (Ex, Hy) [P pol.] + e_far = 0.5 * current_amplitude + h_far = 0.5 * current_amplitude + + phase = np.exp(1j * kz * rz) + e_far *= phase + h_far *= phase + + return e_far, h_far + + +def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, float]: + """Converts a point in spherical to Cartesian coordinates.""" + + x = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.cos(azimuthal_rad) + y = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.sin(azimuthal_rad) + z = FARFIELD_RADIUS_UM * np.cos(polar_rad) + + return x, y, z + + +if __name__ == "__main__": + # Jx --> (Ex, Hy) [P pol.] + 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 * np.pi, NUM_POLAR) + azimuthal_rad = np.linspace(0, 2 * np.pi, NUM_AZIMUTHAL) + + current_components = [Current.Jx, Current.Jy, Current.Kx, Current.Ky] + farfield_components = ["Ex", "Ey", "Hx", "Hy"] + + total_farfields = {} + for component in farfield_components: + total_farfields[component] = np.zeros( + (NUM_POLAR, NUM_AZIMUTHAL), dtype=np.complex128 + ) + + for kx in kxs: + for ky in kys: + + # Skip wavevectors which are outside the light cone. + if np.sqrt(kx**2 + ky**2) >= 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 is +x. + if kx: + rotation_rad = -math.atan(ky / kx) + else: + if ky: + rotation_rad = -0.5 * np.pi + else: + rotation_rad = 0 + + 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)], + ] + ) + + electric_current, magnetic_current = equivalent_currents( + ex_dft, ey_dft, hx_dft, hy_dft + ) + + if DEBUG_OUTPUT: + print(f"electric_current:, {electric_current}") + print(f"magnetic_current:, {magnetic_current}") + + 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_rotated:, {electric_current_rotated}") + print(f"magnetic_current_rotated:, {magnetic_current_rotated}") + + # Verify that ky of the rotated wavevector is 0. + k_rotated = rotation_matrix @ np.transpose(np.array([kx, ky])) + print(f"k_rotated = ({k_rotated[0]:.2f}, {k_rotated[1]:.2f})") + if k_rotated[1] == 0: + print("rotated ky is zero.") + + current_amplitudes = [ + electric_current_rotated[0], + electric_current_rotated[1], + magnetic_current_rotated[0], + magnetic_current_rotated[1], + ] + + farfields = {} + for component in farfield_components: + farfields[component] = np.zeros( + (NUM_POLAR, NUM_AZIMUTHAL), dtype=np.complex128 + ) + + for i, polar in enumerate(polar_rad): + for j, azimuthal in enumerate(azimuthal_rad): + rx, ry, rz = spherical_to_cartesian(polar, azimuthal) + + for current_component, current_amplitude in zip( + current_components, current_amplitudes + ): + if abs(current_amplitude) == 0: + continue + + e_far, h_far = far_fields( + kz, + rz, + current_amplitude, + current_component, + ) + + if ( + current_component.name == "Jx" + or current_component.name == "Ky" + ): + # P polarization + farfields["Ex"][i, j] += e_far + farfields["Hy"][i, j] += h_far + elif ( + current_component.name == "Jy" + or current_component.name == "Kx" + ): + # S polarization + farfields["Ey"][i, j] += e_far + farfields["Hx"][i, j] += h_far + + inverse_rotation_matrix = np.linalg.inv(rotation_matrix) + + for i in range(NUM_POLAR): + for j in range(NUM_AZIMUTHAL): + total_farfields["Ex"][i, j] = ( + inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] + + inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] + ) + total_farfields["Ey"][i, j] = ( + inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] + + inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] + ) + total_farfields["Hx"][i, j] = ( + inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] + + inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] + ) + total_farfields["Hy"][i, j] = ( + inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] + + inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] + ) + + if mp.am_master(): + np.savez( + "dipole_farfields.npz", + **total_farfields, + RESOLUTION_UM=RESOLUTION_UM, + WAVELENGTH_UM=WAVELENGTH_UM, + FARFIELD_RADIUS_UM=FARFIELD_RADIUS_UM, + NUM_POLAR=NUM_POLAR, + NUM_AZIMUTHAL=NUM_AZIMUTHAL, + polar_rad=polar_rad, + azimuthal_rad=azimuthal_rad, + kxs=kxs, + kys=kys, + ) From 73090bd2b173a9fe1b3d2bebf30ba95d80c7169c Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 10 Oct 2024 13:59:14 -0700 Subject: [PATCH 2/7] update script to include integration over kx and ky --- python/examples/dipole_in_vacuum_1D.py | 170 ++++++++++++++----------- 1 file changed, 96 insertions(+), 74 deletions(-) diff --git a/python/examples/dipole_in_vacuum_1D.py b/python/examples/dipole_in_vacuum_1D.py index 1c76a1c5f..bf104caf3 100644 --- a/python/examples/dipole_in_vacuum_1D.py +++ b/python/examples/dipole_in_vacuum_1D.py @@ -25,15 +25,15 @@ 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. + """ + 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. """ pml_um = 1.0 air_um = 10.0 @@ -63,15 +63,15 @@ def dipole_in_vacuum( resolution=RESOLUTION_UM, force_complex_fields=True, cell_size=cell_size, - sources=sources, + 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], + [mp.Ex, mp.Ey, mp.Hx, mp.Hy], frequency, - 0, + 0, 1, center=mp.Vector3(0, 0, 0.5 * size_z_um - pml_um), size=mp.Vector3(), @@ -79,8 +79,8 @@ def dipole_in_vacuum( sim.run( until_after_sources=mp.stop_when_fields_decayed( - 10, src_cmpt, mp.Vector3(0, 0, 0.5423), 1e-6 - ) + 10, src_cmpt, mp.Vector3(0, 0, 0.5423), 1e-6 + ) ) ex_dft = sim.get_dft_array(dft_fields, mp.Ex, 0) @@ -94,13 +94,13 @@ def dipole_in_vacuum( 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 as 2-tuples. + """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 as 2-tuples. """ electric_current = (-hy_dft, hx_dft) @@ -110,45 +110,61 @@ def equivalent_currents( def far_fields( + kx: float, kz: float, + rx: float, rz: float, current_amplitude: complex, current_component: Current, -) -> Tuple[complex, complex]: - """Computes the S- or P-polarized far fields from a sheet current. - - Args: - kz: wavevector of the outgoing planewave in the z direction. - rz: height of the far-field point on the surface of a hemisphere. - current_amplitude: amplitude of the sheet current. - current_component: component of the sheet current. - - Returns: - A 2-tuple of the electric and magnetic far fields. +) -> Tuple[complex, complex, complex]: + """Computes the S- or P-polarized far fields from a sheet current. + + Args: + kx: wavevector of the outgoing planewave in the x direction. + kz: wavevector of the outgoing planewave in the z direction. + rx: x-coordinate of far-field point on the surface of a hemisphere. + rz: z-coordinate of far-field point on the surface of a hemisphere. + current_amplitude: amplitude of the sheet current. + current_component: component of the sheet current. + + Returns: + A 2-tuple of the electric and magnetic far fields. """ if current_component.name == "Jx": - # Jx --> (Ex, Hy) [P pol.] - e_far = -0.5 * current_amplitude - h_far = -0.5 * current_amplitude + # Jx --> (Ex, Hy) [P pol.] + ex0 = frequency * current_amplitude / (2 * kz) + hy0 = frequency * ex0 / kz + ez0 = -kx * hy0 / frequency elif current_component.name == "Jy": - # Jy --> (Hx, Ey) [S pol.] - e_far = -0.5 * current_amplitude - h_far = 0.5 * current_amplitude + # Jy --> (Hx, Ey) [S pol.] + ey0 = -frequency * current_amplitude / (2 * kz) + hx0 = -kz * ey0 / frequency + hz0 = kx * ey0 / frequency elif current_component.name == "Kx": - # Kx --> (Hx, Ey) [S pol.] - e_far = 0.5 * current_amplitude - h_far = -0.5 * current_amplitude + # Kx --> (Hx, Ey) [S pol.] + ey0 = -current_amplitude / 2 + hx0 = -kz * ey0 / frequency + hz0 = kx * ey0 / frequency + elif current_component.name == "Ky": + # Ky --> (Ex, Hy) [P pol.] + ex0 = current_amplitude / 2 + hy0 = frequency * ex0 / kz + ez0 = -kx * hy0 / frequency + + phase = np.exp(1j * (kx * rx + kz * rz)) + if current_component.name == "Jx" or current_component.name == "Ky": + ex = ex0 * phase + hy = hy0 * phase + ez = ez0 * phase + farfields = [ex, hy, ez] else: - # Ky --> (Ex, Hy) [P pol.] - e_far = 0.5 * current_amplitude - h_far = 0.5 * current_amplitude - - phase = np.exp(1j * kz * rz) - e_far *= phase - h_far *= phase + hx = hx0 * phase + ey = ey0 * phase + hz = hz0 * phase + farfields = [hx, ey, hz] - return e_far, h_far + return farfields def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, float]: @@ -162,18 +178,18 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa if __name__ == "__main__": - # Jx --> (Ex, Hy) [P pol.] + # Jx --> (Ex, Hy) [P pol.] 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. + # Far fields are defined on the surface of a hemisphere. polar_rad = np.linspace(0, 0.5 * np.pi, NUM_POLAR) azimuthal_rad = np.linspace(0, 2 * np.pi, NUM_AZIMUTHAL) current_components = [Current.Jx, Current.Jy, Current.Kx, Current.Ky] - farfield_components = ["Ex", "Ey", "Hx", "Hy"] + farfield_components = ["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"] total_farfields = {} for component in farfield_components: @@ -184,7 +200,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa for kx in kxs: for ky in kys: - # Skip wavevectors which are outside the light cone. + # Skip wavevectors which are outside the light cone. if np.sqrt(kx**2 + ky**2) >= frequency: continue @@ -201,7 +217,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa print(f"hx_dft:, {hx_dft}") print(f"hy_dft:, {hy_dft}") - # Rotation angle around z axis to force ky = 0. 0 is +x. + # Rotation angle around z axis to force ky = 0. 0 is +x. if kx: rotation_rad = -math.atan(ky / kx) else: @@ -240,7 +256,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa print(f"electric_current_rotated:, {electric_current_rotated}") print(f"magnetic_current_rotated:, {magnetic_current_rotated}") - # Verify that ky of the rotated wavevector is 0. + # Verify that ky of the rotated wavevector is 0. k_rotated = rotation_matrix @ np.transpose(np.array([kx, ky])) print(f"k_rotated = ({k_rotated[0]:.2f}, {k_rotated[1]:.2f})") if k_rotated[1] == 0: @@ -269,48 +285,54 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa if abs(current_amplitude) == 0: continue - e_far, h_far = far_fields( + farfield_pol = far_fields( + kx, kz, + rx, rz, current_amplitude, current_component, ) if ( - current_component.name == "Jx" - or current_component.name == "Ky" + current_component.name == "Jx" or + current_component.name == "Ky" ): - # P polarization - farfields["Ex"][i, j] += e_far - farfields["Hy"][i, j] += h_far + # 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" + current_component.name == "Jy" or + current_component.name == "Kx" ): - # S polarization - farfields["Ey"][i, j] += e_far - farfields["Hx"][i, j] += h_far + # 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.linalg.inv(rotation_matrix) for i in range(NUM_POLAR): for j in range(NUM_AZIMUTHAL): total_farfields["Ex"][i, j] = ( - inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] - + inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] + inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] + + inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] ) total_farfields["Ey"][i, j] = ( - inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] - + inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] + inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] + + inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] ) + total_farfields["Ez"][i, j] = farfields["Ez"][i, j] total_farfields["Hx"][i, j] = ( - inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] - + inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] + inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] + + inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] ) total_farfields["Hy"][i, j] = ( - inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] - + inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] + inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] + + inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] ) + total_farfields["Hz"][i, j] = farfields["Hz"][i, j] if mp.am_master(): np.savez( From 27b91fdbbce7c337777a363d881ef26a87b0f8cd Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Mon, 14 Oct 2024 19:37:30 -0700 Subject: [PATCH 3/7] fix bug in calculation of phase in farfields due to missing ky term and add missing rotation of farfield points --- python/examples/dipole_in_vacuum_1D.py | 213 +++++++++++++------------ 1 file changed, 108 insertions(+), 105 deletions(-) diff --git a/python/examples/dipole_in_vacuum_1D.py b/python/examples/dipole_in_vacuum_1D.py index bf104caf3..7aa14079d 100644 --- a/python/examples/dipole_in_vacuum_1D.py +++ b/python/examples/dipole_in_vacuum_1D.py @@ -14,9 +14,9 @@ RESOLUTION_UM = 35 WAVELENGTH_UM = 1.0 FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM -NUM_K = 11 -NUM_POLAR = 20 -NUM_AZIMUTHAL = 20 +NUM_K = 17 +NUM_POLAR = 100 +NUM_AZIMUTHAL = 50 DEBUG_OUTPUT = True frequency = 1 / WAVELENGTH_UM @@ -25,30 +25,29 @@ 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. + """ + 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. """ pml_um = 1.0 air_um = 10.0 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 + src_cmpt = mp.Ex elif dipole_component.name == "Jy": - src_cmpt = mp.Ey + src_cmpt = mp.Ey elif dipole_component.name == "Kx": src_cmpt = mp.Hx - else: + elif dipole_component.name == "Ky": src_cmpt = mp.Hy sources = [ @@ -63,24 +62,27 @@ def dipole_in_vacuum( resolution=RESOLUTION_UM, force_complex_fields=True, cell_size=cell_size, - sources=sources, + 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], + [mp.Ex, mp.Ey, mp.Hx, mp.Hy], frequency, - 0, + 0, 1, center=mp.Vector3(0, 0, 0.5 * size_z_um - pml_um), size=mp.Vector3(), ) + # TODO (oskooi): find out why meep.stop_when_fields_decayed fails to + # terminate simulation for certain values of kx and ky. sim.run( - until_after_sources=mp.stop_when_fields_decayed( - 10, src_cmpt, mp.Vector3(0, 0, 0.5423), 1e-6 - ) + until_after_sources=50 + # until_after_sources=mp.stop_when_fields_decayed( + # 10, src_cmpt, mp.Vector3(0, 0, 0.5423), 1e-6 + # ) ) ex_dft = sim.get_dft_array(dft_fields, mp.Ex, 0) @@ -94,22 +96,22 @@ def dipole_in_vacuum( 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 as 2-tuples. + """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 + return (electric_current, magnetic_current) -def far_fields( +def farfield_amplitudes( kx: float, kz: float, rx: float, @@ -117,74 +119,66 @@ def far_fields( current_amplitude: complex, current_component: Current, ) -> Tuple[complex, complex, complex]: - """Computes the S- or P-polarized far fields from a sheet current. - - Args: - kx: wavevector of the outgoing planewave in the x direction. - kz: wavevector of the outgoing planewave in the z direction. - rx: x-coordinate of far-field point on the surface of a hemisphere. - rz: z-coordinate of far-field point on the surface of a hemisphere. - current_amplitude: amplitude of the sheet current. - current_component: component of the sheet current. - - Returns: - A 2-tuple of the electric and magnetic far fields. + """Computes the S- or P-polarized far-field amplitudes from a sheet current. + + Based on the assumption that ky=0 such that the in-plane wavevector is xz. + + Args: + kx, kz: wavevector of the outgoing planewave in the x,z direction. + rx, rz: x,z coordinate of a point on the surface of a hemicircle. + current_amplitude: amplitude of the sheet current. + current_component: component of the sheet current. + + Returns: + A 3-tuple of the per-polarization electric and magnetic far fields. """ - if current_component.name == "Jx": - # Jx --> (Ex, Hy) [P pol.] + # Jx --> (Ex, Hy) [P pol.] ex0 = frequency * current_amplitude / (2 * kz) - hy0 = frequency * ex0 / kz - ez0 = -kx * hy0 / frequency elif current_component.name == "Jy": - # Jy --> (Hx, Ey) [S pol.] + # Jy --> (Hx, Ey) [S pol.] ey0 = -frequency * current_amplitude / (2 * kz) - hx0 = -kz * ey0 / frequency - hz0 = kx * ey0 / frequency elif current_component.name == "Kx": - # Kx --> (Hx, Ey) [S pol.] + # Kx --> (Hx, Ey) [S pol.] ey0 = -current_amplitude / 2 - hx0 = -kz * ey0 / frequency - hz0 = kx * ey0 / frequency elif current_component.name == "Ky": - # Ky --> (Ex, Hy) [P pol.] + # Ky --> (Ex, Hy) [P pol.] ex0 = current_amplitude / 2 - hy0 = frequency * ex0 / kz - ez0 = -kx * hy0 / frequency - phase = np.exp(1j * (kx * rx + kz * rz)) if current_component.name == "Jx" or current_component.name == "Ky": - ex = ex0 * phase - hy = hy0 * phase - ez = ez0 * phase - farfields = [ex, hy, ez] - else: - hx = hx0 * phase - ey = ey0 * phase - hz = hz0 * phase - farfields = [hx, ey, hz] - - return farfields + hy0 = frequency * ex0 / kz + ez0 = -kx * hy0 / frequency + return (ex0, hy0, ez0) + elif current_component.name == "Jy" or current_component.name == "Kx": + hx0 = -kz * ey0 / frequency + hz0 = kx * ey0 / frequency + return (hx0, ey0, hz0) def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, float]: - """Converts a point in spherical to Cartesian coordinates.""" - + """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 = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.cos(azimuthal_rad) y = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.sin(azimuthal_rad) z = FARFIELD_RADIUS_UM * np.cos(polar_rad) - return x, y, z + return (x, y, z) if __name__ == "__main__": - # Jx --> (Ex, Hy) [P pol.] 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. + # Far fields are defined on the surface of a hemisphere. polar_rad = np.linspace(0, 0.5 * np.pi, NUM_POLAR) azimuthal_rad = np.linspace(0, 2 * np.pi, NUM_AZIMUTHAL) @@ -194,13 +188,14 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa total_farfields = {} for component in farfield_components: total_farfields[component] = np.zeros( - (NUM_POLAR, NUM_AZIMUTHAL), dtype=np.complex128 + (NUM_POLAR, NUM_AZIMUTHAL), + dtype=np.complex128 ) for kx in kxs: for ky in kys: - # Skip wavevectors which are outside the light cone. + # Skip wavevectors which are outside the light cone. if np.sqrt(kx**2 + ky**2) >= frequency: continue @@ -217,7 +212,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa print(f"hx_dft:, {hx_dft}") print(f"hy_dft:, {hy_dft}") - # Rotation angle around z axis to force ky = 0. 0 is +x. + # Rotation angle around z axis to force ky=0. 0 radians is +x. if kx: rotation_rad = -math.atan(ky / kx) else: @@ -226,9 +221,6 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa else: rotation_rad = 0 - if DEBUG_OUTPUT: - print(f"rotation angle:, {math.degrees(rotation_rad):.2f}°") - rotation_matrix = np.array( [ [np.cos(rotation_rad), -np.sin(rotation_rad)], @@ -236,14 +228,18 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa ] ) + 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]}") + + if DEBUG_OUTPUT: + print(f"rotation angle:, {math.degrees(rotation_rad):.2f}") + print(f"k_rotated = ({k_rotated[0]:.2f}, {k_rotated[1]:.2f})") + electric_current, magnetic_current = equivalent_currents( ex_dft, ey_dft, hx_dft, hy_dft ) - if DEBUG_OUTPUT: - print(f"electric_current:, {electric_current}") - print(f"magnetic_current:, {magnetic_current}") - electric_current_rotated = rotation_matrix @ np.transpose( np.array([electric_current[0], electric_current[1]]) ) @@ -253,15 +249,11 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa ) 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}") - # Verify that ky of the rotated wavevector is 0. - k_rotated = rotation_matrix @ np.transpose(np.array([kx, ky])) - print(f"k_rotated = ({k_rotated[0]:.2f}, {k_rotated[1]:.2f})") - if k_rotated[1] == 0: - print("rotated ky is zero.") - current_amplitudes = [ electric_current_rotated[0], electric_current_rotated[1], @@ -272,12 +264,15 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa farfields = {} for component in farfield_components: farfields[component] = np.zeros( - (NUM_POLAR, NUM_AZIMUTHAL), dtype=np.complex128 + (NUM_POLAR, NUM_AZIMUTHAL), + dtype=np.complex128 ) for i, polar in enumerate(polar_rad): for j, azimuthal in enumerate(azimuthal_rad): rx, ry, rz = spherical_to_cartesian(polar, azimuthal) + r_rotated = rotation_matrix @ np.transpose(np.array([rx, ry])) + phase = np.exp(1j * (np.dot(k_rotated, r_rotated) + kz * rz)) for current_component, current_amplitude in zip( current_components, current_amplitudes @@ -285,20 +280,22 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa if abs(current_amplitude) == 0: continue - farfield_pol = far_fields( - kx, + farfield_amplitudes_pol = farfield_amplitudes( + k_rotated[0], kz, - rx, + r_rotated[0], rz, current_amplitude, current_component, ) + farfield_pol = np.array(farfield_amplitudes_pol) * phase + if ( current_component.name == "Jx" or current_component.name == "Ky" ): - # P polarization + # P polarization farfields["Ex"][i, j] += farfield_pol[0] farfields["Hy"][i, j] += farfield_pol[1] farfields["Ez"][i, j] += farfield_pol[2] @@ -306,33 +303,39 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jy" or current_component.name == "Kx" ): - # S polarization + # 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.linalg.inv(rotation_matrix) + # Brillouin zone integration. for i in range(NUM_POLAR): for j in range(NUM_AZIMUTHAL): - total_farfields["Ex"][i, j] = ( + total_farfields["Ex"][i, j] += ( inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] + - inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] + + farfields["Ez"][i, j] ) - total_farfields["Ey"][i, j] = ( + total_farfields["Ey"][i, j] += ( inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] + - inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] + inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] + + farfields["Ez"][i, j] ) - total_farfields["Ez"][i, j] = farfields["Ez"][i, j] - total_farfields["Hx"][i, j] = ( + total_farfields["Ez"][i, j] += farfields["Ez"][i, j] + + total_farfields["Hx"][i, j] += ( inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] + - inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] + + farfields["Hz"][i, j] ) - total_farfields["Hy"][i, j] = ( + total_farfields["Hy"][i, j] += ( inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] + - inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] + inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] + + farfields["Hz"][i, j] ) - total_farfields["Hz"][i, j] = farfields["Hz"][i, j] + total_farfields["Hz"][i, j] += farfields["Hz"][i, j] if mp.am_master(): np.savez( From 4fb1ca906c43d4096c530dc6d8eb0e8dc67c4b61 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Fri, 18 Oct 2024 12:07:58 -0700 Subject: [PATCH 4/7] perform BZ integration over wavegrid grid just inside light cone to ensure fields decay away --- python/examples/dipole_in_vacuum_1D.py | 129 ++++++++++++------------- 1 file changed, 62 insertions(+), 67 deletions(-) diff --git a/python/examples/dipole_in_vacuum_1D.py b/python/examples/dipole_in_vacuum_1D.py index 7aa14079d..3d8ffa197 100644 --- a/python/examples/dipole_in_vacuum_1D.py +++ b/python/examples/dipole_in_vacuum_1D.py @@ -11,11 +11,11 @@ Current = Enum("Current", "Jx Jy Kx Ky") -RESOLUTION_UM = 35 +RESOLUTION_UM = 25 WAVELENGTH_UM = 1.0 FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM -NUM_K = 17 -NUM_POLAR = 100 +NUM_K = 20 +NUM_POLAR = 80 NUM_AZIMUTHAL = 50 DEBUG_OUTPUT = True @@ -25,15 +25,15 @@ 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. + """ + 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. """ pml_um = 1.0 air_um = 10.0 @@ -42,9 +42,9 @@ def dipole_in_vacuum( pml_layers = [mp.PML(thickness=pml_um)] if dipole_component.name == "Jx": - src_cmpt = mp.Ex + src_cmpt = mp.Ex elif dipole_component.name == "Jy": - src_cmpt = mp.Ey + src_cmpt = mp.Ey elif dipole_component.name == "Kx": src_cmpt = mp.Hx elif dipole_component.name == "Ky": @@ -76,13 +76,10 @@ def dipole_in_vacuum( size=mp.Vector3(), ) - # TODO (oskooi): find out why meep.stop_when_fields_decayed fails to - # terminate simulation for certain values of kx and ky. sim.run( - until_after_sources=50 - # until_after_sources=mp.stop_when_fields_decayed( - # 10, src_cmpt, mp.Vector3(0, 0, 0.5423), 1e-6 - # ) + until_after_sources=mp.stop_when_fields_decayed( + 10, src_cmpt, mp.Vector3(0, 0, 0.5423), 1e-6 + ) ) ex_dft = sim.get_dft_array(dft_fields, mp.Ex, 0) @@ -96,13 +93,13 @@ def dipole_in_vacuum( 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. + """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) @@ -119,30 +116,30 @@ def farfield_amplitudes( current_amplitude: complex, current_component: Current, ) -> Tuple[complex, complex, complex]: - """Computes the S- or P-polarized far-field amplitudes from a sheet current. - - Based on the assumption that ky=0 such that the in-plane wavevector is xz. - - Args: - kx, kz: wavevector of the outgoing planewave in the x,z direction. - rx, rz: x,z coordinate of a point on the surface of a hemicircle. - current_amplitude: amplitude of the sheet current. - current_component: component of the sheet current. - - Returns: - A 3-tuple of the per-polarization electric and magnetic far fields. + """Computes the S- or P-polarized far-field amplitudes from a sheet current. + + Based on the assumption that ky=0. In-plane wavevector is xz. + + Args: + kx, kz: wavevector of the outgoing planewave in the x,z direction. + rx, rz: x,z coordinate of a point on the surface of a hemicircle. + current_amplitude: amplitude of the sheet current. + current_component: component of the sheet current. + + Returns: + A 3-tuple of the per-polarization electric and magnetic far fields. """ if current_component.name == "Jx": - # Jx --> (Ex, Hy) [P pol.] + # Jx --> (Ex, Hy) [P pol.] ex0 = frequency * current_amplitude / (2 * kz) elif current_component.name == "Jy": - # Jy --> (Hx, Ey) [S pol.] + # Jy --> (Hx, Ey) [S pol.] ey0 = -frequency * current_amplitude / (2 * kz) elif current_component.name == "Kx": - # Kx --> (Hx, Ey) [S pol.] + # Kx --> (Hx, Ey) [S pol.] ey0 = -current_amplitude / 2 elif current_component.name == "Ky": - # Ky --> (Ex, Hy) [P pol.] + # Ky --> (Ex, Hy) [P pol.] ex0 = current_amplitude / 2 if current_component.name == "Jx" or current_component.name == "Ky": @@ -156,14 +153,14 @@ def farfield_amplitudes( 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. + """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 = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.cos(azimuthal_rad) y = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.sin(azimuthal_rad) @@ -178,7 +175,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa kxs = np.linspace(-frequency, frequency, NUM_K) kys = np.linspace(-frequency, frequency, NUM_K) - # Far fields are defined on the surface of a hemisphere. + # Far fields are defined on the surface of a hemisphere. polar_rad = np.linspace(0, 0.5 * np.pi, NUM_POLAR) azimuthal_rad = np.linspace(0, 2 * np.pi, NUM_AZIMUTHAL) @@ -192,11 +189,12 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa dtype=np.complex128 ) + # 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) >= frequency: + # Skip wavevectors which are outside the light cone. + if np.sqrt(kx**2 + ky**2) >= (0.95 * frequency): continue kz = (frequency**2 - kx**2 - ky**2) ** 0.5 @@ -212,14 +210,8 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa 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. - if kx: - rotation_rad = -math.atan(ky / kx) - else: - if ky: - rotation_rad = -0.5 * np.pi - else: - rotation_rad = 0 + # Rotation angle around z axis to force ky=0. 0 radians is +x. + rotation_rad = -math.atan2(ky, kx) rotation_matrix = np.array( [ @@ -261,6 +253,9 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa magnetic_current_rotated[1], ] + # Obtain the far fields over a hemisphere given a + # sheet current with in-plane linear polarization. + farfields = {} for component in farfield_components: farfields[component] = np.zeros( @@ -272,7 +267,8 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa for j, azimuthal in enumerate(azimuthal_rad): rx, ry, rz = spherical_to_cartesian(polar, azimuthal) r_rotated = rotation_matrix @ np.transpose(np.array([rx, ry])) - phase = np.exp(1j * (np.dot(k_rotated, r_rotated) + kz * rz)) + # Note: r_rotated[1] (ry) is not used because k_rotated[1]=0 (ky). + phase = np.exp(1j * 2 * np.pi * (k_rotated[0] * rx + kz * rz)) for current_component, current_amplitude in zip( current_components, current_amplitudes @@ -295,7 +291,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jx" or current_component.name == "Ky" ): - # P polarization + # P polarization farfields["Ex"][i, j] += farfield_pol[0] farfields["Hy"][i, j] += farfield_pol[1] farfields["Ez"][i, j] += farfield_pol[2] @@ -303,14 +299,13 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jy" or current_component.name == "Kx" ): - # S polarization + # 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.linalg.inv(rotation_matrix) + inverse_rotation_matrix = np.transpose(rotation_matrix) - # Brillouin zone integration. for i in range(NUM_POLAR): for j in range(NUM_AZIMUTHAL): total_farfields["Ex"][i, j] += ( From 8ce770b37c3663c4b7f79ae6709d3f906fb5c38f Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sat, 19 Oct 2024 16:54:41 -0700 Subject: [PATCH 5/7] correctly rotate the far fields back into the original coordinate frame --- python/examples/dipole_in_vacuum_1D.py | 150 +++++++++++++------------ 1 file changed, 78 insertions(+), 72 deletions(-) diff --git a/python/examples/dipole_in_vacuum_1D.py b/python/examples/dipole_in_vacuum_1D.py index 3d8ffa197..7f6967397 100644 --- a/python/examples/dipole_in_vacuum_1D.py +++ b/python/examples/dipole_in_vacuum_1D.py @@ -14,9 +14,9 @@ RESOLUTION_UM = 25 WAVELENGTH_UM = 1.0 FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM -NUM_K = 20 +NUM_K = 25 NUM_POLAR = 80 -NUM_AZIMUTHAL = 50 +NUM_AZIMUTHAL = 80 DEBUG_OUTPUT = True frequency = 1 / WAVELENGTH_UM @@ -25,15 +25,15 @@ 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. + """ + 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. """ pml_um = 1.0 air_um = 10.0 @@ -54,7 +54,7 @@ def dipole_in_vacuum( mp.Source( src=mp.GaussianSource(frequency, fwidth=0.1 * frequency), component=src_cmpt, - center=mp.Vector3(), + center=mp.Vector3(0, 0, -0.5 * air_um), ) ] @@ -72,7 +72,7 @@ def dipole_in_vacuum( frequency, 0, 1, - center=mp.Vector3(0, 0, 0.5 * size_z_um - pml_um), + center=mp.Vector3(0, 0, 0.5 * air_um), size=mp.Vector3(), ) @@ -93,13 +93,13 @@ def dipole_in_vacuum( 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. + """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) @@ -116,30 +116,30 @@ def farfield_amplitudes( current_amplitude: complex, current_component: Current, ) -> Tuple[complex, complex, complex]: - """Computes the S- or P-polarized far-field amplitudes from a sheet current. - - Based on the assumption that ky=0. In-plane wavevector is xz. - - Args: - kx, kz: wavevector of the outgoing planewave in the x,z direction. - rx, rz: x,z coordinate of a point on the surface of a hemicircle. - current_amplitude: amplitude of the sheet current. - current_component: component of the sheet current. - - Returns: - A 3-tuple of the per-polarization electric and magnetic far fields. + """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π. + rx, rz: x,z coordinate of a point on the surface of a hemicircle. + current_amplitude: amplitude of the sheet current. + current_component: component of the sheet current. + + Returns: + A 3-tuple of the per-polarization electric and magnetic far fields. """ if current_component.name == "Jx": - # Jx --> (Ex, Hy) [P pol.] + # Jx --> (Ex, Hy, Ez) [P pol.] ex0 = frequency * current_amplitude / (2 * kz) elif current_component.name == "Jy": - # Jy --> (Hx, Ey) [S pol.] + # Jy --> (Hx, Ey, Hz) [S pol.] ey0 = -frequency * current_amplitude / (2 * kz) elif current_component.name == "Kx": - # Kx --> (Hx, Ey) [S pol.] + # Kx --> (Hx, Ey, Hz) [S pol.] ey0 = -current_amplitude / 2 elif current_component.name == "Ky": - # Ky --> (Ex, Hy) [P pol.] + # Ky --> (Ex, Hy, Ez) [P pol.] ex0 = current_amplitude / 2 if current_component.name == "Jx" or current_component.name == "Ky": @@ -153,14 +153,14 @@ def farfield_amplitudes( 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. + """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 = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.cos(azimuthal_rad) y = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.sin(azimuthal_rad) @@ -175,9 +175,10 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa kxs = np.linspace(-frequency, frequency, NUM_K) kys = np.linspace(-frequency, frequency, NUM_K) - # Far fields are defined on the surface of a hemisphere. + # Far fields are defined on the surface of a hemisphere. polar_rad = np.linspace(0, 0.5 * np.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"] @@ -189,11 +190,11 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa dtype=np.complex128 ) - # Brillouin zone integration over 2D grid of (kx, ky). + # 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. + # Skip wavevectors which are outside the light cone. if np.sqrt(kx**2 + ky**2) >= (0.95 * frequency): continue @@ -210,9 +211,12 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa 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 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)], @@ -224,10 +228,6 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa if k_rotated[1] > 1e-10: raise ValueError(f"rotated ky is nonzero: {k_rotated[1]}") - if DEBUG_OUTPUT: - print(f"rotation angle:, {math.degrees(rotation_rad):.2f}") - print(f"k_rotated = ({k_rotated[0]:.2f}, {k_rotated[1]:.2f})") - electric_current, magnetic_current = equivalent_currents( ex_dft, ey_dft, hx_dft, hy_dft ) @@ -253,8 +253,11 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa magnetic_current_rotated[1], ] - # Obtain the far fields over a hemisphere given a - # sheet current with in-plane linear polarization. + # Obtain the far fields over a hemisphere given a sheet current + # with linear in-plane polarization. + + azimuthal_index_shift = int(abs(rotation_rad) // delta_azimuthal_rad) + inverse_rotation_matrix = np.transpose(rotation_matrix) farfields = {} for component in farfield_components: @@ -267,9 +270,11 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa for j, azimuthal in enumerate(azimuthal_rad): rx, ry, rz = spherical_to_cartesian(polar, azimuthal) r_rotated = rotation_matrix @ np.transpose(np.array([rx, ry])) - # Note: r_rotated[1] (ry) is not used because k_rotated[1]=0 (ky). - phase = np.exp(1j * 2 * np.pi * (k_rotated[0] * rx + kz * rz)) + # Note: r_rotated[1] (ry) is not used because k_rotated[1]=0 (ky). + phase = np.exp(1j * 2 * np.pi * (k_rotated[0] * r_rotated[0] + kz * rz)) + # 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 ): @@ -291,7 +296,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jx" or current_component.name == "Ky" ): - # P polarization + # P polarization farfields["Ex"][i, j] += farfield_pol[0] farfields["Hy"][i, j] += farfield_pol[1] farfields["Ez"][i, j] += farfield_pol[2] @@ -299,38 +304,38 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jy" or current_component.name == "Kx" ): - # S polarization + # 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): - total_farfields["Ex"][i, j] += ( + # Map the rotated point to the closest azimuthal angle. + azimuthal_index = (j + azimuthal_index_shift) % NUM_AZIMUTHAL + total_farfields["Ex"][i, azimuthal_index] += ( inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] + - inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] + - farfields["Ez"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] ) - total_farfields["Ey"][i, j] += ( + total_farfields["Ey"][i, azimuthal_index] += ( inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] + - inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] + + inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] + ) + total_farfields["Ez"][i, azimuthal_index] += ( farfields["Ez"][i, j] ) - total_farfields["Ez"][i, j] += farfields["Ez"][i, j] - total_farfields["Hx"][i, j] += ( + total_farfields["Hx"][i, azimuthal_index] += ( inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] + - inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] + - farfields["Hz"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] ) - total_farfields["Hy"][i, j] += ( + total_farfields["Hy"][i, azimuthal_index] += ( inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] + - inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] + + inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] + ) + total_farfields["Hz"][i, azimuthal_index] += ( farfields["Hz"][i, j] ) - total_farfields["Hz"][i, j] += farfields["Hz"][i, j] if mp.am_master(): np.savez( @@ -346,3 +351,4 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa kxs=kxs, kys=kys, ) + From af17883c4018d7ee74270cc1b2aaf98672e4bc07 Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Thu, 14 Nov 2024 23:34:24 -0800 Subject: [PATCH 6/7] use corrected formulas for far fields of sheet current and remove FARFIELD_RADIUS_UM --- python/examples/dipole_in_vacuum_1D.py | 158 ++++++++++++------------- 1 file changed, 79 insertions(+), 79 deletions(-) diff --git a/python/examples/dipole_in_vacuum_1D.py b/python/examples/dipole_in_vacuum_1D.py index 7f6967397..be6520302 100644 --- a/python/examples/dipole_in_vacuum_1D.py +++ b/python/examples/dipole_in_vacuum_1D.py @@ -13,11 +13,10 @@ RESOLUTION_UM = 25 WAVELENGTH_UM = 1.0 -FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM -NUM_K = 25 -NUM_POLAR = 80 -NUM_AZIMUTHAL = 80 -DEBUG_OUTPUT = True +NUM_K = 100 +NUM_POLAR = 50 +NUM_AZIMUTHAL = 200 +DEBUG_OUTPUT = False frequency = 1 / WAVELENGTH_UM @@ -25,18 +24,18 @@ 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. + """ + 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. """ pml_um = 1.0 - air_um = 10.0 + air_um = 50.0 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)] @@ -65,6 +64,7 @@ def dipole_in_vacuum( sources=sources, boundary_layers=pml_layers, k_point=mp.Vector3(kx, ky, 0), + dimensions=3 ) dft_fields = sim.add_dft_fields( @@ -73,7 +73,7 @@ def dipole_in_vacuum( 0, 1, center=mp.Vector3(0, 0, 0.5 * air_um), - size=mp.Vector3(), + size=mp.Vector3() ) sim.run( @@ -93,13 +93,13 @@ def dipole_in_vacuum( 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. + """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) @@ -116,55 +116,57 @@ def farfield_amplitudes( 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π. - rx, rz: x,z coordinate of a point on the surface of a hemicircle. - current_amplitude: amplitude of the sheet current. - current_component: component of the sheet current. - - Returns: - A 3-tuple of the per-polarization electric and magnetic far fields. + """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π. + rx, rz: x,z coordinate of a point on the surface of a hemicircle. + 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 = frequency * current_amplitude / (2 * kz) + # 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.] + # 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.] - ey0 = -current_amplitude / 2 + # 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.] - ex0 = current_amplitude / 2 + # 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": - hy0 = frequency * ex0 / kz ez0 = -kx * hy0 / frequency return (ex0, hy0, ez0) elif current_component.name == "Jy" or current_component.name == "Kx": - hx0 = -kz * ey0 / frequency 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. + """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 = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.cos(azimuthal_rad) - y = FARFIELD_RADIUS_UM * np.sin(polar_rad) * np.sin(azimuthal_rad) - z = FARFIELD_RADIUS_UM * np.cos(polar_rad) + 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) @@ -175,7 +177,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa kxs = np.linspace(-frequency, frequency, NUM_K) kys = np.linspace(-frequency, frequency, NUM_K) - # Far fields are defined on the surface of a hemisphere. + # Far fields are defined on the surface of a hemisphere. polar_rad = np.linspace(0, 0.5 * np.pi, NUM_POLAR) azimuthal_rad = np.linspace(0, 2 * np.pi, NUM_AZIMUTHAL) delta_azimuthal_rad = 2 * np.pi / (NUM_AZIMUTHAL - 1) @@ -190,11 +192,11 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa dtype=np.complex128 ) - # Brillouin zone integration over 2D grid of (kx, ky). + # 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. + # Skip wavevectors which are outside the light cone. if np.sqrt(kx**2 + ky**2) >= (0.95 * frequency): continue @@ -211,11 +213,11 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa 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 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}") + print(f"rotation angle:, {math.degrees(rotation_rad):.2f}°") rotation_matrix = np.array( [ @@ -253,11 +255,8 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa magnetic_current_rotated[1], ] - # Obtain the far fields over a hemisphere given a sheet current - # with linear in-plane polarization. - - azimuthal_index_shift = int(abs(rotation_rad) // delta_azimuthal_rad) - inverse_rotation_matrix = np.transpose(rotation_matrix) + # Obtain the far fields over a hemisphere given a sheet current + # with linear in-plane polarization. farfields = {} for component in farfield_components: @@ -270,11 +269,11 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa for j, azimuthal in enumerate(azimuthal_rad): rx, ry, rz = spherical_to_cartesian(polar, azimuthal) r_rotated = rotation_matrix @ np.transpose(np.array([rx, ry])) - # Note: r_rotated[1] (ry) is not used because k_rotated[1]=0 (ky). + # Note: r_rotated[1] (ry) is not used because k_rotated[1]=0 (ky). phase = np.exp(1j * 2 * np.pi * (k_rotated[0] * r_rotated[0] + kz * rz)) - # Obtain the farfields for each type of sheet current - # and keep a running sum. + # 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 ): @@ -296,7 +295,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jx" or current_component.name == "Ky" ): - # P polarization + # P polarization farfields["Ex"][i, j] += farfield_pol[0] farfields["Hy"][i, j] += farfield_pol[1] farfields["Ez"][i, j] += farfield_pol[2] @@ -304,36 +303,37 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jy" or current_component.name == "Kx" ): - # S polarization + # 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): - # Map the rotated point to the closest azimuthal angle. - azimuthal_index = (j + azimuthal_index_shift) % NUM_AZIMUTHAL - total_farfields["Ex"][i, azimuthal_index] += ( + total_farfields["Ex"][i, j] += ( inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] ) - total_farfields["Ey"][i, azimuthal_index] += ( + total_farfields["Ey"][i, j] += ( inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] + inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] ) - total_farfields["Ez"][i, azimuthal_index] += ( + total_farfields["Ez"][i, j] += ( farfields["Ez"][i, j] ) - total_farfields["Hx"][i, azimuthal_index] += ( + total_farfields["Hx"][i, j] += ( inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] ) - total_farfields["Hy"][i, azimuthal_index] += ( + total_farfields["Hy"][i, j] += ( inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] + inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] ) - total_farfields["Hz"][i, azimuthal_index] += ( + total_farfields["Hz"][i, j] += ( farfields["Hz"][i, j] ) @@ -350,5 +350,5 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa azimuthal_rad=azimuthal_rad, kxs=kxs, kys=kys, + dipole_component=dipole_component ) - From eef5dd16b8f73274bbcdc1e6e3a44e7ec4e35acb Mon Sep 17 00:00:00 2001 From: Ardavan Oskooi Date: Sun, 24 Nov 2024 21:07:28 -0800 Subject: [PATCH 7/7] compute the total radial Poynting flux rather than E, H fields at each grid point in the far field via Brillouin-zone integration --- python/examples/dipole_in_vacuum_1D.py | 191 +++++++++++++------------ 1 file changed, 97 insertions(+), 94 deletions(-) diff --git a/python/examples/dipole_in_vacuum_1D.py b/python/examples/dipole_in_vacuum_1D.py index be6520302..a214dde40 100644 --- a/python/examples/dipole_in_vacuum_1D.py +++ b/python/examples/dipole_in_vacuum_1D.py @@ -13,9 +13,11 @@ RESOLUTION_UM = 25 WAVELENGTH_UM = 1.0 -NUM_K = 100 +NUM_K = 40 NUM_POLAR = 50 -NUM_AZIMUTHAL = 200 +NUM_AZIMUTHAL = 50 +PML_UM = 1.0 +AIR_UM = 10.0 DEBUG_OUTPUT = False frequency = 1 / WAVELENGTH_UM @@ -24,21 +26,19 @@ 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. """ - 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. - """ - pml_um = 1.0 - air_um = 50.0 - size_z_um = pml_um + air_um + pml_um + 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)] + pml_layers = [mp.PML(thickness=PML_UM)] if dipole_component.name == "Jx": src_cmpt = mp.Ex @@ -51,9 +51,9 @@ def dipole_in_vacuum( sources = [ mp.Source( - src=mp.GaussianSource(frequency, fwidth=0.1 * frequency), + src=mp.GaussianSource(frequency, fwidth=0.05 * frequency), component=src_cmpt, - center=mp.Vector3(0, 0, -0.5 * air_um), + center=mp.Vector3(0, 0, -0.5 * AIR_UM), ) ] @@ -64,7 +64,6 @@ def dipole_in_vacuum( sources=sources, boundary_layers=pml_layers, k_point=mp.Vector3(kx, ky, 0), - dimensions=3 ) dft_fields = sim.add_dft_fields( @@ -72,13 +71,13 @@ def dipole_in_vacuum( frequency, 0, 1, - center=mp.Vector3(0, 0, 0.5 * air_um), + 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.5423), 1e-6 + 10, src_cmpt, mp.Vector3(0, 0, 0.5 * AIR_UM), 1e-6 ) ) @@ -93,13 +92,13 @@ def dipole_in_vacuum( 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. + """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) @@ -111,38 +110,35 @@ def equivalent_currents( def farfield_amplitudes( kx: float, kz: float, - rx: float, - rz: 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π. - rx, rz: x,z coordinate of a point on the surface of a hemicircle. - 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. + """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.] + # 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.] + # 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.] + # 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.] + # Ky --> (Ex, Hy, Ez) [P pol.] hy0 = frequency * current_amplitude / (2 * kz) ex0 = kz * hy0 / frequency @@ -155,14 +151,14 @@ def farfield_amplitudes( 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. + """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) @@ -177,27 +173,26 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa 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 * np.pi, NUM_POLAR) + # 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"] - total_farfields = {} - for component in farfield_components: - total_farfields[component] = np.zeros( - (NUM_POLAR, NUM_AZIMUTHAL), - dtype=np.complex128 - ) + poynting_flux = np.zeros((NUM_POLAR, NUM_AZIMUTHAL)) - # Brillouin zone integration over 2D grid of (kx, ky). + # 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.95 * frequency): + # 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 @@ -213,7 +208,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa 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 angle around z axis to force ky=0. 0 radians is +x. rotation_rad = -math.atan2(ky, kx) if DEBUG_OUTPUT: @@ -255,8 +250,8 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa magnetic_current_rotated[1], ] - # Obtain the far fields over a hemisphere given a sheet current - # with linear in-plane polarization. + # Obtain the far fields over a hemisphere given a sheet current + # with linear in-plane polarization. farfields = {} for component in farfield_components: @@ -265,37 +260,29 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa dtype=np.complex128 ) - for i, polar in enumerate(polar_rad): - for j, azimuthal in enumerate(azimuthal_rad): - rx, ry, rz = spherical_to_cartesian(polar, azimuthal) - r_rotated = rotation_matrix @ np.transpose(np.array([rx, ry])) - # Note: r_rotated[1] (ry) is not used because k_rotated[1]=0 (ky). - phase = np.exp(1j * 2 * np.pi * (k_rotated[0] * r_rotated[0] + kz * rz)) + 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. + # 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_amplitudes_pol = farfield_amplitudes( + farfield_pol = farfield_amplitudes( k_rotated[0], kz, - r_rotated[0], - rz, current_amplitude, current_component, ) - farfield_pol = np.array(farfield_amplitudes_pol) * phase - if ( current_component.name == "Jx" or current_component.name == "Ky" ): - # P polarization + # P polarization farfields["Ex"][i, j] += farfield_pol[0] farfields["Hy"][i, j] += farfield_pol[1] farfields["Ez"][i, j] += farfield_pol[2] @@ -303,7 +290,7 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa current_component.name == "Jy" or current_component.name == "Kx" ): - # S polarization + # S polarization farfields["Hx"][i, j] += farfield_pol[0] farfields["Ey"][i, j] += farfield_pol[1] farfields["Hz"][i, j] += farfield_pol[2] @@ -313,42 +300,58 @@ def spherical_to_cartesian(polar_rad, azimuthal_rad) -> Tuple[float, float, floa for i in range(NUM_POLAR): for j in range(NUM_AZIMUTHAL): - total_farfields["Ex"][i, j] += ( + ex_rotated = ( inverse_rotation_matrix[0, 0] * farfields["Ex"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Ey"][i, j] ) - total_farfields["Ey"][i, j] += ( + ey_rotated = ( inverse_rotation_matrix[1, 0] * farfields["Ex"][i, j] + inverse_rotation_matrix[1, 1] * farfields["Ey"][i, j] ) - total_farfields["Ez"][i, j] += ( - farfields["Ez"][i, j] - ) - total_farfields["Hx"][i, j] += ( + hx_rotated = ( inverse_rotation_matrix[0, 0] * farfields["Hx"][i, j] + inverse_rotation_matrix[0, 1] * farfields["Hy"][i, j] ) - total_farfields["Hy"][i, j] += ( + hy_rotated = ( inverse_rotation_matrix[1, 0] * farfields["Hx"][i, j] + inverse_rotation_matrix[1, 1] * farfields["Hy"][i, j] ) - total_farfields["Hz"][i, j] += ( - farfields["Hz"][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", - **total_farfields, RESOLUTION_UM=RESOLUTION_UM, WAVELENGTH_UM=WAVELENGTH_UM, - FARFIELD_RADIUS_UM=FARFIELD_RADIUS_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 + dipole_component=dipole_component, + poynting_flux=poynting_flux ) +