Scattering cross-section calculation#

This tutorial will show you how to compute the radar cross section (RCS) for a dielectric sphere by sampling scattered near fields on a closed surface surrounding the sphere, and transforming them to observation points far away.

This example demonstrates the usefulness of the near field to far field transformation for reducing the simulation size needed for structures involving lots of empty space.

To obtain the scattered field, we will run two simulations: * one with the sphere to compute the total near field on a closed surface around the sphere, and * one without the sphere, to compute just the incident field and subtract it from the above to get the scattered field.

Then, we’ll show how to use the Near2Far feature from Tidy3D to compute the RCS for the sphere using the near field data.

[1]:

# standard python imports
import numpy as np
import matplotlib.pyplot as plt

# tidy3d imports
import tidy3d as td
import tidy3d.web as web


Define Simulation Parameters#

We first need to define our simulation parameters and the structure.

[2]:

# radius and location of the sphere
center = [0, 0, 0]

# permittivity of the sphere
epsr = 4

# free space central wavelength
wavelength = (2.0 * radius) / 2.0
f0 = td.C_0 / wavelength

# distance between the surface of the sphere and the start of the PML layers along each cartesian direction
buffer_PML = 3 * wavelength

# distance between the sphere and the near field monitor along each cartesian direction
buffer_mon = 1 * wavelength

# Define material properties
air = td.Medium(permittivity=1)
diel = td.Medium(permittivity=epsr)

# resolution control
min_steps_per_wvl = 24

# create the sphere
geometry = [sphere]

# define PML layers
pml_layers = 3*[td.PML()]

# set the domain size in x, y, and z
domain_size = buffer_PML + 2 * radius + buffer_PML

# construct simulation size array
sim_size = (domain_size, domain_size, domain_size)


Create Source#

For our incident field, we create a plane wave incident from below the sphere polarized in the x direction.

[3]:

# Bandwidth in Hz
fwidth = f0 / 10.0

# Gaussian source offset; the source peak is at time t = offset/fwidth
offset = 4.

# time dependence of source
gaussian = td.GaussianPulse(freq0=f0, fwidth=fwidth, offset=offset)

# place the source below the sphere, propagating in the +z direction
source = td.PlaneWave(
center=(0,0,-(radius + 3 * buffer_PML / 4)),
size=(td.inf, td.inf, 0),
source_time=gaussian,
direction='+',
pol_angle=0)

# Simulation run time past the source decay (around t=2*offset/fwidth)
run_time = 100 / fwidth



Create Monitor#

Next, we define the monitors that will capture the near field data.

First, we create a FieldMonitor completely enclosing the sphere.

The, using the .surfaces() method, we can extract the 6 planar surfaces surrounding the volume and feed these to the simulation.

This cuts down on the data required vs computing the full volume because only the fields on the enclosing surface are required to get the far field information.

[4]:

# create a volume monitor around the sphere
mon_size = 2 * radius + 2 * buffer_mon
#mon_size = wavelength
monitor = td.FieldMonitor(
center=center,
size=[mon_size, mon_size, mon_size],
freqs=[f0],
name='near_field')

# get the surface monitors associated with the volume monitor
monitors = monitor.surfaces()


Create Simulation#

Now we can put everything together and define the two simulation classes: with the sphere to get the total field, and without the sphere to get the incident field. A uniform grid with a fixed grid size is used in both simulations to allow easily subtracting the incident field in the empty simulation from the total field in the actual simulation.

[5]:

sim = td.Simulation(
size=sim_size,
grid_spec = td.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl),
structures=geometry,
sources=[source],
monitors=monitors,
run_time=run_time,
pml_layers=pml_layers
)

# Here, we add the sphere as an override structure for the messhing such
# that the grids of the two simulations match.
sim_empty = td.Simulation(
size=sim_size,
grid_spec = td.GridSpec.auto(
min_steps_per_wvl=min_steps_per_wvl,
override_structures=geometry,
),
structures=[],
sources=[source],
monitors=monitors,
run_time=run_time,
pml_layers=pml_layers
)


Visualize Geometry#

Let’s take a look and make sure everything is defined properly in both simulations.

[6]:

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim.plot_eps(y=0, ax=ax1);
sim_empty.plot_eps(y=0, ax=ax2);

[15:37:05] INFO     Auto meshing using wavelength 0.5000 defined from        grid_spec.py:466
sources.

           INFO     Auto meshing using wavelength 0.5000 defined from        grid_spec.py:466
sources.

<Figure size 648x216 with 4 Axes>


Run Simulations#

Now we can run both simulations over time and measure the results

[7]:

# Run simulation
import tidy3d.web as web

# compute scattered fields by subtracting out the incident fields from the empty simulation
for mon in monitors:
sim_data.monitor_data[mon.name].Ex -= sim_empty_data[mon.name].Ex
sim_data.monitor_data[mon.name].Ey -= sim_empty_data[mon.name].Ey
sim_data.monitor_data[mon.name].Ez -= sim_empty_data[mon.name].Ez

sim_data.monitor_data[mon.name].Hx -= sim_empty_data[mon.name].Hx
sim_data.monitor_data[mon.name].Hy -= sim_empty_data[mon.name].Hy
sim_data.monitor_data[mon.name].Hz -= sim_empty_data[mon.name].Hz

[15:37:06] INFO     Using Tidy3D credentials from stored file                      auth.py:74

[15:37:08] INFO     Uploaded task 'sphereRCS' with task_id                      webapi.py:120
'5024bd2b-53d8-44be-815e-44ce1382553f'.

[15:37:12] INFO     Maximum flex unit cost: 0.46                                webapi.py:253

           INFO     status = queued                                             webapi.py:262

[15:37:21] INFO     status = preprocess                                         webapi.py:274

[15:37:37] INFO     starting up solver                                          webapi.py:278

[15:37:56] INFO     running solver                                              webapi.py:284

[15:39:31] INFO     early shutoff detected, exiting.                            webapi.py:295

           INFO     status = postprocess                                        webapi.py:301

[15:39:45] INFO     status = success                                            webapi.py:307

[15:39:47] INFO     downloading file "monitor_data.hdf5" to                     webapi.py:537
"data/sphereRCS.hdf5"

[15:39:51] INFO     loading SimulationData from data/sphereRCS.hdf5             webapi.py:369

[15:39:53] INFO     Uploaded task 'sphereRCS_empty' with task_id                webapi.py:120
'60d58032-da91-45a3-b7e2-95db883789cb'.

[15:39:58] INFO     Maximum flex unit cost: 0.46                                webapi.py:253

           INFO     status = queued                                             webapi.py:262

[15:40:07] INFO     status = preprocess                                         webapi.py:274

[15:40:21] INFO     starting up solver                                          webapi.py:278

[15:40:41] INFO     running solver                                              webapi.py:284

[15:40:50] INFO     early shutoff detected, exiting.                            webapi.py:295

           INFO     status = postprocess                                        webapi.py:301

[15:41:03] INFO     status = success                                            webapi.py:307

[15:41:05] INFO     downloading file "monitor_data.hdf5" to                     webapi.py:537
"data/sphereRCS_empty.hdf5"

[15:41:07] INFO     loading SimulationData from data/sphereRCS_empty.hdf5       webapi.py:369


Setting Up Near2Far#

To set up the near-to-far transformation, we need to grab the fields on each surface of the near-field FieldMonitor objects.

So, we simply create a Near2Far object and pass in the surface monitors as shown below. Note that we also need to pass the normal directions of each of the monitors in the list.

In addition to storing the near field data, this object will compute he surface currents and provide various methods for projecting the far field quantities.

We can optionally pass in the number of points per wavelength in the background medium with which to sample fields on the monitors. The default is 10 points per wavelength.

One can also pass in coordinates for the local origin of the set of monitors; the far-field observation points will be defined with respect to this origin. By default, the local origin is set to the average of the centers of all surface monitors passed in.

[8]:

from tidy3d.plugins import Near2Far

n2f = Near2Far.from_surface_monitors(
sim_data=sim_data,
monitors=monitors,
normal_dirs=['-','+','-','+','-','+'],
frequency=f0,
pts_per_wavelength=10
)


Getting Far Field Data#

After the Near2Far object is initialized, each of its surface surface currents are computed.

Then, we just need to call one of its methods to get a far field quantity.

For this example, we use Near2Far.radar_cross_section(theta,phi) to get the RCS at a theta,phi point relative to the center of the original monitor.

Note that this can be set or verified using the Near2Far.origin attribute.

[9]:

print(f'origin at {n2f.origin}')

origin at (0.0, 0.0, 0.0)

[10]:

far_distance = 100 * wavelength

num_theta = 300
num_phi = 2
thetas = np.linspace(0, np.pi, num_theta)
phis = np.linspace(0, np.pi/2, num_phi)



Plot Results#

Now we can plot the RCS and compare it to the analytical RCS computed via the Mie series.

The results match quite well, but there are some errors due to approximations inherent to the near field to far field transformation.

[11]:

def to_db(val):
val = val / np.max(np.abs(val))
return 10.0*np.log10(val)

RCS_phi0 = np.squeeze(RCS[:,0])
RCS_phi90 = np.squeeze(RCS[:,1])

# ------ import analytical data from disk ------

mie_file_id = '2lambda_epsr4'
mie_filename_phi0 = "./data/mie_bRCS_phi0_" + mie_file_id + ".txt"
mie_filename_phi90 = "./data/mie_bRCS_phi90_" + mie_file_id + ".txt"

mie_theta_phi0 = np.squeeze(mie_data_phi0[:,[0]])
mie_phi0 = np.squeeze(mie_data_phi0[:,[1]])

mie_theta_phi90 = np.squeeze(mie_data_phi90[:,[0]])
mie_phi90 = np.squeeze(mie_data_phi90[:,[1]])

# ------ plot for phi = 0 ------

fig, ax = plt.subplots(figsize=(7.5, 5))

ax.plot(mie_theta_phi0, to_db(mie_phi0), '-b', label="$\\phi = 0$, Mie")
ax.plot(thetas, to_db(RCS_phi0), '--r', label="$\\phi = 0$, near2far")
ax.set(xlabel = "$\\theta$ (degrees)",
ylabel = "Bistatic RCS (dBsm)",
yscale = 'linear',
xscale = 'linear'
)
ax.grid(visible = True, which = 'both', axis = 'both', linewidth = 0.4)
plt.legend(loc='best', prop={'size': 14})
plt.tight_layout()

# ------ plot for phi = pi/2 ------

fig, ax = plt.subplots(figsize=(7.5, 5))

ax.plot(mie_theta_phi90, to_db(mie_phi90), '-b', label="$\\phi = \\pi/2$, Mie")
ax.plot(thetas, to_db(RCS_phi90), '--r', label="$\\phi = \\pi/2$, near2far")
ax.set(xlabel = "$\\theta$ (degrees)",
ylabel = "Bistatic RCS (dBsm)",
yscale = 'linear',
xscale = 'linear'
)
ax.grid(visible = True, which = 'both', axis = 'both', linewidth = 0.4)
plt.legend(loc='best', prop={'size': 14})
plt.tight_layout()


<Figure size 540x360 with 1 Axes>

<Figure size 540x360 with 1 Axes>

[ ]: