3D Near-to-far-field

[Snens98]

Hello, for my bachelor thesis, I try to implement a Near-to-Far-Field Transformation in 3D using Meep. It works as expected in 2D. My work is based on the "Near to Far Field Spectra" tutorial from the Meep documentation. For testing in 3D, I used a sphere and a plane wave, and I try to visualize the results. However, I am not getting the expected outcomes according to mie-scattering. Even adding or changing geometry doesn't fundamentally change the result, and I'm not sure why.

Do any of you have experience with a Near-to-Far-Field Transformation in 3D and can help me?

In the following I can show you my basic sim-settings:

cell = mp.Vector3(6, 6, 6)
resolution = 10
simTime = 10
pml_thickness = 0.5
pml_layers = [mp.PML(thickness=pml_thickness)]
wavelength = 0.457 
geometry3DSphere = [
    mp.Sphere(material = mp.Medium(index=1.52), radius = 0.5, center = mp.Vector3(0, 0, 0))
]

sim = mp.Simulation(
        cell_size=cell,
        boundary_layers=pml_layers,
        geometry=Geometry.geometry3DSphere,
        resolution=resolution,

I use a ContinuousSource wave with the parameter:
size = mp.Vector3(0, cell.y, cell.z), center=wavePos, end_time=3, width=0.5, slowness=1

For the nearfield detection I define near2far_regions:

sx, sy, sz = cell.x-pml_thickness*2-1, cell.y-pml_thickness*2-1, cell.z-pml_thickness*2-1 

"""+x-Surface"""
near2far_regions = [
    mp.Near2FarRegion(center=mp.Vector3(sx/2,0,0)
                      size=mp.Vector3(0,sy,sz)
                      direction=mp.X
                      weight=+1),
"""-x-Surface"""
    mp.Near2FarRegion(center=mp.Vector3(-sx/2, 0, 0),
                      size=mp.Vector3(0, sy, sz),
                      direction=mp.X,
                      weight=-1),
"""+y-Surface"""
    mp.Near2FarRegion(center=mp.Vector3(0, sy/2, 0),
                      size=mp.Vector3(sx, 0, sz),
                      direction=mp.Y,
                      weight=+1),
"""-y-Surface"""
    mp.Near2FarRegion(center=mp.Vector3(0, -sy/2, 0),
                      size=mp.Vector3(sx, 0, sz),
                      direction=mp.Y,
                      weight=-1),
"""+z-Surface"""
    mp.Near2FarRegion(center=mp.Vector3(0, 0, sz/2),
                      size=mp.Vector3(sx, sy, 0),
                      direction=mp.Z,
                      weight=+1),
"""-z-Surface"""
    mp.Near2FarRegion(center=mp.Vector3(0, 0, -sz/2),
                      size=mp.Vector3(sx, sy, 0),
                      direction=mp.Z,
                     weight=-1)
]

nearfield_box3D = sim.add_near2far(
    fcen,           
    0,              
    1,              
    *near2far_regions  
)

For the understanding I can show you a simple image. This is how I define the nearfield detection box, the geometry and the wave by the blue arrow:

So then I run the simulation:

sim.run(mp.at_every(timestep, mp.output_efield_z), until=simTime)

For the calculation of the Poynting Vector (Px, Py, Pz) for each theta and phi I use get_farfield to calculate the far-field at each defined point on the sphere. I don't know if this is correct:

def calcPoyntingVec(thetaPoints, phiPoints, nearFieldBox, radius):

    print("--- Calc Poynting Vector ---")

    r = radius           
    field_vectors = []

    for theta in thetaPoints:
        print(f"-- {theta} -- ")
        for phi in phiPoints:
            
            x = r * np.sin(theta) * np.cos(phi)
            y = r * np.sin(theta) * np.sin(phi)
            z = r * np.cos(theta)

            ff = sim.get_farfield(
                nearFieldBox,
                mp.Vector3(x, y, z)
            )

            field_vectors.append(ff)

    field_vectors = np.transpose(field_vectors)
    Ex, Ey, Ez, Hx, Hy, Hz = field_vectors
    # Ex, Ey, Ez, Hx, Hy, Hz = np.transpose(field_vectors, axes=(2, 0, 1))
    
    print(f"-- Get Ex, Ey, Ez, Hx, Hy, Hz  -- ")

    # P = E* x H
    Ex, Ey, Ez = np.conj([Ex, Ey, Ez])
    Px = np.real(Ey*Hz-Ez*Hy)
    Py = np.real(Ez*Hx-Ex*Hz)
    Pz = np.real(Ex*Hy-Ey*Hx)

    Pr = np.sqrt(np.square(Px) + np.square(Py) + np.square(Pz))
    print(f"-- Pr calculated -- ")
    
    saveIntensityToFile(thetaPoints, phiPoints, Pr, f"{directoryImg}/poynting_intensity.csv")

    return Pr

Pr corresponds to the absolut of the Poynting-Vector and consists of phi x theta elements. That basically mena for each angel theta and phi, we have a certain value of the Poynting-Vector. In the following we go back drom the spherical coordinates to a 3D cartesian grad X, Y, Z. To calculate the position of the scattered field using the Poyning-Vector in X, Y and Z we use this command: x = Pr[index] * np.sin(theta) * np.cos(phi), y and z equivalent. Here, the radius is replaced by the absolut eof the Poyntin-Vector for each angle theta and phi. At the end every position in X, Y, Z is displayd in a 3D projection.

radius = 100/fcen
thetaPoints = np.linspace(0, np.pi, 100)  # 0 bis Pi für Theta
phiPoints = np.linspace(0, 2*np.pi, 100)  # 0 bis 2*Pi für Phi
Pr = calcPoyntingVec(thetaPoints, phiPoints, nearfield_box3D, radius)

Px, Py, Pz = [], [], []

index = 0
for phi in thetaPoints:
    for theta in phiPoints:
        x = Pr[index] * np.sin(theta) * np.cos(phi)
        y = Pr[index] * np.sin(theta) * np.sin(phi)
        z = Pr[index] * np.cos(theta)
        index += 1

        X.append(x)
        Y.append(y)
        Z.append(z)
        
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.scatter(X, Y, Z, c='r', marker='o', label="Punkte")

ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Poynting-Vektor in 3D")
plt.legend()

ax.view_init(elev=30, azim=30)
plt.savefig(f"{directoryImg}/3d_Data_Spherical_Scatter_png")

Here in the following is an image of my result for scattering on a sphere.

The other image represents the "correct" scatterin behaviour in 3D:

However, my scattered field in 3D does not look correct.

Do you know how to implement this on the correct way to get a 3D Near-to-Far-Field transform as expected for a sphere for example?

Thank you for your help.

[stevengj]

Near-to-far definitely works in 3d, and has been used many times.

If you just want the scattered field, maybe you're forgetting to subtract the incident field? Or you have some other bug … I don't have time to look through your code in detail.

[Snens98]

Thanks for your response. Do I still need to subtract the incident field, even though I’m starting the wave between the PML and the near-field box? If I understood correctly, the wave is only detected in one direction by the near-field measurement surface, right? So if I start the wave before the measurement box and after the PML layer, it should also work, correct? Or did I misunderstand something?

[stevengj]

I use a ContinuousSource wave with the parameter:

This is wrong. Use a Gaussian source. You want the field to have decayed away completely by the time you use the Fourier transform to compute the far fields.

[Snens98]

Thanks, but I still have a question for understanding:

If I use a ContinuousSource with, for example:

end_time=4
width=0.5
slowness=2

don’t I achieve approximately the same effect as with a GaussianSource?

GaussianSource = mp.Source(
    src=mp.GaussianSource(wavelength=wavelength, fwidth=0.2*fcen, is_integrated=True),
    center=mp.Vector3(4.0, 0, 0),
    component=mp.Ez,
    size=mp.Vector3(0, cell.y, cell.z)
)

Do I also need to run the simulation until the wave has decayed sufficiently, for example to 1e-6?

sim.run(mp.at_every(timestep, mp.output_efield_z), until_after_sources=mp.stop_when_dft_decayed(tol=1e-6, maximum_run_time=300))

One more small question:
Do I add the near-field measurement boxes at the start of the simulation to begin measuring immediately, or only after the wave has turned off?"

[stevengj]

No, it's not the same as a GaussianSource. With a CW source, the Fourier spectrum never converges as you keep running and running.

You should add the near-field boxes at the beginning of the simulation and begin measuring immediately.

[Snens98]

Hey Steven thank your for your quick response. I tryed it on the explained way you mentioned in the comment, however it leads to the same result. Do you have further experience with the implementation of the 3D near-to-far-field method using meep? How did you implement this? Because I think my attempt is close to the right solution however there must be a mistake I do not see.
Do you use a similar method by using the Poynting-Vactor or do you use another method to transform the near field into far field?
Thank you and best regards
Snens

[Snens98]

I am performing a near-to-far-field transformation and calculating the Poynting vector from the far-field data Ex, Ey, Ez, Hx, Hy, Hz.
When I visualize the individual field components Ex, Ey, Ez, Hx, Hy, Hz, I always get the same result. I dont why?

def calcPoyntingVec(thetaPoints, phiPoints, nearFieldBox, radius):

    print("--- Calc Poynting Vector ---")

    r = radius           
    field_vectors = []

    for theta in thetaPoints:
        print(f"-- {theta} -- ")
        for phi in phiPoints:
            
            x = r * np.sin(theta) * np.cos(phi)
            y = r * np.sin(theta) * np.sin(phi)
            z = r * np.cos(theta)

            ff = sim.get_farfield(
                nearFieldBox,
                mp.Vector3(x, y, z)
            )

            field_vectors.append(ff)
    field_vectors = np.transpose(field_vectors)
    print(field_vectors.shape)

    Ex, Ey, Ez, Hx, Hy, Hz = field_vectors
    
    Ex = np.array(Ex)
    Ey = np.array(Ey)
    Ez = np.array(Ez)

With geometry:

Without geometry: