This is an IJulia "notebook", in which we can combine code, results, graphics, formatted text, etcetera to describe a computational study.
In this case, the computation will be focused on solving the scalar-wave (Helmholtz) scattering problem $-\nabla^2 u = \omega^2/c^2 n^2(r) u$ for a set of $N$ concentric nested cylinders with (complex) refractive indices $n_j$, given some superposition of incoming cylindrical waves $H^{(2)}(k_{N+1}r) e^{im\theta -i\omega t}$, for $m\in\mathbb{Z}$ and $k_{N+1}=n_{N+1}\omega/c$ for the outermost medium $n_{N+1}$. We will choose natural units where $c=1$, for convenience, or equivalently specify the "wavenumber" $k=\omega/c$ instead of $\omega$.
See also: A. Ishimaru, Electromagnetic Wave Propagation, Radiation, and Scattering (Prentice Hall, Englewood Cliffs NJ, 1991); and J. J. Bowman et al., Electromagnetic and Acoustic Scattering by Simple Shapes (Wiley, New York, 1969).
Like Python, Julia is a dynamic language that allows you to type commands interactively, with a syntax that is superficially somewhat Matlab-like:
x = linspace(0,10,100)
100-element Array{Float64,1}: 0.0 0.10101 0.20202 0.30303 0.40404 0.505051 0.606061 0.707071 0.808081 0.909091 1.0101 1.11111 1.21212 ⋮ 8.88889 8.9899 9.09091 9.19192 9.29293 9.39394 9.49495 9.59596 9.69697 9.79798 9.89899 10.0
It has lots of built-in functions, similar to Matlab. For example, the Bessel functions $J_m(x)$ and $Y_m(x)$ are built-in:
y = besselj(4,x)
100-element Array{Float64,1}: 0.0 2.70961e-7 4.32874e-6 2.18584e-5 6.88368e-5 0.000167287 0.000344941 0.000634806 0.00107466 0.00170643 0.00257559 0.00373036 0.00522097 ⋮ -0.257128 -0.264824 -0.27028 -0.273466 -0.274369 -0.272996 -0.269371 -0.263537 -0.255555 -0.245506 -0.233484 -0.219603
Plotting is accomplished by external plotting packages. There are several to choose from. For example, PyPlot, is a wrapper around Python's Matplotlib.
Pkg.add("PyPlot") # this only needs to be done once on your machine, to install PyPlot
INFO: Nothing to be done INFO: METADATA is out-of-date — you may not have the latest version of PyPlot INFO: Use `Pkg.update()` to get the latest versions of your packages
using PyPlot # analogous to from PyPlot import * in Python.
Loading help data...
title(L"Bessel function $J_4(x)$") # We use L"..." instead of "..." if we have embedded LaTeX equations $...$ in the string
plot(x,y)
xlabel(L"$x$")
PyObject <matplotlib.text.Text object at 0x1207ac6d0>
Note that the IJulia notebook can contain graphics as well as code, results, and formatted text.
String concatenation in Julia is the *
operator, not +
:
"a" * "b"
"ab"
"a" + "b"
no method +(ASCIIString, ASCIIString) while loading In[7], in expression starting on line 1
However, we can define new functions for String object, and we can define a +
if we want. There are already lots of +
definitions:
methods(+)
import Base: + # tell Julia that we want to add new definitions to + from the Base module (all built-in functions)
+(x::String, y::String) = x * " " * y
"hello" + "world"
"hello world"
Let's make a data structure Scatterer
to describe $N$ layers of (complex) indices $n_j$ for $j=1,\ldots,N$, and an outermost (real) index $n_{N+1}$. The outer radius of index $n_j$ will have radius $R_j$.
# The Scatterer type is "parameterized" by a real type T. Typically, we will use T = Float64, for double-precision
# floating-point arithmetic (about 15 decimal digits), but maybe at some point we will try more accurate arithmetic types.
type Scatterer{T<:Real}
n::Vector{Complex{T}} # the complex indices of each layer 1:N+1, with n[N+1] being the outermost medium.
R::Vector{T} # R[j] is the radius of the interface between n[j] and n[j+1], for j=1:N
# define a "constructor" function to create an object, performing some simple type-checking, where n0
# is the outermost medium.
function Scatterer(n0::T, n::Vector{Complex{T}}, R::Vector{T})
if length(n) != length(R)
throw(ArgumentError("n and R must have the same length"))
end
if any(R .<= 0) || any(diff(R) .< 0)
throw(ArgumentError("R must be positive and monotonic nondecreasing"))
end
new([n, n0],R)
end
end
# For convenience, we also define a constructor that accepts arbitrary numeric types, converting
# everything to the same type T as needed. T is selected by "promoting" the constituent types
# to the "lowest common denominator".
function Scatterer{Tn<:Number,Tr<:Real}(n0::Real, n::Vector{Tn}, R::Vector{Tr})
T = promote_type(typeof(real(complex(one(promote_type(typeof(n0),Tn))))), Tr)
Scatterer{T}(convert(T, n0), copy!(Array(Complex{T},length(n)), n), copy!(Array(T,length(R)), R))
end
Scatterer{T<:Real} (constructor with 1 method)
s = Scatterer(3.0, [4,5,7], [8,9,11])
Scatterer{Float64}(Complex{Float64}[4.0+0.0im,5.0+0.0im,7.0+0.0im,3.0+0.0im],[8.0,9.0,11.0])
Scatterer(3.0, [4,5,7], [8,7,11])
ArgumentError("R must be positive and monotonic nondecreasing") while loading In[12], in expression starting on line 1 in Scatterer at In[10]:14 in Scatterer at In[10]:25
Julia includes a standard function length(A)
to get the length of an array and similar types. We will add a new method to this type in order to get the number of layers ($N$) in the scatterer.
import Base: length # the import statement is needed to add new methods to a built-in function (from the Base module)
length(s::Scatterer) = length(s.R)
length (generic function with 38 methods)
It is also convenient to define a show
method for Scatterer
objects, in order to print them in a nice format.
import Base: show
function show{T}(io::IO, s::Scatterer{T})
println(io, "Scatterer{$T} with ", length(s), " layers surrounded by n = ", real(s.n[end]), ":")
println(io, " n: ", join(s.n[1:end-1], ", "))
println(io, " R: ", join(s.R, ", "))
end
show (generic function with 88 methods)
Now, if we look at s
(defined above), it will print more nicely:
s
Scatterer{Float64} with 3 layers surrounded by n = 3.0: n: 4.0 + 0.0im, 5.0 + 0.0im, 7.0 + 0.0im R: 8.0, 9.0, 11.0
Julia defines hankelh1
and hankelh2
for the Hankel functions $H^{(1)}_m(x) = J_m(x) + iY_m(x)$ and $H^{(2)}_m(x) = J_m(x) - iY_m(x)$, respectively, which represent outgoing and incoming cylindrical waves solutions to a constant-coefficient Helmholtz equation, respectively (also denoted $H^\pm$):
x = linspace(2, 20)
plot(x, abs(hankelh1(3, x)), "r--")
plot(x, real(hankelh1(3, x)), "b-")
legend([L"$|H^{(1)}_3(x)|$", L"$\operatorname{Re} \left[ H^{(1)}_3(x) \right]$"])
xlabel(L"$x$")
grid()
We will also define functions hankel1d
, besseljd
, etcetera, for the derivatives of the Bessel functions, using the identity $Z_m'(x) = [Z_{m-1}(x) - Z_{m+1}(x)]/2$ for $Z \in \{J, Y, H^{(1)}, H^{(2)}\}$. Since this is the same formula repeated several times, we will use a clever feature of Julia called metaprogramming: Julia allows us to write a program to write a program, in this case a loop to construct and evaluate (via @eval
) several related function definitions:
for f in (:besselj, :bessely, :hankelh1, :hankelh2)
fd = symbol(string(f, "d"))
@eval $fd(m, x) = ($f(m-1,x) - $f(m+1,x))*0.5
end
Let's do a quick check to make sure our derivative is correct, by comparing it to a finite-difference approximation:
x = 0.3; m = 2; dx=1e-8
hankelh1d(m, x), (hankelh1(m, x+dx) - hankelh1(m, x-dx))/(2dx)
(0.07387973661267101 + 94.24085493796042im,0.0738797618568543 + 94.2408545867579im)
Now, we can write a function to compute the scattering amplitude, the amplitude $1+2V$ of the outgoing wave $(1+2V) H^{(1)}(k_{N+1}r)$ for a unit-amplitude incoming wave $S H^{(2)}(k_{N+1}r)$. (It is critically important to compute $V$ rather than $1+2V$, because $V\to 0$ for $m\to\infty$ and $1+2V$ is therefore subject to extreme roundoff errors.) We iteratively construct the amplitude by adding the effect of one layer at a time, starting with the scattering amplitude of the innermost interface:
function scat_amplitude(s::Scatterer, ω::Number, m::Integer)
# base case for inner layer
k1R = ω*s.n[1]*s.R[1]
k2R = ω*s.n[2]*s.R[1]
# handle the limiting case of infinite n_1, e.g. a perfect-metal interior
if isinf(abs(k1R))
V = -besselj(m,k2R) / hankelh1(m,k2R)
else
α = s.n[1]/s.n[2] * besseljd(m,k1R)/besselj(m,k1R)
V = (α*besselj(m,k2R) - besseljd(m,k2R)) / (hankelh1d(m,k2R) - α*hankelh1(m,k2R))
end
for j = 2:length(s)
k1R = ω*s.n[j]*s.R[j]
k2R = ω*s.n[j+1]*s.R[j]
α = s.n[j]/s.n[j+1] * (besseljd(m,k1R) + V*hankelh1d(m,k1R)) / (besselj(m,k1R) + V*hankelh1(m,k1R))
V = (α*besselj(m,k2R) - besseljd(m,k2R)) / (hankelh1d(m,k2R) - α*hankelh1(m,k2R))
end
return V
end
scat_amplitude (generic function with 1 method)
(Note that the alternative "transfer-matrix" approach, which constructs the transfer matrix M relating the amplitudes of the Bessel functions in the innermost layer to those in the outermost layer, is problematic because it can yield extremely ill-conditioned matrices.)
Let's do a quick check to make sure it gives some kind of result without encountering an error:
scat_amplitude(s, 0.3, 3)
-0.24349548255735684 - 0.4291916035193559im
For real indices $n$, we must have $|1+V_m| = 1$ to conserve energy (the scattering matrix is unitary), up to roundoff errors:
abs(1+2*scat_amplitude(s, 0.3, 3))
0.9999999999999984
As a consequence of symmetry, which makes $H_m(x) = (-1)^m H_{-m}(x)$, the scattering amplitudes are equal for $\pm m$:
scat_amplitude(s, 0.3, -3)
-0.24349548255735684 - 0.4291916035193559im
Another sanity check is that the scattering amplitudes should go to 0 as $m \to \infty$:
scat_amplitude(s, 0.3, 30)
2.4707060062781136e-37 + 2.164961838917851e-23im
To solve for the scattering of a planewave $e^{i(kx-\omega t)}$ ($k = n_{N+1}\omega/c$) off of the cylinder, we use the Jacobi-Anger identity to expand the planewave in cylindrical waves:
and hence the total (incident + scattered) wave solution (outside the scatterer) is:
where $V_m$ is the scattering amplitude for $H^{(2)}_m(kr)$ from above.
We also need only compute the series for non-negative $m$ because of the $V_{-m} = V_m$ symmetry mentioned above. Assuming $k$ is real, this gives:
Although the exact solution involves an infinite series, the terms converge exponentially to zero for $m \gg kR$, where $R$ is the outermost radius, so we can truncate the series once the terms fall below some desired tolerance tol
.
The exponential convergence results from the asymptotic expansion $J_m(x) \approx \frac{1}{\sqrt{2\pi m}} \left(\frac{e z}{2m}\right)^{m}$, and can be demonstrated with a simple plot:
semilogy(0:20, abs(besselj(0:20, 1)), "ro")
ms = linspace(1,20,200)
semilogy(ms, (e./(2ms)).^ms ./ sqrt(2π*ms), "b-")
legend([L"$|J_m(1)|$", "asymptotic form"])
xlabel(L"Bessel order $m$")
PyObject <matplotlib.text.Text object at 0x101224c50>
# Return (M,V), where M is an array of the m >= 0 and V is an array of the corresponding scattering amplitudes,
# computing enough amplitudes until the remaining terms contribute less than tol up to a radius r
function scat_amplitudes{T}(s::Scatterer{T}, ω::Number, r::Real=s.R[end], tol=sqrt(eps(T)))
M = [0]
V = [scat_amplitude(s, ω, 0)]
abs_kr = abs(s.n[end]*ω*r)
m = 1
while (e*abs_kr/(2m))^m ./ sqrt(2π*m) > tol
push!(M, m)
push!(V, scat_amplitude(s, ω, m))
m += 1
end
return (M,V)
end
scat_amplitudes (generic function with 4 methods)
Let's do a quick check to see how many $m$'s are needed for a typical structure, and plot the phase $\arg(1+2V_m)$ versus $m$ as well as $|V_m|$:
M,V = scat_amplitudes(s, 0.3)
subplot(1,2,1)
plot(M, angle(1+2V)/π, "ro-")
xlabel(L"$m$")
title(L"$\arg (1+2V_m) / \pi$")
grid()
subplot(1,2,2)
semilogy(M, abs(V), "ro-")
xlabel(L"$m$")
title(L"$|V_m|$")
grid()
Now, let's define a function for plotting the total wave solution outside the scatterer, returning a 2d array of the wave amplitudes on a grid of x
and y
coordinates:
function wavesolution{T}(X::AbstractVector, Y::AbstractVector, s::Scatterer{T}, ω::Real)
U = Array(Complex{T}, length(X), length(Y))
M,V = scat_amplitudes(s, ω, hypot(maximum(abs(X)), maximum(abs(Y))))
k = real(s.n[end])*ω
R = s.R[end]
for ix = 1:length(X)
x = X[ix]
for iy = 1:length(Y)
y = Y[iy]
r = hypot(x, y)
if r >= R
θ = atan2(y, x)
kr = k*r
u = besselj(0, kr) + V[1] * hankelh1(0, kr)
for i in 2:length(M) # m > 0 terms
m = M[i]
phase = im^m * complex(cos(m*θ), sin(m*θ))
u += 2*(besselj(m, kr) + V[i] * hankelh1(m, kr)) * (iseven(m) ? real(phase) : im*imag(phase))
end
U[ix,iy] = u
else
U[ix,iy] = 0.0 # don't try to compute interior solution
end
end
end
return U
end
wavesolution (generic function with 1 method)
Let's try a simple example scattering off a perfect-metal ($n = i \infty$) cylinder of radius 1 in air ($n=1$) with a frequency $\omega/c = 2\pi/0.5$ (wavelength $\lambda = 2\pi c / \omega = 0.5$):
s0 = Scatterer(1, [Inf*im], [1])
x = linspace(-4, 4, 200)
y = x
U = wavesolution(x, y, s0, 2π/0.5)
pcolormesh(y, x, real(U)', cmap="RdBu")
axis("equal")
colorbar()
PyObject <matplotlib.colorbar.Colorbar instance at 0x139dd9098>