Adds default parameter values to Bovy2014

cadojo · Dec 2, 2024 · 5fdf88e · 5fdf88e
1 parent 3aacb9f
commit 5fdf88e
Show file tree

Hide file tree

Showing 2 changed files with 36 additions and 11 deletions.
diff --git a/src/potentials.jl b/src/potentials.jl
@@ -257,14 +257,40 @@ A potential field for the Milky Way galaxy, based off of Dr. Bovy's 2015 paper.
 """
 function Bovy2014(; name = :BovyMilkyWayPotential, kwargs...)
 
+    # gala 
     # default_disk = dict(m=68193902782.346756 * u.Msun, a=3.0 * u.kpc, b=280 * u.pc)
     # default_bulge = dict(m=4501365375.06545 * u.Msun, alpha=1.8, r_c=1.9 * u.kpc)
     # default_halo = dict(m=4.3683325e11 * u.Msun, r_s=16 * u.kpc)
 
+    # galpy
+    # bp= PowerSphericalPotentialwCutoff(alpha=1.8,rc=1.9/8.,normalize=0.05)
+    # mp= MiyamotoNagaiPotential(a=3./8.,b=0.28/8.,normalize=.6)
+    # np= NFWPotential(a=16/8.,normalize=.35)
+    # MWPotential2014= bp+mp+np
+
     disk = MiyamotoNagaiPotential(; name = :Disk)
     bulge = PowerLawCutoffPotential(; name = :Bulge)
     halo = NFWPotential(; name = :Halo)
 
+    G = 6.6743e-11
+
+    defaults = [
+        disk.a => 3,
+        disk.b => 280,
+        disk.G => G,
+        disk.m => 68193902782.346756,
+        bulge.c => 1.9,
+        bulge.m => 4501365375.06545,
+        bulge.α => 1.8,
+        bulge.G => G,
+        halo.a => 1.0,
+        halo.b => 1.0,
+        halo.c => 1.0,
+        halo.r => 16,
+        halo.m => 4.3683325e11,
+        halo.G => G
+    ]
+
     aliases = [
         disk.x => x,
         disk.y => y,
@@ -280,7 +306,7 @@ function Bovy2014(; name = :BovyMilkyWayPotential, kwargs...)
     u = [x, y, z]
     du = [ẋ, ẏ, ż]
 
-    grad(sys) = Symbolics.gradient(sys.eqs[1].rhs, [x, y, z]) # TODO remove manual indexing
+    grad(sys) = calculate_jacobian(sys; simplify = true)[(begin + 1):end] # TODO remove manual indexing
 
     eqs = vcat(
         Φ ~ disk.Φ + bulge.Φ + halo.Φ,
@@ -291,7 +317,8 @@ function Bovy2014(; name = :BovyMilkyWayPotential, kwargs...)
 
     return compose(
         ODESystem(
-            eqs, t; name = name, defaults = Dict(aliases)),
+            eqs, t;
+            name = name, defaults = Dict(vcat(defaults, aliases))),
         disk, bulge, halo
     )
 end

diff --git a/test/composition.jl b/test/composition.jl
@@ -18,7 +18,7 @@ function PlummerPotential(; name = :PlummerPotential, kwargs...)
     return ScalarField(value, t, u, p; name = name, kwargs...)
 end
 
-function MilkyWay(; name = :MilkyWay, alias_equations = false, kwargs...)
+function MilkyWay(; name = :MilkyWay, kwargs...)
     @named disk = PlummerPotential()
     @named bulge = PlummerPotential()
     @named halo = PlummerPotential()
@@ -41,23 +41,21 @@ function MilkyWay(; name = :MilkyWay, alias_equations = false, kwargs...)
         halo.z => z
     ]
 
-    grad(sys) = ModelingToolkit.gradient(sys.Φ, [sys.x, sys.y, sys.z])
+    grad(sys) = calculate_jacobian(sys; simplify = true)[(begin + 1):end] # TODO remove manual indexing
 
     eqs = vcat(
         Φ ~ disk.Φ + bulge.Φ + halo.Φ,
         D.(u) .~ du,
-        D.(du) .~ -(grad(disk) .+ grad(bulge) .+ grad(halo))
+        D.(du) .~ -(grad(disk) .+ grad(bulge) .+ grad(halo)),
+        [alias.first ~ alias.second for alias in aliases]
     )
 
-    if alias_equations
-        append!(eqs, [alias.first ~ alias.second for alias in aliases])
-    end
-
     return compose(
         ODESystem(
-            eqs, t; name = :MilkyWay, defaults = Dict(aliases)),
+            eqs, t; name = name, defaults = Dict(aliases)),
         disk, bulge, halo
     )
 end
 
-mw = MilkyWay(alias_equations = true)
+mw = MilkyWay()
+mws = structural_simplify(mw)