Skip to content

Commit 95fa4f2

Browse files
islentislent
committed
init example 08-TrajectoryLookback
1 parent 6417c35 commit 95fa4f2

File tree

1 file changed

+218
-0
lines changed

1 file changed

+218
-0
lines changed

‎examples/08-TrajectoryLookback/08-TrajectoryLookback.jl‎

Lines changed: 218 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,218 @@
1+
# In this example, we lookback the trajectories of MW's satellites
2+
# The original src/data/MW_satellites_origin.csv are extracted from battaglia2022gaia_Astron. Astrophys. - Gaia early DR3 systemic motions of local group dwarf galaxies and orbital properties with a massive Large Magellanic Cloud
3+
# and converted to Galactocentric coordinates using astropy.
4+
# The parameters of LMC are from vasiliev2021tango_MNRAS - Tango for three sagittarius, LMC, and the milky way
5+
6+
using Unitful, UnitfulAstro
7+
using StructArrays
8+
using DataFrames, CSV
9+
using GLMakie
10+
using ProgressMeter
11+
12+
using AstroIC
13+
using AstroNbodySim
14+
using WaveDM
15+
16+
using PhysicalParticles
17+
18+
19+
df_MW_satellites = AstroIC.load_data_MW_satellites() # load from src/data/MW_satellites.csv
20+
21+
22+
##### General parameters
23+
24+
Np = 1000
25+
# Np = 5000
26+
# Np = 10000
27+
# Np = 50000
28+
GravitySolver = Tree()
29+
# GravitySolver = DirectSum()
30+
pids = [1]
31+
32+
TimeEnd = 6.0u"Gyr"
33+
TimeStep = 0.0u"Gyr" # adaptive if zero.
34+
# TimeStep = 1e-5u"Gyr"
35+
TimeBetweenSnapshots = 0.1u"Gyr"
36+
37+
file_LMC = joinpath(@__DIR__, "output/LMC_traj_lookback_Vasiliev2021", "analysis.csv")
38+
flag_load_LMC = false # If false, re-compute the trajectory of LMC
39+
# flag_load_LMC = true # If true and the file exists, skip trajectory lookback of LMC
40+
Np_LMC = 200
41+
# Np_LMC = 1000
42+
# Np_LMC = 10000
43+
44+
##### Initialize MW model
45+
@info "Initializing MW's baryons"
46+
#TODO: WaveDM.jl
47+
baryon_particles = test_MW_MOND(;
48+
V = (x,y,z,ψ)->0.0,
49+
Nx = 32,
50+
Np,
51+
baryon_mode = :particles_static,
52+
export_particles = true,
53+
);
54+
55+
sim_force_baryon = Simulation(baryon_particles;
56+
GravitySolver,
57+
pids,
58+
);
59+
60+
61+
@info "Initializing MW's halo"
62+
MW_model_halo = Zhao(1.55e7u"Msun/kpc^3", 11.75u"kpc", 1.19, 2.95, 0.95)
63+
MW_ρ_halo_func = (r)->GalacticDynamics.density(MW_model_halo, r)
64+
MW_table_r = collect(0.1:0.1:100) * u"kpc";
65+
MW_ρ_halo = MW_ρ_halo_func.(MW_table_r);
66+
MW_table_mass_halo = 4π * cumul_integrate(MW_table_r, MW_table_r.^2 .* MW_ρ_halo);
67+
MW_table_acc_halo = -C.G * MW_table_mass_halo ./ MW_table_r.^2;
68+
MW_table_pot_halo = -C.G * MW_table_mass_halo ./ MW_table_r;
69+
MW_table_pot_halo[1:end-2] .-= [C.G * 4π * NumericalIntegration.integrate(MW_table_r[i:end], MW_ρ_halo[i:end] .* MW_table_r[i:end]) for i in eachindex(MW_table_r)[1:end-2]];
70+
spl_acc = Spline1D(ustrip.(u"kpc", MW_table_r), ustrip.(u"kpc/Gyr^2", MW_table_acc_halo); k=1)
71+
spl_pot = Spline1D(ustrip.(u"kpc", MW_table_r), ustrip.(u"kpc^2/Gyr^2", MW_table_pot_halo); k=1)
72+
73+
74+
### MW background force without LMC
75+
function MW_bg_force(p::AbstractParticle, t;
76+
SofteningLength = 0.01u"kpc",
77+
GravitySolver = Tree(),
78+
# GravitySolver = DirectSum(),
79+
sim_force_baryon = sim_force_baryon,
80+
spl_acc = spl_acc,
81+
)
82+
acc_DM = spl_acc(ustrip(u"kpc", norm(p.Pos))) * ustrip(normalize(p.Pos)) * u"kpc/Gyr^2"
83+
acc_b = compute_force(sim_force_baryon, p.Pos, SofteningLength, GravitySolver, CPU())[1]
84+
return acc = acc_b + acc_DM
85+
end
86+
87+
p1x(sim::Simulation) = find_particle(sim, 1).Pos.x
88+
p1y(sim::Simulation) = find_particle(sim, 1).Pos.y
89+
p1z(sim::Simulation) = find_particle(sim, 1).Pos.z
90+
91+
analysers = Dict(
92+
"x" => p1x,
93+
"y" => p1y,
94+
"z" => p1z,
95+
)
96+
97+
98+
99+
##### LMC
100+
if flag_load_LMC && isfile(file_LMC)
101+
@info "Loading LMC trajectory from $(file_LMC)"
102+
df_traj_LMC = DataFrame(CSV.File(file_LMC));
103+
else
104+
@info "Simulating LMC..."
105+
### Vasiliev 2021
106+
model_LMC = tNFW(5.2629e6u"Msun/kpc^3", 8.5u"kpc" * 1.38^0.6, 85u"kpc" * 1.38^0.6)
107+
config_LMC = SphericalSystem(STAR, Np_LMC,
108+
x -> GalacticDynamics.density(model_LMC, x) * 4π * x^2,
109+
)
110+
particles_LMC = generate(config_LMC;
111+
MaxRadius = 85.0u"kpc" * 1.38^0.6 * 1.5,
112+
);
113+
particles_LMC.Mass .= 1.38e11u"Msun"/Np_LMC; # make sure the total mass is 1.38e11 Msun
114+
115+
116+
tidal_initial_pos = PVector(-0.6, -41.3, -27.1, u"kpc") # Galactocentric
117+
tidal_initial_vel = -PVector(-63.9, -213.8, 206.6, u"km/s") #? Negative for lookback
118+
119+
test_particles = StructArray([
120+
Star(uAstro; id = 1),
121+
])
122+
test_particles.Pos[1] = tidal_initial_pos
123+
test_particles.Vel[1] = tidal_initial_vel
124+
125+
sim_traj_lookback = Simulation(
126+
deepcopy(test_particles);
127+
bgforce = Function[MW_bg_force],
128+
analysers,
129+
TimeEnd = 6.0u"Gyr",
130+
TimeStep = 0.0u"Gyr", # adaptive if zero.
131+
# TimeStep = 1e-5u"Gyr",
132+
TimeBetweenSnapshots = 0.1u"Gyr",
133+
Realtime = false,
134+
OutputDir = "output/LMC_traj_lookback_Vasiliev2021",
135+
);
136+
run(sim_traj_lookback);
137+
138+
139+
df_traj_LMC = DataFrame(CSV.File(joinpath("output/LMC_traj_lookback_Vasiliev2021", "analysis.csv")));
140+
end
141+
142+
143+
function MW_bg_force_LMC(p::AbstractParticle, t;
144+
SofteningLength = 0.01u"kpc",
145+
GravitySolver = Tree(),
146+
# GravitySolver = DirectSum(),
147+
sim_force_baryon = sim_force_baryon,
148+
spl_acc = spl_acc,
149+
df_traj = df_traj_LMC,
150+
particles_LMC = particles_LMC,
151+
)
152+
acc_DM = spl_acc(ustrip(u"kpc", norm(p.Pos))) * ustrip(normalize(p.Pos)) * u"kpc/Gyr^2"
153+
acc_b = compute_force(sim_force_baryon, p.Pos, SofteningLength, GravitySolver, CPU())[1]
154+
155+
if ustrip(u"Gyr", t) <= minimum(df_traj.time)
156+
id = 1
157+
else
158+
id = findfirstvalue(df_traj.time, ustrip(u"Gyr", t))
159+
end
160+
161+
particles_LMC_t = deepcopy(particles_LMC)
162+
particles_LMC_t.Pos .+= PVector(df_traj.x[id], df_traj.y[id], df_traj.z[id]) * u"kpc"
163+
acc_LMC = AstroNbodySim.compute_force_at_point(p.Pos, particles_LMC_t, C.G, SofteningLength) #? DirectSum for LMC particles
164+
return acc = acc_b + acc_DM + acc_LMC
165+
end
166+
167+
168+
@showprogress for i in eachindex(df_MW_satellites.Galaxy)
169+
@info "Simulating galaxy $(df_MW_satellites.Galaxy[i]) without LMC"
170+
tidal_initial_pos = PVector(df_MW_satellites.X[i], df_MW_satellites.Y[i], df_MW_satellites.Z[i]) # Galactocentric
171+
tidal_initial_vel = -PVector(df_MW_satellites.v_X[i], df_MW_satellites.v_Y[i], df_MW_satellites.v_Z[i]) #? Negative for lookback
172+
@info tidal_initial_pos
173+
@info tidal_initial_vel
174+
175+
test_particles = StructArray([
176+
Star(uAstro; id = 1),
177+
])
178+
test_particles.Pos[1] = tidal_initial_pos
179+
test_particles.Vel[1] = tidal_initial_vel
180+
181+
sim_traj_lookback = Simulation(
182+
deepcopy(test_particles);
183+
bgforce = Function[MW_bg_force],
184+
analysers,
185+
TimeEnd,
186+
TimeStep,
187+
TimeBetweenSnapshots,
188+
Realtime = false,
189+
OutputDir = joinpath(@__DIR__, "output/traj_lookback_$(df_MW_satellites.Galaxy[i])_no_LMC"),
190+
);
191+
run(sim_traj_lookback);
192+
193+
194+
@info "Simulating galaxy $(df_MW_satellites.Galaxy[i]) with LMC"
195+
tidal_initial_pos = PVector(df_MW_satellites.X[i], df_MW_satellites.Y[i], df_MW_satellites.Z[i]) # Galactocentric
196+
tidal_initial_vel = -PVector(df_MW_satellites.v_X[i], df_MW_satellites.v_Y[i], df_MW_satellites.v_Z[i]) #? Negative for lookback
197+
198+
test_particles = StructArray([
199+
Star(uAstro; id = 1),
200+
])
201+
test_particles.Pos[1] = tidal_initial_pos
202+
test_particles.Vel[1] = tidal_initial_vel
203+
204+
sim_traj_lookback = Simulation(
205+
deepcopy(test_particles);
206+
bgforce = Function[MW_bg_force_LMC],
207+
analysers,
208+
TimeEnd,
209+
TimeStep,
210+
TimeBetweenSnapshots,
211+
Realtime = false,
212+
OutputDir = joinpath(@__DIR__, "output/traj_lookback_$(df_MW_satellites.Galaxy[i])_with_LMC"),
213+
);
214+
run(sim_traj_lookback);
215+
216+
#TODO release memory of simulations
217+
end
218+

0 commit comments

Comments
 (0)