fixed depolarisation factor, hopefully

.vscode/settings.json

@@ -1,4 +1,3 @@
 {
-    "C_Cpp.errorSquiggles": "Disabled",
-    "julia.environmentPath": "/Users/baptiste/Documents/nano-optics/CoupledDipole.jl"
+    "C_Cpp.errorSquiggles": "Disabled"
 }

Project.toml

@@ -4,10 +4,13 @@ authors = ["baptiste"]
 version = "0.1.0"
 
 [deps]
+Arrow = "69666777-d1a9-59fb-9406-91d4454c9d45"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"

dev/book/dimer_linear.jl

@@ -10,13 +10,6 @@ using FastGaussQuadrature
 using DataFrames
 using VegaLite
 
-## materials
-wavelength = collect(450:2:850.0)
-media = Dict([("Au", epsilon_Au), ("medium", x -> 1.33)])
-mat = Material(wavelength, media)
-
-## reference geometry
-cl0 = cluster_single(20, 20, 40)
 
 
 # dimer_model <- function(d, orientation = c("head-to-tail", "side-by-side"), 
@@ -36,75 +29,130 @@ cl0 = cluster_single(20, 20, 40)
 #   }
 # }
 
-UnitQuaternion(0, 0, 1, 0) # rotation π/2 about y
-UnitQuaternion(√2/2, √2/2, 0, 0)
-rotation_angle(UnitQuaternion(√2/2, √2/2, 0, 0))
-rotation_axis(UnitQuaternion(√2/2, √2/2, 0, 0))
-RotZYZ(π/2,0,0)
-RotMatrix(UnitQuaternion(√2/2, √2/2, 0, 0))
-
-RotY(π/2)
-
-
-function model(; d=80, orientation="head-to-tail")
+# UnitQuaternion(0, 0, 1, 0) # rotation π/2 about y
+# UnitQuaternion(√2/2, √2/2, 0, 0)
+# rotation_angle(UnitQuaternion(√2/2, √2/2, 0, 0))
+# rotation_axis(UnitQuaternion(√2/2, √2/2, 0, 0))
+# RotZYZ(π/2,0,0)
+# RotMatrix(UnitQuaternion(√2/2, √2/2, 0, 0))
+
+# RotY(π/2)
+
+# initial version, but now prefer low level cluster
+# function model(; d=80, orientation="head-to-tail")
+#     if orientation == "head-to-tail"
+#         # dimer axis along y, particle axis along y
+#         cl = cluster_dimer(d, 20., 40, 20, 0)
+#         ## incidence: along z (no rotation)
+#         Incidence = [RotZ(0.)]
+#         res = spectrum_dispersion(cl, mat, Incidence)
+#         d = dispersion_df(res, mat.wavelengths)
+#         return(filter(:polarisation => ==("p"), d))
+#     elseif orientation == "side-by-side"
+#         # dimer axis along y, particle axis along z
+#         cl = cluster_dimer(d, 20., 20, 40, 0)
+#         ## incidence: along x so rotate z around y axis by π/2
+#         Incidence = [RotY(π/2)]
+#         res = spectrum_dispersion(cl, mat, Incidence)
+#         d = dispersion_df(res, mat.wavelengths)
+#         return(filter(:polarisation => ==("s"), d))
+#     end
+
+#     error("unknown orientation given")
+# end
+
+function model(; d=100, orientation="head-to-tail")
     if orientation == "head-to-tail"
-        # dimer axis along y, particle axis along y
-        cl = cluster_dimer(d, 20., 40, 20, 0)
-        ## incidence: along z (no rotation)
-        Incidence = [RotZ(0.)]
-        res = spectrum_dispersion(cl, mat, Incidence)
-        d = dispersion_df(res, mat.wavelengths)
-        return(filter(:polarisation => ==("p"), d))
+        # dimer axis along y
+        positions = [SVector(0.0, y, 0.0) for y in (-d / 2.0, d / 2.0)]
     elseif orientation == "side-by-side"
-        # dimer axis along y, particle axis along z
-        cl = cluster_dimer(d, 20., 20, 40, 0)
-        ## incidence: along x so rotate z around y axis by π/2
-        Incidence = [RotY(π/2)]
-        res = spectrum_dispersion(cl, mat, Incidence)
-        d = dispersion_df(res, mat.wavelengths)
-        return(filter(:polarisation => ==("s"), d))
+        # dimer axis along x
+        positions = [SVector(x, 0.0, 0.0) for x in (-d / 2.0, d / 2.0)]
     end
-
-    error("unknown orientation given")
-end
 
+    sizes = [SVector(20.0, 50.0, 20.0), SVector(20.0, 50.0, 20.0)]
+    rotations = repeat([UnitQuaternion(1.0, 0.0, 0.0, 0.0)], 2)
+    materials = repeat(["Au"], 2)
+    cl = Cluster(positions, rotations, sizes, materials, "particle")
+
+    Incidence = [RotZ(0.0)] ## incidence: along z (no rotation)
+    res = spectrum_dispersion(cl, mat, Incidence)
+    d = dispersion_df(res, mat.wavelengths)
+end
 
-wavelength = collect(400:2:800.0)
+## materials
+wavelength = collect(450:2:850.0)
 media = Dict([("Au", epsilon_Au), ("medium", x -> 1.33)])
 mat = Material(wavelength, media)
-# cl = cluster_dimer(80, 20., 40, 20, 0)
-# Incidence = [RotZ(0.),RotZ(π/2.),RotZ(0.),RotZ(π/2.),RotZ(0.)]
-# res = spectrum_dispersion(cl, mat, Incidence)
-# d = dispersion_df(res, mat.wavelengths)
-# print(d)
 
-params = expand_grid(d=range(80, 200, step=20), orientation=("head-to-tail", "side-by-side"))
 
+params = expand_grid(d=range(100, 500, step=50), orientation=("head-to-tail", "side-by-side"))
 all = pmap_df(params, p -> model(; p...))
 
 
-s = spectrum_dispersion(cl0, mat, [RotZ(0.)])
-single = dispersion_df(s, mat.wavelengths)
-
+## reference geometry
+wavelength = collect(450:2:850.0)
+media = Dict([("Au", epsilon_Au), ("medium", x -> 1.33)])
+mat = Material(wavelength, media)
 
-s = spectrum_oa(cl0, mat)
-single = oa_df(s, mat.wavelengths)
+cl0 = cluster_single(20.0, 50.0, 20.0)
+# cl0 = cluster_single(50.0, 20.0, 20.0)
+# s = spectrum_dispersion(cl0, mat, [UnitQuaternion(RotZ(0.0))])
+s = spectrum_dispersion(cl0, mat, [SVector(0.0, 0.0, 0.0)])
+single = dispersion_df(s, mat.wavelengths)
 
 
-xy = data(single) * mapping(:wavelength, :value, color=:type, row=:crosstype,col=:type)
+xy = data(single) * mapping(:wavelength, :value, row=:crosstype, col=:polarisation)
 layer = visual(Lines)
 draw(layer * xy, facet=(; linkyaxes=:none))
 
 
+# d = [insertcols(all, :cluster => "dimer"),
+#     insertcols(single, :cluster => "single", :d => missing, :orientation => missing)]
 
+# @vlplot(data = d,
+#     width = 400,
+#     height = 300,
+#     mark = {:line},
+#     row = "crosstype",
+#     resolve = {scale = {y = "independent"}},
+#     encoding = {x = "wavelength:q", y = "value:q", color = "d:n", strokeDash = "cluster:n"}
+# )
 
 using AlgebraOfGraphics, CairoMakie
 
 # set_aog_theme!()
 
 
-xy = data(all) * mapping(:wavelength, :value, color=:d => nonnumeric, col=:orientation, row=:crosstype, linestyle=:crosstype => nonnumeric)
+xy = data(filter(:polarisation => ==("p"), all)) * mapping(:wavelength, :value, color=:d => nonnumeric, col=:orientation, row=:crosstype)
 layer = visual(Lines)
 draw(layer * xy, facet=(; linkyaxes=:none))
 
 
+
+
+
+
+# s = spectrum_oa(cl0, mat)
+# single = oa_df(s, mat.wavelengths)
+
+
+# xy = data(single) * mapping(:wavelength, :value, color=:type, row=:crosstype,col=:type)
+# layer = visual(Lines)
+# draw(layer * xy, facet=(; linkyaxes=:none))
+
+using Arrow
+Arrow.write("test.arrow", single)
+Arrow.write("testc.arrow", single, compress=:lz4)
+
+
+y = Arrow.Table("test.arrow") |> DataFrame
+x = Arrow.Table("testr.arrow") |> DataFrame
+eltype.(eachcol(x))
+
+# using CSV
+
+# CSV.write("test.csv", single)
+
+
+

dev/debug.R

@@ -0,0 +1,65 @@
+library(cda)
+# library(rgl)
+library(ggplot2)
+library(dielectric)
+library(cubs)
+library(tidyr)
+library(dplyr)
+library(purrr)
+library(progress)
+
+
+wavelength = c(400, 500)
+gold = epsAu(wavelength)
+
+cl = structure(list(positions = structure(c(-50, 0, 0, 50, 0, 00), dim = c(3,2)),
+                    sizes = structure(c(50, 20, 20, 50, 20, 20), dim = c(3,2)),
+                    angles = structure(c(0, 0, 0, 0, 0, 0), dim = c(3,2))),
+               class = "cluster")
+
+cl0 = structure(list(positions = structure(c(0, 0, 0), dim = c(3,1)),
+                    sizes = structure(c(50, 20, 20), dim = c(3,1)),
+                    angles = structure(c(0, 0, 0), dim = c(3,1))),
+               class = "cluster")
+
+cl0
+cl$positions
+
+
+ii = 1;
+
+lambda = gold$wavelength[ii]
+n_medium = 1.33
+kn = n_medium * 2*pi / lambda
+
+# gold[ii,]
+# -1.6496568840720354 + 5.7717630808981655im
+
+# alpha_kuwata(500, -10+1im, 1.33^2, SVector(30, 30, 50))
+# 3-element SVector{3, ComplexF64} with indices SOneTo(3):
+#   77076.04648078185 + 26235.664281642246im
+# 77076.04648078185 + 26235.664281642246im
+# -98187.15974124733 + 205835.30299929058im
+# V <- 4 * pi/3 * 30*30*50
+# chi <- cda:::depolarisation(30,30,50)
+# cda:::alpha_kuwata(500, -10+1i, V, 30, chi[1], 1.33)
+# cda:::alpha_kuwata(500, -10+1i, V, 50, chi[3], 1.33)
+
+Alpha = cda::alpha_ellipsoid(sizes = cl0$sizes, material = gold[ii,], medium = n_medium)
+
+# Alpha = alpha_spheroids(λ, ε, n_medium^2, cl.sizes)
+# 2-element Vector{SVector{3, ComplexF64}}:
+# [10898.262981422746 + 18516.11377073886im, -903.6988080340118 + 24897.841185369187im, 10898.262981422744 + 18516.113770738863im]
+# [10898.262981422746 + 18516.11377073886im, -903.6988080340118 + 24897.841185369187im, 10898.262981422744 + 18516.113770738863im]
+
+# alpha_kuwata(400, ε, 1.33^2, SVector(20, 20, 50))
+# 3-element SVector{3, ComplexF64} with indices SOneTo(3):
+#   12403.616901231882 + 13704.077043611223im
+# 12403.616901231882 + 13704.077043611223im
+# -6869.361725855006 + 22867.056384367785im
+
+ParticleRotations = cl.rotations
+AlphaBlocks = map((R, A) -> R' * diagm(A) * R, ParticleRotations, Alpha)
+
+
+

dev/read.R

@@ -0,0 +1,14 @@
+# install.packages("arrow")
+
+library(arrow)
+
+a <- arrow::read_feather('test.arrow')
+str(a)
+
+b <- arrow::read_feather('testc.arrow')
+str(b)
+
+c <- arrow::read_csv_arrow('test.csv')
+str(c)
+
+arrow::write_feather(a, 'testr.arrow', compression='lz4')