Skip to content

Commit 4065669

Browse files
GijsVermarienGijsVermarien
committed
Small updates
1 parent a2f8460 commit 4065669

File tree

6 files changed

+167
-155
lines changed

6 files changed

+167
-155
lines changed

‎examples/HTPHermite.py‎

Lines changed: 25 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -46,28 +46,30 @@ def get_trajectories(initial, grav, t_end, pert=None):
4646

4747
return x1, y1, z1, x2, y2, z2, time, E
4848

49-
m1 = 1.441 | units.MSun
50-
m2 = 1.397 | units.MSun
51-
a = 1950100 | units.km
52-
e = 0.6171334
53-
initial = HTpulsar(m1, m2, a, e)
54-
bodies = initial[0]
55-
Porb = 7.751938773864 | units.hour
56-
57-
converter = initial[1]
58-
Nbody_code = Hermite
59-
grav = Nbody_code(converter)
60-
grav.parameters.dt_param = 0.1
6149

62-
x1, y1, z1, x2, y2, z2, time, E = get_trajectories(initial,
63-
grav,
64-
10*Porb)
65-
66-
print("n=", len(x1))
67-
from matplotlib import pyplot
68-
pyplot.plot(x1.value_in(units.au), y1.value_in(units.au), c='r', lw=1)
69-
pyplot.plot(x2.value_in(units.au), y2.value_in(units.au), c='g', lw=1)
70-
pyplot.show()
50+
if __name__=="__main__":
51+
m1 = 1.441 | units.MSun
52+
m2 = 1.397 | units.MSun
53+
a = 1950100 | units.km
54+
e = 0.6171334
55+
initial = HTpulsar(m1, m2, a, e)
56+
bodies = initial[0]
57+
Porb = 7.751938773864 | units.hour
58+
59+
converter = initial[1]
60+
Nbody_code = Hermite
61+
grav = Nbody_code(converter)
62+
grav.parameters.dt_param = 0.1
63+
64+
x1, y1, z1, x2, y2, z2, time, E = get_trajectories(initial,
65+
grav,
66+
10*Porb)
67+
68+
print("n=", len(x1))
69+
from matplotlib import pyplot
70+
pyplot.plot(x1.value_in(units.au), y1.value_in(units.au), c='r', lw=1)
71+
pyplot.plot(x2.value_in(units.au), y2.value_in(units.au), c='g', lw=1)
72+
pyplot.show()
7173

72-
pyplot.scatter(time.value_in(units.yr), E/E[0], c='g')
73-
pyplot.sh
74+
pyplot.scatter(time.value_in(units.yr), E/E[0], c='g')
75+
pyplot.sh

‎examples/HulseTaylor.py‎

Lines changed: 42 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -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)

‎interface.cc‎

Lines changed: 21 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,7 @@ string current_integrator = "Hermite";
2525
string current_perturbation = "None";
2626

2727
Real time_step_parameter = 0.03;
28-
Real speed_of_light = 1;
28+
Real speed_of_light = 1.0;
2929
size_t num_threads = 1;
3030
Real epsilon2 = 0;
3131

@@ -35,18 +35,18 @@ int initialize_code()
3535
{
3636
if (model)
3737
return 0;
38-
38+
3939
model = new Model();
40-
40+
4141
integrators["Hermite"] = new HermiteIntegrator(false);
4242
integrators["SymmetrizedHermite"] = new HermiteIntegrator(true);
4343
integrators["RegularizedHermite"] = new RegularizedHermiteIntegrator(false);
4444
integrators["SymmetrizedRegularizedHermite"] = new RegularizedHermiteIntegrator(true);
45-
45+
4646
perturbations["None"] = new Perturbation();
4747
perturbations["1PN_Pairwise"] = new Pairwise1PNPerturbation();
4848
perturbations["1PN_EIH"] = new EIHPerturbation();
49-
49+
5050
return 0;
5151
}
5252

@@ -56,7 +56,7 @@ int commit_parameters()
5656
{
5757
integrator = integrators.at(current_integrator);
5858
perturbation = perturbations.at(current_perturbation);
59-
59+
6060
integrator->SetPerturbation(perturbation);
6161
integrator->SetNumThreads(num_threads);
6262
integrator->SetTimeStepParameter(time_step_parameter);
@@ -68,7 +68,7 @@ int commit_parameters()
6868
perturbation = nullptr;
6969
return -1;
7070
}
71-
71+
7272
return 0;
7373
}
7474

@@ -138,7 +138,7 @@ int set_light_speed(double c)
138138
int get_integrator(char **i)
139139
{
140140
// Casting away constness
141-
*i = (char *) current_integrator.c_str();
141+
*i = (char *)current_integrator.c_str();
142142
return 0;
143143
}
144144

@@ -151,7 +151,7 @@ int set_integrator(char *i)
151151
int get_perturbation(char **i)
152152
{
153153
// Casting away constness
154-
*i = (char *) current_perturbation.c_str();
154+
*i = (char *)current_perturbation.c_str();
155155
return 0;
156156
}
157157

@@ -163,15 +163,15 @@ int set_perturbation(char *i)
163163

164164
int get_num_threads(int *thread_number)
165165
{
166-
*thread_number = (int) num_threads;
166+
*thread_number = (int)num_threads;
167167
return 0;
168168
}
169169

170170
int set_num_threads(int thread_number)
171171
{
172172
if (thread_number <= 0)
173173
return -1;
174-
174+
175175
num_threads = thread_number;
176176
return 0;
177177
}
@@ -181,17 +181,17 @@ int cleanup_code()
181181
if (model)
182182
delete model;
183183
model = nullptr;
184-
184+
185185
for (auto &p : perturbations)
186186
delete p.second;
187187
perturbations.clear();
188188
perturbation = nullptr;
189-
189+
190190
for (auto &i : integrators)
191191
delete i.second;
192192
integrators.clear();
193193
integrator = nullptr;
194-
194+
195195
return 0;
196196
}
197197

@@ -402,7 +402,7 @@ int evolve_model(double t)
402402
{
403403
if (num_threads > model->GetNumParticles())
404404
return -1;
405-
405+
406406
integrator->Evolve(model, t);
407407
return 0;
408408
}
@@ -425,15 +425,15 @@ int get_total_energy_with(char *type, double *etot)
425425
{
426426
Perturbation *pert = perturbations.at(string(type));
427427
pert->SetSpeedOfLight(speed_of_light);
428-
428+
429429
*etot = model->GetNewtonianEnergy();
430430
*etot += pert->GetEnergy(model);
431431
}
432432
catch (...)
433433
{
434434
return -1;
435435
}
436-
436+
437437
return 0;
438438
}
439439

@@ -446,7 +446,7 @@ int get_total_linear_momentum_with(char *type, double *px, double *py, double *p
446446

447447
Vec p = model->GetNewtonianLinearMomentum();
448448
p += pert->GetLinearMomentum(model);
449-
449+
450450
*px = p[0];
451451
*py = p[1];
452452
*pz = p[2];
@@ -455,7 +455,7 @@ int get_total_linear_momentum_with(char *type, double *px, double *py, double *p
455455
{
456456
return -1;
457457
}
458-
458+
459459
return 0;
460460
}
461461

@@ -483,10 +483,10 @@ int get_total_mass(double *mass)
483483
{
484484
ParticleSetView &all = model->GetAllParticles();
485485
Real m = 0;
486-
486+
487487
for (Particle &p : all)
488488
m += p.mass;
489-
489+
490490
*mass = m;
491491
return 0;
492492
}

0 commit comments

Comments
 (0)