Scattering cross-section calculation
Contents
Scattering cross-section calculation#
To run this notebook from your browser, click this link.
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
radius = 0.5
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
sphere = td.Structure(geometry=td.Sphere(center=center, radius=radius), medium=diel)
geometry = [sphere]
# define PML layers on all sides
boundary_spec = td.BoundarySpec.all_sides(boundary=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
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()
Let’s also create another monitor where the fields are automatically downsampled based on given sampling rates. This is achieved by supplying the interval_space field, and can be useful for further reducing the amount of data that needs to be downloaded from the server, and can still lead to accurate far fields, as shown below.
Here, we downsample by a factor of 2 along the x and y dimensions, and a factor of 3 along z.
[5]:
# create a volume monitor around the sphere
monitor_downsampled = td.FieldMonitor(
center=center,
size=[mon_size, mon_size, mon_size],
freqs=[f0],
name='near_field_downsampled',
interval_space=(2, 2, 3))
# append the surface monitors associated with the downsampled volume monitor
monitors += monitor_downsampled.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.
[6]:
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,
boundary_spec=boundary_spec
)
# 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,
boundary_spec=boundary_spec
)
Visualize Geometry#
Let’s take a look and make sure everything is defined properly in both simulations.
[7]:
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
sim.plot(y=0, ax=ax1);
sim_empty.plot(y=0, ax=ax2);
[16:58:34] INFO Auto meshing using wavelength 0.5000 defined from grid_spec.py:472 sources.
INFO Auto meshing using wavelength 0.5000 defined from grid_spec.py:472 sources.
<AxesSubplot:title={'center':'cross section at y=0.00'}, xlabel='x', ylabel='z'>
<Figure size 648x216 with 2 Axes>
Run Simulations#
Now we can run both simulations over time and measure the results
[8]:
# Run simulation
import tidy3d.web as web
sim_data = web.run(sim, task_name='sphereRCS', path='data/sphereRCS.hdf5')
sim_empty_data = web.run(sim_empty, task_name='sphereRCS_empty', path='data/sphereRCS_empty.hdf5')
# 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
INFO Auto meshing using wavelength 0.5000 defined from grid_spec.py:472 sources.
INFO Auto meshing using wavelength 0.5000 defined from grid_spec.py:472 sources.
INFO Auto meshing using wavelength 0.5000 defined from grid_spec.py:472 sources.
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.
To see the usefulness of downsampling the fields recorded on monitors, we’ll also run the near-to-far transformation with downsampled fields to compare the RCS.
[9]:
from tidy3d.plugins import Near2Far
n2f = Near2Far.from_surface_monitors(
sim_data=sim_data,
monitors=monitors[0:6], # only supply the non-downsampled surface monitors as sources
normal_dirs=['-','+','-','+','-','+'],
frequency=f0,
pts_per_wavelength=10
)
n2f_downsampled = Near2Far.from_surface_monitors(
sim_data=sim_data,
monitors=monitors[6:], # only supply the downsampled surface monitors as sources
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.
[10]:
print(f'origin at {n2f.origin}')
origin at (0.0, 0.0, 0.0)
[11]:
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)
# get the RCS for the original and downsampled case
RCS = n2f.radar_cross_section(thetas, phis).values
RCS_downsampled = n2f_downsampled.radar_cross_section(thetas, phis).values
[17:01:38] INFO Auto meshing using wavelength 0.5000 defined from grid_spec.py:472 sources.
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.
Notice that the downsampled monitors also yield fairly accurate results with less than an eighth of the data.
[12]:
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])
RCS_downsampled_phi0 = np.squeeze(RCS_downsampled[:,0])
RCS_downsampled_phi90 = np.squeeze(RCS_downsampled[:,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_data_phi0 = np.loadtxt(mie_filename_phi0, delimiter='\t', skiprows=2)
mie_theta_phi0 = np.squeeze(mie_data_phi0[:,[0]])
mie_phi0 = np.squeeze(mie_data_phi0[:,[1]])
mie_data_phi90 = np.loadtxt(mie_filename_phi90, delimiter='\t', skiprows=2)
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), '-k', label="$\\phi = 0$, Mie")
ax.plot(thetas, to_db(RCS_phi0), '--.b', label="$\\phi = 0$, near2far")
ax.plot(thetas, to_db(RCS_downsampled_phi0), '--r', label="$\\phi = 0$, near2far downsampled")
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), '-k', label="$\\phi = \\pi/2$, Mie")
ax.plot(thetas, to_db(RCS_phi90), '--.b', label="$\\phi = \\pi/2$, near2far")
ax.plot(thetas, to_db(RCS_downsampled_phi90), '--r', label="$\\phi = \\pi/2$, near2far downsampled")
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>
[ ]: