use SparseArrays

Author: arturgower

Project.toml

@@ -12,6 +12,7 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

src/MultipleScattering.jl

@@ -48,6 +48,7 @@ export scattering_matrix
 
 import Printf: @printf
 
+using SparseArrays: sparse
 using StaticArrays: SVector
 using OffsetArrays: OffsetArray
 

src/physics/special_functions.jl

@@ -183,7 +183,7 @@ A version of the Gaunt coefficients which are used to write the product of two s
 
     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.
+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:
 

src/scattering_matrix.jl

@@ -16,7 +16,7 @@ function scattering_matrix(medium::PhysicalMedium, particles::AbstractParticles,
     # Faire: this could potentially return an MMatrix
     function S_block(j,l)
         if j == l
-            return zeros(Complex{T}, N, N)
+            return sparse(zeros(Complex{T}, N, N))
         else
             x_lj = origin(particles[j]) .- origin(particles[l])
             U = outgoing_translation_matrix(medium, basis_order, basis_order, ω, x_lj)
@@ -28,7 +28,7 @@ function scattering_matrix(medium::PhysicalMedium, particles::AbstractParticles,
     S_blocks = [S_block(j,l) for j in 1:P, l in 1:P]
 
     # Reshape S_blocks into big matrix
-    S = zeros(Complex{T}, P*N, P*N)
+    S = sparse(zeros(Complex{T}, P*N, P*N))
     for i in 1:P
         for j in 1:P
             S[((i-1)*N+1):(i*N), ((j-1)*N+1):(j*N)] .= S_blocks[i,j]

src/source.jl

@@ -37,12 +37,12 @@ struct PlaneSource{T<:Real,Dim,FieldDim,P<:PhysicalMedium} <: AbstractSource{P}
 end
 
 field(s::PlaneSource{T}) where T = function (x::AbstractArray{T}, ω::T)
-        s.amplitude * exp(im * (ω / s.medium.c) * dot(x - s.position,s.direction))
-    end
+    s.amplitude * exp(im * (ω / s.medium.c) * dot(x - s.position,s.direction))
+end
 
 field(s::PlaneSource{T,Dim,1}) where {T, Dim} = function (x::AbstractArray{T}, ω::T)
-        s.amplitude[1] * exp(im * (ω / s.medium.c) * dot(x - s.position,s.direction))
-    end
+    s.amplitude[1] * exp(im * (ω / s.medium.c) * dot(x - s.position,s.direction))
+end
 
 field(s::PlaneSource{T}, x::AbstractArray{T}, ω::T) where T = field(s)(x,ω)