janssen.models¶

Optical beam and material models for testing and validation.

Extended Summary¶

Models for generating datasets for testing and validation.

Routine Listings¶

azimuthally_polarized_beam(): Generate an azimuthally polarized beam.
bessel_beam(): Creates a Bessel beam with specified cone angle.
biological_cell(): Creates a biological cell model with nucleus.
calculate_usaf_group_range(): Calculates the viable USAF group range for given parameters.
circular_polarized_beam(): Generate a circularly polarized beam.
collimated_gaussian(): Creates a collimated Gaussian beam with flat phase.
converging_gaussian(): Creates a Gaussian beam converging to a focus.
create_bar_triplet(): Creates 3 parallel bars (horizontal or vertical).
create_element_pattern(): Creates a single element pattern (horizontal + vertical bars).
create_group_pattern(): Creates a group pattern with multiple elements.
diverging_gaussian(): Creates a Gaussian beam diverging from a virtual source.
gaussian_beam(): Creates a Gaussian beam from complex beam parameter q.
generalized_cylindrical_vector_beam(): Generate a generalized cylindrical vector beam.
generate_usaf_pattern(): Generates USAF 1951 resolution test pattern.
get_bar_width_pixels(): Calculate bar width in pixels for a given group and element.
gradient_index_material(): Creates a gradient-index (GRIN) material with radial profile.
hermite_gaussian(): Creates Hermite-Gaussian modes.
layered_material(): Creates alternating layers of materials.
laguerre_gaussian(): Creates Laguerre-Gaussian modes.
linear_polarized_beam(): Generate a linearly polarized beam with arbitrary angle.
plane_wave(): Creates a uniform plane wave with optional tilt.
propagate_beam(): Generates a beam at multiple z positions as a PropagatingWavefront.
radially_polarized_beam(): Generate a radially polarized beam.
sinusoidal_wave(): Creates a sinusoidal interference pattern.
spherical_inclusion(): Creates a material with spherical inclusion.
uniform_material(): Creates a uniform 3D material with constant refractive index.
x_polarized_beam(): Generate an x-polarized beam.
y_polarized_beam(): Generate a y-polarized beam.

Notes

All functions are JAX-compatible and support automatic differentiation.

Create a Bessel beam with specified cone angle.

Parameters:

wavelength (float) – Wavelength of light in meters.
cone_angle (float) – Cone half-angle in radians. Determines the transverse wave vector component: k_r = k * sin(cone_angle).
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Bessel beam OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

An ideal Bessel beam has a transverse profile given by the zeroth-order Bessel function:

\[E(r) = A J_0(k_r r)\]

where \(k_r = k \\sin(\\theta)\) is the transverse wave vector and \(\\theta\) is the cone half-angle.

Bessel beams are “non-diffracting” - their transverse profile remains constant upon propagation (in the ideal infinite case). In practice, physical Bessel beams have finite extent and eventually diffract.

The central lobe radius (first zero) is approximately \(r_0 \\approx 2.405 / k_r\).

Create a collimated Gaussian beam with flat phase.

Parameters:

wavelength (float) – Wavelength of light in meters.
waist (float) – Beam waist (1/e² intensity radius) in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]],) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
center (Tuple[ScalarFloat, ScalarFloat], optional) – Center position (x0, y0) in meters, by default (0.0, 0.0).
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Collimated Gaussian beam OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

A collimated Gaussian beam has a Gaussian intensity profile and flat (constant) phase across the beam:

\[\begin{split}E(x, y) = A \\exp\\left(-\\frac{(x-x_0)^2 + (y-y_0)^2}{w^2}\\right)\end{split}\]

where \(w\) is the beam waist (1/e² intensity radius).

This represents a beam at its waist position where the wavefront is planar. Upon propagation, the beam will expand due to diffraction.

Create a Gaussian beam converging to a focus.

Parameters:

wavelength (float) – Wavelength of light in meters.
waist (float) – Current beam waist (1/e² intensity radius) in meters. This is the waist at the current plane, not at the focus.
focus_distance (float) – Distance to the focus in meters (positive = focus downstream).
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[int, Tuple[int, int]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
center (Tuple[ScalarFloat, ScalarFloat], optional) – Center position (x0, y0) in meters, by default (0.0, 0.0).
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Converging Gaussian beam OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

A converging Gaussian beam has a Gaussian intensity profile with a spherical (quadratic) converging phase:

\[\begin{split}E(x, y) = A \\exp\\left(-\\frac{r^2}{w^2}\\right) \\exp\\left(-i \\frac{k r^2}{2 f}\\right)\end{split}\]

where \(f\) is the focus distance and the negative sign indicates convergence (wavefront curving inward toward the optical axis).

The radius of curvature R equals the focus distance f for a beam that will come to a focus at distance f downstream.

Upon propagation, this beam will decrease in size until reaching the focus, then expand as a diverging beam.

Create a Gaussian beam diverging from a virtual source.

Parameters:

wavelength (float) – Wavelength of light in meters.
waist (float) – Current beam waist (1/e² intensity radius) in meters. This is the waist at the current plane.
source_distance (float) – Distance from the virtual source point in meters (positive = source was upstream).
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
center (Tuple[ScalarFloat, ScalarFloat], optional) – Center position (x0, y0) in meters, by default (0.0, 0.0).
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Diverging Gaussian beam OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

A diverging Gaussian beam has a Gaussian intensity profile with a spherical (quadratic) diverging phase:

\[\begin{split}E(x, y) = A \\exp\\left(-\\frac{r^2}{w^2}\\right) \\exp\\left(+i \\frac{k r^2}{2 R}\\right)\end{split}\]

where \(R\) is the radius of curvature (positive for diverging, equal to the source distance).

This represents a beam that originated from a point source at distance source_distance upstream. Upon propagation, this beam will continue to expand.

Create a Gaussian beam at arbitrary position from waist.

This is the most general Gaussian beam generator, using the full Gaussian beam propagation formulas including Gouy phase.

Parameters:

wavelength (float) – Wavelength of light in meters.
waist_0 (float) – Beam waist at the waist position (minimum spot size) in meters.
z_from_waist (float) – Distance from the beam waist in meters. Positive = downstream from waist (diverging). Negative = upstream from waist (converging toward waist).
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
center (Tuple[ScalarFloat, ScalarFloat], optional) – Center position (x0, y0) in meters, by default (0.0, 0.0).
amplitude (float, optional) – Peak amplitude at beam waist, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.
include_gouy_phase (bool, optional) – Whether to include the Gouy phase shift, by default True.

Returns:

wavefront – Gaussian beam OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

The complete Gaussian beam field is:

\[\begin{split}E(r, z) = A \\frac{w_0}{w(z)} \\exp\\left(-\\frac{r^2}{w(z)^2}\\right) \\exp\\left(-ikz - i\\frac{kr^2}{2R(z)} + i\\zeta(z)\\right)\end{split}\]

where:

\(w(z) = w_0 \\sqrt{1 + (z/z_R)^2}\) is the beam radius
\(R(z) = z (1 + (z_R/z)^2)\) is the radius of curvature
\(\\zeta(z) = \\arctan(z/z_R)\) is the Gouy phase
\(z_R = \\pi w_0^2 / \\lambda\) is the Rayleigh range

At the waist (z=0), the beam has minimum size and flat phase. The Gouy phase represents an additional phase shift accumulated through the focus.

Create a Hermite-Gaussian mode at the beam waist.

Parameters:

wavelength (float) – Wavelength of light in meters.
waist (float) – Beam waist (1/e² intensity radius of fundamental mode) in meters.
n (int) – Mode index in x direction (number of nodes along x), n >= 0.
m (int) – Mode index in y direction (number of nodes along y), m >= 0.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
center (Tuple[ScalarFloat, ScalarFloat], optional) – Center position (x0, y0) in meters, by default (0.0, 0.0).
amplitude (float, optional) – Peak amplitude normalization, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Hermite-Gaussian mode OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

The Hermite-Gaussian modes at the waist are:

\[\begin{split}E_{n,m}(x, y) = A H_n\\left(\\frac{x\\sqrt{2}}{w}\\right) H_m\\left(\\frac{y\\sqrt{2}}{w}\\right) \\exp\\left(-\\frac{x^2 + y^2}{w^2}\\right)\end{split}\]

where \(H_n\) is the physicist’s Hermite polynomial.

Special cases:

(n=0, m=0): Fundamental Gaussian mode (TEM00)
(n=1, m=0): TEM10 mode with one node along x
(n=0, m=1): TEM01 mode with one node along y

Create a Laguerre-Gaussian mode at the beam waist.

Parameters:

wavelength (float) – Wavelength of light in meters.
waist (float) – Beam waist (1/e² intensity radius of fundamental mode) in meters.
p (int) – Radial mode index (number of radial nodes), p >= 0.
l (int) – Azimuthal mode index (topological charge for vortex beams). Can be positive or negative.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
center (Tuple[ScalarFloat, ScalarFloat], optional) – Center position (x0, y0) in meters, by default (0.0, 0.0).
amplitude (float, optional) – Peak amplitude normalization, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Laguerre-Gaussian mode OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

The Laguerre-Gaussian modes at the waist are:

\[\begin{split}E_{p,l}(r, \\phi) = A \\left(\\frac{r\\sqrt{2}}{w}\\right)^{|l|} L_p^{|l|}\\left(\\frac{2r^2}{w^2}\\right) \\exp\\left(-\\frac{r^2}{w^2}\\right) \\exp(i l \\phi)\end{split}\]

where \(L_p^{|l|}\) is the generalized Laguerre polynomial.

Special cases:

(p=0, l=0): Fundamental Gaussian mode
(p=0, l≠0): Optical vortex beams with topological charge l
(p>0, l=0): Radial modes with p rings

Create a uniform plane wave with optional tilt.

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
amplitude (float, optional) – Amplitude of the plane wave, by default 1.0.
tilt_x (float, optional) – Tilt angle along x-axis in radians (small angle), by default 0.0.
tilt_y (float, optional) – Tilt angle along y-axis in radians (small angle), by default 0.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Plane wave OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

A plane wave has uniform amplitude and linear phase (flat if no tilt). The tilt angles introduce a linear phase ramp corresponding to propagation at an angle to the optical axis:

\[\begin{split}E(x, y) = A \\exp(i k (x \\sin\\theta_x + y \\sin\\theta_y))\end{split}\]

For small angles, \(\\sin\\theta \\approx \\theta\).

janssen.models.propagate_beam(beam_type: str, z_positions: Float[Array, 'zz'], wavelength: float | Float[Array, ''], dx: float | Float[Array, ''], grid_size: Int[Array, '2'] | tuple[int | Int[Array, ''], int | Int[Array, '']], waist: float | Float[Array, ''] = 0.001, amplitude: float | Float[Array, ''] = 1.0, center: tuple[float | Float[Array, ''], float | Float[Array, '']] = (0.0, 0.0), tilt_x: float | Float[Array, ''] = 0.0, tilt_y: float | Float[Array, ''] = 0.0, focus_distance: float | Float[Array, ''] = 1.0, source_distance: float | Float[Array, ''] = 1.0, waist_0: float | Float[Array, ''] = 0.001, include_gouy_phase: bool = True, cone_angle: float | Float[Array, ''] = 0.01, period: float | Float[Array, ''] = 0.0001, direction: float | Float[Array, ''] = 0.0, p: int | Int[Array, ''] = 0, l: int | Int[Array, ''] = 0, n: int | Int[Array, ''] = 0, m: int | Int[Array, ''] = 0) → PropagatingWavefront[source]¶

Generate a beam at multiple z positions as a PropagatingWavefront.

Parameters:

beam_type (str) – Type of beam to generate. One of: - “plane_wave”: Uniform plane wave with optional tilt - “sinusoidal_wave”: Sinusoidal interference pattern - “collimated_gaussian”: Collimated Gaussian beam with flat phase - “converging_gaussian”: Gaussian beam converging to a focus - “diverging_gaussian”: Gaussian beam diverging from a source - “gaussian_beam”: General Gaussian beam at arbitrary z from waist - “bessel_beam”: Bessel beam with specified cone angle - “laguerre_gaussian”: Laguerre-Gaussian modes - “hermite_gaussian”: Hermite-Gaussian modes
z_positions (Float[Array, " zz"]) – Array of z positions at which to evaluate the beam (meters). For “gaussian_beam”, these are distances from the waist. For other beam types, these set the z_position attribute.
wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width).
waist (float) – Beam waist (1/e² intensity radius) in meters. Used by collimated_gaussian, converging_gaussian, diverging_gaussian, laguerre_gaussian, hermite_gaussian. Default is 1e-3.
amplitude (float) – Peak amplitude, by default 1.0.
center (Tuple[ScalarFloat, ScalarFloat]) – Center position (x0, y0) in meters, by default (0.0, 0.0).
tilt_x (float) – Tilt angle along x-axis in radians for plane_wave, by default 0.0.
tilt_y (float) – Tilt angle along y-axis in radians for plane_wave, by default 0.0.
focus_distance (float) – Distance to focus for converging_gaussian (meters). Default is 1.0.
source_distance (float) – Distance from source for diverging_gaussian (meters). Default is 1.0.
waist_0 (float) – Beam waist at the waist position for gaussian_beam (meters). Default is 1e-3.
include_gouy_phase (bool) – Whether to include Gouy phase for gaussian_beam, by default True.
cone_angle (float) – Cone half-angle in radians for bessel_beam. Default is 0.01.
period (float) – Spatial period for sinusoidal_wave (meters). Default is 1e-4.
direction (float) – Direction angle for sinusoidal_wave (radians). Default is 0.0.
p (int) – Radial mode index for laguerre_gaussian, by default 0.
l (int) – Azimuthal mode index for laguerre_gaussian, by default 0.
n (int) – Mode index in x direction for hermite_gaussian, by default 0.
m (int) – Mode index in y direction for hermite_gaussian, by default 0.

Returns:

propagating_wavefront – A PropagatingWavefront containing the beam at all specified z positions.

Return type:

PropagatingWavefront

Raises:

ValueError – If beam_type is not recognized.

Notes

This function uses jax.vmap to efficiently generate the beam at all z positions in parallel. The resulting PropagatingWavefront can be used to visualize beam evolution or as input to propagation algorithms.

For “gaussian_beam”, the z_positions represent distances from the beam waist, allowing visualization of beam evolution through focus. For other beam types, z_positions sets the z_position attribute but the field profile remains constant (as these are evaluated at a single plane).

Create a sinusoidal interference pattern.

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid. If int, creates square grid. If tuple, specifies (height, width).
period (float) – Spatial period of the sinusoidal pattern in meters.
direction (float) – Direction angle of the sinusoidal pattern in radians. 0 = horizontal stripes, pi/2 = vertical stripes. By default 0.0.
amplitude (float) – Peak amplitude of the wave, by default 1.0.
z_position (float) – Initial z position of the wavefront in meters, by default 0.0.

Returns:

wavefront – Sinusoidal wave OpticalWavefront PyTree.

Return type:

OpticalWavefront

Notes

A sinusoidal wave has an intensity profile that varies as:

\[E(x, y) = A \cos\left(\frac{2\pi}{T}(x \cos\theta + y \sin\theta) \right)\]

where \(T\) is the spatial period and \(\theta\) is the direction angle.

This pattern represents the interference of two plane waves and is useful for testing optical systems, creating gratings, and studying diffraction phenomena.

Create biological cell model with nucleus.

Parameters:

shape (Tuple[int, int, int]) – Shape of material (height, width, num_slices)
cell_radius (float) – Outer radius of cell in meters
nucleus_radius (float) – Radius of nucleus in meters
n_cytoplasm (complex) – Complex refractive index of cytoplasm
n_nucleus (complex) – Complex refractive index of nucleus
n_medium (complex) – Complex refractive index of surrounding medium
dx (float) – Pixel spacing in x-y plane in meters
tz (float) – Slice spacing in z direction in meters
center (Tuple[ScalarFloat, float, ScalarFloat], optional) – Center position (y, x, z) in meters. If None, centered in volume.

Returns:

material – Biological cell material PyTree

Return type:

SlicedMaterialFunction

Notes

Algorithm: - Create coordinate grids for x, y, z - Compute radial distance from center - Use nested jnp.where to assign refractive index:

If r < nucleus_radius: n_nucleus

Elif r < cell_radius: n_cytoplasm

Else: n_medium

Convert to SlicedMaterialFunction PyTree

Implementation Details: - Nested jnp.where for three regions - Inner to outer: nucleus, cytoplasm, medium - All operations use JAX arrays - Typical values:

n_medium: 1.337 (water)

n_cytoplasm: 1.360-1.380

n_nucleus: 1.380-1.400

κ increases with density

Examples

Create simple cell model:

>>> cell = biological_cell(
...     shape=(128, 128, 20),
...     cell_radius=25e-6,
...     nucleus_radius=10e-6,
...     n_cytoplasm=1.370 + 0.002j,
...     n_nucleus=1.390 + 0.005j,
...     n_medium=1.337 + 0.0j,
...     dx=1e-6,
...     tz=5e-6
... )

janssen.models.gradient_index_material(shape: tuple[int, int, int], n_center: float | Float[Array, ''], gradient_constant: float | Float[Array, ''], dx: float | Float[Array, ''], tz: float | Float[Array, ''], n_min: float | Float[Array, ''] | None = None) → SlicedMaterialFunction[source]¶

Create gradient-index (GRIN) material with radial profile.

Parameters:

shape (Tuple[int, int, int]) – Shape of material (height, width, num_slices)
n_center (float) – Refractive index at center
gradient_constant (float) – GRIN constant A in n(r) = n0(1 - Ar²/2), units: m^-2
dx (float) – Pixel spacing in x-y plane in meters
tz (float) – Slice spacing in z direction in meters
n_min (float, optional) – Minimum refractive index (for clipping), default is 1.0

Returns:

material – GRIN material PyTree

Return type:

SlicedMaterialFunction

Notes

Algorithm: - Create coordinate grids for x, y - Compute radial distance in x-y plane: r_xy = sqrt(x² + y²) - Apply GRIN formula: n(r) = n_center × (1 - A × r_xy² / 2) - Clip to ensure n >= n_min - Broadcast to all z slices (z-invariant) - Convert to SlicedMaterialFunction PyTree

Implementation Details: - GRIN profile: n(r) = n₀(1 - Ar²/2) - Only depends on x-y position (constant along z) - Uses jnp.clip to enforce physical bounds - Typical A values: 1e8 to 1e12 m^-2 - All operations are differentiable - Real refractive index only (no absorption)

Examples

Create GRIN lens:

>>> grin = gradient_index_material(
...     shape=(128, 128, 20),
...     n_center=1.5,
...     gradient_constant=1e10,
...     dx=1e-6,
...     tz=5e-6
... )

janssen.models.layered_material(shape: tuple[int, int, int], n_layers: tuple[complex | Complex[Array, ''], complex | Complex[Array, '']], dx: float | Float[Array, ''], tz: float | Float[Array, ''], slices_per_layer: int | Int[Array, ''] = 1) → SlicedMaterialFunction[source]¶

Create alternating layers of materials.

Parameters:

shape (Tuple[int, int, int]) – Shape of material (height, width, num_slices)
n_layers (Tuple[ScalarComplex, ScalarComplex]) – Refractive indices of two alternating layers (n1, n2)
dx (float) – Pixel spacing in x-y plane in meters
tz (float) – Slice spacing in z direction in meters
slices_per_layer (int, optional) – Number of slices per layer, by default 1

Returns:

material – Layered material PyTree

Return type:

SlicedMaterialFunction

Notes

Algorithm: - Create array to hold material - Use vmap over z slices to assign refractive index - Layer index = z_idx // slices_per_layer - Use modulo to alternate between n1 and n2 - Convert to SlicedMaterialFunction PyTree

Implementation Details: - Uses jax.lax.cond for layer selection (JAX-compatible) - Layer pattern: n1, n2, n1, n2, … - Each layer is slices_per_layer thick - All operations are differentiable - Uses vmap for efficient parallel processing

Examples

Create optical coating with alternating layers:

>>> coating = layered_material(
...     shape=(128, 128, 20),
...     n_layers=(1.38 + 0.0j, 2.35 + 0.0j),
...     dx=1e-6,
...     tz=50e-9,
...     slices_per_layer=2
... )

Create material with spherical inclusion.

Parameters:

shape (Tuple[int, int, int]) – Shape of material (height, width, num_slices)
radius (float) – Radius of sphere in meters
n_sphere (complex) – Complex refractive index of sphere
n_background (complex) – Complex refractive index of background
dx (float) – Pixel spacing in x-y plane in meters
tz (float) – Slice spacing in z direction in meters
center (Tuple[ScalarFloat, float, ScalarFloat], optional) – Center position (y, x, z) in meters. If None, centered in volume.

Returns:

material – Material with spherical inclusion

Return type:

SlicedMaterialFunction

Notes

Algorithm: - Create coordinate grids for x, y, z - Compute radial distance from center - Use jnp.where to assign refractive index based on distance - Convert to SlicedMaterialFunction PyTree

Implementation Details: - Uses meshgrid with indexing=’ij’ for consistent ordering - If center is None, defaults to volume center - Radial distance: r = sqrt((x-x0)^2 + (y-y0)^2 + (z-z0)^2) - Refractive index: n = n_sphere if r < radius else n_background - All operations use JAX arrays for automatic differentiation

Examples

Create material with absorbing sphere:

>>> sphere = spherical_inclusion(
...     shape=(128, 128, 20),
...     radius=20e-6,
...     n_sphere=1.5 + 0.01j,
...     n_background=1.0 + 0.0j,
...     dx=1e-6,
...     tz=5e-6
... )

janssen.models.uniform_material(shape: tuple[int, int, int], n: complex | Complex[Array, ''], dx: float | Float[Array, ''], tz: float | Float[Array, '']) → SlicedMaterialFunction[source]¶

Create uniform 3D material with constant refractive index.

Parameters:

shape (Tuple[int, int, int]) – Shape of material (height, width, num_slices)
n (complex) – Complex refractive index (n_real + i*kappa)
dx (float) – Pixel spacing in x-y plane in meters
tz (float) – Slice spacing in z direction in meters

Returns:

material – Uniform material PyTree

Return type:

SlicedMaterialFunction

Notes

Algorithm: - Create array filled with constant refractive index - Convert to SlicedMaterialFunction PyTree - Return material

Implementation Details: - Uses jnp.ones for efficient array creation - Multiplies by complex scalar to set refractive index - All operations are JAX-compatible and differentiable

Examples

Create uniform glass material:

>>> glass = uniform_material(
...     shape=(128, 128, 20),
...     n=1.52 + 0.0j,
...     dx=1e-6,
...     tz=5e-6
... )

janssen.models.azimuthally_polarized_beam(wavelength: float | Float[Array, ''], dx: float | Float[Array, ''], grid_size: Int[Array, '2'] | tuple[int | Int[Array, ''], int | Int[Array, '']], beam_radius: float | Float[Array, ''] = None, amplitude: float | Float[Array, ''] = 1.0, z_position: float | Float[Array, ''] = 0.0, apodization: str = 'gaussian') → OpticalWavefront[source]¶

Generate an azimuthally polarized beam.

Creates a cylindrical vector beam where the electric field points in the azimuthal direction (tangent to circles centered on the axis). The polarization direction is:

\[\hat{e}_\phi = -\sin(\phi)\hat{x} + \cos(\phi)\hat{y}\]

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width) or (ny, nx).
beam_radius (float, optional) – Characteristic beam radius in meters. If None, defaults to 1/4 of the grid extent.
amplitude (float, optional) – Peak amplitude, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.
apodization (str, optional) – Amplitude envelope type: “gaussian”, “uniform”, or “bessel”.

Returns:

Azimuthally polarized beam as OpticalWavefront PyTree. Field shape is (H, W, 2) with [Ex, Ey] components.

Return type:

OpticalWavefront

Notes

Azimuthally polarized beams have complementary properties to radial polarization:

Donut focus: When focused by a high-NA lens, azimuthally polarized beams produce NO longitudinal (Ez) field component. The focal spot has a dark center (“donut” shape).
Singular at center: Like radial polarization, the polarization is undefined at r=0, and the amplitude naturally goes to zero.
Orthogonal to radial: At every point, the azimuthal polarization is perpendicular to the radial polarization.

The field components are:

\[ \begin{align}\begin{aligned}E_x(r, \phi) = -A(r) \sin(\phi)\\E_y(r, \phi) = A(r) \cos(\phi)\end{aligned}\end{align} \]

References

janssen.models.circular_polarized_beam(wavelength: float | Float[Array, ''], dx: float | Float[Array, ''], grid_size: Int[Array, '2'] | tuple[int | Int[Array, ''], int | Int[Array, '']], handedness: str = 'right', beam_radius: float | Float[Array, ''] = None, amplitude: float | Float[Array, ''] = 1.0, z_position: float | Float[Array, ''] = 0.0, apodization: str = 'gaussian') → OpticalWavefront[source]¶

Generate a circularly polarized Gaussian beam.

Creates a beam with uniform circular polarization (left or right handed) across the aperture.

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width).
handedness (str, optional) – Polarization handedness: “right” (default) or “left”. - “right” : Right-handed (clockwise when viewed from receiver) - “left” : Left-handed (counter-clockwise when viewed from receiver)
beam_radius (float, optional) – Beam waist (1/e² intensity radius) in meters.
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position in meters, by default 0.0.
apodization (str, optional) – Amplitude envelope type.

Returns:

Circularly polarized beam with field shape (H, W, 2).

Return type:

OpticalWavefront

Notes

The Jones vectors for circular polarization are:

Right-handed (RCP):

\[\begin{split}\vec{E}_{RCP} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ -i \end{pmatrix}\end{split}\]

Left-handed (LCP):

\[\begin{split}\vec{E}_{LCP} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 \\ i \end{pmatrix}\end{split}\]

Note: Sign conventions vary in the literature. We use the physics convention where RCP has E_y = -i*E_x (clockwise rotation for an observer looking into the beam).

Generate a generalized cylindrical vector beam of arbitrary order.

Creates a cylindrical vector beam with polarization pattern determined by the topological order m. Standard radial (m=1, φ₀=0) and azimuthal (m=1, φ₀=π/2) polarizations are special cases.

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width).
order (int, optional) – Topological order m of the polarization pattern, by default 1. - m=1 : Standard radial/azimuthal (1 polarization singularity) - m=2 : Higher-order with 2 singularities - etc.
phase_offset (float, optional) – Phase offset φ₀ in radians, by default 0.0. - φ₀=0 : Radial-like - φ₀=π/2 : Azimuthal-like
beam_radius (float, optional) – Beam waist in meters.
amplitude (float, optional) – Peak amplitude, by default 1.0.
z_position (float, optional) – Initial z position in meters, by default 0.0.

Returns:

Generalized cylindrical vector beam with field shape (H, W, 2).

Return type:

OpticalWavefront

Notes

The generalized cylindrical vector beam has polarization:

\[ \begin{align}\begin{aligned}E_x = A(r) \cos(m\phi + \phi_0)\\E_y = A(r) \sin(m\phi + \phi_0)\end{aligned}\end{align} \]

For m=1: - φ₀=0 gives radial polarization - φ₀=π/2 gives azimuthal polarization

Higher-order beams (m>1) have multiple polarization singularities and create complex focal field distributions.

janssen.models.linear_polarized_beam(wavelength: float | Float[Array, ''], dx: float | Float[Array, ''], grid_size: Int[Array, '2'] | tuple[int | Int[Array, ''], int | Int[Array, '']], polarization_angle: float | Float[Array, ''] = 0.0, beam_radius: float | Float[Array, ''] = None, amplitude: float | Float[Array, ''] = 1.0, z_position: float | Float[Array, ''] = 0.0, apodization: str = 'gaussian') → OpticalWavefront[source]¶

Generate a linearly polarized Gaussian beam.

Creates a beam with uniform linear polarization across the aperture. The polarization direction is specified by an angle from the x-axis.

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width) or (ny, nx).
polarization_angle (float, optional) – Angle of polarization direction measured CCW from +x axis, in radians. Default is 0.0 (x-polarized). - 0 : x-polarized - π/2 : y-polarized - π/4 : 45° polarized
beam_radius (float, optional) – Beam waist (1/e² intensity radius) in meters. If None, defaults to 1/4 of grid extent.
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.
apodization (str, optional) – Amplitude envelope: “gaussian” (default) or “uniform”.

Returns:

Linearly polarized beam as OpticalWavefront PyTree. Field shape is (H, W, 2) with [Ex, Ey] components.

Return type:

OpticalWavefront

Notes

For a linearly polarized beam, the Jones vector is constant across the aperture:

\[\begin{split}\vec{E} = A(r) \begin{pmatrix} \cos(\theta) \\ \sin(\theta) \end{pmatrix}\end{split}\]

where θ is the polarization angle.

When focused by a high-NA lens: - The focal spot is slightly elliptical (elongated along the

polarization direction)

A weak Ez component appears due to depolarization effects
The effect increases with NA

janssen.models.radially_polarized_beam(wavelength: float | Float[Array, ''], dx: float | Float[Array, ''], grid_size: Int[Array, '2'] | tuple[int | Int[Array, ''], int | Int[Array, '']], beam_radius: float | Float[Array, ''] = None, amplitude: float | Float[Array, ''] = 1.0, z_position: float | Float[Array, ''] = 0.0, apodization: str = 'gaussian') → OpticalWavefront[source]¶

Generate a radially polarized beam.

Creates a cylindrical vector beam where the electric field points radially outward from the optical axis at every point. The polarization direction is:

\[\hat{e}_r = \cos(\phi)\hat{x} + \sin(\phi)\hat{y}\]

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width) or (ny, nx).
beam_radius (float, optional) – Characteristic beam radius in meters. If None, defaults to 1/4 of the grid extent. For Gaussian apodization, this is the 1/e² intensity radius.
amplitude (float, optional) – Peak amplitude at beam edge, by default 1.0.
z_position (float, optional) – Initial z position of the wavefront in meters, by default 0.0.
apodization (str, optional) – Amplitude envelope type. Options are: - “gaussian” : Gaussian envelope exp(-r²/w²) (default) - “uniform” : Uniform amplitude within beam_radius - “bessel” : J₁(r) profile (natural for radial pol.)

Returns:

Radially polarized beam as OpticalWavefront PyTree. Field shape is (H, W, 2) with [Ex, Ey] components.

Return type:

OpticalWavefront

Notes

Radially polarized beams have several unique properties:

Tight focusing: When focused by a high-NA lens, radially polarized beams produce a strong longitudinal (Ez) field component at the focus, resulting in a tighter focal spot than linearly polarized beams.
Singular at center: The polarization is undefined at r=0 (the optical axis), creating a phase singularity. The amplitude naturally goes to zero at the center.
Zero orbital angular momentum: Unlike Laguerre-Gaussian vortex beams, radially polarized beams carry no orbital angular momentum.

The field components are:

\[ \begin{align}\begin{aligned}E_x(r, \phi) = A(r) \cos(\phi)\\E_y(r, \phi) = A(r) \sin(\phi)\end{aligned}\end{align} \]

where A(r) is the amplitude envelope (Gaussian, uniform, or Bessel).

References

Generate an x-polarized Gaussian beam.

Convenience wrapper for linear_polarized_beam with angle = 0.

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width).
beam_radius (float, optional) – Beam waist (1/e² intensity radius) in meters.
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position in meters, by default 0.0.
apodization (str, optional) – Amplitude envelope type.

Returns:

X-polarized beam with field shape (H, W, 2).

Return type:

OpticalWavefront

Generate a y-polarized Gaussian beam.

Convenience wrapper for linear_polarized_beam with angle = π/2.

Parameters:

wavelength (float) – Wavelength of light in meters.
dx (float) – Spatial sampling interval (pixel size) in meters.
grid_size (Union[Int[Array, " 2"], Tuple[ScalarInteger, ScalarInteger]]) – Size of the computational grid as (height, width).
beam_radius (float, optional) – Beam waist (1/e² intensity radius) in meters.
amplitude (float, optional) – Peak amplitude at beam center, by default 1.0.
z_position (float, optional) – Initial z position in meters, by default 0.0.
apodization (str, optional) – Amplitude envelope type.

Returns:

Y-polarized beam with field shape (H, W, 2).

Return type:

OpticalWavefront

janssen.models.calculate_usaf_group_range(image_size: int, pixel_size: float, min_bar_pixels: int = 2, grid_fill_fraction: float = 0.95) → dict[source]¶

Calculate the viable USAF group range for given parameters.

Parameters:

image_size (int) – Image size in pixels (square)
pixel_size (float) – Pixel size in meters
min_bar_pixels (int, optional) – Minimum bar width in pixels for visibility, by default 2
grid_fill_fraction (float, optional) – Scale factor for fitting largest group, by default 0.95

Returns:

result – Dictionary containing: - max_group: finest group where bars are still >= min_bar_pixels - min_group: coarsest group that fits in image - recommended_range: suggested range() for generate_usaf_pattern - num_groups: how many groups in recommended range - pixels_per_mm: calculated pixel density - group_info: dict with bar width info for each group

Return type:

dict

Notes

This function maximizes the number of groups that can fit in the image using variable-density row packing (smaller groups pack more per row).

Examples

>>> result = calculate_usaf_group_range(
...     image_size=8192,
...     pixel_size=0.5e-6,
... )
>>> min_g, max_g = result['min_group'], result['max_group']
>>> print(f"Use: groups=range({min_g}, {max_g + 1})")

janssen.models.create_bar_triplet(width: int, length: int, horizontal: bool = True) → Float[Array, '...'][source]¶

Create 3 parallel bars (horizontal or vertical).

Parameters:

width (int) – Width of each bar in pixels (minimum 1)
length (int) – Length of each bar in pixels (minimum 1)
horizontal (bool, optional) – Whether to create horizontal bars, by default True

Returns:

pattern – The bar triplet pattern. Shape depends on orientation: - Horizontal: (5*width, length) - Vertical: (length, 5*width)

Return type:

Float[Array, "..."]

Notes

Creates three bars following USAF specification where bar spacing equals bar width. Total extent is 5 × bar_width (3 bars + 2 spaces).

Pattern structure (for horizontal): - Bar 1: rows [0, width) - Space: rows [width, 2*width) - Bar 2: rows [2*width, 3*width) - Space: rows [3*width, 4*width) - Bar 3: rows [4*width, 5*width)

Both horizontal and vertical patterns are computed, then selected based on the horizontal flag. This is JAX-safe since the flag is a static Python bool that doesn’t change during tracing.

janssen.models.create_element_pattern(bar_width_px: int, gap_factor: float = 0.5) → Float[Array, '...'][source]¶

Create a single USAF element (horizontal + vertical bar triplets).

Parameters:

bar_width_px (int) – Width of each bar in pixels (minimum 1)
gap_factor (float, optional) – Gap between triplets as fraction of bar_width, by default 0.5

Returns:

element – The complete element pattern with both triplets

Return type:

Float[Array, "..."]

Notes

Each USAF element consists of: - 3 horizontal bars (triplet) on the left - 3 vertical bars (triplet) on the right Bar length is 5× the bar width per USAF specification.

The element is composed as: [horizontal triplet] [gap] [vertical triplet]

Triplets are centered vertically within the element canvas.

janssen.models.create_group_pattern(group: int, pixels_per_mm: float) → tuple[Float[Array, '...'], int][source]¶

Create a complete group with 6 elements in 2×3 layout.

Parameters:

group (int) – Group number
pixels_per_mm (float) – Pixel density in pixels per millimeter

Return type:

tuple['], int]

Returns:

group_pattern (Float[Array, "..."]) – The complete group pattern
max_dimension (int) – Maximum dimension of the group

Notes

Elements are arranged in 2 columns × 3 rows: - Column 1: Elements 1, 2, 3 (top to bottom) - Column 2: Elements 4, 5, 6 (top to bottom)

Elements within a group progressively decrease in size following the 2^((element-1)/6) scaling.

The loop over 6 elements is unrolled at trace time since the element count is fixed.

janssen.models.generate_usaf_pattern(image_size: int = 1024, groups: range | None = None, pixel_size: float | Float[Array, ''] = 1e-06, background: float = 0.0, foreground: float = 1.0, max_phase: float = 0.0, auto: bool = False, min_bar_pixels: int = 2) → SampleFunction[source]¶

Generate USAF 1951 resolution test pattern.

Parameters:

image_size (int, optional) – Size of the output image (square), by default 1024
groups (range, optional) – Range of groups to include, by default range(-2, 8). Ignored if auto=True.
pixel_size (float, optional) – Physical size of each pixel in meters, by default 1.0e-6 (1 µm)
background (float, optional) – Background value, by default 0.0 (black)
foreground (float, optional) – Foreground (bar) value, by default 1.0 (white)
max_phase (float, optional) – Maximum phase shift in radians applied to the bars, by default 0.0. The phase pattern follows the same structure as the amplitude, scaling from 0 (at background) to max_phase (at foreground).
auto (bool, optional) – If True, automatically calculate the optimal group range to fill the image based on image_size and pixel_size. Overrides the groups parameter. By default False.
min_bar_pixels (int, optional) – Minimum bar width in pixels for visibility when auto=True, by default 2. Ignored if auto=False.

Returns:

pattern – SampleFunction PyTree containing the USAF test pattern as a complex array with both amplitude and phase information.

Return type:

SampleFunction

Notes

The USAF 1951 test pattern consists of groups arranged in a grid. Each group contains 6 elements of progressively higher resolution.

Resolution formula per MIL-STD-150A:: Resolution = 2^(group + (element-1)/6) line pairs per mm

Standard groups range from -2 (coarsest) to 7 (finest).

Each element consists of: - 3 horizontal bars (bar triplet) - 3 vertical bars (bar triplet) with bar length = 5 × bar width.

The output is a complex field: amplitude * exp(i * phase), where the phase follows the same spatial pattern as the amplitude.

The loop over groups is unrolled at Python trace time since groups_list is known before tracing. Python-level conditionals for bounds checking, scaling, and phase normalization are evaluated at trace time since all controlling values are Python scalars.

A global scale factor is computed from the largest (coarsest) group to ensure all groups fit within their grid cells while preserving the correct relative size ratios between groups.

Examples

>>> from janssen.models import generate_usaf_pattern
>>> pattern = generate_usaf_pattern(image_size=1024, pixel_size=1e-6)
>>> pattern.sample.shape
(1024, 1024)

>>> # Auto mode: fill the image optimally
>>> pattern = generate_usaf_pattern(
...     image_size=8192, pixel_size=0.5e-6, auto=True)

>>> # Camera with 6.5 µm pixels
>>> pattern = generate_usaf_pattern(
...     pixel_size=6.5e-6)

>>> # White background with black bars (typical)
>>> pattern = generate_usaf_pattern(background=1.0, foreground=0.0)

>>> # Phase object with π phase shift on bars
>>> pattern = generate_usaf_pattern(max_phase=jnp.pi)

>>> # Specific group range
>>> pattern = generate_usaf_pattern(groups=range(0, 5))

janssen.models.get_bar_width_pixels(group: int, element: int, pixels_per_mm: float) → int[source]¶

Calculate bar width in pixels for a given group and element.

Parameters:

group (int) – Group number (typically -2 to 7)
element (int) – Element number (1 to 6)
pixels_per_mm (float) – Pixel density in pixels per millimeter

Returns:

bar_width – Bar width in pixels (minimum 1)

Return type:

int

Notes

Resolution formula per MIL-STD-150A:: R = 2^(group + (element-1)/6) line pairs per mm

One line pair = one bar + one space = 2 × bar_width Therefore: bar_width_mm = 1 / (2 × R)