Scattering from an infinite cylinder - Revision history

AndrejKP: Fix the boundary condition arising from H fields (no epsilon there) and the incident wave form. Add solution reference.

2026-01-22T16:24:37Z

Fix the boundary condition arising from H fields (no epsilon there) and the incident wave form. Add solution reference.

AndrejKP: Fix the boundary condition arising from H fields (no epsilon there) and the incident wave form. Add solution reference.

2026-01-22T16:24:10Z

Fix the boundary condition arising from H fields (no epsilon there) and the incident wave form. Add solution reference.

AndrejKP at 16:22, 22 January 2026

2026-01-22T16:22:54Z

VagifO: Corrected expression for gamma

2025-10-03T15:28:50Z

Corrected expression for gamma

BlazS: /* Case */

2020-09-24T08:19:59Z

‎Case

BlazS: /* Case */

2020-09-24T08:17:52Z

‎Case

BlazS: /* Case */

2020-09-24T07:54:00Z

‎Case

BlazS at 07:53, 24 September 2020

2020-09-24T07:53:06Z

BlazS at 07:40, 24 September 2020

2020-09-24T07:40:25Z

BlazS at 07:38, 24 September 2020

2020-09-24T07:38:23Z

@@ Line 16: / Line 16: @@
 \begin{equation}
 E_z^{int} - E_z^{s} =E_z^{inc}  \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
-\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n}  \qquad \text{on} \quad  \partial D.
+\dpar{E_z^{int}}{n} - \dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n}  \qquad \text{on} \quad  \partial D.
 \end{equation}
-The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.
+The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{-ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.
 <syntaxhighlight lang="c++" line>
@@ Line 68: / Line 68: @@
 ==Analytical solution==
-The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
+The solution to this problem is sought in terms of sums of cylindrical harmonics<ref name="analytical solution">M. Kerker , The Scattering of Light and Other Electromagnetic Radiation, 1969.</ref>. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
 \begin{equation}
 \mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
@@ Line 115: / Line 115: @@
 [[File:FDFD_2D_dx_conv2.png|800px]]

@@ Line 16: / Line 16: @@
 \begin{equation}
 E_z^{int} - E_z^{s} =E_z^{inc}  \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
-\dpar{E_z^{int}}{n} - \dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n}  \qquad \text{on} \quad  \partial D.
+\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n}  \qquad \text{on} \quad  \partial D.
 \end{equation}
-The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{-ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.
+The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.
 <syntaxhighlight lang="c++" line>
@@ Line 68: / Line 68: @@
 ==Analytical solution==
-The solution to this problem is sought in terms of sums of cylindrical harmonics<ref name="analytical solution">M. Kerker , The Scattering of Light and Other Electromagnetic Radiation, 1969.</ref>. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
+The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
 \begin{equation}
 \mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
@@ Line 115: / Line 115: @@
 [[File:FDFD_2D_dx_conv2.png|800px]]

@@ Line 16: / Line 16: @@
 \begin{equation}
 E_z^{int} - E_z^{s} =E_z^{inc}  \qquad \text{on} \quad \partial D \label{eq:BC1}, \qquad \text{and} \qquad
-\dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n}  \qquad \text{on} \quad  \partial D.
+\dpar{E_z^{int}}{n} - \dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n}  \qquad \text{on} \quad  \partial D.
 \end{equation}
-The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.
+The incident field $E_z^{inc}$ has an analytical form in our case a plane wave $e^{-ikx}$, since it has a nice closed form we know its derivatives and the boundary conditions are nicely expressed. Around the computational domain we use a PML to simulate an infinite domain.
 <syntaxhighlight lang="c++" line>
@@ Line 68: / Line 68: @@
 ==Analytical solution==
-The solution to this problem is sought in terms of sums of cylindrical harmonics. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
+The solution to this problem is sought in terms of sums of cylindrical harmonics<ref name="analytical solution">M. Kerker , The Scattering of Light and Other Electromagnetic Radiation, 1969.</ref>. The incident wave is a plane wave $\b E_i = \hat e_z E_0 e^{-ikx}$ and can be decomposed into Bessel functions of the first kind as
 \begin{equation}
 \mathbf{E}_{i}=\widehat{z} E_{0} e^{-i k_{0} x}= \widehat{z} E_{0} \sum_{n=-\infty}^{\infty} i^{-n} J_{n}\left(k_{0} \rho\right) e^{i n \phi}.
@@ Line 115: / Line 115: @@
 [[File:FDFD_2D_dx_conv2.png|800px]]

@@ Line 94: / Line 94: @@
 where
 \begin{equation}
-\Gamma_{h}=\frac{Z_{0}}{Z_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
+\Gamma_{h}=\frac{k_{0}}{k_{1}} \frac{J_{n}\left(k_{1} a\right)}{J_{n}^{\prime}\left(k_{1} a\right)}
 \end{equation}
 and

@@ Line 22: / Line 22: @@
 <syntaxhighlight lang="c++" line>
 	// Construct system
-	double l_in = epsr*std::pow(omega/c0, 2.0);
+double l_in = epsr*std::pow(omega/c0, 2.0);
-	double l_out = std::pow(omega/c0, 2.0);
+double l_out = std::pow(omega/c0, 2.0);
-	for (int i : cyl_dom.interior()){
+for (int i : cyl_dom.interior()){
-		op_inner.lap(i) + l_in*op_inner.value(i) = 0.0;
+	op_inner.lap(i) + l_in*op_inner.value(i) = 0.0;
-−
+}
-	for (int i : outer_inter){
+for (int i : outer_inter){
-		op_outer.lap(i) + l_out*op_outer.value(i) = 0.0;
+	op_outer.lap(i) + l_out*op_outer.value(i) = 0.0;
-−
+}
-	for (int c = 0; c < interface_c_idx.size(); ++c){
+for (int c = 0; c < interface_c_idx.size(); ++c){
-		int i = interface_c_idx[c]; // cyl index
+	int i = interface_c_idx[c]; // cyl index
-		int j = interface_o_idx[c]; // outer index
+	int j = interface_o_idx[c]; // outer index
-		Vec2d pos = cyl_dom.pos(i);
+	Vec2d pos = cyl_dom.pos(i);
-		double x = pos[0];
+	double x = pos[0];
-		double y = pos[1];
+	double y = pos[1];
-		// calculate normal derivative of the source funciton at the i-th node
+	// calculate normal derivative of the source funciton at the i-th node
-		// get normal, must point out of the cylinder
+	// get normal, must point out of the cylinder
-		Vec2d normal = cyl_dom.normal(i);
+	Vec2d normal = cyl_dom.normal(i);
-		std::complex<double> incident = std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
+	std::complex<double> incident = std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
-		std::complex<double> din_dx = 1.0_i*k0*std::cos(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
+	std::complex<double> din_dx =
-		std::complex<double> din_dy = 1.0_i*k0*std::sin(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
+.0_i*k0*std::cos(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0)));
-		std::complex<double> dui_dn = normal[0]*din_dx + normal[1]*din_dy;
+	std::complex<double> din_dy =
-		// continuity of the fields, and the incident field
+	// continuity of the fields, and the incident field
-		op_inner.value(i) + (-1)*op_outer.value(j, Nouter + i) = incident;
+	op_inner.value(i) + (-1)*op_outer.value(j, Nouter + i) = incident;
-		// continuity of derivatives
+	// continuity of derivatives
-		op_outer.neumann(j, outer_dom.normal(j)) + (1/epsr)*op_inner.neumann(i, cyl_dom.normal(i), j - Nouter)
+	op_outer.neumann(j, outer_dom.normal(j))
-				= dui_dn;
+              + (1/epsr)*op_inner.neumann(i, cyl_dom.normal(i), j - Nouter)
+			= dui_dn; // todo, change this from Neumanns to derivatives!
-	// SC-PML - where either x or y is constant
-	for (int i : PML){
+// SC-PML - where either x or y is constant
-.0/(sw[i]*sw[i])*op_outer.lap(i) + l_out*op_outer.value(i) +
+for (int i : PML){
-		((-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 0, i)*op_outer.der1(i, 0)
+.0/(sw[i]*sw[i])*op_outer.lap(i) + l_out*op_outer.value(i) +
-			+(-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 1, i)*op_outer.der1(i, 1)) = 0.0;
+	((-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 0, i)*op_outer.der1(i, 0)
+		+(-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 1, i)*op_outer.der1(i, 1)) = 0.0;
-	for (int i: outer_bnd){
-		op_outer.value(i) = 0.0;
+for (int i: outer_bnd){
+	op_outer.value(i) = 0.0;
 </syntaxhighlight>

@@ Line 9: / Line 9: @@
 As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
 \begin{align}
-\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{and} \quad \Omega \label{eq:inner} \\
+\nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{in} \quad D \label{eq:inner} \\
-\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s  = 0 \qquad & \text{outside} \quad \Omega
+\nabla^2 E_z^s + \frac{\omega^2}{c_0^2}E_z^s  = 0 \qquad & \text{outside} \quad D
- \label{eq:outer}
+\label{eq:outer}
 \end{align}
-The field $v$ is the total field inside, $u^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
+The field $E_z^{int}$ is the total field inside, $E_z^s$ is only the scattered part of $E_z$ Outside the cylinder. The boundary conditions are
 \begin{equation}
-E_z^{int} - E_z^{s} =u^{inc}  \qquad \text{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{and} \qquad
+E_z^{int} - E_z^{s} =E_z^{inc}  \qquad \text{on} \quad \partial \Omega \label{eq:BC1}, \qquad \text{and} \qquad
 \dpar{E_z^{int}}{n} - \frac{1}{\eps_r}\dpar{E_z^{s}}{n} = \dpar{E_z^{inc}}{n}  \qquad \text{on} \quad  \partial \Omega.
 \end{equation}

@@ Line 7: / Line 7: @@
 \newcommand{\ddpar}[2]{\frac{\partial^2 #1}{\partial #2^2}}
 \def\doubleunderline#1{\underline{\underline{#1}}}</math>
-As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. Zgled je pomemben iz dveh razlogov, prvi\v c, ogledali si bomo kako vpadno polje zapi\v semo v rob med materialom in prostim prostorom in lahko kot take vire analiti\v cno predstavimo\footnote{Tu to ne pride posebno prav, ker imamo opravka z ravnim valom, v primerih ko imamo to\v ckaste izvire pa je zelo prakti\v cno.}. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
+As a simple two dimensional case let us consider a plane wave scattering on a dielectric cylinder. This case is important for two reasons, firstly, it shows how to ''encode'' the incident wave into the boundary between the cylinder and free space and thus calculate only with the scattered field outside. It shows how to use analytic expressions to represent sources with known analytical solutions, this is especially useful when dealing with point sources. This case is important for two reasons. Firstly we will see how we can describe the solution in terms of scattered fields only, by including the incident wave in the boundary conditions. Secondly the problem has a relatively simple analytical solution in terms of Hankel funitons and can be used as a benchmark case.The problem is defined in two domains, inside and outside the cylinder, we are solving the following system of two PDE that are coupled on the inside boundary
 \begin{align}
 \nabla^2 E_z^{int} +\eps_r \frac{\omega^2}{c_0^2}\thinspace E_z^{int} = 0 \qquad &\text{and} \quad \Omega \label{eq:inner} \\