fix random vel, rejection sampling

src/Tools.jl

@@ -92,7 +92,8 @@ function rotational_velocity_acc(x, y, z, a, ratio = 1.0)
     r = sqrt(x^2 + y^2 + z^2)
     v = sqrt(a * r)
     v_vec = -normalize(PVector(ustrip(u,x), ustrip(u,y), 0.0) × PVector(0.0, 0.0, 1.0)) * v
-    v_vec = ratio * v_vec + (1-ratio) * randn(PVector{Float64}) * v
+    v_vec = ratio * v_vec + (1-ratio) * randn(PVector{Float64}) * v # this will change the average velocity
+    v_vec = v_vec + (1-ratio) * randn(PVector{Float64}) * v
 end
 
 """
@@ -103,5 +104,6 @@ Rotate along the positive z-axis.
 function rotational_velocity(x, y, v, ratio = 1.0)
     u = unit(x)
     v_vec = -normalize(PVector(ustrip(u,x), ustrip(u,y), 0.0) × PVector(0.0, 0.0, 1.0)) * v
-    v_vec = ratio * v_vec + (1-ratio) * randn(PVector{Float64}) * v
+    # v_vec = ratio * v_vec + (1-ratio) * randn(PVector{Float64}) * v # this will change the average velocity
+    v_vec = v_vec + (1-ratio) * randn(PVector{Float64}) * v
 end

src/bulge.jl

@@ -56,8 +56,7 @@ $(TYPEDSIGNATURES)
 """
 function generate(config::Bulge, units = uAstro;
     RotationCurve = nothing,
-    MaxRadius = 5 * config.ScaleRadius,
-    # MaxHeight = 5 * config.ScaleRadius * config.q,
+    MaxRadius = 20 * config.ScaleRadius,
     MaxHeight = MaxRadius,
 )
     uLen = getuLength(units)
@@ -70,7 +69,9 @@ function generate(config::Bulge, units = uAstro;
     end
 
     R = eltype(config.ScaleRadius)[]
+    sizehint!(R, NumSamples)
     z = eltype(config.ScaleRadius)[]
+    sizehint!(z, NumSamples)
 
     target(xy) = -pdf(config,xy[1],xy[2])
     pdf_maximum = -minimum_func(target, [ustrip(config.ScaleRadius), ustrip(config.ScaleRadius)*config.q])[1]

src/disk.jl

@@ -46,7 +46,7 @@ function pdf(config::ExponentialDisc, R, z)
     if iszero(config.HoleRadius)
         return exp(-abs(z)/ustrip(config.ScaleHeight) - R/ustrip(config.ScaleRadius)) / ustrip(config.ScaleHeight) * 2π * R
     else
-        return exp(-ustrip(config.HoleRadius)/R - R/ustrip(config.ScaleRadius)) * sech(z/ustrip(config.ScaleHeight)/2)^2 / 4 / ustrip(config.ScaleHeight) * 2π * R
+        return exp(-ustrip(config.HoleRadius)/R - R/ustrip(config.ScaleRadius)) * sech(z/ustrip(config.ScaleHeight)/2)^2 * 2π * R
     end
 end
 
@@ -56,7 +56,7 @@ $(TYPEDSIGNATURES)
 """
 function generate(config::ExponentialDisc, units = uAstro;
     RotationCurve = nothing,
-    MaxRadius = 5 * config.ScaleRadius,
+    MaxRadius = 20 * config.ScaleRadius,
     MaxHeight = MaxRadius,
     k = 2, # degree of rotation curve interpolation/extrapolation spline (1 = linear, 2 = quadratic, 3 = cubic, up to 5)
     bc = "nearest", # behavior when evaluating the spline outside the support domain, which is (minimum(x), maximum(x)). The allowed values are "nearest", "zero", "extrapolate", "error"
@@ -110,10 +110,17 @@ function generate(config::ExponentialDisc, units = uAstro;
     # end
 
     R = eltype(config.ScaleRadius)[]
+    sizehint!(R, NumSamples)
     z = eltype(config.ScaleRadius)[]
+    sizehint!(z, NumSamples)
 
     target(xy) = -pdf(config,xy[1],xy[2])
-    pdf_maximum = -minimum_func(target, [ustrip(config.ScaleRadius), ustrip(config.ScaleRadius)])[1]
+    if iszero(config.HoleRadius)
+        pdf_maximum = -minimum_func(target, [ustrip(config.ScaleRadius), ustrip(config.ScaleRadius)])[1]
+    else
+        pdf_maximum = pdf(config, ustrip(sqrt(config.ScaleRadius * config.HoleRadius)), 0.0)
+    end
+    # @show pdf_maximum
 
     # rejection sampling
     while length(R) < NumSamples

src/distribution.jl

@@ -5,9 +5,10 @@ Use `Optim.jl` to find the minimum value of `target` function.
 """
 function minimum_func(target::Function, ic)
     result = Optim.optimize(target, ic)
-    maximum_value = Optim.minimum(result)
+    # @show result
+    minimum_value = Optim.minimum(result)
     optimal_xy = Optim.minimizer(result)
-    return maximum_value, optimal_xy
+    return minimum_value, optimal_xy
 end
 
 """
@@ -82,4 +83,31 @@ $(TYPEDSIGNATURES)
 """
 function sech2_cdf_inv(u)
     return atanh(2*u-1)
+end
+
+"""
+$(TYPEDSIGNATURES)
+
+Sample from the discrete pdf using rejection method
+
+!!! attentions
+    - Data must be unitless, and have no NaN
+    - `x` should start from 0
+    - the maximum of `pdf` would be the width of sampling box
+"""
+function rejection_sampling(x, pdf, rMax, NumSamples)
+    spl = Spline1D(x, pdf)
+    pdf_maximum = maximum(pdf)
+
+    R = eltype(pdf)[]
+    sizehint!(R, NumSamples)
+
+    while length(R) < NumSamples
+        R_rand = rand() * rMax
+        if rand() < spl(R_rand) / pdf_maximum
+            push!(R, R_rand)
+        end
+    end
+
+    return R
 end