something broke

Project.toml

@@ -18,6 +18,7 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Query = "1a8c2f83-1ff3-5112-b086-8aa67b057ba1"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 VegaDatasets = "0ae4a718-28b7-58ec-9efb-cded64d6d5b4"
 VegaLite = "112f6efa-9a02-5b7d-90c0-432ed331239a"

dev/book/alpha_comparison.jl

@@ -0,0 +1,42 @@
+# include("../src/CoupledDipole.jl")
+push!(LOAD_PATH, expanduser("~/Documents/nano-optics/CoupledDipole.jl/"))
+using Revise
+using CoupledDipole
+using Rotations
+using LinearAlgebra
+using StaticArrays
+using FastGaussQuadrature
+using DataFrames
+using DataFramesMeta
+using VegaLite
+using AlgebraOfGraphics, CairoMakie
+using ColorSchemes
+set_aog_theme!()
+
+## this example compares various prescriptions of polarisability
+# Mie vs Kuwata vs Majic
+
+## materials
+wavelength = collect(450:2:850.0)
+media = Dict([("Au", epsilon_Au), ("medium", x -> 1.33)])
+mat = Material(wavelength, media)
+
+cl0 = cluster_single(20.0, 20.0, 35.0)
+cl1 = cluster_dimer(100.0, 20.0, 20.0, 35.0, π / 4)
+
+oa0 = spectrum_oa(cl0, mat)
+oa1 = spectrum_oa(cl1, mat)
+oa2 = spectrum_oa(cl1, mat; prescription="majic")
+oa3 = spectrum_oa(cl1, mat; prescription="mie")
+
+
+# d1 = data(@rsubset(all, :polarisation == "pol2"))
+# d2 = data(@rsubset(single, :polarisation == "pol2"))
+
+# m1 = d1 * mapping(:wavelength, :value, color=:d => nonnumeric, col=:orientation, row=:crosstype)
+# m2 = d2 * mapping(:wavelength, :value, row=:crosstype)
+# layer1 = m1 * visual(Lines)
+# layer2 = m2 * visual(Lines, linestyle=:dash)
+# draw(layer1 + layer2, facet=(; linkyaxes=:none),
+#     palettes=(; color=cgrad(ColorSchemes.phase.colors, 12, categorical=true)))
+

src/CoupledDipole.jl

@@ -13,16 +13,15 @@ using DataFrames
 
 import Rotations: RotZYZ
 
+include("Rotations.jl")
+include("Utils.jl")
 include("Clusters.jl")
+include("Mie.jl")
+include("Materials.jl")
 include("CrossSections.jl")
 include("HighLevel.jl")
-include("Materials.jl")
-include("Mie.jl")
 include("PostProcessing.jl")
-include("Rotations.jl")
 include("Visual.jl")
-include("Utils.jl")
-
 
 
 # CoupledDipole

src/HighLevel.jl

@@ -12,6 +12,7 @@ Simulating far-field cross-sections for multiple wavelengths and directions of i
 - `Incidence`: N_inc vector of quaternions describing incidence directions
 - `polarisations`: incident field consists of 2 orthogonal "linear" or "circular" polarisations
 - `N_sca`: number of scattering angles for spherical cubature estimate of σ_sca
+- `prescription`: polarisability prescription for particles
 - `method`: direct or iterative solver
 
 """
@@ -21,6 +22,7 @@ function spectrum_dispersion(
     Incidence;
     N_sca::Int=36,
     polarisations="linear",
+    prescription="kuwata",
     method="direct"
 )
 
@@ -96,11 +98,14 @@ function spectrum_dispersion(
 
         elseif cl.type == "particle"
 
-            Alpha = map(
-                (m, s) -> alpha_kuwata(λ, mat.media[m](λ), n_medium^2, s),
-                cl.materials,
-                cl.sizes,
-            )
+            # Alpha = map(
+            #     (m, s) -> alpha_kuwata(λ, mat.media[m](λ), n_medium^2, s),
+            #     cl.materials,
+            #     cl.sizes,
+            # )
+
+            Epsilon = map(m -> mat.media[m](λ), cl.materials) # evaluate materials at wavelength
+            Alpha = alpha_particles(Epsilon, Sizes, n_medium^2, λ; prescription=prescription)
         end
 
         # update the rotated blocks
@@ -171,6 +176,8 @@ Orientation-averaged far-field cross-sections for multiple wavelengths
 - `Cubature`: spherical cubature method
 - `N_inc`: number of incident angles for spherical cubature
 - `N_sca`: number of scattering angles for spherical cubature estimate of σ_sca
+- `prescription`: polarisability prescription for particles
+- `method`: solution of linear system (direct or iterative)
 
 """
 function spectrum_oa(
@@ -179,6 +186,7 @@ function spectrum_oa(
     cubature="gl",
     N_inc=36,
     N_sca=36,
+    prescription="kuwata",
     method="direct"
 )
 
@@ -266,11 +274,13 @@ function spectrum_oa(
             # ε_name = cl.material # e.g. "Au" to refer to epsilon_Au in mat Dict
             # ε = mat.media[ε_name](λ)
             # Alpha = alpha_particles(λ, ε, n_medium^2, cl.sizes)
-            Alpha = map(
-                (m, s) -> alpha_kuwata(λ, mat.media[m](λ), n_medium^2, s),
-                cl.materials,
-                cl.sizes,
-            )
+            # Alpha = map(
+            #     (m, s) -> alpha_kuwata(λ, mat.media[m](λ), n_medium^2, s),
+            #     cl.materials,
+            #     cl.sizes,
+            # )
+            Epsilon = map(m -> mat.media[m](λ), cl.materials) # evaluate materials at wavelength
+            Alpha = alpha_particles(Epsilon, Sizes, n_medium^2, λ; prescription=prescription)
         end
 
         # update the rotated blocks

src/Materials.jl

@@ -199,16 +199,17 @@ julia> alpha_particles(500, -10+1im, 1.33^3, [SVector(30, 30, 50) for i in 1:4])
 ```
 
 """
-function alpha_particles(λ, ε, ε_m, Sizes; method="kuwata")
-    if method == "kuwata"
-        return (map(s -> alpha_kuwata(λ, ε, ε_m, s), Sizes))
-    elseif method == "majic"
-        return (map(s -> alpha_majic(λ, ε, ε_m, s), Sizes))
-    elseif method == "mie"
-        return (map(s -> alpha_mie(λ, ε, ε_m, s), Sizes))
+function alpha_particles(Epsilon, Sizes, ε_m, λ; prescription="kuwata")
+
+    if prescription == "kuwata"
+        return (map((e, s) -> alpha_kuwata(λ, e, ε_m, s), Epsilon, Sizes))
+    elseif prescription == "majic"
+        return (map((e, s) -> alpha_majic(λ, e, ε_m, s), Epsilon, Sizes))
+    elseif prescription == "mie"
+        return (map((e, s) -> alpha_mie(λ, e, ε_m, s), Epsilon, Sizes))
     else
-        @warning "unknown method $method"
-        return (map(s -> alpha_kuwata(λ, ε, ε_m, s), Sizes))
+        @warning "unknown prescription $prescription"
+        return (map((e, s) -> alpha_kuwata(λ, e, ε_m, s), Epsilon, Sizes))
     end
 end