diff --git a/doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md b/doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md index f4c7ed213..c34916868 100644 --- a/doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md +++ b/doc/docs/Python_Tutorials/Near_to_Far_Field_Spectra.md @@ -664,7 +664,7 @@ The example uses the same setup as the [previous tutorial](#radiation-pattern-of $$P(\theta) = \frac{1}{R} \int_0^R P(r,\theta) s(r) 2\pi rdr = \frac{1}{N} \sum_{n=0}^{N-1} P(r_n,\theta) s(r_n) 2\pi r_n \Delta r$$, -where $s(r)$ is a weighting function necessary for ensuring equal contribution from all dipoles relative to the dipole at $r = 0$. Note: an $E_r$ dipole at $r = 0$ must be placed at $r_0 = 1.5\Delta r$ due to an [interpolation bug](https://github.com/NanoComp/meep/issues/2704). $s(r)$ can be determined empirically by computing the emission in vacuum for a set of dipoles at different radial positions. The emission profile of a dipole in vacuum is a constant independent of its position. This criteria is used to obtain: $s(r) = \frac{1}{2(r/r_0)^2}$. This weighting function is also used to average the flux emitted by each dipole (obtained using using the LDOS feature). This quantity is the denominator in the expression for the extraction efficiency. +where $s(r)$ is a weighting function necessary for ensuring equal contribution from all dipoles relative to the dipole at $r = 0$. Note: an $E_r$ dipole at $r = 0$ must be placed at $r_0 = 1.5\Delta r$ due to an [interpolation bug](https://github.com/NanoComp/meep/issues/2704). $s(r)$ can be determined empirically by computing the emission in vacuum for a set of dipoles at different radial positions. The emission profile of a dipole in vacuum is a constant independent of its position. This criteria is used to obtain: $s(r) = (r_0/r)^2$. This weighting function is also used to average the flux emitted by each dipole (obtained using using the LDOS feature). This quantity is the denominator in the expression for the extraction efficiency. This figure shows the radiation pattern from $N=11$ dipoles with $\lambda$ of 1.0 $\mu$m in the middle of a disc of height 0.29 $\mu$m, radius 1.2 $\mu$m, and refractive index 2.4. @@ -677,36 +677,30 @@ The simulation script is in [examples/disc_extraction_efficiency.py](https://git import math from typing import Tuple -import matplotlib +import matplotlib.pyplot as plt import meep as mp import numpy as np -matplotlib.use("agg") -import matplotlib.pyplot as plt - - -N_DISC = 2.4 # refractive index of disc -DISC_RADIUS_UM = 1.2 # radius of disc -WAVELENGTH_UM = 1.0 # wavelength (in vacuum) -# radius of quarter circle for computing flux in far field -FARFIELD_RADIUS_UM = 1000 * WAVELENGTH_UM - -RESOLUTION = 50 # pixels/μm +RESOLUTION_UM = 50 +N_DISC = 2.4 +DISC_RADIUS_UM = 1.2 +WAVELENGTH_UM = 1.0 +NUM_FARFIELD_PTS = 200 +FARFIELD_ANGLES = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS) +FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM +NUM_DIPOLES = 11 -def plot_radiation_pattern_polar(theta_rad: np.ndarray, radial_flux: np.ndarray): +def plot_radiation_pattern_polar(radial_flux: np.ndarray): """Plots the radiation pattern in polar coordinates. Args: - theta_rad: angles of the radiation pattern. The angles increase clockwise - with zero at the pole (+z direction) and π/2 at the equator (+r - direction). radial_flux: radial flux of the far fields in polar coordinates. """ fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6)) ax.plot( - theta_rad, + FARFIELD_ANGLES, radial_flux, "b-", ) @@ -726,20 +720,19 @@ def plot_radiation_pattern_polar(theta_rad: np.ndarray, radial_flux: np.ndarray) ) -def plot_radiation_pattern_3d(theta_rad: np.ndarray, radial_flux: np.ndarray): +def plot_radiation_pattern_3d(radial_flux: np.ndarray): """Plots the radiation pattern in 3d Cartesian coordinates. Args: - theta_rad: angles of the radiation pattern. radial_flux: radial flux of the far fields in polar coordinates. """ - phis = np.linspace(0, 2 * np.pi, num_angles) + phis = np.linspace(0, 2 * np.pi, NUM_FARFIELD_PTS) - xs = np.zeros((len(theta_rad), len(phis))) - ys = np.zeros((len(theta_rad), len(phis))) - zs = np.zeros((len(theta_rad), len(phis))) + xs = np.zeros((NUM_FARFIELD_PTS, NUM_FARFIELD_PTS)) + ys = np.zeros((NUM_FARFIELD_PTS, NUM_FARFIELD_PTS)) + zs = np.zeros((NUM_FARFIELD_PTS, NUM_FARFIELD_PTS)) - for i, theta in enumerate(theta_rad): + for i, theta in enumerate(FARFIELD_ANGLES): for j, phi in enumerate(phis): xs[i, j] = radial_flux[i] * np.sin(theta) * np.cos(phi) ys[i, j] = radial_flux[i] * np.sin(theta) * np.sin(phi) @@ -760,54 +753,55 @@ def plot_radiation_pattern_3d(theta_rad: np.ndarray, radial_flux: np.ndarray): ) -def radiation_pattern( - theta_rad: np.ndarray, sim: mp.Simulation, n2f_mon: mp.DftNear2Far -) -> np.ndarray: +def radiation_pattern(sim: mp.Simulation, n2f_mon: mp.DftNear2Far) -> np.ndarray: """Computes the radiation pattern from the near fields. Args: - theta_rad: angles of the radiation pattern. sim: a `Simulation` object. n2f_mon: a `DftNear2Far` object returned by `Simulation.add_near2far`. + + Returns: + The radiation pattern (radial flux at each angle) as a 1d array. """ - e_field = np.zeros((theta_rad.shape[0], 3), dtype=np.complex128) - h_field = np.zeros((theta_rad.shape[0], 3), dtype=np.complex128) - for n in range(num_angles): + e_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128) + h_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128) + for n in range(NUM_FARFIELD_PTS): far_field = sim.get_farfield( n2f_mon, mp.Vector3( - FARFIELD_RADIUS_UM * math.sin(theta_rad[n]), + FARFIELD_RADIUS_UM * math.sin(FARFIELD_ANGLES[n]), 0, - FARFIELD_RADIUS_UM * math.cos(theta_rad[n]), + FARFIELD_RADIUS_UM * math.cos(FARFIELD_ANGLES[n]), ), ) - e_field[n, :] = [np.conj(far_field[j]) for j in range(3)] + e_field[n, :] = [far_field[j] for j in range(3)] h_field[n, :] = [far_field[j + 3] for j in range(3)] - flux_r = np.real(e_field[:, 1] * h_field[:, 2] - e_field[:, 2] * h_field[:, 1]) - flux_z = np.real(e_field[:, 0] * h_field[:, 1] - e_field[:, 1] * h_field[:, 0]) - flux_rz = np.sqrt(np.square(flux_r) + np.square(flux_z)) + flux_x = np.real(np.conj(e_field[:, 1]) * h_field[:, 2] - + np.conj(e_field[:, 2]) * h_field[:, 1]) + flux_z = np.real(np.conj(e_field[:, 0]) * h_field[:, 1] - + np.conj(e_field[:, 1]) * h_field[:, 0]) + flux_r = np.sqrt(np.square(flux_x) + np.square(flux_z)) - return flux_rz + return flux_r -def radiation_pattern_flux(theta_rad: np.ndarray, radial_flux: np.ndarray) -> float: +def radiation_pattern_flux(radial_flux: np.ndarray) -> float: """Computes the total flux from the radiation pattern. Based on integrating the radiation pattern over solid angles spanned by polar angles in the range of [0, π/2]. Args: - theta_rad: angles of the radiation pattern. radial_flux: radial flux of the far fields in polar coordinates. """ dphi = 2 * math.pi - dtheta_rad = theta_rad[1] - theta_rad[0] + dtheta = FARFIELD_ANGLES[1] - FARFIELD_ANGLES[0] total_flux = ( - np.sum(radial_flux * np.sin(theta_rad)) + np.sum(radial_flux * np.sin(FARFIELD_ANGLES)) * FARFIELD_RADIUS_UM**2 - * dtheta_rad + * dtheta * dphi ) @@ -815,40 +809,39 @@ def radiation_pattern_flux(theta_rad: np.ndarray, radial_flux: np.ndarray) -> fl def dipole_in_disc( - t_disc_um: float, h_disc: float, rpos: float, m: int, theta_rad: np.ndarray + disc_um: float, zpos: float, rpos_um: float, m: int ) -> Tuple[float, np.ndarray]: """Computes the total flux and radiation pattern of a dipole in a disc. Args: - t_disc_um: thickness of disc. - h_disc: height of dipole above ground plane as fraction of t_disc_um. - rpos: radial position of dipole. + disc_um: thickness of disc. + zpos: height of dipole above ground plane as fraction of disc_um. + rpos_um: radial position of dipole. m: angular φ dependence of the fields exp(imφ). - theta_rad: angles of the radiation pattern. Returns: A 2-tuple of the total flux and the radiation pattern. """ - t_pml_um = 0.5 # thickness of PML - t_air_um = 1.0 # thickness of air padding above disc - length_r_um = 5.0 # length of cell in r + pml_um = 1.0 # thickness of PML + padding_um = 1.0 # thickness of air padding above disc + r_um = 5.0 # length of cell in r frequency = 1 / WAVELENGTH_UM # center frequency of source/monitor # field decay threshold for runtime termination criteria - decay_tol = 1e-6 + field_decay_threshold = 1e-6 - size_r = length_r_um + t_pml_um - size_z = t_disc_um + t_air_um + t_pml_um + size_r = r_um + pml_um + size_z = disc_um + padding_um + pml_um cell_size = mp.Vector3(size_r, 0, size_z) boundary_layers = [ - mp.PML(t_pml_um, direction=mp.R), - mp.PML(t_pml_um, direction=mp.Z, side=mp.High), + mp.PML(pml_um, direction=mp.R), + mp.PML(pml_um, direction=mp.Z, side=mp.High), ] src_cmpt = mp.Er - src_pt = mp.Vector3(rpos, 0, -0.5 * size_z + h_disc * t_disc_um) + src_pt = mp.Vector3(rpos_um, 0, -0.5 * size_z + zpos * disc_um) sources = [ mp.Source( src=mp.GaussianSource(frequency, fwidth=0.1 * frequency), @@ -860,13 +853,15 @@ def dipole_in_disc( geometry = [ mp.Block( material=mp.Medium(index=N_DISC), - center=mp.Vector3(0.5 * DISC_RADIUS_UM, 0, -0.5 * size_z + 0.5 * t_disc_um), - size=mp.Vector3(DISC_RADIUS_UM, mp.inf, t_disc_um), + center=mp.Vector3( + 0.5 * DISC_RADIUS_UM, 0, -0.5 * size_z + 0.5 * disc_um + ), + size=mp.Vector3(DISC_RADIUS_UM, mp.inf, disc_um), ) ] sim = mp.Simulation( - resolution=RESOLUTION, + resolution=RESOLUTION_UM, cell_size=cell_size, dimensions=mp.CYLINDRICAL, m=m, @@ -880,14 +875,16 @@ def dipole_in_disc( 0, 1, mp.FluxRegion( - center=mp.Vector3(0.5 * length_r_um, 0, 0.5 * size_z - t_pml_um), - size=mp.Vector3(length_r_um, 0, 0), + center=mp.Vector3(0.5 * r_um, 0, 0.5 * size_z - pml_um), + size=mp.Vector3(r_um, 0, 0), ), mp.FluxRegion( center=mp.Vector3( - length_r_um, 0, 0.5 * size_z - t_pml_um - 0.5 * (t_air_um + t_disc_um) + r_um, + 0, + 0.5 * size_z - pml_um - 0.5 * (padding_um + disc_um) ), - size=mp.Vector3(0, 0, t_air_um + t_disc_um), + size=mp.Vector3(0, 0, padding_um + disc_um), ), ) @@ -897,14 +894,15 @@ def dipole_in_disc( 50, src_cmpt, src_pt, - decay_tol, + field_decay_threshold, ), ) - delta_vol = 2 * np.pi * rpos / (RESOLUTION**2) - dipole_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * delta_vol + delta_vol = 2 * np.pi * rpos_um / (RESOLUTION_UM**2) + dipole_flux = (-np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * + delta_vol) - dipole_radiation_pattern = radiation_pattern(theta_rad, sim, n2f_mon) + dipole_radiation_pattern = radiation_pattern(sim, n2f_mon) return dipole_flux, dipole_radiation_pattern @@ -912,92 +910,79 @@ def dipole_in_disc( if __name__ == "__main__": disc_thickness = 0.7 * WAVELENGTH_UM / N_DISC dipole_height = 0.5 - num_dipoles = 11 - dipole_rpos = np.linspace(0, DISC_RADIUS_UM, num_dipoles) + dipole_rpos_um = np.linspace(0, DISC_RADIUS_UM, NUM_DIPOLES) - # An Er source at r=0 needs to be slighty offset. + # An Er source at r = 0 needs to be slighty offset. # https://github.com/NanoComp/meep/issues/2704 - dipole_rpos[0] = 1.5 / RESOLUTION - - delta_rpos = dipole_rpos[2] - dipole_rpos[1] + dipole_rpos_um[0] = 1.5 / RESOLUTION_UM - # number of angular grid points in [0, π/2] - num_angles = 100 + delta_rpos_um = dipole_rpos_um[2] - dipole_rpos_um[1] - # grid of polar angles for computing radiated flux in far field - theta_rad = np.linspace(0, 0.5 * math.pi, num_angles) - - # r = 0 source requires a single simulation with m = ±1. + # Er source at r = 0 requires a single simulation with m = ±1. m = -1 dipole_flux, dipole_radiation_pattern = dipole_in_disc( - disc_thickness, dipole_height, dipole_rpos[0], m, theta_rad + disc_thickness, dipole_height, dipole_rpos_um[0], m, ) - flux_total = dipole_flux * dipole_rpos[0] * delta_rpos - radiation_pattern_total = dipole_radiation_pattern * dipole_rpos[0] * delta_rpos + flux_total = dipole_flux * dipole_rpos_um[0] * delta_rpos_um + radiation_pattern_total = dipole_radiation_pattern * dipole_rpos_um[0] * delta_rpos_um print( - f"dipole:, {dipole_rpos[0]:.4f}, " - f"{radiation_pattern_flux(theta_rad, dipole_radiation_pattern):.6f}" + f"dipole:, {dipole_rpos_um[0]:.4f}, " + f"{radiation_pattern_flux(dipole_radiation_pattern):.6f}" ) - # r > 0 source requires Fourier-series expansion of φ. - flux_tol = 1e-3 # threshold flux to determine when to truncate expansion - for rpos in dipole_rpos[1:]: + # Er source at r > 0 requires Fourier-series expansion of φ. + + # Threshold flux to determine when to truncate expansion. + flux_decay_threshold = 1e-2 + + for rpos_um in dipole_rpos_um[1:]: dipole_flux_total = 0 - dipole_radiation_pattern_total = np.zeros((num_angles,)) + dipole_radiation_pattern_total = np.zeros((NUM_FARFIELD_PTS,)) dipole_radiation_pattern_flux_max = 0 m = 0 while True: dipole_flux, dipole_radiation_pattern = dipole_in_disc( - disc_thickness, dipole_height, rpos, m, theta_rad - ) - dipole_flux_total += dipole_flux if m == 0 else 2 * dipole_flux - dipole_radiation_pattern_total += ( - dipole_radiation_pattern if m == 0 else 2 * dipole_radiation_pattern + disc_thickness, dipole_height, rpos_um, m ) + dipole_flux_total += dipole_flux * (2 if m == 0 else 1) + dipole_radiation_pattern_total += (dipole_radiation_pattern * + (2 if m == 0 else 1)) dipole_radiation_pattern_flux = radiation_pattern_flux( - theta_rad, dipole_radiation_pattern + dipole_radiation_pattern ) - if dipole_radiation_pattern_flux > dipole_radiation_pattern_flux_max: + if (dipole_radiation_pattern_flux > + dipole_radiation_pattern_flux_max): dipole_radiation_pattern_flux_max = dipole_radiation_pattern_flux - if m > 0 and ( - (dipole_radiation_pattern_flux / dipole_radiation_pattern_flux_max) - < flux_tol - ): + if (m > 0 and + (dipole_radiation_pattern_flux / + dipole_radiation_pattern_flux_max) < flux_decay_threshold): break print( - f"dipole-m:, {rpos:.4f}, {m}, " f"{dipole_radiation_pattern_flux:.6f}" + f"dipole:, {rpos_um:.4f}, {m}, {dipole_radiation_pattern_flux:.6f}" ) m += 1 - scale_factor = 1 / (2 * (rpos / dipole_rpos[0]) ** 2) - flux_total += dipole_flux_total * scale_factor * rpos * delta_rpos + scale_factor = (dipole_rpos_um[0] / rpos_um)**2 + flux_total += dipole_flux_total * scale_factor * rpos_um * delta_rpos_um radiation_pattern_total += ( - dipole_radiation_pattern_total * scale_factor * rpos * delta_rpos + dipole_radiation_pattern_total * scale_factor * rpos_um * delta_rpos_um ) - dipole_radiation_pattern_total_flux = radiation_pattern_flux( - theta_rad, dipole_radiation_pattern_total - ) - print( - f"dipole:, {rpos:.4f}, {m}, " f"{dipole_radiation_pattern_total_flux:.6f}" - ) - flux_total /= num_dipoles - radiation_pattern_total /= num_dipoles - - radiation_pattern_total_flux = radiation_pattern_flux( - theta_rad, radiation_pattern_total - ) + flux_total /= NUM_DIPOLES + radiation_pattern_total /= NUM_DIPOLES + radiation_pattern_total_flux = radiation_pattern_flux(radiation_pattern_total) extraction_efficiency = radiation_pattern_total_flux / flux_total print(f"exteff:, {extraction_efficiency:.6f}") - plot_radiation_pattern_polar(radiation_pattern_total * FARFIELD_RADIUS_UM**2) - plot_radiation_pattern_3d(radiation_pattern_total * FARFIELD_RADIUS_UM**2) + radiation_pattern_scaled = radiation_pattern_total * FARFIELD_RADIUS_UM**2 + plot_radiation_pattern_polar(radiation_pattern_scaled) + plot_radiation_pattern_3d(radiation_pattern_scaled) ``` Focusing Properties of a Metasurface Lens diff --git a/python/examples/disc_extraction_efficiency.py b/python/examples/disc_extraction_efficiency.py index b6e5446d4..5dde21108 100644 --- a/python/examples/disc_extraction_efficiency.py +++ b/python/examples/disc_extraction_efficiency.py @@ -1,41 +1,36 @@ """Computes the extraction efficiency of a collection of dipoles in a disc. -tutorial reference: https://meep.readthedocs.io/en/latest/Python_Tutorials/Near_to_Far_Field_Spectra/#extraction-efficiency-of-a-disc-in-cylindrical-coordinates +tutorial reference: +https://meep.readthedocs.io/en/latest/Python_Tutorials/Near_to_Far_Field_Spectra/#extraction-efficiency-of-a-disc-in-cylindrical-coordinates """ import math from typing import Tuple -import matplotlib +import matplotlib.pyplot as plt import meep as mp import numpy as np -matplotlib.use("agg") -import matplotlib.pyplot as plt - - -N_DISC = 2.4 # refractive index of disc -DISC_RADIUS_UM = 1.2 # radius of disc -WAVELENGTH_UM = 1.0 # wavelength (in vacuum) - -# radius of quarter circle for computing flux in far field -FARFIELD_RADIUS_UM = 1000 * WAVELENGTH_UM -RESOLUTION = 50 # pixels/μm +RESOLUTION_UM = 50 +N_DISC = 2.4 +DISC_RADIUS_UM = 1.2 +WAVELENGTH_UM = 1.0 +NUM_FARFIELD_PTS = 200 +FARFIELD_ANGLES = np.linspace(0, 0.5 * math.pi, NUM_FARFIELD_PTS) +FARFIELD_RADIUS_UM = 1e6 * WAVELENGTH_UM +NUM_DIPOLES = 11 -def plot_radiation_pattern_polar(theta_rad: np.ndarray, radial_flux: np.ndarray): +def plot_radiation_pattern_polar(radial_flux: np.ndarray): """Plots the radiation pattern in polar coordinates. Args: - theta_rad: angles of the radiation pattern. The angles increase clockwise - with zero at the pole (+z direction) and π/2 at the equator (+r - direction). radial_flux: radial flux of the far fields in polar coordinates. """ fig, ax = plt.subplots(subplot_kw={"projection": "polar"}, figsize=(6, 6)) ax.plot( - theta_rad, + FARFIELD_ANGLES, radial_flux, "b-", ) @@ -55,20 +50,19 @@ def plot_radiation_pattern_polar(theta_rad: np.ndarray, radial_flux: np.ndarray) ) -def plot_radiation_pattern_3d(theta_rad: np.ndarray, radial_flux: np.ndarray): +def plot_radiation_pattern_3d(radial_flux: np.ndarray): """Plots the radiation pattern in 3d Cartesian coordinates. Args: - theta_rad: angles of the radiation pattern. radial_flux: radial flux of the far fields in polar coordinates. """ - phis = np.linspace(0, 2 * np.pi, num_angles) + phis = np.linspace(0, 2 * np.pi, NUM_FARFIELD_PTS) - xs = np.zeros((len(theta_rad), len(phis))) - ys = np.zeros((len(theta_rad), len(phis))) - zs = np.zeros((len(theta_rad), len(phis))) + xs = np.zeros((NUM_FARFIELD_PTS, NUM_FARFIELD_PTS)) + ys = np.zeros((NUM_FARFIELD_PTS, NUM_FARFIELD_PTS)) + zs = np.zeros((NUM_FARFIELD_PTS, NUM_FARFIELD_PTS)) - for i, theta in enumerate(theta_rad): + for i, theta in enumerate(FARFIELD_ANGLES): for j, phi in enumerate(phis): xs[i, j] = radial_flux[i] * np.sin(theta) * np.cos(phi) ys[i, j] = radial_flux[i] * np.sin(theta) * np.sin(phi) @@ -89,54 +83,57 @@ def plot_radiation_pattern_3d(theta_rad: np.ndarray, radial_flux: np.ndarray): ) -def radiation_pattern( - theta_rad: np.ndarray, sim: mp.Simulation, n2f_mon: mp.DftNear2Far -) -> np.ndarray: +def radiation_pattern(sim: mp.Simulation, n2f_mon: mp.DftNear2Far) -> np.ndarray: """Computes the radiation pattern from the near fields. Args: - theta_rad: angles of the radiation pattern. sim: a `Simulation` object. n2f_mon: a `DftNear2Far` object returned by `Simulation.add_near2far`. + + Returns: + The radiation pattern (radial flux at each angle) as a 1d array. """ - e_field = np.zeros((theta_rad.shape[0], 3), dtype=np.complex128) - h_field = np.zeros((theta_rad.shape[0], 3), dtype=np.complex128) - for n in range(num_angles): + e_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128) + h_field = np.zeros((NUM_FARFIELD_PTS, 3), dtype=np.complex128) + for n in range(NUM_FARFIELD_PTS): far_field = sim.get_farfield( n2f_mon, mp.Vector3( - FARFIELD_RADIUS_UM * math.sin(theta_rad[n]), + FARFIELD_RADIUS_UM * math.sin(FARFIELD_ANGLES[n]), 0, - FARFIELD_RADIUS_UM * math.cos(theta_rad[n]), + FARFIELD_RADIUS_UM * math.cos(FARFIELD_ANGLES[n]), ), ) - e_field[n, :] = [np.conj(far_field[j]) for j in range(3)] + e_field[n, :] = [far_field[j] for j in range(3)] h_field[n, :] = [far_field[j + 3] for j in range(3)] - flux_r = np.real(e_field[:, 1] * h_field[:, 2] - e_field[:, 2] * h_field[:, 1]) - flux_z = np.real(e_field[:, 0] * h_field[:, 1] - e_field[:, 1] * h_field[:, 0]) - flux_rz = np.sqrt(np.square(flux_r) + np.square(flux_z)) + flux_x = np.real( + np.conj(e_field[:, 1]) * h_field[:, 2] - np.conj(e_field[:, 2]) * h_field[:, 1] + ) + flux_z = np.real( + np.conj(e_field[:, 0]) * h_field[:, 1] - np.conj(e_field[:, 1]) * h_field[:, 0] + ) + flux_r = np.sqrt(np.square(flux_x) + np.square(flux_z)) - return flux_rz + return flux_r -def radiation_pattern_flux(theta_rad: np.ndarray, radial_flux: np.ndarray) -> float: +def radiation_pattern_flux(radial_flux: np.ndarray) -> float: """Computes the total flux from the radiation pattern. Based on integrating the radiation pattern over solid angles spanned by polar angles in the range of [0, π/2]. Args: - theta_rad: angles of the radiation pattern. radial_flux: radial flux of the far fields in polar coordinates. """ dphi = 2 * math.pi - dtheta_rad = theta_rad[1] - theta_rad[0] + dtheta = FARFIELD_ANGLES[1] - FARFIELD_ANGLES[0] total_flux = ( - np.sum(radial_flux * np.sin(theta_rad)) + np.sum(radial_flux * np.sin(FARFIELD_ANGLES)) * FARFIELD_RADIUS_UM**2 - * dtheta_rad + * dtheta * dphi ) @@ -144,40 +141,39 @@ def radiation_pattern_flux(theta_rad: np.ndarray, radial_flux: np.ndarray) -> fl def dipole_in_disc( - t_disc_um: float, h_disc: float, rpos: float, m: int, theta_rad: np.ndarray + disc_um: float, zpos: float, rpos_um: float, m: int ) -> Tuple[float, np.ndarray]: """Computes the total flux and radiation pattern of a dipole in a disc. Args: - t_disc_um: thickness of disc. - h_disc: height of dipole above ground plane as fraction of t_disc_um. - rpos: radial position of dipole. + disc_um: thickness of disc. + zpos: height of dipole above ground plane as fraction of disc_um. + rpos_um: radial position of dipole. m: angular φ dependence of the fields exp(imφ). - theta_rad: angles of the radiation pattern. Returns: A 2-tuple of the total flux and the radiation pattern. """ - t_pml_um = 0.5 # thickness of PML - t_air_um = 1.0 # thickness of air padding above disc - length_r_um = 5.0 # length of cell in r + pml_um = 1.0 # thickness of PML + padding_um = 1.0 # thickness of air padding above disc + r_um = 5.0 # length of cell in r frequency = 1 / WAVELENGTH_UM # center frequency of source/monitor # field decay threshold for runtime termination criteria - decay_tol = 1e-6 + field_decay_threshold = 1e-6 - size_r = length_r_um + t_pml_um - size_z = t_disc_um + t_air_um + t_pml_um + size_r = r_um + pml_um + size_z = disc_um + padding_um + pml_um cell_size = mp.Vector3(size_r, 0, size_z) boundary_layers = [ - mp.PML(t_pml_um, direction=mp.R), - mp.PML(t_pml_um, direction=mp.Z, side=mp.High), + mp.PML(pml_um, direction=mp.R), + mp.PML(pml_um, direction=mp.Z, side=mp.High), ] src_cmpt = mp.Er - src_pt = mp.Vector3(rpos, 0, -0.5 * size_z + h_disc * t_disc_um) + src_pt = mp.Vector3(rpos_um, 0, -0.5 * size_z + zpos * disc_um) sources = [ mp.Source( src=mp.GaussianSource(frequency, fwidth=0.1 * frequency), @@ -189,13 +185,13 @@ def dipole_in_disc( geometry = [ mp.Block( material=mp.Medium(index=N_DISC), - center=mp.Vector3(0.5 * DISC_RADIUS_UM, 0, -0.5 * size_z + 0.5 * t_disc_um), - size=mp.Vector3(DISC_RADIUS_UM, mp.inf, t_disc_um), + center=mp.Vector3(0.5 * DISC_RADIUS_UM, 0, -0.5 * size_z + 0.5 * disc_um), + size=mp.Vector3(DISC_RADIUS_UM, mp.inf, disc_um), ) ] sim = mp.Simulation( - resolution=RESOLUTION, + resolution=RESOLUTION_UM, cell_size=cell_size, dimensions=mp.CYLINDRICAL, m=m, @@ -209,14 +205,14 @@ def dipole_in_disc( 0, 1, mp.FluxRegion( - center=mp.Vector3(0.5 * length_r_um, 0, 0.5 * size_z - t_pml_um), - size=mp.Vector3(length_r_um, 0, 0), + center=mp.Vector3(0.5 * r_um, 0, 0.5 * size_z - pml_um), + size=mp.Vector3(r_um, 0, 0), ), mp.FluxRegion( center=mp.Vector3( - length_r_um, 0, 0.5 * size_z - t_pml_um - 0.5 * (t_air_um + t_disc_um) + r_um, 0, 0.5 * size_z - pml_um - 0.5 * (padding_um + disc_um) ), - size=mp.Vector3(0, 0, t_air_um + t_disc_um), + size=mp.Vector3(0, 0, padding_um + disc_um), ), ) @@ -226,14 +222,14 @@ def dipole_in_disc( 50, src_cmpt, src_pt, - decay_tol, + field_decay_threshold, ), ) - delta_vol = 2 * np.pi * rpos / (RESOLUTION**2) + delta_vol = 2 * np.pi * rpos_um / (RESOLUTION_UM**2) dipole_flux = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * delta_vol - dipole_radiation_pattern = radiation_pattern(theta_rad, sim, n2f_mon) + dipole_radiation_pattern = radiation_pattern(sim, n2f_mon) return dipole_flux, dipole_radiation_pattern @@ -241,89 +237,80 @@ def dipole_in_disc( if __name__ == "__main__": disc_thickness = 0.7 * WAVELENGTH_UM / N_DISC dipole_height = 0.5 - num_dipoles = 11 - dipole_rpos = np.linspace(0, DISC_RADIUS_UM, num_dipoles) + dipole_rpos_um = np.linspace(0, DISC_RADIUS_UM, NUM_DIPOLES) - # An Er source at r=0 needs to be slighty offset. + # An Er source at r = 0 needs to be slighty offset. # https://github.com/NanoComp/meep/issues/2704 - dipole_rpos[0] = 1.5 / RESOLUTION - - delta_rpos = dipole_rpos[2] - dipole_rpos[1] + dipole_rpos_um[0] = 1.5 / RESOLUTION_UM - # number of angular grid points in [0, π/2] - num_angles = 100 + delta_rpos_um = dipole_rpos_um[2] - dipole_rpos_um[1] - # grid of polar angles for computing radiated flux in far field - theta_rad = np.linspace(0, 0.5 * math.pi, num_angles) - - # r = 0 source requires a single simulation with m = ±1. + # Er source at r = 0 requires a single simulation with m = ±1. m = -1 dipole_flux, dipole_radiation_pattern = dipole_in_disc( - disc_thickness, dipole_height, dipole_rpos[0], m, theta_rad + disc_thickness, + dipole_height, + dipole_rpos_um[0], + m, ) - flux_total = dipole_flux * dipole_rpos[0] * delta_rpos - radiation_pattern_total = dipole_radiation_pattern * dipole_rpos[0] * delta_rpos + flux_total = dipole_flux * dipole_rpos_um[0] * delta_rpos_um + radiation_pattern_total = ( + dipole_radiation_pattern * dipole_rpos_um[0] * delta_rpos_um + ) print( - f"dipole:, {dipole_rpos[0]:.4f}, " - f"{radiation_pattern_flux(theta_rad, dipole_radiation_pattern):.6f}" + f"dipole:, {dipole_rpos_um[0]:.4f}, " + f"{radiation_pattern_flux(dipole_radiation_pattern):.6f}" ) - # r > 0 source requires Fourier-series expansion of φ. - flux_tol = 1e-3 # threshold flux to determine when to truncate expansion - for rpos in dipole_rpos[1:]: + # Er source at r > 0 requires Fourier-series expansion of φ. + + # Threshold flux to determine when to truncate expansion. + flux_decay_threshold = 1e-2 + + for rpos_um in dipole_rpos_um[1:]: dipole_flux_total = 0 - dipole_radiation_pattern_total = np.zeros((num_angles,)) + dipole_radiation_pattern_total = np.zeros((NUM_FARFIELD_PTS,)) dipole_radiation_pattern_flux_max = 0 m = 0 while True: dipole_flux, dipole_radiation_pattern = dipole_in_disc( - disc_thickness, dipole_height, rpos, m, theta_rad + disc_thickness, dipole_height, rpos_um, m ) - dipole_flux_total += dipole_flux if m == 0 else 2 * dipole_flux - dipole_radiation_pattern_total += ( - dipole_radiation_pattern if m == 0 else 2 * dipole_radiation_pattern + dipole_flux_total += dipole_flux * (2 if m == 0 else 1) + dipole_radiation_pattern_total += dipole_radiation_pattern * ( + 2 if m == 0 else 1 ) dipole_radiation_pattern_flux = radiation_pattern_flux( - theta_rad, dipole_radiation_pattern + dipole_radiation_pattern ) if dipole_radiation_pattern_flux > dipole_radiation_pattern_flux_max: dipole_radiation_pattern_flux_max = dipole_radiation_pattern_flux - if m > 0 and ( - (dipole_radiation_pattern_flux / dipole_radiation_pattern_flux_max) - < flux_tol + if ( + m > 0 + and (dipole_radiation_pattern_flux / dipole_radiation_pattern_flux_max) + < flux_decay_threshold ): break - print( - f"dipole-m:, {rpos:.4f}, {m}, " f"{dipole_radiation_pattern_flux:.6f}" - ) + print(f"dipole:, {rpos_um:.4f}, {m}, {dipole_radiation_pattern_flux:.6f}") m += 1 - scale_factor = 1 / (2 * (rpos / dipole_rpos[0]) ** 2) - flux_total += dipole_flux_total * scale_factor * rpos * delta_rpos + scale_factor = (dipole_rpos_um[0] / rpos_um) ** 2 + flux_total += dipole_flux_total * scale_factor * rpos_um * delta_rpos_um radiation_pattern_total += ( - dipole_radiation_pattern_total * scale_factor * rpos * delta_rpos + dipole_radiation_pattern_total * scale_factor * rpos_um * delta_rpos_um ) - dipole_radiation_pattern_total_flux = radiation_pattern_flux( - theta_rad, dipole_radiation_pattern_total - ) - print( - f"dipole:, {rpos:.4f}, {m}, " f"{dipole_radiation_pattern_total_flux:.6f}" - ) - - flux_total /= num_dipoles - radiation_pattern_total /= num_dipoles - - radiation_pattern_total_flux = radiation_pattern_flux( - theta_rad, radiation_pattern_total - ) + flux_total /= NUM_DIPOLES + radiation_pattern_total /= NUM_DIPOLES + radiation_pattern_total_flux = radiation_pattern_flux(radiation_pattern_total) extraction_efficiency = radiation_pattern_total_flux / flux_total print(f"exteff:, {extraction_efficiency:.6f}") - plot_radiation_pattern_polar(radiation_pattern_total * FARFIELD_RADIUS_UM**2) - plot_radiation_pattern_3d(radiation_pattern_total * FARFIELD_RADIUS_UM**2) + radiation_pattern_scaled = radiation_pattern_total * FARFIELD_RADIUS_UM**2 + plot_radiation_pattern_polar(radiation_pattern_scaled) + plot_radiation_pattern_3d(radiation_pattern_scaled)