updated Makefile

hermitepn/Makefile → Makefile

@@ -17,29 +17,29 @@ LDFLAGS  += -lm $(MUSE_LD_FLAGS)
 
 OBJS = interface.o
 
-CODELIB = src/libhermitepn.a
+CODELIB = src/libhermite_grx.a
 
 CODE_GENERATOR = $(AMUSE_DIR)/build.py
 
-all: hermitepn_worker
+all: hermite_grx_worker
 
 clean:
 	$(RM) -f *.so *.o *.pyc worker_code.cc worker_code.h 
-	$(RM) *~ hermitepn_worker worker_code.cc
+	$(RM) *~ hermite_grx_worker worker_code.cc
 	$(RM) -f *~
 	make -C src clean
 
 $(CODELIB):
 	make -C src all
 
 worker_code.cc: interface.py
-	$(CODE_GENERATOR) --type=c interface.py HermitePNInterface -o $@
+	$(CODE_GENERATOR) --type=c interface.py HermiteGRXInterface -o $@
 
 worker_code.h: interface.py
-	$(CODE_GENERATOR) --type=H -i amuse.community.interface.stopping_conditions.StoppingConditionInterface interface.py HermitePNInterface -o $@
+	$(CODE_GENERATOR) --type=H -i amuse.community.interface.stopping_conditions.StoppingConditionInterface interface.py HermiteGRXInterface -o $@
 
-hermitepn_worker: worker_code.cc worker_code.h $(CODELIB) $(OBJS)
-	$(MPICXX) $(CXXFLAGS) $< -o $@ $(OBJS) $(CODELIB) -L./src -L$(AMUSE_DIR)/lib/stopcond -lstopcond -lhermitepn
+hermite_grx_worker: worker_code.cc worker_code.h $(CODELIB) $(OBJS)
+	$(MPICXX) $(CXXFLAGS) $< -o $@ $(OBJS) $(CODELIB) -L./src -L$(AMUSE_DIR)/lib/stopcond -lstopcond -lhermite_grx
 
 interface.o: interface.cc
 	$(CXX) $(CXXFLAGS) -c -o $@ $<

README.md

@@ -1,2 +1,7 @@
-# Hermite_GRX
-location: https://github.com/amusecode/Hermite_GRX
+# Hermite Post-Newtonian
+
+This code expands the original hermite code with post-newtonian terms.
+
+## Installation
+From the AMUSE source directory run the following git command `git submodule add [email protected]:amusecode/Hermite_GRX.git src/amuse/community/hermite_grx`. This will add the code as a submodule such that the respective git structures do not interfere. 
+The code can then be compiled using `make hermite_grx.code` or all codes can be built together using `python setup.py develop_build`.

__init__.py

@@ -0,0 +1 @@
+# generated file

GalacticCenter.py → examples/GalacticCenter.py

@@ -3,7 +3,7 @@
 from amuse.units import constants
 
 import orbital_elements as orb_elem
-from hermitepn.interface import *
+from amuse.community.hermite_grx.interface import *
 
 import numpy as np
 import time
@@ -93,7 +93,7 @@ def setup_sim(initial, num_threads, dt_param, integrator, perturbation):
 	mbh, bhs, converter = initial
 
 	if perturbation != None:
-		grav = HermitePN(converter)
+		grav = HermiteGRX(converter)
 		grav.parameters.light_speed = constants.c
 		grav.parameters.perturbation = perturbation
 		grav.parameters.integrator = integrator

HTPHermite.py → examples/HTPHermite.py

@@ -1,7 +1,7 @@
 from amuse.lab import *
 from amuse.units.units import *
 from amuse.units import constants
-from hermitepn.interface import *
+from amuse.community.hermite_grx.interface import *
 
 def HTpulsar(m1, m2, a, e):
     converter = nbody_system.nbody_to_si(m1+m2, a)
@@ -46,28 +46,30 @@ def get_trajectories(initial, grav, t_end, pert=None):
 
     return x1, y1, z1, x2, y2, z2, time, E
 
-m1 = 1.441 | units.MSun
-m2 = 1.397 | units.MSun
-a = 1950100 | units.km
-e = 0.6171334
-initial = HTpulsar(m1, m2, a, e)
-bodies = initial[0]
-Porb = 7.751938773864 | units.hour
-
-converter = initial[1]
-Nbody_code = Hermite
-grav = Nbody_code(converter)
-grav.parameters.dt_param = 0.1
 
-x1, y1, z1, x2, y2, z2, time, E = get_trajectories(initial,
-                                                   grav,
-                                                   10*Porb)
-
-print("n=", len(x1))
-from matplotlib import pyplot
-pyplot.plot(x1.value_in(units.au), y1.value_in(units.au), c='r', lw=1)
-pyplot.plot(x2.value_in(units.au), y2.value_in(units.au), c='g', lw=1)
-pyplot.show()
+if __name__=="__main__":
+    m1 = 1.441 | units.MSun
+    m2 = 1.397 | units.MSun
+    a = 1950100 | units.km
+    e = 0.6171334
+    initial = HTpulsar(m1, m2, a, e)
+    bodies = initial[0]
+    Porb = 7.751938773864 | units.hour
+
+    converter = initial[1]
+    Nbody_code = Hermite
+    grav = Nbody_code(converter)
+    grav.parameters.dt_param = 0.1
+
+    x1, y1, z1, x2, y2, z2, time, E = get_trajectories(initial,
+                                                    grav,
+                                                    10*Porb)
+
+    print("n=", len(x1))
+    from matplotlib import pyplot
+    pyplot.plot(x1.value_in(units.au), y1.value_in(units.au), c='r', lw=1)
+    pyplot.plot(x2.value_in(units.au), y2.value_in(units.au), c='g', lw=1)
+    pyplot.show()
 
-pyplot.scatter(time.value_in(units.yr), E/E[0], c='g')
-pyplot.sh
+    pyplot.scatter(time.value_in(units.yr), E/E[0], c='g')
+    pyplot.sh

HulseTaylor.py → examples/HulseTaylor.py

@@ -3,7 +3,7 @@
 from matplotlib import pyplot
 from amuse.units import constants
 import numpy as np
-from hermitepn.interface import *
+from amuse.community.hermite_grx.interface import *
 from amuse.ext.orbital_elements import new_binary_from_orbital_elements
 from amuse.ext.orbital_elements import orbital_elements_from_binary
 
@@ -74,7 +74,7 @@ def get_trajectories(initial, grav, t_end, dt, pert=None):
 
 def run_nbody_code(Nbody_code, label, dt, tend):
     grav = Nbody_code(converter)
-    if Nbody_code==HermitePN:
+    if Nbody_code==HermiteGRX:
         pert = '1PN_Pairwise'
         grav.parameters.integrator = 'RegularizedHermite'
         if "EIH" in label:
@@ -95,50 +95,52 @@ def run_nbody_code(Nbody_code, label, dt, tend):
                                                                  dt,
                                                                  pert)
     return x1, y1, z1, x2, y2, z2, time, E, sma, ecc
-m1 = 1.441 | units.MSun
-m2 = 1.397 | units.MSun
-a = 1950100 | units.km
-e = 0.6171334
-initial = HTpulsar(m1, m2, a, e)
-bodies = initial[0]
-Porb = 7.751938773864 | units.hour
-converter = initial[1]
-Nbody_codes = [Hermite, HermitePN, HermitePN]
-colors = ['r', 'b', 'g']
-labels = ["Hermite", "hermitePN", "HermitePN\_EIH"]
-dt = 0.1*Porb
-t_end = 10*Porb
 
-figure = pyplot.figure(figsize = (10, 10))
+if __name__=="__main__":
+    m1 = 1.441 | units.MSun
+    m2 = 1.397 | units.MSun
+    a = 1950100 | units.km
+    e = 0.6171334
+    initial = HTpulsar(m1, m2, a, e)
+    bodies = initial[0]
+    Porb = 7.751938773864 | units.hour
+    converter = initial[1]
+    Nbody_codes = [Hermite, HermiteGRX, HermiteGRX]
+    colors = ['r', 'b', 'g']
+    labels = ["Hermite", "hermitePN", "HermitePN\_EIH"]
+    dt = 0.1*Porb
+    t_end = 10*Porb
 
-for Nbody_code, ci, li in zip(Nbody_codes, colors, labels):
-    x1, y1, z1, x2, y2, z2, time, E, sma, ecc = run_nbody_code(Nbody_code, li, dt, t_end)
-    print("n=", len(x1))
-    print("t=", time/Porb)
+    figure = pyplot.figure(figsize = (10, 10))
 
-    subplot = figure.add_subplot(2, 2, 1)
-    pyplot.plot(x1.value_in(units.au), y1.value_in(units.au), c=ci, lw=1)
-    pyplot.plot(x2.value_in(units.au), y2.value_in(units.au), c=ci, lw=1)
-    pyplot.xlabel('x')
-    pyplot.ylabel('y')
+    for Nbody_code, ci, li in zip(Nbody_codes, colors, labels):
+        x1, y1, z1, x2, y2, z2, time, E, sma, ecc = run_nbody_code(Nbody_code, li, dt, t_end)
+        print("n=", len(x1))
+        print("t=", time/Porb)
 
-    subplot = figure.add_subplot(2, 2, 2)
-    pyplot.plot(sma.value_in(units.RSun), ecc, c=ci, label=li)
-    pyplot.xlabel('a')
-    pyplot.ylabel('e')
+        subplot = figure.add_subplot(2, 2, 1)
+        pyplot.plot(x1.value_in(units.au), y1.value_in(units.au), c=ci, lw=1)
+        pyplot.plot(x2.value_in(units.au), y2.value_in(units.au), c=ci, lw=1)
+        pyplot.xlabel('x')
+        pyplot.ylabel('y')
 
-    subplot = figure.add_subplot(2, 2, 3)
-    pyplot.plot(time.value_in(units.yr), sma.value_in(units.RSun), c=ci)
-    pyplot.xlabel('t')
-    pyplot.ylabel('a')
+        subplot = figure.add_subplot(2, 2, 2)
+        pyplot.plot(sma.value_in(units.RSun), ecc, c=ci, label=li)
+        pyplot.xlabel('a')
+        pyplot.ylabel('e')
 
-    subplot = figure.add_subplot(2, 2, 4)
-    pyplot.plot(time.value_in(units.yr), ecc, c=ci)
-    pyplot.xlabel('t')
-    pyplot.ylabel('e')
-
-subplot = figure.add_subplot(2, 2, 2)
-pyplot.legend()
+        subplot = figure.add_subplot(2, 2, 3)
+        pyplot.plot(time.value_in(units.yr), sma.value_in(units.RSun), c=ci)
+        pyplot.xlabel('t')
+        pyplot.ylabel('a')
+
+        subplot = figure.add_subplot(2, 2, 4)
+        pyplot.plot(time.value_in(units.yr), ecc, c=ci)
+        pyplot.xlabel('t')
+        pyplot.ylabel('e')
+
+    subplot = figure.add_subplot(2, 2, 2)
+    pyplot.legend()
 
-#pyplot.show()
-pyplot.savefig("HTpulsar.pdf", fontsize=10)
+    #pyplot.show()
+    pyplot.savefig("HTpulsar.pdf", fontsize=10)

Montgomery.py → examples/Montgomery.py

@@ -2,7 +2,7 @@
 from amuse.units.units import *
 from amuse.units import constants
 import numpy as np
-from hermitepn.interface import *
+from amuse.community.hermite_grx.interface import *
 
 def get_initial_conditions_figure_eight(unit_mass, unit_length):
 	converter = nbody_system.nbody_to_si(unit_mass, unit_length)
@@ -41,7 +41,7 @@ def get_trajectories(initial, grav, t_end, pert=None):
     E = [] | units.J 
 
     t = 0 * t_end
-    dt = t_end / 10000.0
+    dt = t_end / 1000.0
     i = 0.
 
     while t < t_end:
@@ -69,19 +69,25 @@ def get_trajectories(initial, grav, t_end, pert=None):
 #initial[0].velocity += 10 |units.kms
 bodies = initial[0]
 converter = initial[1]
-Nbody_code = HermitePN
+Nbody_code = HermiteGRX
 #Nbody_code = Hermite
-#grav = HermitePN(converter)
 grav = Nbody_code(converter)
 pert = '1PN_Pairwise'
-if Nbody_code=="HermitePN":
+print(Nbody_code)
+if "HermiteGRX" in str(Nbody_code):
     print("perturbation=", grav.parameters.perturbation)
     grav.parameters.perturbation = pert
     grav.parameters.integrator = 'RegularizedHermite'
     grav.parameters.dt_param = 0.1
-    grav.parameters.light_speed = 0.001*constants.c
+    grav.parameters.light_speed = 0.01*constants.c
+
+    #grav.parameters.integrator = 'Hermite'
+    #grav.parameters.perturbation = '1PN_EIH'
+    #grav.parameters.integrator = 'SymmetrizedRegularizedHermite'
+    #grav.parameters.light_speed = relative_lightspeed # * constants.c
 else:
     grav.parameters.dt_param = 0.1
+print(grav.parameters)
 print("v=", bodies[0].velocity.length()/grav.parameters.light_speed)
 
 x1, y1, z1, x2, y2, z2, x3, y3, z3, time, E = get_trajectories(initial,

Pythagoras.py → examples/Pythagoras.py

@@ -3,7 +3,7 @@
 from amuse.units import constants
 import numpy as np
 import matplotlib.pyplot as plt
-from hermitepn.interface import *
+from amuse.community.hermite_grx.interface import *
 
 import orbital_elements as orb_elem
 
@@ -96,7 +96,7 @@ def get_trajectories(initial, grav, t_end, pert=None):
 
 initial = get_initial_conditions_pythagorean(1 | units.MSun, 1 | AU)
 converter = initial[1]
-grav = HermitePN(converter)
+grav = HermiteGRX(converter)
 grav.parameters.perturbation = 'None'
 grav.parameters.integrator = 'RegularizedHermite'
 grav.parameters.dt_param = 0.001

examples/README.md

@@ -0,0 +1,2 @@
+# Hermite_GRX
+location: https://github.com/amusecode/Hermite_GRX