33from matplotlib import pyplot
44from amuse .units import constants
55import numpy as np
6- from hermitepn .interface import *
6+ from amuse . community . hermite_grx .interface import *
77from amuse .ext .orbital_elements import new_binary_from_orbital_elements
88from amuse .ext .orbital_elements import orbital_elements_from_binary
99
@@ -74,7 +74,7 @@ def get_trajectories(initial, grav, t_end, dt, pert=None):
7474
7575def run_nbody_code (Nbody_code , label , dt , tend ):
7676 grav = Nbody_code (converter )
77- if Nbody_code == HermitePN :
77+ if Nbody_code == HermiteGRX :
7878 pert = '1PN_Pairwise'
7979 grav .parameters .integrator = 'RegularizedHermite'
8080 if "EIH" in label :
@@ -95,50 +95,52 @@ def run_nbody_code(Nbody_code, label, dt, tend):
9595 dt ,
9696 pert )
9797 return x1 , y1 , z1 , x2 , y2 , z2 , time , E , sma , ecc
98- m1 = 1.441 | units .MSun
99- m2 = 1.397 | units .MSun
100- a = 1950100 | units .km
101- e = 0.6171334
102- initial = HTpulsar (m1 , m2 , a , e )
103- bodies = initial [0 ]
104- Porb = 7.751938773864 | units .hour
105- converter = initial [1 ]
106- Nbody_codes = [Hermite , HermitePN , HermitePN ]
107- colors = ['r' , 'b' , 'g' ]
108- labels = ["Hermite" , "hermitePN" , "HermitePN\_EIH" ]
109- dt = 0.1 * Porb
110- t_end = 10 * Porb
11198
112- figure = pyplot .figure (figsize = (10 , 10 ))
99+ if __name__ == "__main__" :
100+ m1 = 1.441 | units .MSun
101+ m2 = 1.397 | units .MSun
102+ a = 1950100 | units .km
103+ e = 0.6171334
104+ initial = HTpulsar (m1 , m2 , a , e )
105+ bodies = initial [0 ]
106+ Porb = 7.751938773864 | units .hour
107+ converter = initial [1 ]
108+ Nbody_codes = [Hermite , HermiteGRX , HermiteGRX ]
109+ colors = ['r' , 'b' , 'g' ]
110+ labels = ["Hermite" , "hermitePN" , "HermitePN\_EIH" ]
111+ dt = 0.1 * Porb
112+ t_end = 10 * Porb
113113
114- for Nbody_code , ci , li in zip (Nbody_codes , colors , labels ):
115- x1 , y1 , z1 , x2 , y2 , z2 , time , E , sma , ecc = run_nbody_code (Nbody_code , li , dt , t_end )
116- print ("n=" , len (x1 ))
117- print ("t=" , time / Porb )
114+ figure = pyplot .figure (figsize = (10 , 10 ))
118115
119- subplot = figure .add_subplot (2 , 2 , 1 )
120- pyplot .plot (x1 .value_in (units .au ), y1 .value_in (units .au ), c = ci , lw = 1 )
121- pyplot .plot (x2 .value_in (units .au ), y2 .value_in (units .au ), c = ci , lw = 1 )
122- pyplot .xlabel ('x' )
123- pyplot .ylabel ('y' )
116+ for Nbody_code , ci , li in zip (Nbody_codes , colors , labels ):
117+ x1 , y1 , z1 , x2 , y2 , z2 , time , E , sma , ecc = run_nbody_code (Nbody_code , li , dt , t_end )
118+ print ("n=" , len (x1 ))
119+ print ("t=" , time / Porb )
124120
125- subplot = figure .add_subplot (2 , 2 , 2 )
126- pyplot .plot (sma .value_in (units .RSun ), ecc , c = ci , label = li )
127- pyplot .xlabel ('a' )
128- pyplot .ylabel ('e' )
121+ subplot = figure .add_subplot (2 , 2 , 1 )
122+ pyplot .plot (x1 .value_in (units .au ), y1 .value_in (units .au ), c = ci , lw = 1 )
123+ pyplot .plot (x2 .value_in (units .au ), y2 .value_in (units .au ), c = ci , lw = 1 )
124+ pyplot .xlabel ('x' )
125+ pyplot .ylabel ('y' )
129126
130- subplot = figure .add_subplot (2 , 2 , 3 )
131- pyplot .plot (time . value_in ( units . yr ), sma .value_in (units .RSun ), c = ci )
132- pyplot .xlabel ('t ' )
133- pyplot .ylabel ('a ' )
127+ subplot = figure .add_subplot (2 , 2 , 2 )
128+ pyplot .plot (sma .value_in (units .RSun ), ecc , c = ci , label = li )
129+ pyplot .xlabel ('a ' )
130+ pyplot .ylabel ('e ' )
134131
135- subplot = figure .add_subplot (2 , 2 , 4 )
136- pyplot .plot (time .value_in (units .yr ), ecc , c = ci )
137- pyplot .xlabel ('t' )
138- pyplot .ylabel ('e' )
139-
140- subplot = figure .add_subplot (2 , 2 , 2 )
141- pyplot .legend ()
132+ subplot = figure .add_subplot (2 , 2 , 3 )
133+ pyplot .plot (time .value_in (units .yr ), sma .value_in (units .RSun ), c = ci )
134+ pyplot .xlabel ('t' )
135+ pyplot .ylabel ('a' )
136+
137+ subplot = figure .add_subplot (2 , 2 , 4 )
138+ pyplot .plot (time .value_in (units .yr ), ecc , c = ci )
139+ pyplot .xlabel ('t' )
140+ pyplot .ylabel ('e' )
141+
142+ subplot = figure .add_subplot (2 , 2 , 2 )
143+ pyplot .legend ()
142144
143- #pyplot.show()
144- pyplot .savefig ("HTpulsar.pdf" , fontsize = 10 )
145+ #pyplot.show()
146+ pyplot .savefig ("HTpulsar.pdf" , fontsize = 10 )
0 commit comments