fix cylindrical 2*pi*R pdf

Author: islent

.github/workflows/CI.yml

@@ -11,7 +11,7 @@ jobs:
     strategy:
       matrix:
         version:
-          - '1.10.4'
+          - '1.11.1'
         os:
           - ubuntu-latest
           - macOS-latest

.gitignore

@@ -1,2 +1,3 @@
 Manifest.toml
-docs/build
+docs/build
+test/output

Project.toml

@@ -1,7 +1,7 @@
 name = "AstroIC"
 uuid = "5b3ca070-48c3-41ea-b469-89825dd01fde"
 authors = ["islent <[email protected]>"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 AstroIO = "c85a633c-0c3f-44a2-bffe-7f9d0681b3e7"
@@ -12,29 +12,40 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab"
 PhysicalParticles = "78724250-4eb2-11e9-0f34-13a7aab864c5"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
 
 [compat]
 AstroIO = "0.1"
 AstroLib = "0.4"
 AstroSimBase = "0.1"
-BangBang = "0.3, 0.4"
+BangBang = "0.3, 0.4, 0.5"
 Dierckx = "0.5"
 Distributions = "0.25"
 DocStringExtensions = "0.8, 0.9"
+Optim = "1.9.4"
 PhysicalConstants = "0.2"
 PhysicalParticles = "1"
+Plots = "1.40.8"
 PrecompileTools = "1"
+QuadGK = "2"
 Reexport = "1"
 StructArrays = "0.6"
 Unitful = "1"
 UnitfulAstro = "1"
 julia = "1.6"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+test = ["Test"]

Readme.md

@@ -51,7 +51,7 @@ config = PlummerStarCluster(
     NumSamples = 100,
     VirialRadius = 0.010u"kpc",
     TotalMass = 1.0e5u"Msun",
-    model = AstroIC.Newton(),
+    model = AstroSimBase.Newton(),
 )
 
 ## Now generate particles. MaxRadius restricts the sampling region.

docs/src/manual/plummer.md

@@ -15,7 +15,7 @@ config = PlummerStarCluster(
     NumSamples = 1000,
     VirialRadius = 0.010u"kpc",
     TotalMass = 1.0e5u"Msun",
-    model = AstroIC.Newton(),
+    model = AstroSimBase.Newton(),
 )
 
 # generate

src/AstroIC.jl

@@ -12,15 +12,18 @@ using BangBang
 using Dates
 using StructArrays
 using Dierckx
+using QuadGK
+using Optim
 
 @reexport using PhysicalParticles
 @reexport using AstroIO
 
 import Base: show
 import Unitful: Units
+import Distributions: pdf
 import PhysicalConstants: CODATA2018
 import PhysicalParticles: rotate, rotate_x, rotate_y, rotate_z
-using AstroSimBase
+@reexport using AstroSimBase
 
 using AstroLib
 import AstroLib: planets
@@ -56,6 +59,7 @@ _common_keywords = """
 
 include("Tools.jl")
 include("physics.jl")
+include("distribution.jl")
 
 include("plummer.jl")
 include("disk.jl")

src/Tools.jl

@@ -80,4 +80,28 @@ function rand_pos_3d(r::Number)
     y = r * sin(acos(2.0eta-1.0)) * sin(2.0pi*gamma)
     z = r * (2.0eta-1.0)
     return PVector(x, y, z)
+end
+
+"""
+$(TYPEDSIGNATURES)
+Rotate along the positive z-axis.
+- `ratio`: rotational motion v.s. random motion. If equals `1`, only rotational component
+"""
+function rotational_velocity_acc(x, y, z, a, ratio = 1.0)
+    u = unit(x)
+    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
+end
+
+"""
+$(TYPEDSIGNATURES)
+Rotate along the positive z-axis.
+- `ratio`: rotational motion v.s. random motion. If equals `1`, only rotational component
+"""
+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
 end

src/bulge.jl

@@ -45,13 +45,19 @@ function Base.show(io::IO, config::Bulge)
     )
 end
 
+function pdf(config::Bulge, R, z)
+    r_prime = sqrt(R^2 + (z / config.q)^2)
+    return exp(-r_prime^2/ustrip(config.CutRadius)^2) / (1 + r_prime/ustrip(config.ScaleRadius))^config.α * 2π*R
+end
+
 """
 $(TYPEDSIGNATURES)
 
 """
 function generate(config::Bulge, units = uAstro;
     RotationCurve = nothing,
     MaxRadius = 5 * config.ScaleRadius,
+    # MaxHeight = 5 * config.ScaleRadius * config.q,
     MaxHeight = MaxRadius,
 )
     uLen = getuLength(units)
@@ -65,17 +71,30 @@ function generate(config::Bulge, units = uAstro;
 
     R = eltype(config.ScaleRadius)[]
     z = eltype(config.ScaleRadius)[]
-    
+
+    target(xy) = -pdf(config,xy[1],xy[2])
+    pdf_maximum = -minimum_func(target, [ustrip(config.ScaleRadius), ustrip(config.ScaleRadius)*config.q])[1]
+
+    # rejection sampling
     while length(R) < NumSamples
-        R_rand = randn() * config.CutRadius / sqrt(2)
-        z_rand = randn() * config.q * config.CutRadius / sqrt(2)
-        r_prime = sqrt(R_rand^2 + (z_rand/config.q)^2)
-        if rand() <= (1 + r_prime / config.ScaleRadius)^(-config.α)
-            push!(R, abs(R_rand))
-            push!(z, z_rand)
+        R_rand = rand() * ustrip(MaxRadius)
+        z_rand = rand() * 2*ustrip(MaxHeight) - ustrip(MaxHeight)
+        if rand() < pdf(config, R_rand, z_rand) / pdf_maximum
+            push!(R, R_rand * uLen)
+            push!(z, z_rand * uLen)
         end
     end
 
+    # while length(R) < NumSamples
+    #     R_rand = randn() * config.CutRadius / sqrt(2)
+    #     z_rand = randn() * config.q * config.CutRadius / sqrt(2)
+    #     r_prime = sqrt(R_rand^2 + (z_rand/config.q)^2)
+    #     if rand() <= (1 + r_prime / config.ScaleRadius)^(-config.α)
+    #         push!(R, abs(R_rand))
+    #         push!(z, z_rand)
+    #     end
+    # end
+
     pos2d = StructArray(rand_pos_2d.(R))
     pos = StructArray(PVector.(pos2d.x, pos2d.y, z))