How to perform near field to far field projections in Tidy3D FDTD¶

sim = td.Simulation(
    size=sim_size,
    center=[0, 0, 0],
    grid_spec=td.GridSpec.uniform(dl=wavelength / min_cells_per_wvl),
    structures=geometry,
    sources=[source],
    monitors=[monitor_near],
    run_time=run_time,
    boundary_spec=boundary_spec,
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim.plot(x=0, ax=ax1)
sim.plot(y=0, ax=ax2)
plt.show()

sim_data = web.run(
    sim, task_name="aperture_1", path="data/aperture_1.hdf5", verbose=True
)

# far field projection monitor
monitor_far = td.FieldProjectionAngleMonitor(
    center=[
        0,
        offset_mon,
        0,
    ],  # the monitor's center defined the local origin - the projection distance
    # and angles will all be measured with respect to this local origin
    size=[td.inf, 0, td.inf],
    # the size and center of any far field monitor should indicate where the *near* fields are recorded
    freqs=[f0],
    name="far_field",
    phi=list(phi_proj),
    theta=list(theta_proj),
    proj_distance=r_proj,
    far_field_approx=True,  # we leave this to its default value of 'True' because we are interested in fields sufficiently
    # far away that geometric far field approximations can be invoked to speed up the calculation
)

# helper functin to call the projector
def get_proj_fields(sim_data, monitor_near, monitor_far, pts_per_wavelength=10):
    # object that does projections is constructed using the near-field monitor, because those are the fields to be projected
    projector = td.FieldProjector.from_near_field_monitors(
        sim_data=sim_data,
        near_monitors=[monitor_near],
        normal_dirs=["+"],  # we are projecting along the + direction
        pts_per_wavelength=pts_per_wavelength,  # to speed up calculations, the fields on the near-field monitor can be downsampled to these
        # many points per wavelength (default is already 10)
    )
    return projector.project_fields(monitor_far)


# execute the projector, with the far field monitor as input, to do the projection
# let's also time this, for later use
import time

t0 = time.perf_counter()
projected_field_data = get_proj_fields(sim_data, monitor_near, monitor_far)
t1 = time.perf_counter()
proj_time = t1 - t0

def make_field_plot(phi, theta, vals1, vals2):
    n_plots = 2
    fig, ax = plt.subplots(1, n_plots, tight_layout=True, figsize=(8, 3.8))
    im1 = ax[0].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals1),
        cmap="RdBu",
        shading="auto",
    )
    im2 = ax[1].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals2),
        cmap="RdBu",
        shading="auto",
    )
    fig.colorbar(im1, ax=ax[0])
    fig.colorbar(im2, ax=ax[1])
    ax[0].set_title("Analytic")
    ax[1].set_title("Field projection")
    for _ax in ax:
        _ax.set_xlabel("$\phi$ (deg)")
        _ax.set_ylabel("$\\theta$ (deg)")


# RMSE
def rmse(array_ref, array_test):
    error = array_test - array_ref
    rmse = np.sqrt(np.mean(np.abs(error.flatten()) ** 2))
    nrmse = rmse / np.abs(np.max(array_ref.flatten()) - np.min(array_ref.flatten()))
    return nrmse


# plot Etheta
Etheta_analytic = analytic_field_data.Etheta.isel(f=0, r=0)
Etheta_proj = projected_field_data.Etheta.isel(f=0, r=0)
make_field_plot(phi_proj, theta_proj, Etheta_analytic, Etheta_proj)

# print the normalized RMSE
print(
    f"Normalized root mean squared error: {rmse(Etheta_analytic.values, Etheta_proj.values) * 100:.2f} %"
)

plt.show()

sim_data2 = web.run(
    sim2, task_name="aperture_2", path="data/aperture_2.hdf5", verbose=True
)

# extract the computed projected fields
projected_field_data_server = sim_data2[monitor_far.name]

# plot Etheta
Etheta_proj_server = projected_field_data_server.Etheta.isel(f=0, r=0)
make_field_plot(phi_proj, theta_proj, Etheta_analytic, Etheta_proj_server)

# print the normalized RMSE
print(
    f"Normalized root mean squared error: {rmse(Etheta_analytic.values, Etheta_proj_server.values) * 100:.2f} %"
)

plt.show()

# use the simulation log to find the time taken for server-side computations
server_time = float(
    sim_data2.log.split("Field projection time (s):    ", 1)[1].split("\n", 1)[0]
)
print(f"Client-side field projection took {proj_time:.2f} s")
print(f"Server-side field projection took {server_time:.2f} s")

def make_cart_plot(phi, theta, vals1, vals2, vals3):
    n_plots = 3
    fig, ax = plt.subplots(1, n_plots, tight_layout=True, figsize=(12.6, 3.8))
    im1 = ax[0].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals1),
        cmap="RdBu",
        shading="auto",
    )
    im2 = ax[1].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals2),
        cmap="RdBu",
        shading="auto",
    )
    im3 = ax[2].pcolormesh(
        phi * 180 / np.pi,
        theta * 180 / np.pi,
        np.real(vals3),
        cmap="RdBu",
        shading="auto",
    )
    fig.colorbar(im1, ax=ax[0])
    fig.colorbar(im2, ax=ax[1])
    fig.colorbar(im3, ax=ax[2])
    ax[0].set_title("Ex")
    ax[1].set_title("Ey")
    ax[2].set_title("Ez")
    for _ax in ax:
        _ax.set_xlabel("$\phi$ (deg)")
        _ax.set_ylabel("$\\theta$ (deg)")


# get the fields in Cartesian coordinates from the projected data we already computed above
fields_cartesian = projected_field_data.fields_cartesian.isel(f=0, r=0)

# plot Ex, Ey, Ez
make_cart_plot(
    phi_proj, theta_proj, fields_cartesian.Ex, fields_cartesian.Ey, fields_cartesian.Ez
)

# get the power
power = projected_field_data.power.isel(f=0, r=0)

# plot the power
fig, ax = plt.subplots(1, 1, tight_layout=True, figsize=(4.3, 3.8))
im = ax.pcolormesh(
    phi_proj * 180 / np.pi,
    theta_proj * 180 / np.pi,
    power,
    cmap="inferno",
    shading="auto",
)
fig.colorbar(im, ax=ax)
_ = ax.set_title("Power")
_ = ax.set_xlabel("$\phi$ (deg)")
_ = ax.set_ylabel("$\\theta$ (deg)")

plt.show()

# new projection distance
r_proj_new = 20 * wavelength

# re-project our far field data above to this new distance
reprojected_field_data = projected_field_data.renormalize_fields(r_proj_new)

# now all the fields stored in 'projected_field_data' correspond to this new distance
# compare to the analytical fields at this new distance
analytic_field_data_new = analytic_fields_aperture(
    monitor_far, sim_size, height, width, r_proj_new
)

# plot Etheta
Etheta_analytic = analytic_field_data_new.Etheta.isel(f=0, r=0)
Etheta_proj = reprojected_field_data.Etheta.isel(f=0, r=0)
make_field_plot(phi_proj, theta_proj, Etheta_analytic, Etheta_proj)

# print the normalized RMSE
print(
    f"Normalized root mean squared error: {rmse(Etheta_analytic.values, Etheta_proj.values) * 100:.2f} %"
)

plt.show()

# project fields only a few wavelengths away from the aperture
r_proj_intermediate = 4 * wavelength

# create a field monitor to measure these fields at the intermediate projection distance,
# so that we have something to which we can compare the 'FieldProjector' results
monitor_intermediate = td.FieldMonitor(
    center=[0, offset_mon + r_proj_intermediate, 0],
    size=[td.inf, 0, td.inf],
    freqs=[f0],
    name="inter_field",
)

# make a larger simulation along y to accommodate the plane at which the intermediate fields need to be measured
shift = 1.2 * r_proj_intermediate
sim_size3 = [sim_size[0], sim_size[1] + shift, sim_size[2]]
# move the sim center
sim_center = [0, (sim_size[1] + shift) / 2 - sim_size[1] / 2, 0]
sim3 = td.Simulation(
    size=sim_size3,
    center=sim_center,
    grid_spec=td.GridSpec.uniform(dl=wavelength / min_cells_per_wvl),
    structures=geometry,
    sources=[source],
    monitors=[
        monitor_near,
        monitor_intermediate,
    ],  # provide both near field and intermediate field monitors
    run_time=run_time,
    boundary_spec=boundary_spec,
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim3.plot(x=0, ax=ax1)
sim3.plot(y=0, ax=ax2)
plt.show()

sim_data3 = web.run(
    sim3, task_name="aperture_3", path="data/aperture_3.hdf5", verbose=True
)

# make the projection monitor which projects fields without approximations
xs = np.linspace(-sim_size[0] / 2, sim_size[0] / 2, 100)
ys = np.linspace(-sim_size[1] / 2, sim_size[1] / 2, 100)
monitor_intermediate_proj = td.FieldProjectionCartesianMonitor(
    center=[
        0,
        offset_mon,
        0,
    ],  # the monitor's center defined the local origin - the projection distance
    # and angles will all be measured with respect to this local origin
    size=[td.inf, 0, td.inf],
    freqs=[f0],
    name="inter_field_proj",
    proj_axis=1,  # because we're projecting along the +y axis
    x=list(xs),  # local x coordinates - corresponds to global x in this case
    y=list(ys),  # local y coordinates - corresponds to global z in this case
    proj_distance=r_proj_intermediate,
    far_field_approx=False,  # turn off the far-field approximation (is 'True' by default)
)

# execute the projector, with the projection field monitor as input, to do the projection
# let's also time this, for later use
import time

t0 = time.perf_counter()
projected_field_data_noapprox = get_proj_fields(
    sim_data3, monitor_near, monitor_intermediate_proj
)
t1 = time.perf_counter()
proj_time_new = t1 - t0

# Compute projected fields *with* far field approximations, to facilitate an accuracy comparison
monitor_intermediate_proj_approx = td.FieldProjectionCartesianMonitor(
    center=[
        0,
        offset_mon,
        0,
    ],  # the monitor's center defined the local origin - the projection distance
    # and angles will all be measured with respect to this local origin
    size=[td.inf, 0, td.inf],
    freqs=[f0],
    name="inter_field_proj_approx",
    proj_axis=1,  # because we're projecting along the +y axis
    x=list(xs),  # local x coordinates - corresponds to global x in this case
    y=list(ys),  # local y coordinates - corresponds to global z in this case
    proj_distance=r_proj_intermediate,
    far_field_approx=True,  # turn on the far-field approximation
)

# execute the projector, with the projection field monitor as input, to do the projection
# let's also time this, for later use
import time

t0 = time.perf_counter()
projected_field_data_approx = get_proj_fields(
    sim_data3, monitor_near, monitor_intermediate_proj_approx
)
t1 = time.perf_counter()
proj_time_new_approx = t1 - t0

# let's see how long this took compared to the previous case when the approximations were turned on
print(
    f"Client-side field projection *with approximations on* took {proj_time_new_approx:.2f} s"
)
print(
    f"Client-side field projection *with approximations off* took {proj_time_new:.2f} s"
)

# plot the actual measured fields
fields_meas = sim_data3[monitor_intermediate.name].colocate(x=xs, z=ys)
make_cart_plot(
    ys,
    xs,
    fields_meas.Ex.isel(f=0, y=0),
    fields_meas.Ey.isel(f=0, y=0),
    fields_meas.Ez.isel(f=0, y=0),
)
plt.suptitle("Measured fields")

# projected field without approximations - get them in Cartesian coords
fields_proj_noapprox = projected_field_data_noapprox.fields_cartesian
make_cart_plot(
    ys,
    xs,
    fields_proj_noapprox.Ex.isel(f=0, y=0),
    fields_proj_noapprox.Ey.isel(f=0, y=0),
    fields_proj_noapprox.Ez.isel(f=0, y=0),
)
plt.suptitle("Projected, no approximations")

# projected field with approximations - get them in Cartesian coords
fields_proj_approx = projected_field_data_approx.fields_cartesian
make_cart_plot(
    ys,
    xs,
    fields_proj_approx.Ex.isel(f=0, y=0),
    fields_proj_approx.Ey.isel(f=0, y=0),
    fields_proj_approx.Ez.isel(f=0, y=0),
)
_ = plt.suptitle("Projected, with far field approximations")

# RMSE
Emag_meas = np.sqrt(
    np.abs(fields_meas.Ex) ** 2
    + np.abs(fields_meas.Ey) ** 2
    + np.abs(fields_meas.Ez) ** 2
)
Emag_proj_noapprox = np.sqrt(
    np.abs(fields_proj_noapprox.Ex) ** 2
    + np.abs(fields_proj_noapprox.Ey) ** 2
    + np.abs(fields_proj_noapprox.Ez) ** 2
)
Emag_proj_approx = np.sqrt(
    np.abs(fields_proj_approx.Ex) ** 2
    + np.abs(fields_proj_approx.Ey) ** 2
    + np.abs(fields_proj_approx.Ez) ** 2
)
print(
    f"Normalized RMSE for |E|, no far field approximation: {rmse(Emag_meas.values, Emag_proj_noapprox.values) * 100:.2f} %"
)
print(
    f"Normalized RMSE for |E|, with far field approximation: {rmse(Emag_meas.values, Emag_proj_approx.values) * 100:.2f} %"
)

plt.show()

sim4 = td.Simulation(
    size=sim_size,
    center=[0, 0, 0],
    grid_spec=td.GridSpec.uniform(dl=wavelength / min_cells_per_wvl),
    structures=geometry,
    sources=[source],
    monitors=[monitor_intermediate_proj],  # only need to supply the projection monitor
    run_time=run_time,
    boundary_spec=boundary_spec,
)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim4.plot(x=0, ax=ax1)
sim4.plot(y=0, ax=ax2)
plt.show()

# run the simulation
sim_data4 = web.run(
    sim4, task_name="aperture_4", path="data/aperture_4.hdf5", verbose=True
)

# extract the projected fields as before and plot them
projected_field_data_server = sim_data4[monitor_intermediate_proj.name]

# plot the actual measured fields from the previous simulation
fields_meas = sim_data3[monitor_intermediate.name].colocate(x=xs, z=ys)
make_cart_plot(
    ys,
    xs,
    fields_meas.Ex.isel(f=0, y=0),
    fields_meas.Ey.isel(f=0, y=0),
    fields_meas.Ez.isel(f=0, y=0),
)
plt.suptitle("Measured fields")

# projected field without approximations computed on the server
fields_proj_noapprox = projected_field_data_server.fields_cartesian
make_cart_plot(
    ys,
    xs,
    fields_proj_noapprox.Ex.isel(f=0, y=0),
    fields_proj_noapprox.Ey.isel(f=0, y=0),
    fields_proj_noapprox.Ez.isel(f=0, y=0),
)
plt.suptitle("Projected, no approximations, computed on the server")

# RMSE
Emag_proj_server = np.sqrt(
    np.abs(fields_proj_noapprox.Ex) ** 2
    + np.abs(fields_proj_noapprox.Ey) ** 2
    + np.abs(fields_proj_noapprox.Ez) ** 2
)
print(
    f"Normalized RMSE for |E|, no far field approximation, computed on the server: {rmse(Emag_meas.values, Emag_proj_server.values) * 100:.2f} %\n"
)

# use the simulation log to find the time taken for server-side computations
server_time = float(
    sim_data4.log.split("Field projection time (s):    ", 1)[1].split("\n", 1)[0]
)
print(
    f"Client-side field projection *without approximations* took {proj_time_new:.2f} s"
)
print(f"Server-side field projection *without approximations* took {server_time:.2f} s")

plt.show()