-
Notifications
You must be signed in to change notification settings - Fork 115
Expand file tree
/
Copy pathgravity_kepler_disks.py
More file actions
162 lines (134 loc) · 6.2 KB
/
Copy pathgravity_kepler_disks.py
File metadata and controls
162 lines (134 loc) · 6.2 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
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
from __future__ import print_function
from amuse.lab import *
#from amuse.io import store
#from amuse.community.seba.interface import SeBa
from amuse.community.fractalcluster.interface import new_fractal_cluster_model
# #BOOKLISTSTART2# #
def resolve_close_encounter(time, bodies, Johannes):
Johannes.initialize_from_particles(bodies)
rcom = bodies.center_of_mass()
vcom = bodies.center_of_mass_velocity()
a, e = Johannes.get_elements()
p = Johannes.get_periastron()
print("Close encounter at t=", time.in_(units.Myr), "a=", a.in_(units.AU), "e=", e, "p=", p.in_(units.AU), "M=", bodies.mass.max().in_(units.MSun), bodies.mass.min().in_(units.MSun), "at d=", rcom.in_(units.parsec), "with v=", vcom.in_(units.kms))
truncate_disks_due_to_encounter(bodies, p)
# #BOOKLISTSTOP2# #
# #BOOKLISTSTART3# #
def truncate_disks_due_to_encounter(bodies, p):
q = bodies[1].mass/bodies[0].mass
rtr_prim = 0.28*p / q**(0.32)
rtr_sec = 0.28*p * q**(0.32)
dm0 = truncate_disk_due_to_encounter(bodies[0], rtr_prim)
dm1 = truncate_disk_due_to_encounter(bodies[1], rtr_sec)
mtot = bodies.mass.sum()
bodies[0].accreted_mass += dm1 * bodies[0].mass/mtot
bodies[1].accreted_mass += dm0 * bodies[1].mass/mtot
bodies[0].radius = min(bodies[0].radius, 0.5*p)
bodies[1].radius = min(bodies[1].radius, 0.5*p)
bodies[0].mass = bodies[0].stellar_mass + \
bodies[0].disk_mass + bodies[0].accreted_mass
bodies[1].mass = bodies[1].stellar_mass + \
bodies[1].disk_mass + bodies[1].accreted_mass
# #BOOKLISTSTOP3# #
def stripped_disk_mass(body, dr):
rold = body.disk_radius
rnew = rold - dr
dm = body.disk_mass * (rold**0.5-rnew**0.5)/rold**0.5
return max(0|units.MSun, dm)
def truncate_disk_due_to_encounter(body, r_tr):
dr = max(0|units.AU, body.disk_radius-r_tr)
dm = stripped_disk_mass(body, dr)
body.disk_radius -= dr
body.disk_mass -= dm
return dm
def main(N, Rvir, Qvir, Fd):
filename= 'Cl_N%g_R%gpc_Q%g_F%g.h5'%(N, Rvir.value_in(units.parsec), Qvir, Fd)
t_end = 1.0 | units.Myr
dt = 0.1 | units.Myr
Mmax = 100 | units.MSun
masses = new_kroupa_mass_distribution(N, Mmax)
Mtot_init = masses.sum()
converter=nbody_system.nbody_to_si(Mtot_init,Rvir)
bodies = new_fractal_cluster_model(N=N, fractal_dimension=Fd,
convert_nbody=converter)
bodies.scale_to_standard(converter, virial_ratio=Qvir)
bodies.stellar_mass = masses
bodies.disk_mass = 0.1*bodies.stellar_mass
bodies.mass = bodies.stellar_mass + bodies.disk_mass
bodies.accreted_mass = 0 | units.MSun
bodies.disk_radius = 400 | units.AU
bodies.radius = 10 * bodies.disk_radius
gravity = ph4(converter)
gravity.parameters.epsilon_squared = (100|units.AU)**2
gravity.particles.add_particles(bodies)
channel_from_gd_to_framework = gravity.particles.new_channel_to(bodies)
channel_from_framework_to_gd = bodies.new_channel_to(gravity.particles)
stopping_condition = gravity.stopping_conditions.collision_detection
stopping_condition.enable()
Johannes = Kepler(converter)
Johannes.initialize_code()
write_set_to_file(bodies.savepoint(0|units.Myr), filename, 'hdf5', append_to_file=False)
Etot_init = gravity.kinetic_energy + gravity.potential_energy
Etot_prev = Etot_init
Nenc = 0
dEk_enc = zero
dEp_enc = zero
time = 0.0 | t_end.unit
# #BOOKLISTSTART1# #
while time < t_end:
time += dt
gravity.evolve_model(time)
Etot_prev_se = gravity.kinetic_energy + gravity.potential_energy
while stopping_condition.is_set():
channel_from_gd_to_framework.copy()
Ek_enc = gravity.kinetic_energy
Ep_enc = gravity.potential_energy
for ci in range(len(stopping_condition.particles(0))):
particles_in_encounter = Particles(
particles=[stopping_condition.particles(0)[ci],
stopping_condition.particles(1)[ci]])
local_particles_in_encounter = particles_in_encounter.get_intersecting_subset_in(bodies)
resolve_close_encounter(gravity.model_time,
local_particles_in_encounter, Johannes)
Nenc+=1
print(f"At time= {gravity.model_time} Nenc= {Nenc} ",
(f"Rdisk= {local_particles_in_encounter.disk_radius}")
channel_from_framework_to_gd.copy_attributes(["radius"])
dEk_enc += Ek_enc - gravity.kinetic_energy
dEp_enc += Ep_enc - gravity.potential_energy
gravity.evolve_model(time)
channel_from_framework_to_gd.copy_attributes(["mass"])
# #BOOKLISTSTOP1# #
write_set_to_file(bodies.savepoint(time), filename, 'hdf5')
Ekin = gravity.kinetic_energy
Epot = gravity.potential_energy
Etot = Ekin + Epot
dE = Etot_prev-Etot
dE_se = Etot_prev_se-Etot
Mtot = bodies.mass.sum()
print("T=", time, end=' ')
print("M=", Mtot, "(dM[SE]=", Mtot/Mtot_init, ")", end=' ')
print("E= ", Etot, "Q= ", Ekin/Epot, end=' ')
print("dE=", (Etot_init-Etot)/Etot, "ddE=", (Etot_prev-Etot)/Etot, end=' ')
print("(dE[SE]=", dE_se/Etot, ")")
print("dE(enc)=", dEk_enc, dEp_enc)
Etot_init -= dE
Etot_prev = Etot
gravity.stop()
Johannes.stop()
def new_option_parser():
from amuse.units.optparse import OptionParser
result = OptionParser()
result.add_option("-N", dest="N", type="int",default = 2000,
help="number of stars [%default]")
result.add_option("-R", dest="Rvir", type="float",
unit=units.parsec, default = 0.5|units.parsec,
help="cluser virial radius [%default]")
result.add_option("-Q", dest="Qvir", type="float",default = 0.5,
help="virial ratio [%default]")
result.add_option("-F", dest="Fd", type="float",default = 1.6,
help="fractal dimension [%default]")
return result
if __name__ in ('__main__', '__plot__'):
o, arguments = new_option_parser().parse_args()
main(**o.__dict__)