FieldProjectionAngleMonitor weirdness

yaugenst commented 10 months ago

We have come across some potential issues with the td.FieldProjectionAngleMonitor. In particular, we seem to be unable to reproduce far field radiation patterns when compared to either Lumerical or our own implementation of the Fraunhofer approximation, which both match each other (courtesy of @bhnzhang). Maybe we are not using the monitor correctly? But what is especially weird is that the monitor does not seem to converge with higher resolutions. Is there documentation on what exactly the td.FieldProjectionAngleMonitor implements? Since there is a proj_distance parameter, I'm guessing it's a different calculation.

This is what our comparison looks like: gc_convergence

It is strange that the two plots (very roughly) match only for the lowest resolution. If we zoom in a little on the higher resolutions the difference becomes even more pronounced: gc_convergence

This is the simulation setup: screenshot_20231201_105453

Here is the code to reproduce these results. Note that we originally came across this issue in a different setting (grating coupler simulations), but I reduced it to the MWE below:

#!/usr/bin/env python3

import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
import xrft
from tidy3d.web import Batch, BatchData

def near2far(field, angles=(-15, 5), pad_mult=5, downsample=1):
    nx = int(np.floor(len(field.x) / downsample))
    x_interp = np.linspace(field.x[0], field.x[-1], nx)

    field_interp = field.interp(x=x_interp)

    npad = pad_mult * len(x_interp)
    field_pad = xrft.padding.pad(field_interp, {"x": npad})
    field_k = xrft.fft(field_pad)

    lambda0 = np.squeeze(td.C_0 / field.f).values
    deg = np.rad2deg(field_k.freq_x * lambda0)
    field_k = field_k.assign_coords(thetax=deg)
    return field_k.where(
        (field_k.thetax >= angles[0]) & (field_k.thetax <= angles[1]), drop=True
    )

def make_gc_sim(resolution, angles=(-15, 5)):
    rng = np.random.default_rng(123456)

    sim_size = (40, 0, 2.0)
    lcen = 1.0
    fcen = td.C_0 / lcen

    grid_spec = td.GridSpec.auto(wavelength=lcen, min_steps_per_wvl=resolution)
    boundary_spec = td.BoundarySpec(
        x=td.Boundary.pml(),
        y=td.Boundary.periodic(),
        z=td.Boundary.pml(),
    )

    structures = [
        td.Structure(
            geometry=td.Box(
                center=(cx, 0, 0), size=(rng.uniform(0.2, 0.8, 1), td.inf, 0.5)
            ),
            medium=td.Medium(name="test", permittivity=8),
        )
        for cx in np.linspace(-10, 10, 20)
    ]

    nf_monitor = td.FieldMonitor(
        center=(0, 0, -0.5),
        size=(td.inf, td.inf, 0),
        freqs=[fcen],
        colocate=True,
        name="near_fields",
    )

    ff_monitor = td.FieldProjectionAngleMonitor(
        center=nf_monitor.center,
        size=nf_monitor.size,
        freqs=nf_monitor.freqs,
        normal_dir="-",
        phi=[0.0],
        theta=np.linspace(*map(np.deg2rad, angles), 501),
        colocate=True,
        name="far_fields",
    )

    source = td.GaussianBeam(
        center=(0, 0, 0.5),
        size=(td.inf, td.inf, 0),
        source_time=td.GaussianPulse(freq0=fcen, fwidth=fcen / 10),
        direction="-",
        waist_radius=20,
        pol_angle=np.pi / 2,
        angle_theta=np.deg2rad(5),
    )

    sim = td.Simulation(
        size=sim_size,
        sources=[source],
        structures=structures,
        monitors=[nf_monitor, ff_monitor],
        grid_spec=grid_spec,
        boundary_spec=boundary_spec,
        shutoff=1e-6,
        run_time=5e-12,
    )

    return sim

def main():
    run_dir = "ff_test"

    batch = Batch(simulations={res: make_gc_sim(res) for res in (10, 20, 30, 40, 50)})
    batch.run(path_dir=run_dir)

    batch_data = BatchData.load(run_dir)

    fig, ax = plt.subplots(2, 1, sharex=True, tight_layout=True, dpi=100)

    for resolution in batch_data.task_ids:
        sim_data = batch_data.load_sim_data(resolution)

        power = np.squeeze(sim_data["far_fields"].power)

        ek = near2far(np.squeeze(sim_data["near_fields"].Ey))
        abs2k = np.abs(ek) ** 2

        ax[0].plot(np.rad2deg(power.theta), power, label=resolution)
        ax[1].plot(abs2k.thetax, abs2k, label=resolution)

    ax[0].set_title("FieldProjectionAngleMonitor")
    ax[0].legend(title="Resolution")
    ax[1].set_title("Fraunhofer approximation")
    ax[1].set_xlabel("Angle (deg)")
    fig.savefig("gc_convergence.png", bbox_inches="tight", pad_inches=0.1, dpi=300)

if __name__ == "__main__":
    main()

The match does become better when the near field is not perturbed as much.

momchil-flex commented 10 months ago

Interesting! We have this example showcasing various filed pojrection results, as well as various tests vs analytic results (e.g. far field emission of a dipole, scattering from nanoparticles,...). We know that our projection works as expected in many cases, so it's important to figure out what is going on here. The dependence on resolution is indeed very weird. The projection distance should not matter when you have the far field approximation on, which is true by default (more precisely, it should just be an overall scaling factor).

@shashwat-sh could you take a look?