From 8ee8de281e150f5ccda737f8e05f4a6d71d10bea Mon Sep 17 00:00:00 2001
From: baptiste <auguieba@gmail.com>
Date: Mon, 5 Dec 2022 21:57:00 +1300
Subject: [PATCH] fix mie

---
 dev/book/alpha_comparison.jl | 13 +++++++++++++
 dev/test-alpha.jl            | 22 ----------------------
 src/HighLevel.jl             |  4 ++--
 src/Materials.jl             |  9 +++++----
 4 files changed, 20 insertions(+), 28 deletions(-)

diff --git a/dev/book/alpha_comparison.jl b/dev/book/alpha_comparison.jl
index 267a72a..25fc9c7 100644
--- a/dev/book/alpha_comparison.jl
+++ b/dev/book/alpha_comparison.jl
@@ -29,6 +29,19 @@ oa1 = spectrum_oa(cl1, mat)
 oa2 = spectrum_oa(cl1, mat; prescription="majic")
 oa3 = spectrum_oa(cl1, mat; prescription="mie")
 
+# λ = 633.0
+# e = -10.0 + 1im
+# ε_m = 1.5^2
+# s = SVector(1.0, 1.0, 1.0)
+# CoupledDipole.alpha_majic(λ, e, ε_m, s)
+# CoupledDipole.alpha_mie(λ, e, ε_m, s)
+# Epsilon = [e, e]
+# Sizes = [s, s]
+# CoupledDipole.alpha_particles(Epsilon, Sizes, ε_m, λ; prescription="kuwata")
+# CoupledDipole.alpha_particles(Epsilon, Sizes, ε_m, λ; prescription="majic")
+# CoupledDipole.alpha_particles(Epsilon, Sizes, ε_m, λ; prescription="mie")
+
+
 
 # d1 = data(@rsubset(all, :polarisation == "pol2"))
 # d2 = data(@rsubset(single, :polarisation == "pol2"))
diff --git a/dev/test-alpha.jl b/dev/test-alpha.jl
index b313829..e69de29 100644
--- a/dev/test-alpha.jl
+++ b/dev/test-alpha.jl
@@ -1,22 +0,0 @@
-
-function alpha_majic(λ, ε, ε_m, a, c)
-
-    ε_r = (ε / ε_m)
-    k = 2π / λ * √ε_m
-
-    V = 4 / 3 * π * a * a * c
-    e = sqrt(c^2 - a^2) / c
-    La = -1 / (2 * e^2) * ((1 - e^2) / e * atanh(e) - 1)
-    Lc = (1 - e^2) / e^2 * (atanh(e) / e - 1)
-
-    AlphaS_x = V / 4 / π * (ε_r - 1) / (1 + (ε_r - 1) * La)
-    AlphaS_z = V / 4 / π * (ε_r - 1) / (1 + (ε_r - 1) * Lc)
-
-    Omega_x = 1 / 5 * (ε_r - 2 + 3 * e^2) / (1 + (ε_r - 1) * La) - 12 / 25 * e^2
-    Omega_z = 1 / 5 * (ε_r - 2 - ε_r * e^2) / (1 + (ε_r - 1) * Lc) + 9 / 25 * e^2
-
-    Ax = AlphaS_x / (1 - Omega_x * k^2 * c^2 - 1im * 2 / 3 * k^3 * AlphaS_x)
-    Az = AlphaS_z / (1 - Omega_z * k^2 * c^2 - 1im * 2 / 3 * k^3 * AlphaS_z)
-
-    SVector(Ax, Az)
-end
diff --git a/src/HighLevel.jl b/src/HighLevel.jl
index b38935e..560cf2b 100644
--- a/src/HighLevel.jl
+++ b/src/HighLevel.jl
@@ -105,7 +105,7 @@ function spectrum_dispersion(
             # )
 
             Epsilon = map(m -> mat.media[m](λ), cl.materials) # evaluate materials at wavelength
-            Alpha = alpha_particles(Epsilon, Sizes, n_medium^2, λ; prescription=prescription)
+            Alpha = alpha_particles(Epsilon, cl.sizes, n_medium^2, λ; prescription=prescription)
         end
 
         # update the rotated blocks
@@ -280,7 +280,7 @@ function spectrum_oa(
             #     cl.sizes,
             # )
             Epsilon = map(m -> mat.media[m](λ), cl.materials) # evaluate materials at wavelength
-            Alpha = alpha_particles(Epsilon, Sizes, n_medium^2, λ; prescription=prescription)
+            Alpha = alpha_particles(Epsilon, cl.sizes, n_medium^2, λ; prescription=prescription)
         end
 
         # update the rotated blocks
diff --git a/src/Materials.jl b/src/Materials.jl
index 48200c6..e31f431 100644
--- a/src/Materials.jl
+++ b/src/Materials.jl
@@ -179,7 +179,7 @@ end
 
 
 """
-    alpha_particles(λ, ε, ε_m, Sizes)
+alpha_particles(Epsilon, Sizes, ε_m, λ; prescription="kuwata")
 
 Principal polarisability components of N particles
 - `λ`: wavelength
@@ -208,7 +208,7 @@ function alpha_particles(Epsilon, Sizes, ε_m, λ; prescription="kuwata")
     elseif prescription == "mie"
         return (map((e, s) -> alpha_mie(λ, e, ε_m, s), Epsilon, Sizes))
     else
-        @warning "unknown prescription $prescription"
+        @warn "unknown prescription $prescription"
         return (map((e, s) -> alpha_kuwata(λ, e, ε_m, s), Epsilon, Sizes))
     end
 end
@@ -238,10 +238,11 @@ function alpha_mie(λ, ε, ε_m, Size)
     s = sqrt(ε) / n_medium
     k = n_medium * 2π / λ
     x = k * Size[1]
-    conversion = 2 / 3 * 1i * k^3
+    conversion = 2 / 3 * 1im * k^3
     Γ, Δ, A, B = mie_susceptibility(x, s, 1)
 
-    return Δ / conversion
+    A1 = Δ[1] / conversion
+    return SVector(A1, A1, A1)
 end
 
 """