Curl_curl function shoud be k * k' - I abs2(k)

Project.toml

@@ -1,7 +1,7 @@
 name = "LinearMaxwellVlasov"
 uuid = "b1137d70-0135-468f-bb3e-ace6f597c457"
 authors = ["James Cook <[email protected]>"]
-version = "0.1.8"
+version = "0.1.9"
 
 [deps]
 CommonSubexpressions = "bbf7d656-a473-5ed7-a52c-81e309532950"

src/Tensors.jl

@@ -189,11 +189,12 @@ end
 Calculate the tensor representing the linear Maxwell-Vlasov set of equations.
 The determinant is zero when the wavenumber and frequency represent a solution
 to the linear Maxwell-Vlasov system of equations for these species.
+Note that the curlcurl operator keeps track of the factor of -1 from im^2.
 """
 function tensor(plasma::AbstractPlasma, config::Configuration,
     cache::Cache=Cache())
   ϵᵢⱼ = dielectric(plasma, config, cache)
-  return ϵᵢⱼ + curl_curl(config.wavenumber) * (c₀ / config.frequency)^2
+  return ϵᵢⱼ - curlcurl(config.wavenumber) * (c₀ / config.frequency)^2
 end
 
 """

src/Wavenumbers.jl

@@ -47,9 +47,24 @@ Base.angle(K::Wavenumber) = propagationangle(K)
 
 @inline cartesian_vector(K::Wavenumber) = @SArray [perp(K), 0.0, para(K)]
 
-@inline function curl_curl(K::Wavenumber)
-  k = cartesian_vector(K)
-  return k*k' - dot(k, k) * I
+"""
+julia> using Symbolics
+julia> @syms kx ky kz;
+julia> k = [kx, ky, kz];
+julia> curl = im * [0 -k[3] k[2]; k[3] 0 -k[1]; -k[2] k[1] 0];
+julia> curl * curl
+3×3 Matrix{Any}:
+ ky^2 + kz^2  -kx*ky       -kx*kz
+ -kx*ky       kx^2 + kz^2  -ky*kz
+ -kx*kz       -ky*kz       kx^2 + ky^2
+
+"""
+@inline function curlcurl(K::Wavenumber)
+  kx = perp(K) # ky = 0
+  kz = para(K)
+  return @SArray [ kz^2    0           -kx * kz;
+                   0       kx^2 + kz^2  0      ;
+                  -kx * kz 0            kx^2   ]
 end
 
 Base.abs(K::Wavenumber) = sqrt(abs2(K))

test/Wavenumbers.jl

@@ -48,12 +48,14 @@ const k0 = -1.234
   end
 
   @testset "curl curl operator" begin
-    K = Wavenumber(wavenumber=k0, propagationangle=π/4)
-    k_x_k = LMV.curl_curl(K)
-    c = LMV.cartesian_vector(K)
+    for k1 in (k0, k0 + im)
+      K = Wavenumber(wavenumber=k0, propagationangle=π/4)
+      ∇x∇x = LMV.curlcurl(K)
+      c = LMV.cartesian_vector(K)
 
-    x = [0.0 -c[3] c[2]; c[3] 0.0 -c[1]; [-c[2] c[1] 0.0]]
-    @test all(x*x .== k_x_k)
+      curl = im * [0.0 -c[3] c[2]; c[3] 0.0 -c[1]; [-c[2] c[1] 0.0]]
+      @test all(curl * curl .== ∇x∇x)
+    end
   end
 
   @testset "Costructor wavenumber propagationangle" begin