Add functionality for quiet start

Project.toml

@@ -4,13 +4,15 @@ authors = ["Luke Adams <[email protected]> and contributors"]
 version = "1.0.0"
 
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructEquality = "6ec83bb0-ed9f-11e9-3b4c-2b04cb4e219c"
 
 [compat]
+Distributions = "0.25.115"
 FFTW = "1"
 LinearSolve = "2.9"
 SparseArrays = "1.9"

src/ParticleInCell.jl

@@ -12,6 +12,7 @@ using StaticArrays
 using StructEquality
 using SparseArrays
 using LinearSolve
+using Distributions
 
 import Base: eachindex
 

src/species/species_test.jl

@@ -1,5 +1,6 @@
 @testitem "Species utilities" tags = [:unit] begin
     using StaticArrays
+    using Distributions
 
     @test electrons([1.0 2.0 3.0], [4.0 5.0 6.0], [1.0, 2.0, 3.0]) ==
           VariableWeightSpecies{1,1,Float64}(
@@ -36,4 +37,13 @@
         1.602176634e-19,
         9.1093837015e-31,
     )
+
+    dist = Normal(0.0, 1.0)
+    @test ParticleInCell.quiet_velocities(4, dist) ≈
+          [-1.1503493803760087, -0.3186393639643751, 0.3186393639643751, 1.1503493803760087]
+
+    # See Birdsall and Langdon, Table 16.5a
+    @test ParticleInCell.bit_reversed_sequence(6) == [0, 1 / 2, 1 / 4, 3 / 4, 1 / 8, 5 / 8]
+    @test ParticleInCell.bit_reversed_sequence(6, 3) ==
+          [0, 1 / 3, 2 / 3, 1 / 9, 4 / 9, 7 / 9]
 end

src/species/species_utils.jl

@@ -36,3 +36,28 @@ function electrons(positions::Vector{T}, momentums::Vector{T}, weight::T) where
     weights = fill(weight, length(positions))
     return electrons(positions, momentums, weights)
 end
+
+function quiet_velocities(num_macroparticles, distribution, T = Float64)
+    velocities = Vector{T}(undef, num_macroparticles)
+    for i = 1:num_macroparticles
+        velocities[i] = quantile(distribution, (i - 0.5) / num_macroparticles)
+    end
+    return velocities
+end
+
+function bit_reversed_sequence(n, base = 2, T = Float64)
+    seq = Vector{T}(undef, n)
+
+    m = round(Int, log(base, n), RoundUp)
+
+    for i = 0:n-1
+        j = 0
+        for (k, digit) in enumerate(digits(i, base = base, pad = m))
+            j += digit * base^(m - k)
+        end
+
+        seq[i+1] = j / base^m
+    end
+
+    return seq
+end

test/Project.toml

@@ -1,6 +1,7 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"