A Julia package for solving large-scale electromagnetic scattering problems in two dimensions; specifically, those containing a large number of penetrable smooth particles. Provides the ability to optimize over the particle parameters for various design problems.

R_multipole = 1.1 is the constant ratio between scattering disks and the maximal radius of their particles, and thus half the minimal distance between neighboring particles. While mathematically this can be reduced to 1 + eps(), that will increase the necessary P.

AbstractShapeParams

Abstract type which all shape types inherit from.

CircleParams(R)

Returns object for a circular shape, containing its radius in the field R (which is also the radius of the scattering disk).

See also: ShapeParams,R_multipole.

CurrentSource(x1, y1, x2, y2, σ)

Constructor for the CurrentSource type, where (x1,y1) and (x2,y2) denote the start and end points of the source, and σ contains the potential density.

FMMoptions(FMM; nx = 0, dx = 0.0, acc = 0, tol = 0.0, method = "pre")

Constructor for struct containing all FMM options. FMM decides if FMM is used, and the following keyword arguments dictate its behavior:

• nx::Integer: number of groups in x direction (for division)
• dx::Real: group height/width (alternative division)
• acc::Integer: accuracy digits for translation truncation, and also for gmres if tol is not given
• tol::Real: gmres tolerance
• method::String: method used: for now can be "pre". Mainly used for development.

LineSource(x, y)

Constructor for the LineSource type, where (x,y) is the coordinate of the current filament.

OptimBuffer(Ns::Integer, P::Integer, Npoints::Integer, [J::Integer])

Constructor for the OptimBuffer type, which stores some of the buffers and shared variables necessary for optimization. Includes the cylindrical harmonics coefficient vector β, field values at points of interest (f), the partial derivatives ∂β, and storage for the various right-hand side vectors used while solving for ∂β.

If the number of optimization variables J is not supplied, it is assumed that J = Ns.

PlaneWave(θi)

Constructor for the PlaneWave incident field type, where θi is the angle between the wavevector and the x-axis.

ScatteringProblem(shapes, ids, centers, φs)

Constructor for the ScatteringProblem type, including particle shape information for multiple-scattering problems.

ShapeParams(t,ft,dft)

Returns ShapeParams object containing the parametrization of a two-dimensional shape. t is a uniform sampling of [0,2π), ft = [x(t) y(t)], and dft = [x'(t) y'(t)]. The field R contains the radius of the shape's scattering disk.

See also: CircleParams,R_multipole.

calc_near_field(k0, kin, P, sp::ScatteringProblem, points, ui::Einc;
                        opt::FMMoptions = FMMoptions(), method = "multipole",
                        verbose = true)

Calculates the total electric field as a result of a plane wave with incident field ui scattering from the ScatteringProblem sp, at points. Uses the FMM options given by opt (default behaviour is disabled FMM); method = "multipole" dictates whether electric field is calculated using the multipole/cylindrical harmonics, uses a faster but less accurate Hankel recurrence formula ("recurrence"), or falls back on potential densities ("density"). Either way, the multiple-scattering system is solved in the cylindrical harmonics space, and the field by a particular scatterer inside its own scattering discs is calculated by potential densities, as the cylindrical harmonics approximation is not valid there.

draw_shapes(shapes, ids, centers, φs; ax = gca(), normalize = 1.0)
draw_shapes(sp; ax = gca(), normalize = 1.0)

Draws all of the shapes in a given scattering problem in the PyPlot axis 'ax'. Parametrized shapes are drawn as polygons while circles are drawn using matplotlib's patch.Circle. Divides all lengths by 'normalize'.

ellipse(r1, r2, N)

Return a ShapeParams object containing the shape parametrized by (x/r1)^2 + (y/r2)^2 = 1 with 2N nodes.

find_border(sp::ScatteringProblem, points::Array{Float64,2}) -> [x_min; x_max; y_min; y_max]

Returns bounding box that contains all of the shapes in sp as well as specified points.

find_border(sp::ScatteringProblem) -> [x_min; x_max; y_min; y_max]

Returns bounding box that contains all of the shapes in sp.

get_potential(kout, kin, P, s::ShapeParams) -> sigma_mu

Given a shape s with 2N discretization nodes, outer and inner wavenumbers kout,kin, and the cylindrical harmonics parameter P, returns the potential densities sigma_mu. Each column contains the response to a different harmonic, where the first 2N entries contain the single-layer potential density ($\sigma$), and the lower entries contain the double-layer density ($\mu$).

get_potential(kout, kin, P, t, ft, dft) -> sigma_mu

Same, but with the ShapeParams supplied directly.

get_potentialPW(kout, kin, s::ShapeParams, θ_i) -> sigma_mu

Given a shape s with 2N discretization nodes, outer and inner wavenumbers kout,kin, and an incident plane-wave angle, returns the potential densities vector sigma_mu. The first 2N entries contain the single-layer potential density ($\sigma$), and the lower entries contain the double-layer density ($\mu$).

hex_grid(a::Integer, rows::Integer, d; minus1 = false)

Return centers, an (M,2) array containing points on a hexagonal lattice with horizontal rows, with a points distanced d in each row and rows rows. If minus1 is true, the last point in every odd row is omitted.

luneburg_grid(R_lens, N_cells, er; levels = 0, TM = true) -> centers, ids, rs

Returns the coordinates and radii of the circular inclusions in a Luneburg lens device of radius R_lens with N_cells unit cells across its diameter. Radii are determined by averaging over cell permittivity, assuming air outside and relative permittivity er in the rods, and depends on incident field polarization (TM/TE with respect to z-axis). If levels == 0, groups identical radii together, such that rs[ids[n]] is the radius of the rod centered at (center[n,1],center[n,2]). Otherwise quantizes the radii to uniformly spaced levels.

minimumN(kout, kin, shape_function; tol = 1e-9, N_points = 10_000,
    N_start = 400, N_min = 100, N_max = 1_000) -> N, err

Return the minimum N necessary (i.e. 2N nodes) to achieve error of at most tol in the electric field for a ShapeParams inclusion created by shape_function(N) which is filled with material of wavenumber kin and surrounded by free space with wavenumber k0. Error is calculated on N_points points on the scattering disk (s.R), by assuming a fictitious line source and comparing its field to that produced by the resulting potential densities.

Since the error scales with $N^{-3}$ for moderate wavelengths and errors, we estimate N using the error of N_start, then binary search based on that guess between N_min and N_max.

minimumP(k0, kin, s::ShapeParams; tol = 1e-9, N_points = 10_000, P_min = 1,
    P_max = 60, dist = 2) -> P, errP

Return the minimum P necessary to achieve error of at most tol in the electric field, when compared to that obtained with 2N discretization, for a ShapeParams inclusion filled with material of wavenumber kin and surrounded by free space with wavenumber k0. Error is calculated on N_points points on a disk of radius dist*s.R.

Uses binary search between P_min and P_max.

optimize_radius(rs0, r_min, r_max, points, ids, P, ui, k0, kin, centers,
    fmmopts, optimopts::Optim.Options; minimize = true, method = "BFGS")

Optimize the radii of circular particles for minimization or maximization of the field intensity at points, depending on minimize. Uses Optim's Fminbox box-contrained optimization to contain radii in feasible rangle, given in scalar or vector form by r_min and r_max.

Here, ids allows for grouping particles - for example, to maintain symmetry of the optimized device. optimopts defines the convergence criteria and other optimization parameters for both the inner and outer iterations. method can be either "BFGS" or "LBFGS". See the Optim.Fminbox documentation for more details. adjoint dictates whether the gradient is calculated using the adjoint method (faster) or directly.

Returns an object of type Optim.MultivariateOptimizationResults.

optimize_φ(φs0, points, P, ui::Einc, k0, kin, shapes, centers, ids, fmmopts,
    optimopts::Optim.Options, minimize = true)

Optimize the rotation angles of a particle collection for minimization or maximization (depending on minimize) of the field intensity at points. optimopts and method define the optimization method, convergence criteria, and other optimization parameters. adjoint dictates whether the gradient is calculated using the adjoint method (faster) or directly. Returns an object of type Optim.MultivariateOptimizationResults.

plot_far_field(k0, kin, P, sp::ScatteringProblem, pw::PlaneWave;
                    opt::FMMoptions = FMMoptions(), method = "multipole",
                    plot_points = 200)

Plots and returns the echo width (radar cross section in two dimensions) for a given scattering problem. opt, method are as in plot_near_field.

plot_near_field(k0, kin, P, sp::ScatteringProblem, ui::Einc;
                    opt::FMMoptions = FMMoptions(), method = "multipole",
                    x_points = 201, y_points = 201, border = find_border(sp),
                    normalize = 1.0)

Plots the total electric field as a result of a plane wave with incident TM field ui scattering from the ScatteringProblem sp, using matplotlib's pcolormesh. Can accept number of sampling points in each direction plus bounding box or calculate automatically.

Uses the FMM options given by opt (FMM is disabled by default); method = "multipole" dictates whether electric field is calculated using the multipole/cylindrical harmonics, uses a faster but less accurate Hankel recurrence formula ("recurrence"), or falls back on potential densities ("density"). Either way, the multiple-scattering system is solved in the cylindrical harmonics space. Normalizes all distances and sizes in plot (but not output) by normalize.

Returns the calculated field in two formats:

1. (points, Ez) where Ez[i] is the total electric field at points[i,:], and
2. (xgrid,ygrid,zgrid), the format suitable for pcolormesh, where zgrid[i,j]

contains the field at (mean(xgrid[i, j:j+1]), mean(ygrid[i:i+1, j])).

randpoints(M, dmin, width, height; failures = 100)

Return centers, an (M,2) array containing M points distanced at least dmin in a width by height box. Fails failures times successively before giving up.

randpoints(M, dmin, width, height, points; failures = 100)

Same as randpoints(M, dmin, width, height; failures = 100) but also requires centers to be distanced at least dmin from points.

rect_grid(a::Integer, b::Integer, dx, dy)

Return centers, an (a*b,2) array containing the points spanned by a points distanced dx and b points distanced dy, in the x and y directions, respectively.

rounded_star(r, d, num, N)

Return a ShapeParams object containing the shape parametrized by $(x(θ),y(θ)) = (r + d*cos(θ*num))*(cos(θ),sin(θ))$ with 2N nodes.

scatteredfield(sigma_mu, k, t, ft, dft, p) -> u_s

Same, but with the ShapeParams supplied directly. Useful for computing u_s for rotated shapes.

scatteredfield(sigma_mu, k, s::ShapeParams, p) -> u_s

Computes field scattered by the particle s with pre-computed potential densities sigma_mu at points p. All points must either be inside k = kin or outside k = kout the particle.

solve_particle_scattering(k0, kin, P, sp::ScatteringProblem, u::Einc; get_inner = true, verbose = true) -> beta, inner

Solve the scattering problem sp with outer wavenumber k0, inner wavenumber kin, 2P+1 cylindrical harmonics per inclusion and incident TM field u. Solves multiple-scattering equation directly. Returns the cylindrical harmonics basis beta along with potential densities (in case of arbitrary inclusion) or inner cylindrical coefficients (in case of circular). By default, incident wave propagates left->right.

Inner coefficients are only calculated if get_inner is true, and timing is printed if verbose is true.

solve_particle_scattering_FMM(k0, kin, P, sp::ScatteringProblem, u::Einc, opt::FMMoptions; plot_res = false, get_inner = true, verbose = true) -> result, inner

Solve the scattering problem sp with outer wavenumber k0, inner wavenumber kin, 2P+1 cylindrical harmonics per inclusion and incident TM field u. Utilizes FMM with options opt to solve multiple-scattering equation. Returns the cylindrical harmonics basis beta along with convergence data in result. inner contains potential densities (in case of arbitrary inclusion) or inner cylindrical coefficients (in case of circular).

plot_res controls plotting of the residual. Inner coefficients are calculated only if get_inner is true, and timing is printed if verbose is true.

square_grid(a::Integer, d)

Return centers, an (a^2,2) array containing the points on an a by a grid of points distanced d.

squircle(r, N)

Return a ShapeParams object containing the shape parametrized by $x(θ)^4 + y(θ)^4 = r^4$ with 2N nodes.

uinc(k, points, u::Einc)

Computes incident field produced by u with outer wavenumber k at points.

uniqueind(v::Vector{T}) where T <: Number -> inds,u

Given a vector of numbers v of length n, returns the unique subset u as well as a vector of indices inds of length n such that v == u[inds].

verify_min_distance(shapes, centers::Array{Float64,2}, ids, points::Array{Float64,2})
verify_min_distance(sp::ScatteringProblem, points)

Returns true if the shapes placed at centers are properly distanced (non-intersecting scattering disks), and all points are outside the scattering disks.

verify_min_distance(shapes, centers::Array{Float64,2}, ids)
verify_min_distance(sp::ScatteringProblem)

Returns true if the shapes placed at centers are properly distanced (non-intersecting scattering disks).