Base

PhysicalMedium{Dim,FieldDim}

An abstract type used to represent the physical medium, the dimension of the field, and the number of spatial dimensions. We expect every sub-type of PhysicalMedium{Dim,1} to have a field c which is the complex wave speed.

ScalarMedium{Dim}

An type used to represent a scalar wave which satisfies a Helmoholtz equations. Dim represents the number of spatial dimensions.

basislength_to_basisorder(P::Type{<:PhysicalMedium}, len::Int)

The inverse of basisorder_to_basislength.

basisorder_to_basislength(P::Type{<:PhysicalMedium}, order::Int)

The number of basis functions, of the type of medium P, with basis order up to order. For scalar waves this is 2*order + 1 in 2D and (order + 1)^2 in 3D.

A tuples of vectors of the field close to the boundary of the shape. The field is calculated from sim::FrequencySimulation, but the PhysicalMedium inside and outside of the shape are assumed to be given by insidemedium and outsidemedium.

Extract the dimension of the field of this physical property

Extract the dimension of the field of this type of physical property

the field inside an AbstractParticle a some given point x.

Basis of outgoing wave. A series expansion in this basis should converge to any scattered field outside of a ball which contains the scatterer.

Basis of outgoing wave when the dependance on the angles has been removed.

outgoing_translation_matrix(medium::PhysicalMedium, in_order, out_order, ω, d)

Returns the matrix U which translates the outgoing basis into the regular basis. That is, U * vs(in_order, x) ≈ us(out_order, x + d) for |x| < |d|, where vs = regular_basis_function(medium, ω) and us = outgoing_basis_function(medium, ω).

A basis for regular functions, that is, smooth functions. A series expansion in this basis should converge to any regular function within a ball.

regular_translation_matrix(medium::PhysicalMedium, in_order, out_order, ω, d)

Returns the matrix V which translates the regular basis. That is, V * vs(in_order, x) ≈ vs(out_order, x + d), where vs = regular_basis_function(medium, ω).

Extract the dimension of the space that this physical property lives in

Extract the dimension of the space that this type of physical property lives in

complex_legendre_array(l_max::Int, z::Complex{T})

returns a vector of all associated legendre functions with degree l <= l_max. The degree and order (indices) of the elements of the vector are given by spherical_harmonics_indices(l_max::Int).

The associated legendre polynomials are taken from the package AssociatedLegendrePolynomials.jl.

diffbessel(f::Function,m,x,n::Int)

Differentiates 'n' times any bessel function 'f' of order 'm' and at the argument 'x'.

m-th Derivative of Bessel function of the first kind

Derivative of Bessel function of first kind

Derivative of Hankel function of the first kind

diffhankelh1(m,x,n::Int)

Differentiates 'n' times the hankelh1 function of order 'm' and at the argument 'x'.

diffsbessel(f::Function,m,x)

Differentiates the spherical bessel function 'f' (for any spherical bessel). The order is 'm' and the argument is 'x'. Note 'x' can be a complex number.

gaunt_coefficient(l1,m1,l2,m2,l3,m3)

A version of the Gaunt coefficients which are used to write the product of two spherical harmonics. If Y_{l,m} is a complex spherical harmonic, with the typical phase conventions from quantum mechanics, then:

gaunt_coefficient(l1,m1,l2,m2,l3,m3) = 4*π*im^{l2+l3-l1} Integral[Y_{l1,m1}*conj(Y_{l2,m2})*conj(Y_{l3,m3})]

where the integral is over the solid angle. In the papers by Art Gower, this symbol is denoted by c_{(l1,m1),(l2,m2)(l3,m3)}.

The most standard gaunt coefficients G(l1,m1;l2,m2;l3) are related through the identity:

4pi * G(l1,m1;l2,m2;l3) = im^(l1-l2-l3) * (-1)^m2 * gaunt_coefficient(l1,m1,l2,-m2,l3,m1+m2)

gaunt_cross_coefficient(l1,m1,l2,m2,l3,m3)

Similar to the Gaunt coefficient, but use for vector spherical harmonics, and specifically how one shear wave coverts to the other.

sbesselj(m,x)

Returns the spherical besselj function. The order is 'm' and the argument is 'x'. Note 'x' can be a complex number.

shankelh1(m,x)

Returns the spherical hankel function of the first kind. The order is 'm' and the argument is 'x'. Note 'x' can be a complex number.

spherical_harmonics(l_max::Int, θ::Complex{T}, φ::Union{T,Complex{T}})

returns a vector of all spherical harmonics with degree l <= l_max with complex arguments. The degree and order (indices) of the elements of the vector are given by spherical_harmonics_indices(l_max::Int).

spherical_harmonics(l_max::Int, θ::T, φ::T)

Returns a vector of all the spherial harmonics $Y_{(l,m)}(\theta,\phi)$ for all the degrees l and orders m. The degree and order (indices) of the elements of the vector are given by spherical_harmonics_indices(l_max::Int).

spherical_harmonics_dθ(l_max::Int, θ::T, φ::T)

Returns a vector of all $\partial Y_{(l,m)}(\theta,\phi) / \partial \theta$ for all the degrees l and orders m.

Electromagnetic

Represents the physical properties for a homogenous isotropic electromagnetic medium. Produces a three dimensional vector field.

Object we can scatter waves off

Subtypes will contain information about shape and material properties. Most crucially, they will implement the t_matrix function

CapsuleParticle(outer::Particle, inner::Particle)

A particle within another particle, both with the same shape type and origin.

Particle(medium::PhysicalMedium, shape::Shape)

Create particle with inner medium and shape (dimensions must agree).

Particle(medium, radius)

Returns a particle shaped like a sphere or circle, when the particle shape is not given and with the specified radius.

in(vector, particle)::Bool

Returns true if vector is in interior of shape of particle, false otherwise.

issubset(shape1, shape2)::Bool

Returns true if shape1 is entirely contained within shape2, false otherwise (also works with particles).

issubset(shape1, shape2)::Bool

Returns true if shape1 is entirely contained within shape2, false otherwise (also works with particles).

issubset(shape1, shape2)::Bool

Returns true if shape1 is entirely contained within shape2, false otherwise (also works with particles).

iscongruent(p1::AbstractParticle, p2::AbstractParticle)::Bool
≅(p1::AbstractParticle, p2::AbstractParticle)::Bool

Returns true if medium and shape of particles are the same, ignoring origin, false otherwise.

estimate_outgoing_basisorder(medium::PhysicalMedium, p::Particle, ω::Number; tol = 1e-6)

Estimate the basis order needed to represent the field scattered by the particle p. This is the smallest basis order such that increasing it changes the (normalised) diagonal of the T-matrix by less than tol.

iscongruent(p1::AbstractParticle, p2::AbstractParticle)::Bool
≅(p1::AbstractParticle, p2::AbstractParticle)::Bool

Returns true if medium and shape of particles are the same, ignoring origin, false otherwise.

Symmetry(T)
Symmetry(T1,T2)

returns the symmetry of the object T, or the combined symmetry of the two objects T1 and T2. Some possible symmetries include PlanarSymmetry, RadialSymmetry, and WithoutSymmetry.

AbstractSymmetry

An abstract types which dictates the symmetry of the setup. That is, the symmetry shared between the incident wave and the shape of the material.

WithoutSymmetry

A type used to describe materials and incident waves which share no common symmetry. This will lead to the most general regular wave expansion for the eignvectors.

An incident plane-wave and halfspace material will result in all fields being plane-waves.

PlanarSymmetry

A type used to describe materials and incident waves that both share a planar symmetry.

PlanarAzimuthalSymmetry

A type used to describe a materail and incident wave which have both PlanarSymmetry and AzimuthalSymmetry.

For spatial dimension > 2, we can consider problems that have azimuthal symmetry. For example, a plane-wave incident on a sphere.

AzimuthalSymmetry

A type used to describe materials and incident waves in 3 spatial dimensions that share a symmetry about the azimuth.

RadialSymmetry

A type used to describe materials and incident waves in that are both radially symmetric. That is, the material is a sphere, and the incident wave is radially symmetric around the center of this spherical material.

An incident plane-wave and infinite cylinder material will result in all fields being cylindrical Bessel waves with a phase.

TranslationSymmetry

A type used to describe materials and incident waves that both share a translation symmetry.

Represent any source (incident) wave

Subtypes may have a symmetry (such as PlaneSource) and will contain information about physical medium.

PlaneSource(medium::P, amplitude::SVector, direction::SVector)

Is a struct type which describes a plane-wave source that drives/forces the whole system. It has three fields: a physical medium, an amplitude of the source, and the direction the propagate in direction.

For any given angular frequency ω, the PlaneSource has the value $e^{i ω/c \mathbf v \cdot \mathbf x }$ at the point $\mathbf x$, where $c$ is the medium wavespeed and $\mathbf v$ is the direction.

RegularSource(medium::P, field::Function, coef::Function)

Is a struct type which describes the source field that drives/forces the whole system. It is also known as an incident wave. It has three fields RegularSource.medium, RegularSource.field, and RegularSource.coefficients.

The source field at the position 'x' and angular frequency 'ω' is given by

x = [1.0,0.0]
ω = 1.0
RegularSource.field(x,ω)

The field RegularSource.coefficients(order,centre,ω) returns the coefficients of the source field when expanded in the regular basis regular_basis_function centred at centre, up to the basis order order.

regular_spherical_coefficients(source::RegularSource)

return a function which can calculate the coefficients of a regular spherical wave basis.

regular_spherical_source(PhysicalMedium,regular_coefficients; amplitude = one(T), position = zeros(T,Dim))

$c_n =$regular_coefficients[n], $x_o=$position, and let $v_n(kx)$ represent the regular spherical basis with wavenumber $k$ at position $x$. The above returns a source where source.field(x,ω) =$\sum_n c_n v_n(kx -k x_0)$

Check that the source functions return the correct types

source_expand(RegularSource, centre; basis_order = 4)

Returns a function of (x,ω) which approximates the value of the source at (x,ω). That is, the source is written in terms of a regular basis expansion centred at centre.

See also: DiscreteImpulse, frequency_to_time

ContinuousImpulse{T<:AbstractFloat}

A struct used to represent an analytic impulse function. Has two fields: in_time a function of time t, and in_freq a function of the angular frequency ω. in_freq should be the Fourier transform of in_time, though this is not enforced.

We use the Fourier transform convention: F(ω) = ∫ f(t)exp(imωt) dt, f(t) = (2π)^(-1) * ∫ F(ω)exp(-imωt) dω.

An impluse f(t) is convoluted in time with the field u(t), however we avoid the convolution by working with the fourier transform F(ω) of the impulse f(t), which results in

frequency to time: (2π)^(-1) * ∫ F(ω)U(ω)exp(-imωt) dω

See also: ContinuousImpulse, frequency_to_time, DiscreteGaussianImpulse

DiscreteImpulse{T<:AbstractFloat}

A struct used to represent a numerical impulse. Only the fields: in_freq which is the frequency response vector, and the frequency vector ω are required to use this struct to use in the function frequency_to_time.

See also: ContinuousImpulse, TimeDiracImpulse

DiscreteGaussianImpulse(t_vec[, ω_vec])

Returns a discretised gaussian impulse.

Dirac Delta function of unit area in the frequency domain centred at ω=ω0.

Warning: in frequency space this is a singuarity and so may lead to unexpected behaviour.

See also: ContinuousImpulse, TimeDiracImpulse

GaussianImpulse(maxω[; σ = 3.0/maxω^2])

Returns a gaussian impulse function, which in the frequency domain is exp(-σ*ω^2)*(2sqrt(σ*pi)).

TimeDiracImpulse(t0::T)

Dirac Delta function of unit area in the time domain centred at t=t0.

Warning: in the time domain this is a singuarity and so may lead to unexpected behaviour.

continuous_to_discrete_impulse(impulse::ContinuousImpulse, t_vec, ω_vec = t_to_ω(t_vec); t_shift = 0.0, ω_shift = 0.0)

Returns a DiscreteImpulse by sampling impulse. The signal can be shifted in time and frequency by choosing t_shift and ω_shift.

Check that the continuous impulse functions return the correct types

Check that the discrete impulse vectors are the right sizes

Returns the first element of array which isn't zero (assumes elements are increasing and distinct)

Convert a FrequencySimulationResult into a TimeSimulationResult by using the inverse fourier transform. Assumes only positive frequencies and a real time signal

frequency_to_time(field_mat::AbstractArray, ω_vec::AbstractVector,
    t_vec::AbstractArray = ω_to_t(ω_vec);
    method = :DFT,
    impulse::ContinuousImpulse = TimeDiracImpulse(zero(T)),
    discrete_impulse::DiscreteImpulse = continuous_to_discrete_impulse(impulse, t_vec, ω_vec))

See also: DiscreteImpulse, ContinuousImpulse

Calculates the time response from the frequency response by approximating an inverse Fourier transform. The time signal is assumed to be real and the frequenices ω_vec are assumed to be positive (can include zero) and sorted. The result is convoluted in time with the user specified impulse.

We use the Fourier transform convention: F(ω) = ∫ f(t)exp(imωt) dt f(t) = (2π)^(-1) * ∫ F(ω)exp(-imωt) dt

To easily sample any time, the default is not FFT, but a discrete version of the transform above.

The inverse of ωtot if ω_vec[1] == 0

Convert a TimeSimulationResult into a FrequencySimulationResult by using the fourier transform. Assumes only positive frequencies and a real time signal

function time_to_frequency(field_mat::AbstractArray, t_vec::AbstractVector,
    ω_vec::AbstractArray = t_to_ω(t_vec);
    method = :DFT,
    impulse::ContinuousImpulse = TimeDiracImpulse(zero(T)),
    discrete_impulse::DiscreteImpulse = continuous_to_discrete_impulse(impulse, t_vec, ω_vec)
)

The inverse of the function frequencytotime (only an exact inverse when using :DFT integration). We use the Fourier transform convention: F(ω) = ∫ f(t)exp(imω*t) dt

returns an array of time from the frequency array ωvec. Uses the same convention for sampling the time as the discrete Fourier transfrom. Assumes ωvec is ordered and non-negative.

FrequencySimulation([particles::AbstractParticles=[],]
                    source::AbstractSource)

Build a FrequencySimulation. If particles are not provided, an empty array is used.

After building, you can run the simulation to get a FrequencySimulationResult.

run(sim::FrequencySimulation, x, ω; basis_order=5)

Run the simulation sim for the position x and angular frequency ω.

Position can be an SVector or Vector{SVector} and frequency can be a float or vector of floats.

run(sim::FrequencySimulation, region::Shape;
    res=20, xres=res, yres=res, basis_order=5)

Run the simulation sim for a grid of positions in region and for angular frequency ω.

Frequency can be a float or vector of floats. The resolution of the grid points is defined by xres and yres.

overlapping_pairs(particles::AbstractParticles)
overlapping_pairs(origins, radii)

Returns the pairs of indices [i, j] of the particles which overlap, where two particles overlap if the distance between their origins is less than the sum of their outer radiuses.

Abstract class for Results of Simulations

Struct to hold results of a FrequencySimulation

Struct to hold results of a simulation through time

basis_coefficients(sim::FrequencySimulation, ω::AbstractFloat; basis_order::Int=5)::Matrix{Complex}

Return coefficients for bases around each particle for a given simulation and angular frequency (ω).

field(result::SimulationResult, [i::Integer, j::Integer])

Get field from result, optionally specifying indices.

Returns single value of/matrix of complex SVectors() if vector field, and complex float if scalar field.

Create the matrix S which will be inverted to find the scattering coefficients.

t_matrix(particle, medium, ω, order)

Returns a finite T-matrix, with size depending on order, for a specific particle within a medium with specific physical properties.

t_matrix(Particle{2,Acoustic{T,2},Sphere{T,2}}, Acoustic{T,2}, ω, order)

The T-matrix for a 2D circlular acoustic particle in a 2D acoustic medium.

t_matrix(Particle{2,Acoustic{T,2},IsotropicHelmholtz{T,2}}, Acoustic{T,2}, ω, order)

The T-matrix for a 2D isotropic Helmholtz resonator in a 2D acoustic medium.

t_matrix(Particle{2,Acoustic{T,2},Helmholtz{T,2}}, Acoustic{T,2}, ω, order)

The T-matrix for a 2D non-isotropic Helmholtz resonator in a 2D acoustic medium.

t_matrix(CapsuleParticle{2,Acoustic{T,2},Sphere{T,2}}, Acoustic{T,2}, ω, order)

The T-matrix for a 2D circlular capsule particle in an acoustic medium.

t_matrix(Particle{3,Acoustic{T,3},Sphere{T,3}}, Acoustic{T,3}, ω, order)

The T-matrix for a 3D spherical acoustic particle in a 3D acoustic medium.

get_t_matrices(PhysicalMedium, Vector{Particles}, ω, basis_order::Integer)

Returns vector of T-matrices from a vector of particles in a specific domain. Can save computation if multiple of the same kind of particle are present in the vector.