From 7f9e3dfb17b1070f526b39c398c4c38eff370065 Mon Sep 17 00:00:00 2001
From: baptiste <baptiste.auguie@gmail.com>
Date: Fri, 2 Dec 2022 16:44:36 +1300
Subject: [PATCH] cluster shell

---
 Project.toml               |  2 +
 dev/book/molecule_shell.jl | 77 ++++++++++++++++++++++++++++++++++++++
 dev/book/particle_helix.jl |  2 +
 src/Clusters.jl            |  4 +-
 src/HighLevel.jl           |  6 +--
 5 files changed, 84 insertions(+), 7 deletions(-)
 create mode 100644 dev/book/molecule_shell.jl

diff --git a/Project.toml b/Project.toml
index e75a952..62af057 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,12 +4,14 @@ authors = ["baptiste"]
 version = "0.1.0"
 
 [deps]
+AlgebraOfGraphics = "cbdf2221-f076-402e-a563-3d30da359d67"
 Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
diff --git a/dev/book/molecule_shell.jl b/dev/book/molecule_shell.jl
new file mode 100644
index 0000000..646be56
--- /dev/null
+++ b/dev/book/molecule_shell.jl
@@ -0,0 +1,77 @@
+# include("../src/CoupledDipole.jl")
+push!(LOAD_PATH, expanduser("~/Documents/nano-optics/CoupledDipole.jl/"))
+using Revise
+using CoupledDipole
+using LinearAlgebra
+using StaticArrays
+using FastGaussQuadrature
+using DataFrames
+using VegaLite
+using Rotations
+
+## this example looks at a spherical shell of uniaxial molecules in water
+## contrasting radial and tangential configurations
+## extinction spectra for varying concentrations
+## orientation averaged extinction
+
+function dye_coverage(ρ, R) 
+    area = 4π * R^2
+    return Int(ceil(area * ρ))
+end
+# dye_coverage(1.2,5)
+
+## materials
+wavelength = collect(450:1:650.0)
+media = Dict([("Rhodamine", alpha_bare), ("medium", x -> 1.33)])
+mat = Material(wavelength, media)
+
+function model(; ρ=1, R0=2, d = 0.5, medium=1.33, orientation="radial")
+
+
+    R = R0 + d
+    N = dye_coverage(ρ, R) 
+    @info "$N dipoles"
+    if N > 1e4
+        @warn "$N dipoles is rather a lot, are you sure? "
+    end
+    if orientation == "radial"
+        # cl = cluster_shell(N, a, b, c, R; orientation="radial", material="Rhodamine", type="point")
+        cl = cluster_shell(N, 0, 0, 1, R; orientation="radial", material="Rhodamine", type="point")
+    elseif orientation == "flat"
+        cl = cluster_shell(N, 1, 1, 0, R; orientation="radial", material="Rhodamine", type="point")
+    else
+        @warn "orientation not recognised"
+    end
+    # testing
+    # cl = cluster_shell(5, 1, 1, 0, R; orientation="radial", material="Rhodamine", type="point")
+    res = spectrum_oa(cl, mat)
+    d = oa_df(res, mat.wavelengths)
+end
+
+
+# model(ρ  =100)
+params = expand_grid(ρ=range(0.2, 1.2, step=0.2), orientation=("radial", "flat"))
+
+all = pmap_df(params, p -> model(; p...))
+
+## reference molecule
+cl0 = cluster_single(0, 0,1,  0, 0, 0, "Rhodamine", "point")
+
+
+s = spectrum_oa(cl0, mat)
+single = oa_df(s, mat.wavelengths)
+
+using ColorSchemes
+set_aog_theme!()
+d1 = data(filter(:crosstype => ==("extinction"), all))
+d2 = data(filter(:crosstype => ==("extinction"), single))
+m1 = d1 * mapping(:wavelength, :value, color=:ρ => nonnumeric, col=:orientation, row=:crosstype)
+m2 = d2 * mapping(:wavelength, :value, row=:crosstype)
+layer1 = m1 * visual(Lines)
+layer2 = m2 * visual(Lines, linestyle=:dash)
+fg = draw(layer1 + layer2, facet=(; linkyaxes=:none),
+    palettes=(; color=cgrad(ColorSchemes.phase.colors, 12, categorical=true)))
+# https://docs.juliaplots.org/latest/generated/colorschemes/
+
+fg
+# save("figure.pdf", fg, px_per_unit=3)
\ No newline at end of file
diff --git a/dev/book/particle_helix.jl b/dev/book/particle_helix.jl
index f2220d7..ec49214 100644
--- a/dev/book/particle_helix.jl
+++ b/dev/book/particle_helix.jl
@@ -7,6 +7,8 @@ using StaticArrays
 using FastGaussQuadrature
 using DataFrames
 using VegaLite
+using AlgebraOfGraphics
+using Makie
 using Rotations
 
 ## this example looks at a helical strand of Au nanorods in water
diff --git a/src/Clusters.jl b/src/Clusters.jl
index 2097623..2314b0d 100644
--- a/src/Clusters.jl
+++ b/src/Clusters.jl
@@ -216,9 +216,9 @@ end
 
 function cluster_shell(N, a, b, c, R; orientation="radial", material="Rhodamine", type="point")
 
-    sizes = [SVector(a, b, c) for _ ∈ 1:N] # identical particles
-
     positions = R .* sample_fibonacci(N)
+    N = length(positions) # might be +1
+    sizes = [SVector(a, b, c) for _ ∈ 1:N] # identical particles
 
     if orientation == "radial"
 
diff --git a/src/HighLevel.jl b/src/HighLevel.jl
index f14c932..6631721 100644
--- a/src/HighLevel.jl
+++ b/src/HighLevel.jl
@@ -108,8 +108,6 @@ function spectrum_dispersion(
         # we'd compute rotation matrices on the fly
         # alpha_blocks!(AlphaBlocks, Alpha, cl.angles)
         # but instead we've prestored them, since wavelength-independent
-        # AlphaBlocks =
-        #     map((R, A) -> R * (diagm(A) * R'), ParticleRotations, Alpha)
         AlphaBlocks =
             map((R, A) -> R' * (diagm(A) * R), ParticleRotations, Alpha)
 
@@ -277,10 +275,8 @@ function spectrum_oa(
 
         # update the rotated blocks
         # alpha_blocks!(AlphaBlocks, Alpha, cl.angles)
-        # AlphaBlocks =
-        #     map((R, A) -> R' * (diagm(A) * R), ParticleRotations, Alpha)
         AlphaBlocks =
-            map((R, A) -> R * (diagm(A) * R'), ParticleRotations, Alpha)
+             map((R, A) -> R' * (diagm(A) * R), ParticleRotations, Alpha)
 
         propagator_freespace_labframe!(F, kn, cl.positions, AlphaBlocks)