EM Fields’ ABC: Fundamentals of Electromagnetic Interactions#

In radio propagation and ray tracing, high-frequency electromagnetic fields interact with the physical environment through several distinct physical mechanisms:

  1. Specular Reflection: Plane surfaces reflect rays according to Snell’s law, with reflection coefficients given by Fresnel formulas.

  2. Edge Diffraction: Sharp edges scatter energy into geometric shadow zones governed by the Uniform Theory of Diffraction (UTD).

  3. Transmission: Rays passing through finite-thickness dielectric slabs experience attenuation and phase delay.

  4. Diffuse Scattering: Surface roughness scatters power into a broad angular lobe, following a scattering pattern (e.g., Lambertian, directive, or backscattering).

This tutorial demonstrates each interaction mathematically and visually using DiffeRT.

import differt.plotting as dplt
import jax.numpy as jnp
import matplotlib.pyplot as plt
from differt.em import (
    BackscatteringPattern,
    DirectivePattern,
    InteractionType,
    Material,
    compute_received_fields,
    compute_received_power,
    materials,
)
from differt.geometry import Mesh, Scene, TracedPaths

dplt.set_backend("plotly")
Warp CUDA warning: Could not find or load the NVIDIA CUDA driver. GPU execution will not be available.

1. Specular Reflection: Two-Ray Ground Interference#

When both a direct Line-of-Sight (LoS) path and a ground-reflected path reach the receiver, their complex electric fields superimpose as

\[E_{\text{total}} = E_{\text{LoS}} + E_{\text{refl}}.\]

Because the propagation path lengths differ, the phase difference \(\Delta \phi = \frac{2\pi}{\lambda} \Delta d\) cycles periodically as the receiver moves, creating characteristic constructive and destructive interference fringes (two-ray multipath fading).

Hide code cell source

frequency = 1e9  # 1 GHz carrier frequency

# Large ground plate at z=0 made of concrete
ground = (
    Mesh(
        vertices=jnp.array([
            [-50.0, -50.0, 0.0],
            [200.0, -50.0, 0.0],
            [200.0, 50.0, 0.0],
            [-50.0, 50.0, 0.0],
        ]),
        triangles=jnp.array([[0, 1, 2], [0, 2, 3]]),
    )
    .set_face_materials(0)
    .set_materials("itu_concrete")
)

ht, hr = 10.0, 1.5  # Transmitter height = 10m, Receiver height = 1.5m
d = jnp.linspace(5.0, 100.0, 500)  # Receiver distance along ground
tx = jnp.array([0.0, 0.0, ht])
rx = jnp.stack([d, jnp.zeros_like(d), jnp.full_like(d, hr)], axis=-1)

# Reflection bounce point on the ground plane (z=0)
bounce_x = d * (ht / (ht + hr))
bounce = jnp.stack([bounce_x, jnp.zeros_like(d), jnp.zeros_like(d)], axis=-1)
tx_arr = jnp.broadcast_to(tx, rx.shape)

los_paths = TracedPaths(
    vertices=jnp.stack([tx_arr, rx], axis=-2),
    objects=-jnp.ones((*d.shape, 2), dtype=int),
    mask=jnp.ones(d.shape, dtype=bool),
    interaction_types=jnp.zeros((*d.shape, 0), dtype=int),
)

refl_paths = TracedPaths(
    vertices=jnp.stack([tx_arr, bounce, rx], axis=-2),
    objects=jnp.stack(
        [
            -jnp.ones_like(d, dtype=int),
            jnp.zeros_like(d, dtype=int),
            -jnp.ones_like(d, dtype=int),
        ],
        axis=-1,
    ),
    mask=jnp.ones(d.shape, dtype=bool),
    interaction_types=jnp.full((*d.shape, 1), int(InteractionType.REFLECTION)),
)

los_field = compute_received_fields(los_paths, ground, frequency)
refl_field = compute_received_fields(refl_paths, ground, frequency)
total_field = los_field + refl_field

with dplt.reuse(backend="plotly") as fig:
    ground.plot(opacity=0.5)
    dplt.draw_markers(
        tx[None, :],
        labels=["tx"],
        marker={"color": "red", "size": 5},
        name="tx",
    )
    dplt.draw_markers(
        rx[::20], marker={"color": "blue", "size": 3}, name="sampled rx"
    )
    dplt.draw_markers(
        bounce[::20],
        marker={"color": "orange", "size": 3},
        name="specular bounce point",
    )
    fig.show()

plt.figure(figsize=(8, 4.5))
plt.plot(
    d,
    compute_received_power(los_field),
    linestyle="--",
    label="Direct (LoS)",
)
plt.plot(
    d,
    compute_received_power(refl_field),
    linestyle=":",
    label="Reflected only",
)
plt.plot(
    d,
    compute_received_power(total_field),
    label="Total (two-ray interference)",
)
plt.xlabel("Receiver distance (m)")
plt.ylabel("Received power (dBW)")
plt.title("Two-ray interference creates multipath fading fringes")
plt.legend()
plt.grid(visible=True, alpha=0.3)
plt.show()
../_images/cc28aaa93866648390bb7e0b08a233f8155e4bf3cd06a9831e11cb32667a673c.png

2. Edge Diffraction: The Shadow Region & UTD#

Classical geometrical optics predicts abrupt discontinuities at shadow boundaries:

  • Incident Shadow Boundary (ISB): where the direct LoS ray becomes occluded.

  • Reflection Shadow Boundary (RSB): where the specular reflected ray can no longer exist.

The Uniform Theory of Diffraction (UTD) introduces diffraction coefficients containing transition Fresnel integrals \(F(X)\) that compensate for these geometric shadow discontinuities, ensuring the total field remains continuous across shadow boundaries.

Hide code cell source

# Canonical 90-degree wedge made of metal plates sharing an edge along the
# z-axis. This is Sionna RT's own "simple_wedge" scene, and the TX/RX
# placement below matches its wedge diffraction tutorial exactly
# (https://nvlabs.github.io/sionna/rt/tutorials/Diffraction.html), so the
# two can be compared directly.
wedge = Mesh(
    vertices=jnp.array([
        [0.0, -30.0, -15.0],
        [0.0, -30.0, 15.0],
        [0.0, 0.0, 15.0],
        [0.0, 0.0, -15.0],
        [30.0, 0.0, 15.0],
        [30.0, 0.0, -15.0],
    ]),
    triangles=jnp.array([[0, 1, 2], [0, 2, 3], [3, 2, 4], [3, 4, 5]]),
    face_materials=jnp.array([0, 0, 0, 0]),
    material_names=("Metal",),
)

freq_wedge = 1e9

# TX: 30 deg from the wedge's 0-face, 50 m from the edge.
tx_angle = jnp.radians(30.0)
tx_dist = 50.0
tx_diff = tx_dist * jnp.array([jnp.cos(tx_angle), jnp.sin(tx_angle), 0.0])

# RX: an arc 5 m from the edge, spanning the wedge's *entire* 270-degree
# exterior (free-space) angle, i.e., all 3 quadrants not occupied by the
# wedge itself -- not just the 2 quadrants nearest the ISB/RSB.
angle = jnp.linspace(jnp.radians(0.01), 1.5 * jnp.pi - jnp.radians(0.01), 1000)
angle_deg = jnp.degrees(angle)
radius = 5.0
rx_diff = radius * jnp.stack(
    [jnp.cos(angle), jnp.sin(angle), jnp.zeros_like(angle)], axis=-1
)

diff_scene = Scene(transmitters=tx_diff, receivers=rx_diff, mesh=wedge)
los_paths = diff_scene.trace_paths(order=0)
refl_paths = diff_scene.trace_paths(order=1)

los_field = compute_received_fields(los_paths, wedge, freq_wedge)[..., 0]
refl_field = compute_received_fields(refl_paths, wedge, freq_wedge).sum(axis=-1)

# Half-edge index 5 corresponds to the shared wedge edge along the z-axis
diffraction_paths = TracedPaths(
    vertices=jnp.stack(
        [
            jnp.broadcast_to(tx_diff, rx_diff.shape),
            jnp.broadcast_to(jnp.zeros(3), rx_diff.shape),
            rx_diff,
        ],
        axis=-2,
    ),
    objects=jnp.stack(
        [
            -jnp.ones(rx_diff.shape[0], dtype=int),
            jnp.full(rx_diff.shape[0], 5, dtype=int),
            -jnp.ones(rx_diff.shape[0], dtype=int),
        ],
        axis=-1,
    ),
    mask=jnp.ones(rx_diff.shape[0], dtype=bool),
    interaction_types=jnp.full(
        (rx_diff.shape[0], 1), InteractionType.DIFFRACTION
    ),
)

diff_field = compute_received_fields(diffraction_paths, wedge, freq_wedge)

los_valid = jnp.any(los_paths.mask, axis=-1)
refl_valid = jnp.any(refl_paths.mask, axis=-1)

# RSB/ISB are where a valid path becomes invalid as the angle increases (a
# "falling edge" in the mask); near-grazing incidence right at the start of
# the sweep (angle ~ 0, next to the 0-face) can trigger a spurious
# opposite-direction transition, so only falling edges are considered. RSB
# comes before ISB here, since TX sits close to the 0-face: receivers first
# cross into the reflection shadow, then (much further along) into the
# incident shadow.
refl_falling = jnp.diff(refl_valid.astype(int)) < 0
los_falling = jnp.diff(los_valid.astype(int)) < 0
rsb = float(angle_deg[jnp.nonzero(refl_falling)[0][0] + 1])
isb = float(angle_deg[jnp.nonzero(los_falling)[0][0] + 1])

los_only = jnp.where(los_valid, los_field, 0.0)
refl_only = jnp.where(refl_valid, refl_field, 0.0)
total_field = los_only + refl_only + diff_field

p_los_masked = jnp.where(los_valid, compute_received_power(los_field), jnp.nan)
p_refl_masked = jnp.where(
    refl_valid, compute_received_power(refl_field), jnp.nan
)
p_diff = compute_received_power(diff_field)
p_tot = compute_received_power(total_field)

with dplt.reuse(backend="plotly") as fig:
    wedge.plot(opacity=0.5)
    dplt.draw_markers(
        tx_diff[None, :],
        labels=["tx"],
        marker={"color": "red", "size": 5},
        name="tx",
    )
    dplt.draw_markers(
        jnp.zeros(3)[None, :],
        labels=["edge point"],
        marker={"color": "green", "size": 5},
        name="diffraction point",
    )
    dplt.draw_markers(
        rx_diff[:: rx_diff.shape[0] // 36],
        marker={"color": "blue", "size": 3},
        name="sampled rx",
    )
    fig.show()

plt.figure(figsize=(9, 5))
plt.plot(
    angle_deg,
    p_tot,
    label="Total (LoS + Refl + Diff)",
    color="k",
    linewidth=1.2,
)
plt.plot(angle_deg, p_los_masked, label="LoS", linewidth=0.9)
plt.plot(angle_deg, p_refl_masked, label="Reflected", linewidth=0.9)
plt.plot(angle_deg, p_diff, label="Diffracted", color="#2ca02c", linewidth=0.9)

plt.axvline(
    rsb, color="gray", linestyle="--", linewidth=1, label=f"RSB ({rsb:.1f}°)"
)
plt.axvline(
    isb, color="gray", linestyle=":", linewidth=1, label=f"ISB ({isb:.1f}°)"
)

# Leave headroom above the highest curve peak so the region labels sit in
# clear space instead of overlapping the two-ray interference fringes.
curve_top = float(jnp.nanmax(p_tot))
ymax = curve_top + 16
ymin = ymax - 90
plt.ylim(ymin, ymax)
plt.xlim(float(angle_deg[0]), float(angle_deg[-1]))

for x, label in (
    (rsb / 2, "Region I\n(LoS + reflection)"),
    ((rsb + isb) / 2, "Region II\n(LoS only)"),
    (
        (isb + float(angle_deg[-1])) / 2,
        "Region III\n(diffraction only)",
    ),
):
    plt.text(
        x,
        curve_top + 4,
        label,
        ha="center",
        va="bottom",
        bbox={"facecolor": "none", "edgecolor": "black", "pad": 4.0},
    )

plt.xlabel("Receiver angle around the shared edge (deg)")
plt.ylabel("Received power (dBW)")
plt.title(
    f'f = {freq_wedge / 1e9:g} GHz ("Metal"): UTD keeps the total field '
    "continuous across the RSB and ISB"
)
plt.legend(loc="lower left")
plt.grid(visible=True, alpha=0.3)
plt.show()
../_images/ce180ca3556dfb62eb8264daa2fa84f78eb78c8165db723e569aaa266acae8f4.png

3. Transmission: Attenuation Through a Wall#

Transmission accounts for waves passing through finite-thickness dielectric obstacles. A material requires an explicit thickness parameter \(d\). The ray traverses the wall, attenuated by the slab transmission coefficient \(T(\theta)\).

Hide code cell source

freq_diff = 3.5e9  # 3.5 GHz carrier frequency

wall_mesh = (
    Mesh(
        vertices=jnp.array([
            [0.0, -30.0, -10.0],
            [0.0, -30.0, 10.0],
            [0.0, 30.0, 10.0],
            [0.0, 30.0, -10.0],
        ]),
        triangles=jnp.array([[0, 1, 2], [0, 2, 3]]),
    )
    .set_face_materials(0)
    .set_materials("itu_concrete")
)

concrete_slab = Material(
    name="itu_concrete",
    properties=materials["itu_concrete"].properties,
    thickness=0.2,  # 20 cm thick concrete slab
)

tx_trans = jnp.array([-5.0, 0.0, 0.0])
y_rx = jnp.linspace(-20.0, 20.0, 300)
rx_trans = jnp.stack(
    [jnp.full_like(y_rx, 5.0), y_rx, jnp.zeros_like(y_rx)], axis=-1
)
bounce_trans = jnp.stack(
    [jnp.zeros_like(y_rx), 0.5 * y_rx, jnp.zeros_like(y_rx)], axis=-1
)
tx_arr_trans = jnp.broadcast_to(tx_trans, rx_trans.shape)

trans_paths = TracedPaths(
    vertices=jnp.stack([tx_arr_trans, bounce_trans, rx_trans], axis=-2),
    objects=jnp.stack(
        [
            -jnp.ones_like(y_rx, dtype=int),
            jnp.zeros_like(y_rx, dtype=int),
            -jnp.ones_like(y_rx, dtype=int),
        ],
        axis=-1,
    ),
    mask=jnp.ones_like(y_rx, dtype=bool),
    interaction_types=jnp.full(
        (*y_rx.shape, 1), int(InteractionType.TRANSMISSION)
    ),
)

free_paths = TracedPaths(
    vertices=jnp.stack([tx_arr_trans, rx_trans], axis=-2),
    objects=-jnp.ones((*y_rx.shape, 2), dtype=int),
    mask=jnp.ones_like(y_rx, dtype=bool),
    interaction_types=jnp.zeros((*y_rx.shape, 0), dtype=int),
)

free_space_power = compute_received_power(
    compute_received_fields(free_paths, wall_mesh, freq_diff)
)
trans_power = compute_received_power(
    compute_received_fields(
        trans_paths,
        wall_mesh,
        freq_diff,
        radio_materials={"itu_concrete": concrete_slab},
    )
)
with dplt.reuse(backend="plotly") as fig:
    wall_mesh.plot(opacity=0.5)
    dplt.draw_markers(
        tx_trans[None, :],
        labels=["tx"],
        marker={"color": "red", "size": 5},
        name="tx",
    )
    dplt.draw_markers(
        rx_trans[::15], marker={"color": "blue", "size": 3}, name="sampled rx"
    )
    dplt.draw_markers(
        bounce_trans[::15],
        marker={"color": "orange", "size": 3},
        name="wall crossing point",
    )
    fig.show()

plt.figure(figsize=(8, 4.5))
plt.plot(y_rx, free_space_power, linestyle="--", label="Free space (no wall)")
plt.plot(y_rx, trans_power, label="Through 20cm concrete wall")
plt.xlabel("Receiver position along the wall (m)")
plt.ylabel("Received power (dBW)")
plt.title("Transmission: insertion loss across dielectric slab")
plt.legend()
plt.grid(visible=True, alpha=0.3)
plt.show()
../_images/cf4f8f1c1df02c1b31b719b2550198a41775a6787fc58d2ddd25bc0c7cf2c423.png

4. Diffuse Scattering: Lambertian Rough Surfaces#

Real-world rough surfaces diffuse energy across a broad angular distribution. Using a material’s scattering_coefficient \(S \in [0, 1]\), scattered fields follow a Lambertian cosine pattern centered around the surface normal vector \(\hat{n}\).

Hide code cell source

rough_wall = Material(
    name="itu_concrete",
    properties=materials["itu_concrete"].properties,
    scattering_coefficient=0.6,
)

tx_scat = jnp.array([-8.0, 0.0, 8.0])
bounce_scat = jnp.array([0.0, 0.0, 0.0])
radius_scat = 15.0
elevation_deg = jnp.linspace(-80.0, 80.0, 300)  # measured from surface normal
elevation_rad = jnp.radians(elevation_deg)
rx_scat = bounce_scat + radius_scat * jnp.stack(
    [
        -jnp.cos(elevation_rad),
        jnp.zeros_like(elevation_rad),
        jnp.sin(elevation_rad),
    ],
    axis=-1,
)

tx_arr_scat = jnp.broadcast_to(tx_scat, rx_scat.shape)
bounce_arr_scat = jnp.broadcast_to(bounce_scat, rx_scat.shape)

scat_paths = TracedPaths(
    vertices=jnp.stack([tx_arr_scat, bounce_arr_scat, rx_scat], axis=-2),
    objects=jnp.stack(
        [
            -jnp.ones_like(elevation_deg, dtype=int),
            jnp.zeros_like(elevation_deg, dtype=int),
            -jnp.ones_like(elevation_deg, dtype=int),
        ],
        axis=-1,
    ),
    mask=jnp.ones_like(elevation_deg, dtype=bool),
    interaction_types=jnp.full(
        (*elevation_deg.shape, 1), int(InteractionType.SCATTERING)
    ),
)

scattering_power = compute_received_power(
    compute_received_fields(
        scat_paths,
        wall_mesh,
        freq_diff,
        radio_materials={"itu_concrete": rough_wall},
    )
)

# TX sits at +45 deg from the surface normal, so the specular direction is at -45 deg
specular_elevation = -45.0

with dplt.reuse(backend="plotly") as fig:
    wall_mesh.plot(opacity=0.5)
    dplt.draw_markers(
        tx_scat[None, :],
        labels=["tx"],
        marker={"color": "red", "size": 5},
        name="tx",
    )
    dplt.draw_markers(
        bounce_scat[None, :],
        labels=["bounce point"],
        marker={"color": "green", "size": 5},
        name="scattering point",
    )
    dplt.draw_markers(
        rx_scat[::15], marker={"color": "blue", "size": 3}, name="sampled rx"
    )
    fig.show()

plt.figure(figsize=(8, 4.5))
plt.plot(
    elevation_deg,
    scattering_power,
    label="Diffusely scattered power",
)
plt.axvline(
    specular_elevation,
    color="gray",
    linestyle="--",
    label="Specular direction (-45°)",
)
plt.xlabel("Receiver elevation from the surface normal (deg)")
plt.ylabel("Received power (dBW)")
plt.title(
    "Diffuse scattering produces a broad lobe centered on the surface normal"
)
plt.legend()
plt.grid(visible=True, alpha=0.3)
plt.show()
../_images/d4972cad769b8e528819510247cbc788a0aa9da7f09eb6d7dbbda2af3dab21e4.png

5. Beyond Lambertian: Directive & Backscattering Patterns#

The Lambertian pattern above spreads diffusely-scattered power broadly around the surface normal \(\hat{n}\), independently of the incidence direction. Real rough surfaces, however, often scatter more power near the specular direction \(\hat{k}_{sp}\) (a directive lobe), or even back toward the transmitter (a retroreflective, or backscattering, lobe, as observed e.g. on foliage or building corners).

Material.scattering_pattern accepts any AbstractScatteringPattern subclass, so swapping in DirectivePattern or BackscatteringPattern (both following Degli-Esposti’s directive scattering model) is enough to change a material’s scattering behavior, without touching the field solver itself: this is exactly the point of the AbstractScatteringPattern abstraction. Implementing a new scattering physics only requires subclassing it and providing a __call__ method.

directive_wall = Material(
    name="itu_concrete",
    properties=materials["itu_concrete"].properties,
    scattering_coefficient=0.6,
    scattering_pattern=DirectivePattern(alpha_r=6.0),
)
backscattering_wall = Material(
    name="itu_concrete",
    properties=materials["itu_concrete"].properties,
    scattering_coefficient=0.6,
    scattering_pattern=BackscatteringPattern(
        alpha_r=6.0, alpha_i=6.0, lambda_=0.2
    ),
)

retroreflection_elevation = -specular_elevation  # +45 deg, back toward the TX

plt.figure(figsize=(8, 4.5))
for label, material in (
    ("Lambertian", rough_wall),
    (r"Directive ($\alpha_R = 6$)", directive_wall),
    (
        r"Backscattering ($\alpha_R = \alpha_I = 6$, $\Lambda = 0.2$)",
        backscattering_wall,
    ),
):
    power = compute_received_power(
        compute_received_fields(
            scat_paths,
            wall_mesh,
            freq_diff,
            radio_materials={"itu_concrete": material},
        )
    )
    plt.plot(elevation_deg, power, label=label)

plt.axvline(
    specular_elevation,
    color="gray",
    linestyle="--",
    label="Specular direction (-45°)",
)
plt.axvline(
    retroreflection_elevation,
    color="black",
    linestyle=":",
    label="Retroreflection direction (+45°)",
)
plt.xlabel("Receiver elevation from the surface normal (deg)")
plt.ylabel("Received power (dBW)")
plt.title(
    "Directive and backscattering patterns concentrate power away from the normal"
)
plt.legend()
plt.grid(visible=True, alpha=0.3)
plt.show()
../_images/3ba0bf213c007f8d020d3c98659599aa1863e2c5e418432852713d800bca9e03.png

Summary & Custom Solvers#

DiffeRT provides a modular EM solver architecture: