Total field scatter field modelling using MEEP

[athulshaji]

Hello, I would like to implement a total field/scattered field (TF/SF) source in 2D using MEEP. How can this be done? I am looking to solve an RCS problem where an irregularly shaped object is illuminated by a source coming from a particular azimuthal angle.

Thank you, Athul Shaji

[smartalecH]

If you're after RCS and don't need to use a TF/SF implementation, then you can simply use a near-to-far transform (like the antenna tutorial) in conjunction with a normalization run (like the bend flux tutorial) and a diagonal eigenmode plane source (there's a tutorial for this too). I squashed the three tutorials together and created something like this:

First, we'll generate a calibration run over a simulation domain with plenty of PML, an eigenmode source to launch a plane-wave, and a near-to-far monitor. We'll run the simulation until the fields decay and we'll also save the fluxes to calibrate the same simulation domain, but with a PEC cylinder. This way we can remove the incident field from the total field and obtain the scattered field:

import meep as mp
import math
import numpy as np
import matplotlib.pyplot as plt

resolution = 50  
sxy = 10 # length of simulation domain excluding PMLs
dft_mon = 4 # length of DFT box
dpml = 2.0 # PML thickness 
rot_angle = np.radians(-45) # rotation angle of the source

cell = mp.Vector3(sxy+2*dpml,sxy+2*dpml,0)
pml_layers = [mp.PML(dpml)]

 # ------------------ Source parameters ----------------------- #

fcen = 1.0 # center frequency at 1 micron
df = 0.1 # bandwidth of gaussian source
k_point = mp.Vector3(fcen*1).rotate(mp.Vector3(z=1), rot_angle)
sources = [mp.EigenModeSource(src=mp.GaussianSource(fcen,fwidth=df),
                              center=mp.Vector3(x=-4),
                              size=mp.Vector3(y=(sxy + 2*dpml)),
                              direction=mp.AUTOMATIC if rot_angle == 0 else mp.NO_DIRECTION,
                              eig_kpoint=k_point,
                              eig_band=1,
                              eig_parity=mp.EVEN_Y+mp.ODD_Z if rot_angle == 0 else mp.ODD_Z,
                              eig_match_freq=True)]

 # ------------------ Formulate simulation -------------------- #

sim = mp.Simulation(cell_size=cell,
                    resolution=resolution,
                    sources=sources,
                    boundary_layers=pml_layers)

nearfield_box = sim.add_near2far(fcen, 0, 1,
                                 mp.Near2FarRegion(mp.Vector3(y=0.5*dft_mon), size=mp.Vector3(dft_mon)),
                                 mp.Near2FarRegion(mp.Vector3(y=-0.5*dft_mon), size=mp.Vector3(dft_mon), weight=-1),
                                 mp.Near2FarRegion(mp.Vector3(0.5*dft_mon), size=mp.Vector3(y=dft_mon)),
                                 mp.Near2FarRegion(mp.Vector3(-0.5*dft_mon), size=mp.Vector3(y=dft_mon), weight=-1))

sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(x=-4), 1e-4))

 # ------------------ Save flux data ------------------------- #

# for normalization run, save flux fields data
nearfield_refl_data = sim.get_near2far_data(nearfield_box)

sim.reset_meep()

Now that we have our calibration run done, we'll clear meep, add the infinite cylinder, run another simulation with the loaded calibration fluxes, and compute the near-to-far transform.

r = 1 # radius of cylinder
geometry = [mp.Cylinder(center=mp.Vector3(),radius = r, material=mp.metal, height=mp.inf)]
sim = mp.Simulation(cell_size=cell,
                    resolution=resolution,
                    sources=sources,
                    geometry = geometry,
                    boundary_layers=pml_layers)

nearfield_box = sim.add_near2far(fcen, 0, 1,
                                 mp.Near2FarRegion(mp.Vector3(y=0.5*dft_mon), size=mp.Vector3(dft_mon)),
                                 mp.Near2FarRegion(mp.Vector3(y=-0.5*dft_mon), size=mp.Vector3(dft_mon), weight=-1),
                                 mp.Near2FarRegion(mp.Vector3(0.5*dft_mon), size=mp.Vector3(y=dft_mon)),
                                 mp.Near2FarRegion(mp.Vector3(-0.5*dft_mon), size=mp.Vector3(y=dft_mon), weight=-1))

# for normal run, load negated fields to subtract incident from total fields
sim.load_minus_near2far_data(nearfield_box, nearfield_refl_data)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ez, mp.Vector3(x=-4), 1e-4))

We can use the resultant N2F data to look at the scattering characteristics of the PEC cylinder:

r = 1000/fcen      # half side length of far-field square box OR radius of far-field circle
res_ff = 1         # resolution of far fields (points/μm)
far_flux_box = (nearfield_box.flux(mp.Y, mp.Volume(center=mp.Vector3(y=r), size=mp.Vector3(2*r)), res_ff)[0]
               - nearfield_box.flux(mp.Y, mp.Volume(center=mp.Vector3(y=-r), size=mp.Vector3(2*r)), res_ff)[0]
               + nearfield_box.flux(mp.X, mp.Volume(center=mp.Vector3(r), size=mp.Vector3(y=2*r)), res_ff)[0]
               - nearfield_box.flux(mp.X, mp.Volume(center=mp.Vector3(-r), size=mp.Vector3(y=2*r)), res_ff)[0])

npts = 100         # number of points in [0,2*pi) range of angles
angles = 2*math.pi/npts*np.arange(npts)

E = np.zeros((npts,3),dtype=np.complex128)
H = np.zeros((npts,3),dtype=np.complex128)
for n in range(npts):
    ff = sim.get_farfield(nearfield_box,
                          mp.Vector3(r*math.cos(angles[n]),
                                     r*math.sin(angles[n])))
    E[n,:] = [np.conj(ff[j]) for j in range(3)]
    H[n,:] = [ff[j+3] for j in range(3)]

Px = np.real(np.multiply(E[:,1],H[:,2])-np.multiply(E[:,2],H[:,1]))
Py = np.real(np.multiply(E[:,2],H[:,0])-np.multiply(E[:,0],H[:,2]))
Pr = np.sqrt(np.square(Px)+np.square(Py))

far_flux_circle = np.sum(Pr)*2*np.pi*r/len(Pr)

plt.figure()
plt.plot(angles*180/np.pi,Pr/max(Pr),'b-')
plt.grid(True)
plt.xlabel('Angle (degrees)')
plt.ylabel('Scattering (a.u.)')
plt.savefig('angles.png')
plt.show()

You can run different polarizations by changing the eigenmode that you launch (TE vs TM). You can get fancy with your broadband source and get RCS data for multiple wavelengths using a single simulation.

Once you have a calibration run done, you can even use the adjoint solver to find a geometry that produces a desired RCS profile. This is super efficient, as it calculates the gradient of your cost function (you can define any arbitrary cost function with respect to your geometry) with just one extra simulation.

In practice, you would want to benchmark this with the analytic solution to make sure the calibration was done correctly (there are lots of papers that analytically scatter off of PEC cylinders, for example). As such there's no guarantee that the code above is bug-free...

@oskooi All of this might be a good tutorial/notebook to add to the docs.

[athulshaji]

Thank you.

[stevengj]

(Note that there is no need to compute the far-field at all to get the RCS profile — the scattered field in any direction and polarization is exactly given by relatively simple integrals of the local field pattern.)