Oblique planewaves are discontinuous around the boundary under periodic boundary conditions.

[radillus]

Hello,
When I set the oblique planewave source following this FAQ, I notice that one shoud set the amp_func as exp(i*2*pi*k.dot(x)) rather than exp(i*k.dot(x)) when using the periodic boundary conditions just like this tutorial does. Otherwise it will cause discontinuous wave as below shows. Another file this FAQ metioned does not set the k in its simulation so it will not be a problem at that case.

setting with exp(i*k.dot(x))

exp(i*2*pi*k.dot(x))

Here is my test code

import numpy as np
import meep as mp
mp.verbosity(0)
import matplotlib.pyplot as plt
from matplotlib.axes import Axes
from matplotlib.animation import FuncAnimation


air = mp.Medium(index=1)

default_material = air

# assume propagating along z axis

dpml = 1.0
gh = 2.0
dpad = 0.5
cell_xy_size = 5.0
cell_z_size = dpml + dpad + gh + dpad + dpml

resolution = 20

wavelength = 1
frequency = 1/wavelength

# polar_angle is the angle between k vector and z axis (between 0 and π/2)
# azimuth_angle is the angle between k vector projection on xy plane and x axis (between 0 and 2π)
polar_angle = np.pi/6
azimuth_angle = 0
k_point = mp.Vector3(
    x=np.sin(polar_angle)*np.cos(azimuth_angle),
    y=np.sin(polar_angle)*np.sin(azimuth_angle),
    z=np.cos(polar_angle),
)

# ------------------------------------------------------------------------
# case A
k_point = k_point.unit().scale(frequency)

def pw_amp(k: mp.Vector3, x0: mp.Vector3):
    def _pw_amp(x):
        return np.exp(1j * 2*np.pi * k.dot(x-x0))
    return _pw_amp
# ------------------------------------------------------------------------
# case B
k_point = k_point.unit().scale(2*np.pi*frequency)

def pw_amp(k: mp.Vector3, x0: mp.Vector3):
    def _pw_amp(x):
        return np.exp(1j * k.dot(x-x0))
    return _pw_amp
# ------------------------------------------------------------------------

source_center = mp.Vector3(0, 0, -0.5*cell_z_size+dpml+dpad)
source_size = mp.Vector3(cell_xy_size, cell_xy_size, 0)

cell_size = mp.Vector3(cell_xy_size, cell_xy_size, cell_z_size)

boundary_layers = [mp.PML(thickness=dpml, direction=mp.Z)]

# also applied for gaussian source
sources = [
    mp.Source(
        src=mp.ContinuousSource(
            frequency,
            fwidth=0.1*frequency,
            is_integrated=True,
        ),
        component=mp.Ex,
        center=source_center,
        size=source_size,
        amp_func=pw_amp(k_point, source_center),
    )
]

simulation = mp.Simulation(
    cell_size=cell_size,
    resolution=resolution,
    geometry=[],
    sources=sources,
    force_complex_fields=True,
    boundary_layers=boundary_layers,
    default_material=default_material,
    k_point=k_point,
)

components = [mp.Ex, mp.Ey, mp.Ez, mp.Hx, mp.Hy, mp.Hz]
components_names = ['Ex', 'Ey', 'Ez', 'Hx', 'Hy', 'Hz']


simulation_end_time = 30
simulation_step_time = 0.5

# colunm format: [Ex, Ey, Ez, Hx, Hy, Hz]
# row format: [real, imag, abs, angle]
animation_time = 10  # seconds
frames = int(simulation_end_time/simulation_step_time) + 1
fps = frames//animation_time
fig = plt.figure(figsize=(20, 12))

def draw_component(name: str, field: np.ndarray, axes: list[Axes]):
    plt.colorbar(axes[0].imshow(field.real))
    axes[0].set_title(f'{name} real')
    plt.colorbar(axes[1].imshow(field.imag))
    axes[1].set_title(f'{name} imag')
    plt.colorbar(axes[2].imshow(np.abs(field)))
    axes[2].set_title(f'{name} abs')
    plt.colorbar(axes[3].imshow(np.angle(field), cmap='hsv', vmin=-np.pi, vmax=np.pi))
    axes[3].set_title(f'{name} angle')

def draw_frame(frame_index: int):
    fig.clf()

    print(f'frame {frame_index}/{frames}')
    fig.suptitle(f'{frame_index*simulation_step_time:.1f} us')

    if frame_index == 0:
        simulation.run(until=0)
    else:
        simulation.run(until=simulation_step_time)

    axes = fig.subplots(4, 6)
    for cn, c, a in zip(components_names, components, axes.T):
        c_data = simulation.get_array(
            component=c,
            center=mp.Vector3(0, 0, 0),
            size=cell_size,
        )
        draw_component(
            cn,
            c_data[:, c_data.shape[1]//2, :],
            a,
        )

animation = FuncAnimation(fig, draw_frame, frames=frames, cache_frame_data=False)
animation.save('oblique_source@xz_B.mp4', fps=fps)

[stevengj]

It depends on how you define k. Conventionally, in math, the k vector already includes the 2π. However, when you set the boundary conditions by k_point in Meep, it is divided by 2π, so if you use the k_point vector then you need to multiply by 2π in your amp_fun.

[stevengj]

I updated the FAQ fix the statement to say that k_point = k/2π