66import matplotlib .pyplot as plt
77from amuse .units import nbody_system
88from amuse .ic .plummer import new_plummer_model
9- from amuse .community .bhtree import Bhtree
109from amuse .community .huayno import Huayno
11- from amuse .community .kepler import Kepler
1210from amuse .community .ph4 import Ph4
11+ from amuse .community .bhtree import Bhtree
1312
1413
15- def virial_ratio_evolution (code , bodies , Q_init , t_end ):
16- dt = 0.06125 | t_end .unit
14+ def virial_ratio_evolution (code , bodies , Q_init , time_end ):
15+ dt = 0.06125 | time_end .unit
1716 bodies .scale_to_standard (virial_ratio = Q_init )
1817 bodies .radius = 0 | nbody_system .length
1918 gravity = code ()
@@ -22,9 +21,9 @@ def virial_ratio_evolution(code, bodies, Q_init, t_end):
2221 channel_from_gravity_to_framework = gravity .particles .new_channel_to (bodies )
2322
2423 Etot_prev = Etot_init = gravity .kinetic_energy + gravity .potential_energy
25- time = [0.0 ] | t_end .unit
24+ time = [0.0 ] | time_end .unit
2625 Q = [Q_init ]
27- while time [- 1 ] < t_end :
26+ while time [- 1 ] < time_end :
2827 time .append (time [- 1 ] + dt )
2928 gravity .evolve_model (time [- 1 ])
3029 channel_from_gravity_to_framework .copy ()
@@ -40,36 +39,45 @@ def virial_ratio_evolution(code, bodies, Q_init, t_end):
4039 return time , Q
4140
4241
43- def gravity_to_virial (N , t_end ):
42+ def gravity_to_virial (number_of_stars , time_end ):
4443 Q_init = 0.2
45- particles = new_plummer_model (N )
44+ particles = new_plummer_model (number_of_stars )
4645 codes = [Ph4 , Huayno , Bhtree ]
4746 ci = 0
4847 x_label = "time [N-body units]"
4948 # y_label = "$Q [\equiv -E_{\rm kin}/E_{\rm pot}]$"
5049 y_label = "virial ratio $Q$"
51- figure = plt .figure (figsize = ( 14 , 10 ) )
50+ figure = plt .figure ()
5251 ax1 = figure .add_subplot (1 , 1 , 1 )
5352 ax1 .set_xlabel (x_label )
5453 ax1 .set_ylabel (y_label )
55- ax1 .set_xlim (0 , t_end .value_in (t_end .unit ))
54+ ax1 .set_xlim (0 , time_end .value_in (time_end .unit ))
5655 ax1 .set_ylim (0 , 0.65 )
57- plt .plot ([0 , t_end .value_in (t_end .unit )], [0.5 , 0.5 ], lw = 1 , ls = "--" , c = "k" )
56+ ax1 .plot ([0 , time_end .value_in (time_end .unit )], [0.5 , 0.5 ], lw = 1 , ls = "--" , c = "k" )
5857 for code in codes :
59- time , Q = virial_ratio_evolution (code , particles , Q_init , t_end )
60- plt .plot (time .value_in (t_end .unit ), Q )
58+ time , Q = virial_ratio_evolution (code , particles , Q_init , time_end )
59+ ax1 .plot (
60+ time .value_in (time_end .unit ),
61+ Q ,
62+ )
6163 ci += 1
6264 plt .savefig ("gravity_to_virial" )
6365
6466
6567def new_argument_parser ():
6668 parser = argparse .ArgumentParser (
67- formatter_class = argparse .ArgumentDefaultsHelpFormatter
69+ formatter_class = argparse .ArgumentDefaultsHelpFormatter ,
70+ )
71+ parser .add_argument (
72+ "-N" ,
73+ "--number_of_stars" ,
74+ type = int ,
75+ default = 1000 ,
76+ help = "number of stars" ,
6877 )
69- parser .add_argument ("-N" , type = int , default = 1000 , help = "number of stars" )
7078 parser .add_argument (
7179 "-t" ,
72- "--t_end " ,
80+ "--time_end " ,
7381 type = nbody_system .time ,
7482 default = 2 | nbody_system .time ,
7583 help = "end time of the simulation" ,
0 commit comments