Skip to content

Beams API

beams

Gaussian beams on planar stacks: Goos–Hänchen and Imbert–Fedorov shifts.

A finite beam is a superposition of plane waves. On a planar stack each keeps its in-plane wavevector q, so the beam's angular spectrum passes through component by component, each with the plane-wave S-matrix of :func:hyperbolic_optics.bcm.stack_smatrix at its own q. :func:reflected_beam rebuilds the reflected field on the interface plane, and :func:beam_shift measures its centroid: the Goos–Hänchen shift along the plane of incidence and the Imbert–Fedorov shift across it, which anisotropic and hyperbolic media make large and asymmetric.

Functions:

  • reflected_beam –

    The reflected beam's field on the reference plane, from its angular spectrum.

  • beam_shift –

    Centroid (Δx, Δy) of |field|²: the Goos–Hänchen and Imbert–Fedorov shifts.

reflected_beam

reflected_beam(
    eps,
    thickness_um,
    wavelength_um,
    angle_rad,
    eps_top,
    eps_bottom,
    waist_um,
    polarization=(1.0, 0.0),
    mu=None,
    samples=61,
    span=5.0,
    points=129,
    extent=4.0,
)

The reflected beam's field on the reference plane, from its angular spectrum.

Parameters:

  • polarization –

    The incident beam's Jones vector (p, s) on the central wave.

  • points, extent –

    The real-space grid: points per side over ±extent footprint widths.

  • eps, thickness_um, mu –

    The layers, as in :func:hyperbolic_optics.bcm.stack_smatrix (lab frame; the plane of incidence of the central wave is x-z).

  • angle_rad –

    Angle of incidence of the central wave in the top medium.

  • eps_top, eps_bottom –

    The half-spaces' permittivities.

  • waist_um –

    The beam waist w₀.

  • samples, span –

    The angular-spectrum grid (see :func:_beam_grid).

Returns:

  • –

    (x, y, field): coordinates in µm on the interface plane, [P] each,

  • –

    and the reflected field's p and s amplitudes (in the reflected central

  • –

    wave's basis) [P, P, 2].

Source code in hyperbolic_optics/beams.py
def reflected_beam(
    eps, thickness_um, wavelength_um, angle_rad, eps_top, eps_bottom, waist_um,
    polarization=(1.0, 0.0), mu=None, samples: int = 61, span: float = 5.0, points: int = 129,
    extent: float = 4.0,
):  # fmt: skip
    """The reflected beam's field on the reference plane, from its angular spectrum.

    Args:
        polarization: The incident beam's Jones vector ``(p, s)`` on the central wave.
        points, extent: The real-space grid: ``points`` per side over ``±extent``
            footprint widths.
        eps, thickness_um, mu: The layers, as in :func:`hyperbolic_optics.bcm.stack_smatrix`
            (lab frame; the plane of incidence of the central wave is x-z).
        angle_rad: Angle of incidence of the central wave in the top medium.
        eps_top, eps_bottom: The half-spaces' permittivities.
        waist_um: The beam waist ``w₀``.
        samples, span: The angular-spectrum grid (see :func:`_beam_grid`).

    Returns:
        ``(x, y, field)``: coordinates in µm on the interface plane, ``[P]`` each,
        and the reflected field's p and s amplitudes (in the reflected central
        wave's basis) ``[P, P, 2]``.
    """
    grid = _beam_grid(wavelength_um, waist_um, angle_rad, eps_top, samples, span)
    s = _component_smatrices(eps, thickness_um, wavelength_um, grid, eps_top, eps_bottom, mu)
    c_in = _jones(grid, eps_top, forward=True) @ jnp.asarray(polarization, dtype=complex)
    c_out = _jones(grid, eps_top, forward=False)
    local = s[:, :2, :2] @ c_in[..., None]
    amplitude = (jnp.conj(jnp.swapaxes(c_out, -1, -2)) @ local)[..., 0]  # [K, 2]
    spectrum = jnp.sqrt(grid.weight)[:, None] * amplitude
    k0 = 2 * jnp.pi / jnp.asarray(wavelength_um)
    footprint = jnp.asarray(waist_um) / jnp.cos(angle_rad)
    x = jnp.linspace(-extent, extent, points) * footprint
    y = jnp.linspace(-extent, extent, points) * jnp.asarray(waist_um)
    dqx, dqy = grid.qx - grid.q0[0], grid.qy - grid.q0[1]
    phase = jnp.exp(
        1j * k0 * (dqx[:, None, None] * x[None, :, None] + dqy[:, None, None] * y[None, None, :])
    )
    field = jnp.einsum("kxy,kc->xyc", phase, spectrum)
    return x, y, field

beam_shift

beam_shift(x, y, field)

Centroid (Δx, Δy) of |field|²: the Goos–Hänchen and Imbert–Fedorov shifts.

Source code in hyperbolic_optics/beams.py
def beam_shift(x, y, field) -> tuple[jax.Array, jax.Array]:
    """Centroid ``(Δx, Δy)`` of ``|field|²``: the Goos–Hänchen and Imbert–Fedorov shifts."""
    intensity = jnp.sum(jnp.abs(field) ** 2, -1)
    total = jnp.sum(intensity)
    return (
        jnp.sum(intensity * x[:, None]) / total,
        jnp.sum(intensity * y[None, :]) / total,
    )