tutorial for oblique planewave in cylindrical coordinate

doc/docs/Python_Tutorials/Cylindrical_Coordinates.md

@@ -455,6 +455,158 @@ As shown below, the results for the scattering cross section computed using cyli
 ![](../images/cylinder_cross_section.png#center)
 
 
+Scattering of Sphere with Oblique Planewave
+-------------------------------------------
+
+It is also possible to launch an oblique incident planewave in cylindrical coordinate by decomposing the planewave $A_xe^{ik_xx+ik_yy}\hat{x} + A_ye^{ik_xx+ik_yy}\hat{y}$ into $\sum_m (J_r(r, m)\hat{r} + J_\phi(r, m)\hat{\phi})e^{im\phi}$ through [Jacobi-Anger expansion](https://en.wikipedia.org/wiki/Jacobi%E2%80%93Anger_expansion). The exact expressions of $J_r(r,m)$ and $J_\phi(r,m)$ are given [here](http://github.com/zlin-opt/axisym_meta3d_inverse_design/blob/master/Implementation_of_FDFD_with_Cylindrical_Coordinates.pdf) by Zin Lin. In the simplest case of normal incidence, $J_r(r,m)$ and $J_\phi(r,m)$ are nonzero only when $m = \pm 1$, as shown in the [previous tutorial](https://meep.readthedocs.io/en/latest/Python_Tutorials/Cylindrical_Coordinates/#scattering-cross-section-of-a-finite-dielectric-cylinder).
+
+Given the decomposition of planewave into the sum of different current sources at each $m$, we can run individual simulations at each $m$ with their corresponding source amplitudes and record the relevant physical quantities. For some quantities such as fields, linearity implies that we can simply sum the results from each simulations; for some other quantities such as flux, orthogonality implies cross terms will be zero, and we can again simply sum the results. Moreover, simulations
+at each $m$ values are embarrassingly parallel so they can be run simultaneously.
+
+We present an example below that calculates the scattered flux of a sphere. Because of the spherical symmetry, incidence at different angle should have identical results. We can thus use this feature to check our approach. Note that because of the axial symmetry in the cylindrical coordinates, we cannot distinguish different azimuthal angles but we can distinguish different polar angles. We thus simply choose our incidence to be of form $E_ye^{ik_xx}$, and we can vary the angle of incidence by varying $k_x$.
+
+On the other hand, because the source amplitudes $J_r(r,m)$ and $J_\phi(r,m)$ are generally not constant and extend to infinity, we used the principle of equivalence (for reference, see [Electromagnetic wave source condition](https://arxiv.org/pdf/1301.5366.pdf)) to create equivalent sources that are of finite sizes. Specifically, with the chosen incidence, the E fields in space are $E_ye^{ik_xx+ik_zz}$, and thus H fields can be computed by taking the curl; then Jacobi-Anger expansion can express the dependencies in $x$ and $y$ in terms of $m$ and $r$; afterwards, we created a box of sources surrounding the geometry and specify sources of amplitude $J = n \times H$ and $K = - n \times E$.
+
+Empirically, we found that the Courant factor has to scale as $1/(|m|+0.5)$ in cylindrical coordinate to maintain numerical stability. By default, Meep uses the same Courant factor but instead zeros out fields near axis for $|m| > 1 $. In this tutorial, we choose to scale the Courant factor accordingly and force Meep to use the actual fields near axis via `accurate_fields_near_cylorigin=True`.
+
+```py
+import numpy as np
+from scipy import special
+import meep as mp
+mp.verbosity(0)
+r = 0.6  # size of flux box
+cyl_r = 0.5 # radius of sphere
+h = 2 * r  # height/diameter of sphere
+
+wvl = 2 * np.pi * cyl_r / 4
+frq_cen = 1 / wvl
+dfrq = 0.2
+nfrq = 1
+resolution, mrange = 50, 5
+dpml = 0.5 * wvl
+dair = 1.0 * wvl
+pml_layers = [mp.PML(thickness=dpml)]
+sr = r + dair + dpml
+sz = dpml + dair + h + dair + dpml
+cell_size = mp.Vector3(sr, 0, sz)
+n_cyl = 2.0
+geometry = [mp.Sphere(material=mp.Medium(index=n_cyl), center=mp.Vector3(), radius=cyl_r)]
+
+k_cen = 2 * np.pi * frq_cen
+alpha_list = [0, np.pi/36, np.pi/24, np.pi/18, np.pi/12]
+alpha_range = len(alpha_list)
+
+
+src_size_tb = 2*r
+src_size_side = 3*r
+src_center_top = mp.Vector3(src_size_tb/2, 0, src_size_side/2)
+src_center_bottom = mp.Vector3(src_size_tb/2, 0, -src_size_side/2)
+src_center_side = mp.Vector3(src_size_tb, 0, 0)
+
+scatt_flux_m = np.zeros((alpha_range, mrange+1))
+for alpha_i in range(alpha_range):
+    alpha = alpha_list[alpha_i]
+    kxy, kz = k_cen*np.sin(alpha), k_cen * np.cos(alpha)
+    amp_side = lambda v3: np.exp(1j * kz*(v3.z+src_size_side/2))
+    phase_top = amp_side(src_center_top)
+
+    for cur_m in range(0, mrange+1):
+        if alpha!=0 or cur_m == 1:
+            coeff_p1 = 0.5 * (1j)**(cur_m+1)
+            coeff_m1 = 0.5 * (1j)**(cur_m-1)
+
+            src_cen = src_size_tb/2
+            Jpm = lambda v3: coeff_p1 * special.jv(cur_m+1, kxy * (v3.x+src_cen)) + coeff_m1 * special.jv(cur_m-1, kxy * (v3.x+src_cen))
+            Jrm = lambda v3: 1j * coeff_p1 * special.jv(cur_m+1, kxy * (v3.x+src_cen)) - 1j * coeff_m1 * special.jv(cur_m-1, kxy * (v3.x+src_cen))
+            Jside = (1j)**cur_m * special.jv(cur_m, kxy*src_size_tb) * kxy/k_cen
+
+            src_t  = mp.GaussianSource(frq_cen, fwidth=dfrq)
+            sourcesp = [
+                mp.Source(src_t,component=mp.Er, center=src_center_bottom,size=mp.Vector3(src_size_tb), amplitude = -kz/k_cen, amp_func = Jrm),
+                mp.Source(src_t,component=mp.Ep, center=src_center_bottom,size=mp.Vector3(src_size_tb), amplitude = -kz/k_cen, amp_func = Jpm),
+                mp.Source(src_t,component=mp.Hr, center=src_center_bottom,size=mp.Vector3(src_size_tb), amp_func = Jpm),
+                mp.Source(src_t,component=mp.Hp, center=src_center_bottom,size=mp.Vector3(src_size_tb), amplitude = -1, amp_func = Jrm),
+                mp.Source(src_t,component=mp.Er, center=src_center_top,size=mp.Vector3(src_size_tb), amplitude = phase_top*kz/k_cen, amp_func = Jrm),
+                mp.Source(src_t,component=mp.Ep, center=src_center_top,size=mp.Vector3(src_size_tb), amplitude = phase_top*kz/k_cen, amp_func = Jpm),
+                mp.Source(src_t,component=mp.Hr, center=src_center_top,size=mp.Vector3(src_size_tb), amplitude = -phase_top, amp_func = Jpm),
+                mp.Source(src_t,component=mp.Hp, center=src_center_top,size=mp.Vector3(src_size_tb), amplitude = phase_top, amp_func = Jrm),
+                mp.Source(src_t,component=mp.Ez, center=src_center_side,size=mp.Vector3(z=src_size_side), amplitude = -Jrm(src_center_top)*kz/k_cen, amp_func = amp_side),
+                mp.Source(src_t,component=mp.Hz, center=src_center_side,size=mp.Vector3(z=src_size_side), amplitude = Jpm(src_center_top), amp_func = amp_side),
+                mp.Source(src_t,component=mp.Ep, center=src_center_side,size=mp.Vector3(z=src_size_side), amplitude = Jside, amp_func = amp_side),
+            ]
+
+
+            sim = mp.Simulation(
+                cell_size=cell_size,
+                boundary_layers=pml_layers,
+                resolution=resolution,
+                sources=sourcesp,
+                dimensions=mp.CYLINDRICAL,
+                m=cur_m,
+                force_complex_fields = True,
+                accurate_fields_near_cylorigin=True,
+                Courant=min(0.5, 1/(abs(cur_m)+0.5)))
+
+            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()
+
+            sim = mp.Simulation(
+                cell_size=cell_size,
+                geometry=geometry,
+                boundary_layers=pml_layers,
+                resolution=resolution,
+                sources=sourcesp,
+                dimensions=mp.CYLINDRICAL,
+                m=cur_m,
+                force_complex_fields = True,
+                accurate_fields_near_cylorigin=True,
+                Courant=min(0.5, 1/(abs(cur_m)+0.5)))
+
+            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_m[alpha_i, cur_m] = box_z1_flux[0] - box_z2_flux[0] - box_r_flux[0]
+            sim.reset_meep()
+
+scatt_power_m = np.zeros((alpha_range, mrange+1))
+for i in range(mrange+1):
+    scatt_power_m[:,i] = - 2*np.sum(scatt_flux_m[:,0:(i+1)], axis=1) + scatt_flux_m[:,0]
+
+print(scatt_power_m)
+
+```
 
 Focusing Properties of a Binary-Phase Zone Plate
 ------------------------------------------------