Repository navigation
Expand file tree
/
Copy pathmerge_two_stars.py
More file actions
93 lines (77 loc) · 2.34 KB
/
Copy pathmerge_two_stars.py
File metadata and controls
93 lines (77 loc) · 2.34 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
"""
Initialize two stars to a certain age and merge them using MMAMS
"""
import sys
import argparse
import matplotlib.pyplot as plt
from amuse.datamodel import Particle
from amuse.units import units
from amuse.community.mesa import Mesa
from amuse.community.mmams import Mmams
# #BOOKLISTSTART# #
def merge_two_stars(
mass_primary=5 | units.MSun,
mass_secondary=3 | units.MSun,
time_collision=1.0 | units.Myr,
):
"""
Merges two stars using the Make Me a Massive Star code and returns the
resulting density profile.
"""
primary = Particle(mass=mass_primary)
secondary = Particle(mass=mass_secondary)
stellar = Mesa()
primary = stellar.particles.add_particle(primary)
secondary = stellar.particles.add_particle(secondary)
stellar.evolve_model(time_collision)
stellar.merge_colliding(
primary.copy(), secondary.copy(), Mmams, return_merge_products=["se"]
)
radius = stellar.particles[0].get_radius_profile()
rho = stellar.particles[0].get_density_profile()
stellar.stop()
plot_density_profile(radius, rho)
# #BOOKLISTSTOP# #
def plot_density_profile(radius, rho):
"""
Plots density against radius.
"""
figure = plt.figure()
ax = figure.add_subplot(1, 1, 1)
ax.plot(radius.value_in(units.RSun), rho.value_in(units.g / units.cm**3))
ax.set_xlabel(r"$R$ [$R_\odot$]")
ax.set_ylabel("density [$g/cm^3$]")
ax.set_yscale("log", nonpositive="clip")
plt.savefig("merge_two_stars.pdf")
print("Saved figure in file merge_two_stars.pdf")
def new_argument_parser():
result = argparse.ArgumentParser(
formatter_class=argparse.ArgumentDefaultsHelpFormatter,
)
result.add_argument(
"-M",
"--mass_primary",
type=units.MSun,
default=5 | units.MSun,
help="Mass of the primary star",
)
result.add_argument(
"-m",
"--mass_secondary",
type=units.MSun,
default=3 | units.MSun,
help="Mass of the secondary star",
)
result.add_argument(
"-t",
"--time_collision",
type=units.Myr,
default=1.0 | units.Myr,
help="end time of the simulation",
)
return result
def main(**kwargs):
merge_two_stars(**kwargs)
if __name__ == "__main__":
arguments = new_argument_parser().parse_args()
main(**arguments.__dict__)