It can often be tedious to create easily usable customized simulation objects from scratch, as there are various details to consider. For instance, if one wants to create an easily-usable beam source from scratch, one must create functionality to automatically interpret the source plane into field coordinates, handle conversion for polar and azimuthal angles, handle different direction specifications, etc. This all becomes much more complicated to build than Tidy3D's natively supported Gaussian beam.

Fortunately, we can take advantage of the fact that Tidy3D is open source to create complex Gaussian sources without needing to consider any of these extraneous functionalities. Namely, we use the Tidy3D BeamProfile class to handle all extra usability considerations, so we only need to define the scalar field profile for z being the propagation direction. The BeamProfile object will automatically rotate the field amplitudes according to the geometry of the source plane and specified angles.

This is just one way in which we can take advantage of the open source nature of Tidy3D, and these principles can be applied to any base object class that Tidy3D provides.

[1]:

# Standard python imports
import matplotlib.pyplot as plt
import numpy as np

# Tidy3D import
import tidy3d as td
from tidy3d import web

# Import base object
from tidy3d.components.beam import BeamProfile

Quasi Gaussian source¶

Here we show how to create a Quasi-Gaussian mode in Tidy3D using our CustomFieldSource object, based on J. Tuovinen, 1992. When a small beam waist (on the order of a wavelength or smaller) is required, this source is a better solution to the full-wave equation as opposed to the paraxial wave equation:

$E(r,z) = \frac{w_0}{w(z)}\frac{1}{{F^{\prime\prime}}^2}\exp{\Bigl(-\frac{(r/F^{\prime\prime})^2}{w(z)^2} - jkz - jkR(z)(F^{\prime\prime} - 1) + j\phi_0}\Bigr)$

where:

$\omega_0$ is the beam waist

$z_R=\frac{\pi n_\epsilon\omega_0^2}{\lambda}$ is the Rayleigh length

$\omega(z)=\omega_0\sqrt{1+\frac{z^2}{z_R^2}}$ is the beam radius

$R(z)=z(1+\frac{z_R^2}{z^2})$ is the wavefront curvature

$k=\frac{2\pi n}{\lambda}$ is the wavenumber

$F^{\prime\prime}=[1 + (r/R(z))^2]^{0.5}$

$\phi_0(z)=\arctan{\left(\frac{z}{z_R}\right)}$ is the phase shift

We then validate the field profile along the propagation direction.

[2]:

def QuasiGaussian(
    size,
    source_time,
    center=(0, 0, 0),
    direction="+",
    pol_angle=0,
    angle_phi=0,
    angle_theta=0,
    w0=1,
    waist_distance=0,
):
    freqs = [source_time.freq0]

    class QuasiGaussian_obj(BeamProfile):
        """Component for constructing OAM beam data. The normal direction is implicitly
        defined by the ``size`` parameter.
        """

        def scalar_field(self, points, background_n):
            """Scalar field for OAM beam.
            Scalar field corresponding to the analytic beam in coordinate system such that the
            propagation direction is z and the ``E``-field is entirely ``x``-polarized. The field is
            computed on an unstructured array ``points`` of shape ``(3, ...)``.
            """
            x, y, z = points
            z -= waist_distance
            r = np.sqrt(x**2 + y**2)

            k = 2 * np.pi * np.array(freqs) / td.C_0 * background_n
            zR = np.real(w0**2 * k / 2)
            wz = w0 * np.sqrt(1 + z**2 / zR**2)
            inv_Rz = z / (z**2 + zR**2)
            Fpp = np.sqrt(1 + r**2 * inv_Rz**2)
            phi0 = np.arctan2(z, zR)

            phase = -((r / Fpp) ** 2) / wz**2 - 1j * k * z + 1j * phi0
            if np.all(inv_Rz) != 0:
                phase -= 1j * k * (Fpp - 1) / inv_Rz

            E = w0 / wz * 1 / Fpp**2 * np.exp(phase)

            return E

    QuasiGaussian_obj = QuasiGaussian_obj(
        center=center,
        size=size,
        freqs=freqs,
        angle_theta=angle_theta,
        angle_phi=angle_phi,
        pol_angle=pol_angle,
        direction=direction,
    )
    field_data = QuasiGaussian_obj.field_data
    field_dataset = td.FieldDataset(
        Ex=field_data.Ex,
        Ey=field_data.Ey,
        Ez=field_data.Ez,
        Hx=field_data.Hx,
        Hy=field_data.Hy,
        Hz=field_data.Hz,
    )
    QuasiGaussian_source = td.CustomFieldSource(
        center=center, size=size, source_time=source_time, field_dataset=field_dataset
    )
    return QuasiGaussian_source

[3]:

freq0 = td.C_0 / 0.5

QG_center = (0, 0, -5)
QG_size = (5, 5, 0)
QG_source_time = td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10)
QG_w0 = td.C_0 / freq0 / 5
QG_waist_distance = -5

QG_source = QuasiGaussian(
    center=QG_center,
    size=QG_size,
    source_time=QG_source_time,
    w0=QG_w0,
    waist_distance=QG_waist_distance,
)

G_source = td.GaussianBeam(
    center=QG_center,
    size=QG_size,
    source_time=QG_source_time,
    direction="+",
    waist_radius=QG_w0,
    waist_distance=QG_waist_distance,
)

[4]:

field = td.FieldMonitor(
    size=(td.inf, 0, td.inf), name="field", freqs=[freq0 * 0.9, freq0, freq0 * 1.1]
)

QG_sim = td.Simulation(
    size=(5, 5, 15), sources=[QG_source], structures=[], monitors=[field], run_time=50 / freq0
)

G_sim = QG_sim.updated_copy(sources=[G_source])

[5]:

QG_sim_data = web.run(simulation=QG_sim, task_name="gaussian beam profile test")
G_sim_data = web.run(simulation=G_sim, task_name="gaussian beam profile test")

10:11:53 EDT Created task 'gaussian beam profile test' with task_id             
             'fdve-dae0a7c2-2ca3-47b5-9530-d2b10e10e60e' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-dae0a7c2-2ca
             3-47b5-9530-d2b10e10e60e'.

             Task folder: 'default'.

Output()

10:11:56 EDT Maximum FlexCredit cost: 0.025. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

             status = success

Output()

10:11:58 EDT loading simulation from simulation_data.hdf5

             Created task 'gaussian beam profile test' with task_id             
             'fdve-f63904a9-d6a1-464c-a209-bb091fc16af2' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-f63904a9-d6a
             1-464c-a209-bb091fc16af2'.

             Task folder: 'default'.

Output()

10:12:00 EDT Maximum FlexCredit cost: 0.025. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

             status = success

Output()

10:12:02 EDT loading simulation from simulation_data.hdf5

[6]:

fig, ax = plt.subplots(1, 2)

QG_sim_data.plot_field(field_monitor_name="field", field_name="Ex", f=freq0, ax=ax[0])
G_sim_data.plot_field(field_monitor_name="field", field_name="Ex", f=freq0, ax=ax[1])
ax[0].set_title("Quasi-Gaussian")
ax[1].set_title("Gaussian")
plt.tight_layout()
plt.show()

No description has been provided for this image

Hermite-Gaussian source¶

Here we show how to create a Hermite-Gaussian mode in Tidy3D using our CustomFieldSource object. We then validate the phase change along the propagation direction.

In cartesian coordinates, where $x$ is the propagation direction, the OAM has field profile

$E_{l,m}(x,y,z)=E_0\frac{\omega_0}{\omega(z)}H_l\Bigl(\frac{\sqrt{2}x}{\omega(z)}\Bigr) H_m\Bigl(\frac{\sqrt{2}y}{\omega(z)}\Bigr)\exp\Bigl(-\frac{x^2+y^2}{\omega^2(z)}\Bigr)\exp\Bigl(-i\frac{k(x^2+y^2)}{2R(z)}\Bigr)\exp\Bigl( i\psi(z)\Bigr)\exp( -ikz)$

where:

$\omega_0$ is the beam waist

$z_R=\frac{\pi n_\epsilon\omega_0^2}{\lambda}$ is the Rayleigh length

$\omega(z)=\omega_0\sqrt{1+\frac{z^2}{z_R^2}}$ is the beam radius

$k=\frac{2\pi n}{\lambda}$ is the wavenumber

$R(z)=z(1+\frac{z_R^2}{z^2})$ is the wavefront curvature

$H_l,H_m$ are Hermite polynomials of orders $l$ and $m$

$\psi(z)=(l+m+1)\arctan\Bigl(\frac{z}{z_R}\Bigr)$ is the phase function

[7]:

from scipy.special import hermite


def HGM(
    size,
    m,
    l,
    source_time,
    center=(0, 0, 0),
    direction="+",
    pol_angle=0,
    angle_phi=0,
    angle_theta=0,
    w0=1,
    waist_distance=0,
):
    freqs = [source_time.freq0]

    class HGM_obj(BeamProfile):
        """Component for constructing HGM beam data. The normal direction is implicitly
        defined by the ``size`` parameter.
        """

        def scalar_field(self, points, background_n):
            """Scalar field for HGM beam.
            Scalar field corresponding to the analytic beam in coordinate system such that the
            propagation direction is z and the ``E``-field is entirely ``x``-polarized. The field is
            computed on an unstructured array ``points`` of shape ``(3, ...)``.
            """
            x, y, z = points
            z += waist_distance

            zR = np.pi * background_n * w0**2 / (td.C_0 / np.array(self.freqs))
            wz = w0 * np.sqrt(1 + z**2 / zR**2)
            Rz = td.inf if np.any(z == 0) else z * (1 + zR**2 / z**2)
            k = 2 * np.pi * np.array(self.freqs) / td.C_0 * background_n
            psi_z = (l + m + 1) * np.arctan2(z, np.real(zR))

            E = source_time.amplitude * w0 / wz
            E *= hermite(l)(np.real(np.sqrt(2) * x / wz)) * hermite(m)(np.real(np.sqrt(2) * y / wz))
            E *= np.exp(-(x**2 + y**2) / wz**2) * np.exp(-1j * k * (x**2 + y**2) / (2 * Rz))
            E *= np.exp(1j * psi_z) * np.exp(-1j * k * z)

            return E

    HGM_obj = HGM_obj(
        center=center,
        size=size,
        freqs=freqs,
        angle_theta=angle_theta,
        angle_phi=angle_phi,
        pol_angle=pol_angle,
        direction=direction,
    )
    field_data = HGM_obj.field_data
    field_dataset = td.FieldDataset(
        Ex=field_data.Ex,
        Ey=field_data.Ey,
        Ez=field_data.Ez,
        Hx=field_data.Hx,
        Hy=field_data.Hy,
        Hz=field_data.Hz,
    )
    HGM_source = td.CustomFieldSource(
        center=center, size=size, source_time=source_time, field_dataset=field_dataset
    )
    return HGM_source

Create and run example Hermite-Gaussian beam:

[8]:

freq0 = td.C_0 / 1.55
hgm = HGM(size=(0, 4, 4), m=2, l=1, source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10))

[9]:

sim_size = (5.1, 5, 5)
sim_center = (sim_size[0] / 2 - 0.05, 0, 0)

num_field_monitors = 61
field_monitor_xs = np.linspace(
    sim_center[0] - sim_size[0] / 2 + 0.1, sim_size[0] - 0.1, num_field_monitors
)
field_monitors = []
for x in field_monitor_xs:
    field_monitor = td.FieldMonitor(
        center=(x, 0, 0), size=(0, sim_size[1], sim_size[2]), name=str(x), freqs=[freq0]
    )
    field_monitors.append(field_monitor)

HGM_test_sim = td.Simulation(
    center=sim_center,
    size=sim_size,
    structures=[],
    sources=[hgm],
    monitors=field_monitors,
    run_time=1e-13,
)

[10]:

hgm_sim_data = web.run(simulation=HGM_test_sim, task_name="HGM test")

10:12:03 EDT Created task 'HGM test' with task_id                               
             'fdve-b6465c5c-20b7-429e-aba4-29a00b8b28bf' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-b6465c5c-20b
             7-429e-aba4-29a00b8b28bf'.

             Task folder: 'default'.

Output()

10:12:05 EDT Maximum FlexCredit cost: 0.025. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

10:12:06 EDT status = success

Output()

10:12:08 EDT loading simulation from simulation_data.hdf5

Display field in propagation direction:

[11]:

hgm_sim_data.plot_field(str(field_monitor_xs[0]), "Ey", f=freq0)
plt.show()

To see how the field evolves, run the following code to visualize along the propagation direction.

[12]:

import ipywidgets as widgets

component = "Ey"

field_slider = widgets.FloatSlider(
    value=0,
    min=field_monitor_xs[0],
    max=field_monitor_xs[-1],
    step=None,
    description=component,
    disabled=False,
    continuous_update=True,
    orientation="horizontal",
    readout=True,
    readout_format=".3f",
)


def update_plot(x):
    if x not in field_monitor_xs:
        x = min(field_monitor_xs, key=lambda ll: abs(ll - x))
    hgm_sim_data.plot_field(str(x), component, f=freq0)
    plt.show()


widgets.interactive(update_plot, x=field_slider)

[12]:

interactive(children=(FloatSlider(value=0.05000000000000018, description='Ey', max=5.0, min=0.0500000000000001…

Laguerre-Gaussian source¶

Here we show how to create an Orbital Angular Momentum (OAM) beam in Tidy3D using our CustomFieldSource object. We then validate the phase change along the propagation direction.

In cylindrical coordinates, where $z$ is the propagation direction, the OAM has field profile

$E_{l,p}(r,\phi,z)=E_0\sqrt{\frac{2p!}{\pi(p+|l|)!}} \frac{1}{\omega(z)} \Bigl(\frac{r\sqrt{2}}{\omega(z)}\Bigr)^{|l|} \exp\Bigl(-\frac{r^2}{\omega^2(z)}\Bigr) L_p^{|l|}\Bigl(\frac{2r^2}{\omega^2(z)}\Bigr) \exp\Bigl(-ik\frac{r^2}{2R(z)}\Bigr) \exp(-il\phi) \exp\Bigl( i\psi(z)\Bigr)$

where:

$\omega_0$ is the beam waist

$z_R=\frac{\pi n_\epsilon\omega_0^2}{\lambda}$ is the Rayleigh length

$\omega(z)=\omega_0\sqrt{1+\frac{z^2}{z_R^2}}$ is the beam radius

$R(z)=z(1+\frac{z_R^2}{z^2})$ is the wavefront curvature

$k=\frac{2\pi n}{\lambda}$ is the wavenumber

$L_p^{|l|}$ is the generalized Laguerre polynomial

$\psi(z)=(|l| + 2p + 1)\arctan\Bigl(\frac{z}{z_R}\Bigr)$ is the phase function

[13]:

from scipy.special import factorial, genlaguerre


def OAM(
    size,
    p,
    l,
    source_time,
    center=(0, 0, 0),
    direction="+",
    pol_angle=0,
    angle_phi=0,
    angle_theta=0,
    w0=1,
    waist_distance=0,
):
    freqs = [source_time.freq0]

    class OAM_obj(BeamProfile):
        """Component for constructing OAM beam data. The normal direction is implicitly
        defined by the ``size`` parameter.
        """

        def scalar_field(self, points, background_n):
            """Scalar field for OAM beam.
            Scalar field corresponding to the analytic beam in coordinate system such that the
            propagation direction is z and the ``E``-field is entirely ``x``-polarized. The field is
            computed on an unstructured array ``points`` of shape ``(3, ...)``.
            """
            x, y, z = points
            z += waist_distance
            r = np.sqrt(x**2 + y**2)
            phi = np.arctan2(y, x)
            absl = np.abs(l)

            zR = np.pi * background_n * w0**2 / (td.C_0 / np.array(self.freqs))
            wz = w0 * np.sqrt(1 + z**2 / zR**2)
            Rz = td.inf if np.any(z == 0) else z * (1 + zR**2 / z**2)
            k = 2 * np.pi * np.array(self.freqs) / td.C_0 * background_n
            psi_z = (np.abs(l) + 2 * p + 1) * np.arctan2(z, np.real(zR))

            E = genlaguerre(p, absl)(2 * r**2 / wz**2) if p != 0 else 1
            E *= (
                source_time.amplitude
                * np.sqrt(2 * factorial(p) / (np.pi * factorial(p + absl)))
                / wz
            )
            E *= (r * np.sqrt(2) / wz) ** absl * np.exp(-(r**2) / wz**2)
            E *= np.exp(-1j * k * r**2 / (2 * Rz)) * np.exp(-1j * l * phi) * np.exp(1j * psi_z)

            return E

    OAM_obj = OAM_obj(
        center=center,
        size=size,
        freqs=freqs,
        angle_theta=angle_theta,
        angle_phi=angle_phi,
        pol_angle=pol_angle,
        direction=direction,
    )
    field_data = OAM_obj.field_data
    field_dataset = td.FieldDataset(
        Ex=field_data.Ex,
        Ey=field_data.Ey,
        Ez=field_data.Ez,
        Hx=field_data.Hx,
        Hy=field_data.Hy,
        Hz=field_data.Hz,
    )
    OAM_source = td.CustomFieldSource(
        center=center, size=size, source_time=source_time, field_dataset=field_dataset
    )
    return OAM_source

We will demonstrate this functionality by defining an OAM source that propagates in the $x$ direction.

[14]:

freq0 = td.C_0 / 1.55
oam = OAM(size=(0, 4, 4), p=0, l=1, source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10))

[15]:

sim_size = (5.1, 5, 5)
sim_center = (sim_size[0] / 2 - 0.05, 0, 0)

num_field_monitors = 61
field_monitor_xs = np.linspace(
    sim_center[0] - sim_size[0] / 2 + 0.1, sim_size[0] - 0.1, num_field_monitors
)
field_monitors = []
for x in field_monitor_xs:
    field_monitor = td.FieldMonitor(
        center=(x, 0, 0), size=(0, sim_size[1], sim_size[2]), name=str(x), freqs=[freq0]
    )
    field_monitors.append(field_monitor)

OAM_test_sim = td.Simulation(
    center=sim_center,
    size=sim_size,
    structures=[],
    sources=[oam],
    monitors=field_monitors,
    run_time=1e-13,
)

[16]:

oam_sim_data = web.run(simulation=OAM_test_sim, task_name="OAM test")

10:12:09 EDT Created task 'OAM test' with task_id                               
             'fdve-a562c6dc-5923-4326-b7aa-28f4b4fa5313' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-a562c6dc-592
             3-4326-b7aa-28f4b4fa5313'.

             Task folder: 'default'.

Output()

10:12:11 EDT Maximum FlexCredit cost: 0.025. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

10:12:12 EDT status = success

Output()

10:12:14 EDT loading simulation from simulation_data.hdf5

Display field in propagation direction:

[17]:

oam_sim_data.plot_field(str(field_monitor_xs[30]), "E", f=freq0, val="abs^2")
plt.show()

To see how the field amplitude evolves, run the following code to visualize along the propagation direction.

[18]:

import ipywidgets as widgets

component = "Ey"

field_slider = widgets.FloatSlider(
    value=0,
    min=field_monitor_xs[0],
    max=field_monitor_xs[-1],
    step=None,
    description=component,
    disabled=False,
    continuous_update=True,
    orientation="horizontal",
    readout=True,
    readout_format=".3f",
)


def update_plot(x):
    if x not in field_monitor_xs:
        x = min(field_monitor_xs, key=lambda ll: abs(ll - x))
    oam_sim_data.plot_field(str(x), component, f=freq0)
    plt.show()


widgets.interactive(update_plot, x=field_slider)

[18]:

interactive(children=(FloatSlider(value=0.05000000000000018, description='Ey', max=5.0, min=0.0500000000000001…

Bessel Beam¶

Here we show how to create an Bessel beam in Tidy3D using our CustomFieldSource object. We then validate the phase change along the propagation direction.

In cylindrical coordinates $(\rho, \phi, z)$, where $z$ is the propagation direction, the Bessel beam has field profile

$U_\beta(r)=J_m(k_T\rho)e^{im\phi}e^{-i\beta z}$

where:

$k=\frac{2\pi n}{\lambda}$ is the wavenumber

$J_m$ is the Bessel function of the first kind and $m^{th}$ order

$k_T=\sqrt{k^2-\beta^2}$

Source: Saleh & Teich, Fundamentals of Photonics

[19]:

from scipy.special import jv


def Bessel(
    size,
    m,
    beta,
    source_time,
    center=(0, 0, 0),
    direction="+",
    pol_angle=0,
    angle_phi=0,
    angle_theta=0,
    w0=1,
    waist_distance=0,
):
    if not isinstance(m, int):
        raise Exception("Error: m must be an integer")
    freqs = [source_time.freq0]

    class Bessel_obj(BeamProfile):
        """Component for constructing Bessel beam data. The normal direction is implicitly
        defined by the ``size`` parameter.
        """

        def scalar_field(self, points, background_n):
            """Scalar field for Bessel beam.
            Scalar field corresponding to the analytic beam in coordinate system such that the
            propagation direction is z and the ``E``-field is entirely ``x``-polarized. The field is
            computed on an unstructured array ``points`` of shape ``(3, ...)``.
            """
            x, y, z = points
            z += waist_distance
            rho = np.sqrt(x**2 + y**2)
            phi = np.arctan2(y, x)
            beta % (2 * np.pi)

            k = 2 * np.pi * np.array(self.freqs) / td.C_0 * background_n
            k_T = np.sqrt(k**2 - beta**2)

            E = source_time.amplitude * jv(m, k_T * rho)
            E *= np.exp(1j * m * phi) * np.exp(-1j * beta * z)

            return E

    Bessel_obj = Bessel_obj(
        center=center,
        size=size,
        freqs=freqs,
        angle_theta=angle_theta,
        angle_phi=angle_phi,
        pol_angle=pol_angle,
        direction=direction,
    )
    field_data = Bessel_obj.field_data
    field_dataset = td.FieldDataset(
        Ex=field_data.Ex,
        Ey=field_data.Ey,
        Ez=field_data.Ez,
        Hx=field_data.Hx,
        Hy=field_data.Hy,
        Hz=field_data.Hz,
    )
    Bessel_source = td.CustomFieldSource(
        center=center, size=size, source_time=source_time, field_dataset=field_dataset
    )
    return Bessel_source

Create and run example Bessel beam:

[20]:

freq0 = td.C_0 / 1.55
bessel = Bessel(
    size=(0, 4, 4), m=0, beta=0, source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10)
)

[21]:

sim_size = (5.1, 5, 5)
sim_center = (sim_size[0] / 2 - 0.05, 0, 0)

num_field_monitors = 61
field_monitor_xs = np.linspace(
    sim_center[0] - sim_size[0] / 2 + 0.1, sim_size[0] - 0.1, num_field_monitors
)
field_monitors = []
for x in field_monitor_xs:
    field_monitor = td.FieldMonitor(
        center=(x, 0, 0), size=(0, sim_size[1], sim_size[2]), name=str(x), freqs=[freq0]
    )
    field_monitors.append(field_monitor)

Bessel_test_sim = td.Simulation(
    center=sim_center,
    size=sim_size,
    structures=[],
    sources=[bessel],
    monitors=field_monitors,
    run_time=1e-13,
)

[22]:

bessel_sim_data = web.run(simulation=Bessel_test_sim, task_name="Bessel test")

10:12:15 EDT Created task 'Bessel test' with task_id                            
             'fdve-d9e0774e-934a-4f85-8cf2-9a6f71dbc4ef' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d9e0774e-934
             a-4f85-8cf2-9a6f71dbc4ef'.

             Task folder: 'default'.

Output()

10:12:17 EDT Maximum FlexCredit cost: 0.025. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

10:12:18 EDT status = success

Output()

10:12:19 EDT loading simulation from simulation_data.hdf5

Display field amplitude in propagation direction:

[23]:

bessel_sim_data.plot_field(str(field_monitor_xs[0]), "E", val="abs", f=freq0)
plt.show()

To see how the field amplitude evolves, run the following code to visualize along the propagation direction.

[24]:

import ipywidgets as widgets

component = "E"

field_slider = widgets.FloatSlider(
    value=0,
    min=field_monitor_xs[0],
    max=field_monitor_xs[-1],
    step=None,
    description=component,
    disabled=False,
    continuous_update=True,
    orientation="horizontal",
    readout=True,
    readout_format=".3f",
)


def update_plot(x):
    if x not in field_monitor_xs:
        x = min(field_monitor_xs, key=lambda ll: abs(ll - x))
    bessel_sim_data.plot_field(str(x), component, val="abs", f=freq0)
    plt.show()


widgets.interactive(update_plot, x=field_slider)

[24]:

interactive(children=(FloatSlider(value=0.05000000000000018, description='E', max=5.0, min=0.05000000000000018…

Ince-Gaussian modes¶

Here we show how to create Ince-Gaussian (IG) beams in Tidy3D using our CustomFieldSource object. We then validate the phase change along the propagation direction.

In elliptic-cylindrical coordinates $(\xi,\eta)$, where

$x=\sqrt{\varepsilon/2}\omega(z)\cosh\xi\cos\eta$
$y=\sqrt{\varepsilon/2}\omega(z)\sinh\xi\sin\eta$, the IGB has field profile $u_\varepsilon(\xi,\eta,z)=\frac{\omega_0}{\omega(z)}C_p^m(i\xi,\varepsilon)C_p^m(\eta,\varepsilon)\exp\Bigr[-\frac{r^2}{\omega(z)^2} + ik\frac{r^2}{2R(z)}\Bigr]\exp(ikz)\exp\Bigl(-i\psi(z)\Bigr)$

where:

$k=\frac{2\pi n}{\lambda}$ is the wavenumber

$\omega_0$ is the beam waist

$z_R=\frac{\pi n_\epsilon\omega_0^2}{\lambda}$ is the Rayleigh length

$\omega(z)=\omega_0\sqrt{1+\frac{z^2}{z_R^2}}$ is the beam radius

$R(z)=z(1+\frac{z_R^2}{z^2})$ is the wavefront curvature

$z_R=\frac{\pi n_\epsilon\omega_0^2}{\lambda}$ is the Rayleigh length

$\psi(z)=(p+1)\arctan(\frac{z}{z_R})$ is the phase function

$C_p^m(\cdot,\varepsilon)$ is an Ince-Gaussian polynomial of order $p$ and degree $m$ with ellipticity $\varepsilon$. Note that $m$ and $p$ must have the same parity with $m\leq p$ and $\varepsilon\in(0,\infty)$

We also add a visualization option for the beam to ensure that the correct Ince polynomial is solved for.

[25]:

from scipy.interpolate import RectBivariateSpline
from scipy.special import gamma


def InceGaussian(
    size,
    m,
    p,
    e,
    source_time,
    center=(0, 0, 0),
    direction="+",
    pol_angle=0,
    angle_phi=0,
    angle_theta=0,
    w0=1,
    waist_distance=0,
    num_sampling_points=501,
    normalize=True,
    display_profile=False,
):
    if not isinstance(m, int):
        raise TypeError("Degree m must be an integer.")
    if not isinstance(p, int):
        raise TypeError("Order p must be an integer.")
    if m < 0 or p < 0:
        raise ValueError("m and p must be nonnegative.")
    if (p - m) % 2 == 1:
        raise ValueError("m and p must have the same parity.")
    if m > p:
        raise ValueError("m must be less than or equal to p.")
    if e <= 0:
        raise ValueError("Ellipticity must be greater than 0.")

    freqs = [source_time.freq0]
    L = max(size) / 2
    N = num_sampling_points if num_sampling_points % 2 == 1 else num_sampling_points + 1
    x = np.linspace(-L, L, N)
    y = np.linspace(-L, L, N)
    X, Y = np.meshgrid(x, y)
    R = np.sqrt(X**2 + Y**2)

    class IG_obj(BeamProfile):
        """Component for constructing Ince-Gaussian beam data. The normal direction is implicitly
        defined by the ``size`` parameter.

        from https://www.mathworks.com/matlabcentral/fileexchange/46222-ince-gaussian-beam

        Copyright (c) 2014, Miguel A. Bandres
        All rights reserved.

        Redistribution and use in source and binary forms, with or without
        modification, are permitted provided that the following conditions are met:

        * Redistributions of source code must retain the above copyright notice, this
          list of conditions and the following disclaimer.

        * Redistributions in binary form must reproduce the above copyright notice,
          this list of conditions and the following disclaimer in the documentation
          and/or other materials provided with the distribution
        THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
        AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
        IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
        DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
        FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
        DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
        SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
        CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
        OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
        OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
        """

        def scalar_field(self, points, background_n):
            x_p, y_p, z = points
            z += waist_distance

            wvl0 = td.C_0 / source_time.freq0
            k = 2 * np.pi * background_n / wvl0
            zR = np.real(np.pi * background_n * w0**2 / wvl0)
            wz = w0 * np.sqrt(1 + z[0] ** 2 / zR**2)

            def mesh_elliptic(f_m):
                xi_m, eta_m = np.zeros((N, N)), np.zeros((N, N))

                elliptic_complex = np.arccosh(
                    (
                        X[0 : (N + 1) // 2, (N - 1) // 2 : N]
                        + 1j * Y[0 : (N + 1) // 2, (N - 1) // 2 : N]
                    )
                    / f_m
                )
                xi0 = np.real(elliptic_complex)
                eta0 = np.imag(elliptic_complex)
                eta0 = np.where(eta0 < 0, eta0 + 2 * np.pi, eta0)

                xi_m[: (N + 1) // 2, (N - 1) // 2 :] = xi0
                xi_m[: (N + 1) // 2, : (N - 1) // 2] = xi_m[: (N + 1) // 2, (N + 3) // 2 - 1 :][
                    :, ::-1
                ]
                xi_m[(N + 3) // 2 - 1 :, :] = xi_m[: (N - 1) // 2, :][::-1, :]

                eta_m[: (N + 1) // 2, (N - 1) // 2 :] = eta0
                eta_m[: (N + 1) // 2, : (N - 1) // 2] = (
                    np.pi - eta_m[: (N + 1) // 2, (N + 3) // 2 - 1 :][:, ::-1]
                )
                eta_m[(N + 3) // 2 - 1 :, :] = np.pi + eta_m[: (N - 1) // 2, :][::-1, ::-1]

                return xi_m, eta_m

            def C(e_coord):  # Even Ince polynomials
                if p == 0 and m == 0:
                    return 1 + 0j, [1 + 0j], 0j
                if np.isscalar(e_coord):
                    e_coord = np.asarray(e_coord)
                shape = e_coord.shape
                if len(shape) == 1:
                    largo, ancho = (shape[0], 1)
                elif len(shape) >= 2:
                    largo, ancho = (shape[0], shape[1])  # largo = rows, ancho = cols
                else:  # scalar or empty
                    largo, ancho = (1, 1)
                e_coord = e_coord.flatten(order="F")  # MATLAB uses column-major order
                e_coord = e_coord.reshape(1, -1)

                if p % 2 == 0:
                    j, N_c, n = p / 2, p / 2 + 1, int(m / 2 + 1)
                    M = np.diag(e * (j + np.arange(1, N_c)), k=1)
                    M += np.diag(
                        np.concatenate(([2 * e * j], e * (j - np.arange(1, N_c - 1)))), k=-1
                    )
                    M += np.diag(np.concatenate(([0], 4 * ((np.arange(0, N_c - 1) + 1) ** 2))))

                    ets, A = np.linalg.eig(M)  # ets = eigenvalues, A = eigenvectors
                    index = np.argsort(ets)  # indices to sort eigenvalues
                    ets = ets[index]  # sorted eigenvalues
                    A = A[:, index]

                    mv = np.arange(2, p, 2)
                    N2 = np.sqrt(
                        2 * A[0, n - 1] ** 2 * gamma(p / 2 + 1) ** 2
                        + np.sum(
                            (
                                np.sqrt(gamma((p + mv) / 2 + 1) * gamma((p - mv) / 2 + 1))
                                * A[2 : int(p / 2 + 2), n - 1]
                            )
                            ** 2
                        )
                    )
                    NS = np.sign(np.sum(A[:, n - 1]))
                    A = A / N2 * NS

                    rr = np.arange(0, N_c)
                    Rr, Xr = np.meshgrid(rr, e_coord)
                    IP = np.cos(2 * Xr * Rr) @ A[:, n - 1]
                    IP = np.reshape(IP, (largo, ancho), order="F")

                    dIP = -2 * Rr * np.sin(2 * Xr * Rr) @ A[:, n - 1]
                    dIP = np.reshape(dIP, (largo, ancho), order="F")

                else:  # p is odd
                    j, N_c, n = int((p - 1) / 2), (p - 1) / 2 + 1, int((m + 1) / 2)
                    M = np.diag(e / 2 * (p + (2 * np.arange(0, N_c - 1) + 2)), k=1)
                    M += np.diag(e / 2 * (p - (2 * np.arange(1, N_c) - 1)), k=-1)
                    M += np.diag(
                        np.concatenate(([e / 2 + p * e / 2 + 1], (2 * np.arange(1, N_c) + 1) ** 2))
                    )

                    ets, A = np.linalg.eig(M)  # ets = eigenvalues, A = eigenvectors
                    index = np.argsort(ets)  # indices to sort eigenvalues
                    ets = ets[index]  # sorted eigenvalues
                    A = A[:, index]

                    mv = np.arange(1, p + 1, 2)
                    N2 = np.sqrt(
                        np.sum(
                            (
                                np.sqrt(gamma((p + mv) / 2 + 1) * gamma((p - mv) / 2 + 1))
                                * A[:, n - 1]
                            )
                            ** 2
                        )
                    )
                    NS = np.sign(np.sum(A[:, n - 1]))
                    A = A / N2 * NS

                    rr = 2 * np.arange(0, N_c) + 1
                    Rr, Xr = np.meshgrid(rr, e_coord)
                    IP = np.cos(Xr * Rr) @ A[:, n - 1]
                    IP = np.reshape(IP, (largo, ancho), order="F")

                    dIP = -Rr * np.sin(Xr * Rr) @ A[:, n - 1]
                    dIP = np.reshape(dIP, (largo, ancho), order="F")

                return IP, A[:, n - 1], dIP

            f0 = np.sqrt(e / 2) * w0

            if z.any() == 0:
                xi, eta = mesh_elliptic(f0)
                Ince_eta, _, _ = C(eta)
                Ince_ixi, _, _ = C(1j * xi)
                IGB = Ince_eta * Ince_ixi * np.exp(-((R / w0) ** 2))
            else:
                Rz = z[0] * (1 + (zR / z[0]) ** 2)
                f = f0 * wz / w0
                xi, eta = mesh_elliptic(f)
                Ince_eta, _, _ = C(eta)
                Ince_ixi, _, _ = C(1j * xi)
                IGB = (w0 / wz) * Ince_eta * Ince_ixi * np.exp(-((R / wz) ** 2))
                IGB *= np.exp(
                    1j * (k * z[0] + k * R**2 / (2 * Rz) - (p + 1) * np.arctan2(z[0], zR))
                )

            C0, coef, _ = C(0)
            Cpi2, _, DCpi2 = C(np.pi / 2)

            if normalize:
                if p % 2 == 0:
                    norm = (
                        (-1) ** (m / 2)
                        * np.sqrt(2)
                        * gamma(p / 2 + 1)
                        * coef[0]
                        * np.sqrt(2 / np.pi)
                        / w0
                        / C0
                        / Cpi2
                    )
                else:
                    norm = (
                        (-1) ** ((m + 1) / 2)
                        * gamma((p + 1) / 2 + 1)
                        * np.sqrt(4 * e / np.pi)
                        * coef[0]
                        / w0
                        / C0
                        / DCpi2
                    )

                IGB = IGB * norm

            # interpolate 2D grid to points' coordinates
            IGB = np.real(IGB)

            shifted_points_x = x_p.reshape(-1) - np.min(np.abs(x_p))
            shifted_points_y = y_p.reshape(-1) - np.min(np.abs(y_p))
            shifted_points_X, shifted_points_Y = np.meshgrid(shifted_points_x, shifted_points_y)
            IGB_interp = RectBivariateSpline(x, y, IGB)

            sqrt = int(np.sqrt(len(x_p)))
            square_x = shifted_points_x.reshape(sqrt, sqrt)
            square_y = shifted_points_y.reshape(sqrt, sqrt)

            IGB_recast = IGB_interp(square_x[:, 0], square_y[0, :]).T

            # check profile
            if display_profile:
                plt.imshow(np.abs(IGB_recast).T, cmap="hot")

            return IGB_recast.ravel().reshape(-1, 1)

        plt.show()

    IG_obj = IG_obj(
        center=center,
        size=size,
        freqs=freqs,
        angle_theta=angle_theta,
        angle_phi=angle_phi,
        pol_angle=pol_angle,
        direction=direction,
    )
    field_data = IG_obj.field_data
    field_dataset = td.FieldDataset(
        Ex=field_data.Ex,
        Ey=field_data.Ey,
        Ez=field_data.Ez,
        Hx=field_data.Hx,
        Hy=field_data.Hy,
        Hz=field_data.Hz,
    )
    IG_source = td.CustomFieldSource(
        center=center, size=size, source_time=source_time, field_dataset=field_dataset
    )
    return IG_source

[26]:

freq0 = td.C_0 / 0.532
IGB_x = -1
IGB_99 = InceGaussian(
    center=(IGB_x, 0, 0),
    size=(0, 6, 6),
    m=9,
    p=9,
    e=2,
    source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10),
    waist_distance=IGB_x,
    display_profile=True,
)

[27]:

IGB = InceGaussian(
    center=(IGB_x, 0, 0),
    size=(0, 6, 6),
    m=4,
    p=10,
    e=2,
    source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10),
    waist_distance=IGB_x,
    display_profile=True,
)

[28]:

sim_size = (10.5, 10, 10)
sim_center = (IGB_x + sim_size[0] / 2 - 0.1, 0, 0)

num_field_monitors = 61
field_monitor_xs = np.linspace(IGB_x + 0.1, IGB_x + sim_size[0] / 2 - 0.1, num_field_monitors)
field_monitors = []
for x in field_monitor_xs:
    field_monitor = td.FieldMonitor(
        center=(x, 0, 0), size=(0, sim_size[1], sim_size[2]), name=str(x), freqs=[freq0]
    )
    field_monitors.append(field_monitor)

auto = td.AutoGrid(min_steps_per_wvl=40)

Ince_Gaussian_test_sim = td.Simulation(
    center=sim_center,
    size=sim_size,
    structures=[],
    sources=[IGB],
    monitors=field_monitors,
    run_time=1e-13,
    grid_spec=td.GridSpec(grid_x=auto, grid_y=auto, grid_z=auto),
)

[29]:

ince_gaussian_sim_data = web.run(simulation=Ince_Gaussian_test_sim, task_name="Ince-Gaussian test")

10:13:04 EDT Created task 'Ince-Gaussian test' with task_id                     
             'fdve-5f31a9c6-cc0b-4aa3-ac06-e22797c3657c' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-5f31a9c6-cc0
             b-4aa3-ac06-e22797c3657c'.

             Task folder: 'default'.

Output()

10:13:07 EDT Maximum FlexCredit cost: 0.586. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

10:13:08 EDT status = success

Output()

10:14:19 EDT loading simulation from simulation_data.hdf5

[30]:

ince_gaussian_sim_data.plot_field(str(field_monitor_xs[0]), "E", val="abs", f=freq0)
plt.show()

[31]:

import ipywidgets as widgets

component = "E"

field_slider = widgets.FloatSlider(
    value=0,
    min=field_monitor_xs[0],
    max=field_monitor_xs[-1],
    step=None,
    description=component,
    disabled=False,
    continuous_update=True,
    orientation="horizontal",
    readout=True,
    readout_format=".3f",
)


def update_plot(x):
    if x not in field_monitor_xs:
        x = min(field_monitor_xs, key=lambda ll: abs(ll - x))
    ince_gaussian_sim_data.plot_field(str(x), component, val="abs", f=freq0)
    plt.show()


widgets.interactive(update_plot, x=field_slider)

[31]:

interactive(children=(FloatSlider(value=0.0, description='E', max=4.15, min=-0.9, readout_format='.3f', step=N…

Hyper-geometric Gaussian modes¶

Here we show how to create a Hyper-geometric Gaussian (HGG) beam in Tidy3D using our CustomFieldSource object. We then validate the phase change along the propagation direction.

In cylindrical coordinates, where $\rho=\frac{r}{\omega_0}$ is the normalized radial coordinate and $Z=\frac{z}{z_R}$ is the normalized longitude direction (the direction of propagation), the HGG has field profile

$E_{p,m}\Bigl(\rho,\phi,Z\Bigr)=\sqrt{\frac{2^{p+|m|+1}}{\pi\Gamma(p+|m|+1)}}\frac{\Gamma\Bigl(\frac{p}{2}+|m|+1\Bigr)}{\Gamma(|m|+1)}i^{|m|+1}Z^{\frac{p}{2}}(Z+i)^{-\Bigl(\frac{p}{2}+|m|+1\Bigr)}\rho^{|m|}\exp\Bigl(-\frac{i\rho^2}{Z+i}\Bigr) e^{im\phi} {_1}F_1\Bigl(-\frac{p}{2},|m|+1;\frac{\rho^2}{Z(Z+i)}\Bigr)$

where:

$\omega_0$ is the beam waist

$z_R=\frac{\pi n_\epsilon\omega_0^2}{\lambda}$ is the Rayleigh length

$m$ is the rotational index

$p\geq-|m|$

$\Gamma(x)$ is the gamma function

${_1}F_1(a,b;x)$ is a confluent hypergeometric function

[32]:

from scipy.special import hyp1f1


def HGG(
    size,
    p,
    m,
    source_time,
    center=(0, 0, 0),
    direction="+",
    pol_angle=0,
    angle_phi=0,
    angle_theta=0,
    w0=1,
    waist_distance=0,
):
    if p < -np.abs(m):
        raise Exception("Error: p must be greater than or equal to -|m|")
    freqs = [source_time.freq0]

    class HGG_obj(BeamProfile):
        """Component for constructing HGG beam data. The normal direction is implicitly
        defined by the ``size`` parameter.
        """

        def scalar_field(self, points, background_n):
            """Scalar field for HGG beam.
            Scalar field corresponding to the analytic beam in coordinate system such that the
            propagation direction is z and the ``E``-field is entirely ``x``-polarized. The field is
            computed on an unstructured array ``points`` of shape ``(3, ...)``.
            """
            x, y, z = points
            z += waist_distance

            r = np.sqrt(x**2 + y**2) / w0
            phi = np.arctan2(y, x)
            abs_m = np.abs(m)

            E = source_time.amplitude * np.sqrt(
                2 ** (p + abs_m + 1) / (np.pi * gamma(p + abs_m + 1))
            )

            Z = z / (np.pi * background_n * w0**2 / (td.C_0 / np.array(self.freqs)))
            if np.all(Z) == 0:
                return E * r ** (p + abs_m) * np.exp(-(r**2) + 1j * m * phi)

            E *= gamma(p / 2 + abs_m + 1) / gamma(abs_m + 1) * 1j ** (abs_m + 1) * Z ** (p / 2)
            E *= (Z + 1j) ** (-(p / 2 + abs_m + 1)) * r**abs_m * np.exp(-(1j * r**2) / (Z + 1j))
            E *= np.exp(1j * m * phi) * hyp1f1(-p / 2, abs_m + 1, r**2 / (Z * (Z + 1j)))

            return E

    HGG_obj = HGG_obj(
        center=center,
        size=size,
        freqs=freqs,
        angle_theta=angle_theta,
        angle_phi=angle_phi,
        pol_angle=pol_angle,
        direction=direction,
    )
    field_data = HGG_obj.field_data
    field_dataset = td.FieldDataset(
        Ex=field_data.Ex,
        Ey=field_data.Ey,
        Ez=field_data.Ez,
        Hx=field_data.Hx,
        Hy=field_data.Hy,
        Hz=field_data.Hz,
    )
    HGG_source = td.CustomFieldSource(
        center=center, size=size, source_time=source_time, field_dataset=field_dataset
    )
    return HGG_source

Create and run example HGG beam:

[33]:

freq0 = td.C_0 / 1.55
hgg = HGG(size=(0, 4, 4), p=1, m=2, source_time=td.GaussianPulse(freq0=freq0, fwidth=freq0 / 10))

[34]:

sim_size = (5.1, 5, 5)
sim_center = (sim_size[0] / 2 - 0.05, 0, 0)

num_field_monitors = 61
field_monitor_xs = np.linspace(
    sim_center[0] - sim_size[0] / 2 + 0.1, sim_size[0] - 0.1, num_field_monitors
)
field_monitors = []
for x in field_monitor_xs:
    field_monitor = td.FieldMonitor(
        center=(x, 0, 0), size=(0, sim_size[1], sim_size[2]), name=str(x), freqs=[freq0]
    )
    field_monitors.append(field_monitor)

HGG_test_sim = td.Simulation(
    center=sim_center,
    size=sim_size,
    structures=[],
    sources=[hgg],
    monitors=field_monitors,
    run_time=1e-13,
)

[35]:

hgg_sim_data = web.run(simulation=HGG_test_sim, task_name="HGG test")

10:14:22 EDT Created task 'HGG test' with task_id                               
             'fdve-d412f0b6-b6f9-42c0-bb26-c60d20ef004e' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d412f0b6-b6f
             9-42c0-bb26-c60d20ef004e'.

             Task folder: 'default'.

Output()

10:14:25 EDT Maximum FlexCredit cost: 0.025. Minimum cost depends on task       
             execution details. Use 'web.real_cost(task_id)' to get the billed  
             FlexCredit cost after a simulation run.

10:14:26 EDT status = success

Output()

10:14:31 EDT loading simulation from simulation_data.hdf5

Display field in propagation direction:

[36]:

hgg_sim_data.plot_field(str(field_monitor_xs[30]), "Ey", f=freq0)
plt.show()

To see how the field evolves, run the following code to visualize along the propagation direction.

[37]:

import ipywidgets as widgets

component = "Ey"

field_slider = widgets.FloatSlider(
    value=0,
    min=field_monitor_xs[0],
    max=field_monitor_xs[-1],
    step=None,
    description=component,
    disabled=False,
    continuous_update=True,
    orientation="horizontal",
    readout=True,
    readout_format=".3f",
)


def update_plot(x):
    if x not in field_monitor_xs:
        x = min(field_monitor_xs, key=lambda ll: abs(ll - x))
    hgg_sim_data.plot_field(str(x), component, f=freq0)
    plt.show()


widgets.interactive(update_plot, x=field_slider)

[37]:

interactive(children=(FloatSlider(value=0.05000000000000018, description='Ey', max=5.0, min=0.0500000000000001…

[ ]:

TIDY3D
LEARNING CENTER

Advanced Gaussian sources¶

Quasi Gaussian source¶

Hermite-Gaussian source¶

Laguerre-Gaussian source¶

Bessel Beam¶

Ince-Gaussian modes¶

Hyper-geometric Gaussian modes¶

TIDY3D LEARNING CENTER

Advanced Gaussian sources¶

Quasi Gaussian source¶

Hermite-Gaussian source¶

Laguerre-Gaussian source¶

Bessel Beam¶

Ince-Gaussian modes¶

Hyper-geometric Gaussian modes¶

TIDY3D
LEARNING CENTER