tutorial with cylindrical coordinates and planewaves

[stevengj]

It would be nice to show an example using the cylindrical-coordinates feature for an incident planewave, illustrating that you have to decompose the field into exp(imφ) components. For example, scattering from a finite-height dielectric cylinder (if you google "scattering finite height dielectric cylinder" you can find lots of examples). The cylindrical results can be compared to a corresponding "brute force" 3d simulation.

If the incident planewave is on-axis, it is easy: just a superposition of m=±1. This case alone would be nice to show.

More generally, for an off-axis planewave, you have to use the Jacobi–Anger expansion to express the incident field as an infinite series of m's. If the wavelength is not too small compared to the scatterer diameter, the results will converge quickly, but you will still have to do a set of (embarrassingly parallel) simulations for different m's and sum the results.

[stevengj]

For now I would use an Hp source (H in the φ direction) that extends all the way to the edge of the cell (inside the PML). If you used an Ep source, you would switch to a magnetic boundary condition in the r direction (so that Ep can be nonzero at the edge)

[stevengj]

Note that a linearly polarized source is m=±1 because it is the superposition of left- and right-circularly polarized waves.

It is not m=0, because m=0 corresponds to a field pattern that is invariant under rotations. A linear polarization is not invariant, because if you rotate it by 180°, it flips sign (corresponding to m=±1). An m=0 wave would be something like a TE01 or TM01 mode, where either the electric or magnetic field looks like this:

[villadsegede]

Hello, I would like to implement plane wave propagation in circular coordinates (and would be happy if it will also contribute to your example database), but I am not sure I follow your thoughts completely:

Say, for a plane x-polarised TEM^z wave, it can be expressed in polar coordinates as

And from Faraday's law, we would also have both rho and a phi component for the H-field.

I am therefore a bit unsure about how you would generate a TEM^z wave just by exciting phi-polarised sources for m=+/-1, or is it me who've misunderstood something?

[stevengj]

@villadsegede, you're right, if you want a linearly polarized planewave I guess you need both an Er and Ep source, 90° out of phase (i.e. Ep amplitude multiplied by ±i depending on the sign of m=±1).

[stevengj]

Simple test: put in the above source in vacuum, and check that it produces a (circularly polarized) planewave. (If you add the fields from m=±1 simulations, then you should get the xz plane of a linearly polarized wave.)

To compute things like the scattered power, I think you can just compute the powers for m=±1 separately and add them, because the cross terms in E*×H cancel between different m's when you integrate over φ.

[oskooi]

Two sources for Er and Ep which are 90° out-of-phase for m=±1 based on the formula above does indeed produce a planewave propagating in the z direction as shown below.

Note the artifacts for Ep at r=0.

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

resolution = 25

dpml = 1.0
pml_layers = [mp.PML(thickness=dpml)]

dair_r = 3.0
dair_z = 5.0
sr = dair_r+dpml
sz = dpml+dair_z+dpml

cell_size = mp.Vector3(sr,0,sz)

fcen = 1.0

sources = [mp.Source(mp.ContinuousSource(fcen,fwidth=0.2*fcen,is_integrated=True),
                     component=mp.Er,
                     center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
                     size=mp.Vector3(sr)),
           mp.Source(mp.ContinuousSource(fcen,fwidth=0.2*fcen,is_integrated=True),
                     component=mp.Ep,
                     center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
                     size=mp.Vector3(sr),
                     amplitude=-1j)]

sim = mp.Simulation(cell_size=cell_size,
                    boundary_layers=pml_layers,
                    resolution=resolution,
                    sources=sources,
                    dimensions=mp.CYLINDRICAL,
                    m=-1)

sim.run(until=50)

nonpml_vol = mp.Volume(center=mp.Vector3(0.5*dair_r),size=mp.Vector3(dair_r,0,dair_z))
er_data = sim.get_array(component=mp.Er,vol=nonpml_vol)
ep_data = sim.get_array(component=mp.Ep,vol=nonpml_vol)

## get_array_metadata cannot (yet) be applied to cylindrical coordinates                                                                                   
r = np.linspace(0,dair_r,er_data.shape[1])
z = np.linspace(-0.5*dair_z,0.5*dair_z,er_data.shape[0])

plt.figure(dpi=150)

plt.subplot(1,2,1)
plt.pcolormesh(r,z,np.real(er_data), cmap='RdBu', shading='gouraud')
plt.xlabel(r'$r$')
plt.ylabel(r'$z$')
plt.title(r'$E_r$')
plt.axis('scaled')

plt.subplot(1,2,2)
plt.pcolormesh(r,z,np.real(ep_data), cmap='RdBu', shading='gouraud')
plt.xlabel(r'$r$')
plt.ylabel(r'$z$')
plt.title(r'$E_{\phi}$')
plt.axis('scaled')

plt.subplots_adjust(hspace=0.3)
plt.tight_layout()
plt.show()

Next step is to use this planewave setup to compute the scattering cross section of a finite-height dielectric cylinder and validate the results by comparing to the same simulation in 3d Cartesian coordinates.

[oskooi]

The scattering cross section of a finite-height dielectric cylinder (adapted from Tutorial/Basics/Mie Scattering of a Lossless Dielectric Sphere) computed using cylindrical and 3d Cartesian coordinates is shown below for two different cases. Results show good agreement. As @stevengj mentioned above, the scattered power for m=±1 can be computed separately and simply summed without worrying about the cross terms.

It would be good to generalize this calculation for an input planewave at an arbitrary angle.

Case 1: r=1.1, h=3.0

Case 2: r=0.5, h=2.1

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

r = 1.1  # radius of cylinder                                                                                                                                          
h = 3.0  # height of cylinder                                                                                                                                          

wvl_min = 2*np.pi*r/10
wvl_max = 2*np.pi*r/2

frq_min = 1/wvl_max
frq_max = 1/wvl_min
frq_cen = 0.5*(frq_min+frq_max)
dfrq = frq_max-frq_min
nfrq = 100

## at least 8 pixels per smallest wavelength, i.e. np.floor(8/wvl_min)                                                                                                 
resolution = 20

dpml = 0.5*wvl_max
dair = 1.0*wvl_max

pml_layers = [mp.PML(thickness=dpml)]

n_cyl = 2.0

def cross_sect_cyl():
    sr = r+dair+dpml
    sz = dpml+dair+h+dair+dpml
    cell_size = mp.Vector3(sr,0,sz)

    sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
                         component=mp.Er,
                         center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
                         size=mp.Vector3(sr)),
               mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
                         component=mp.Ep,
                         center=mp.Vector3(0.5*sr,0,-0.5*sz+dpml),
                         size=mp.Vector3(sr),
                         amplitude=-1j)]

    sim = mp.Simulation(cell_size=cell_size,
                        boundary_layers=pml_layers,
                        resolution=resolution,
                        sources=sources,
                        dimensions=mp.CYLINDRICAL,
                        m=-1)

    box_z1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(0.5*r,0,-0.5*h),size=mp.Vector3(r)))
    box_z2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(0.5*r,0,+0.5*h),size=mp.Vector3(r)))
    box_r  = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(r),size=mp.Vector3(z=h)))

    sim.run(until_after_sources=10)

    freqs = mp.get_flux_freqs(box_z1)
    box_z1_data = sim.get_flux_data(box_z1)
    box_z2_data = sim.get_flux_data(box_z2)
    box_r_data  = sim.get_flux_data(box_r)

    box_z1_flux0 = mp.get_fluxes(box_z1)

    sim.reset_meep()

    geometry = [mp.Block(material=mp.Medium(index=n_cyl),
                         center=mp.Vector3(0.5*r),
                         size=mp.Vector3(r,0,h))]

    sim = mp.Simulation(cell_size=cell_size,
                        geometry=geometry,
                        boundary_layers=pml_layers,
                        resolution=resolution,
                        sources=sources,
                        dimensions=mp.CYLINDRICAL,
                        m=-1)

    box_z1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(0.5*r,0,-0.5*h),size=mp.Vector3(r)))
    box_z2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(0.5*r,0,+0.5*h),size=mp.Vector3(r)))
    box_r  = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(r),size=mp.Vector3(z=h)))

    sim.load_minus_flux_data(box_z1, box_z1_data)
    sim.load_minus_flux_data(box_z2, box_z2_data)
    sim.load_minus_flux_data(box_r,  box_r_data)

    sim.run(until_after_sources=100)

    box_z1_flux = mp.get_fluxes(box_z1)
    box_z2_flux = mp.get_fluxes(box_z2)
    box_r_flux  = mp.get_fluxes(box_r)

    scatt_flux = np.asarray(box_z1_flux)-np.asarray(box_z2_flux)-np.asarray(box_r_flux)
    intensity = np.asarray(box_z1_flux0)/(0.5*np.pi*r**2)

    ## multiply by 2 for m=+1 circular polarization                                                                                                                    
    scatt_cross_section = 2*np.divide(-scatt_flux,intensity)

    return freqs, scatt_cross_section

def cross_sect_3d():
    sx = dpml+dair+h+dair+dpml
    syz = 2*(dpml+dair+r)
    cell_size = mp.Vector3(sx,syz,syz)

    # is_integrated=True necessary for any planewave source extending into PML                                                                                         
    sources = [mp.Source(mp.GaussianSource(frq_cen,fwidth=dfrq,is_integrated=True),
                         center=mp.Vector3(-0.5*sx+dpml),
                         size=mp.Vector3(0,syz,syz),
                         component=mp.Ez)]

    symmetries = [mp.Mirror(mp.Y),
                  mp.Mirror(mp.Z,phase=-1)]

    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell_size,
                        boundary_layers=pml_layers,
                        sources=sources,
                        k_point=mp.Vector3(),
                        symmetries=symmetries)

    box_x1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=-0.5*h),size=mp.Vector3(0,2*r,2*r)))
    box_x2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=+0.5*h),size=mp.Vector3(0,2*r,2*r)))
    box_y1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=-r),size=mp.Vector3(h,0,2*r)))
    box_y2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=+r),size=mp.Vector3(h,0,2*r)))
    box_z1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=-r),size=mp.Vector3(h,2*r,0)))
    box_z2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=+r),size=mp.Vector3(h,2*r,0)))

    sim.run(until_after_sources=10)

    freqs = mp.get_flux_freqs(box_x1)
    box_x1_data = sim.get_flux_data(box_x1)
    box_x2_data = sim.get_flux_data(box_x2)
    box_y1_data = sim.get_flux_data(box_y1)
    box_y2_data = sim.get_flux_data(box_y2)
    box_z1_data = sim.get_flux_data(box_z1)
    box_z2_data = sim.get_flux_data(box_z2)

    box_x1_flux0 = mp.get_fluxes(box_x1)

    sim.reset_meep()

    geometry = [mp.Cylinder(material=mp.Medium(index=n_cyl),
                            axis=mp.Vector3(1),
                            center=mp.Vector3(),
                            radius=r,
                            height=h)]

    sim = mp.Simulation(resolution=resolution,
                        cell_size=cell_size,
                        boundary_layers=pml_layers,
                        sources=sources,
                        k_point=mp.Vector3(),
                        symmetries=symmetries,
                        geometry=geometry)

    box_x1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=-0.5*h),size=mp.Vector3(0,2*r,2*r)))
    box_x2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(x=+0.5*h),size=mp.Vector3(0,2*r,2*r)))
    box_y1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=-r),size=mp.Vector3(h,0,2*r)))
    box_y2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(y=+r),size=mp.Vector3(h,0,2*r)))
    box_z1 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=-r),size=mp.Vector3(h,2*r,0)))
    box_z2 = sim.add_flux(frq_cen, dfrq, nfrq, mp.FluxRegion(center=mp.Vector3(z=+r),size=mp.Vector3(h,2*r,0)))

    sim.load_minus_flux_data(box_x1, box_x1_data)
    sim.load_minus_flux_data(box_x2, box_x2_data)
    sim.load_minus_flux_data(box_y1, box_y1_data)
    sim.load_minus_flux_data(box_y2, box_y2_data)
    sim.load_minus_flux_data(box_z1, box_z1_data)
    sim.load_minus_flux_data(box_z2, box_z2_data)

    sim.run(until_after_sources=100)

    box_x1_flux = mp.get_fluxes(box_x1)
    box_x2_flux = mp.get_fluxes(box_x2)
    box_y1_flux = mp.get_fluxes(box_y1)
    box_y2_flux = mp.get_fluxes(box_y2)
    box_z1_flux = mp.get_fluxes(box_z1)
    box_z2_flux = mp.get_fluxes(box_z2)

    scatt_flux = np.asarray(box_x1_flux)-np.asarray(box_x2_flux)+np.asarray(box_y1_flux)-np.asarray(box_y2_flux)+np.asarray(box_z1_flux)-np.asarray(box_z2_flux)
    intensity = np.asarray(box_x1_flux0)/(2*r)**2
    scatt_cross_section = np.divide(-scatt_flux,intensity)

    return freqs, scatt_cross_section

freqs, cs_cyl = cross_sect_cyl()
freqs, cs_3d = cross_sect_3d()

if mp.am_master():
    plt.figure(dpi=150)
    plt.loglog(2*np.pi*r*np.asarray(freqs),cs_cyl,'bo-',label='cyl')
    plt.loglog(2*np.pi*r*np.asarray(freqs),cs_3d,'ro-',label='3d')
    plt.grid(True,which="both",ls="-")
    plt.xlabel('(cylinder circumference)/wavelength, 2πr/λ')
    plt.ylabel('scattering cross section, σ')
    plt.legend(loc='upper right')
    plt.title('Scattering Cross Section a Lossless Dielectric Cylinder')
    plt.tight_layout()
    plt.savefig("cyl_scattering_cross_section_r{}_h{}.png".format(r,h))