Check the sampling again

Kolaru · Nov 23, 2023 · 0a49d9f · 0a49d9f
1 parent 446dcd6
commit 0a49d9f
Showing 1 changed file with 6 additions and 5 deletions.
diff --git a/src/NormalModes.jl b/src/NormalModes.jl
@@ -154,13 +154,14 @@ Return the deviation from the average geometry.
 """
 function StatsBase.sample(nm::NormalDecomposition, n_samples)
     hbar = 1  # Atomic units
-    m = 1  # All modes have mass 1
-    Δz_dist = MvNormal(Diagonal(1/2 * (hbar ./ (m * nm.ωs))))
-    Δp_dist = MvNormal(Diagonal(1/2 * hbar * m * nm.ωs))
+    Δz_dist = MvNormal(Diagonal(1/2 * (hbar ./ nm.ωs)))
+    Δp_dist = MvNormal(Diagonal(1/2 * hbar * nm.ωs))
     MU = nm.M * nm.U
 
-    Δx = MU*rand(Δz_dist, n_samples) 
-    Δp = nm.M.^2 * MU*rand(Δp_dist, n_samples)
+    Δx = MU * rand(Δz_dist, n_samples)
+
+    @warning "I'm really not sure about the sampling in momentum space"
+    Δp = nm.M.^-2 * MU * rand(Δp_dist, n_samples)
     return Δx * aunit(u"m"), Δp * aunit(u"kg*m/s")
 end