@@ -77,7 +77,7 @@ void EIHPerturbation::CalculateAccelerations(Model *model, int thread_id, int nu
77
77
Real r_bc = sqrt (r_bc2);
78
78
Real r_bc3 = r_bc2 * r_bc;
79
79
Real r_ac = x_ac.Norm ();
80
- Vec n_bc = x_bc / r_bc;
80
+ // Vec n_bc = x_bc / r_bc;
81
81
82
82
Real m_c_over_r_bc3 = c.mass / r_bc3;
83
83
// Add the 3 cross terms
@@ -163,7 +163,7 @@ Real EIHPerturbation::GetEnergy(Model *model)
163
163
Real v_a_dot_v_a = v_a.SquaredNorm ();
164
164
Real square_v_a_dot_v_a = v_a_dot_v_a * v_a_dot_v_a; // PN1.0
165
165
// 1PN kinetic term
166
- energy += c2_recipr * a.mass * 0 .375 * v_a_dot_v_a * v_a_dot_v_a;
166
+ energy += c2_recipr * a.mass * 0.375 * v_a_dot_v_a * v_a_dot_v_a;
167
167
// 2PN kinetic term
168
168
energy += c4_recipr * a.mass * 0.3125 * v_a_dot_v_a * v_a_dot_v_a * v_a_dot_v_a;
169
169
for (Particle &b : all)
@@ -173,32 +173,32 @@ Real EIHPerturbation::GetEnergy(Model *model)
173
173
174
174
// Compute the differences in position and velocity
175
175
Vec x_ab = a.pos - b.pos ;
176
- Vec v_ab = a.vel - b.vel ;
176
+ // Vec v_ab = a.vel - b.vel;
177
177
// compute the distances and its powers
178
178
Real r_ab2 = x_ab.SquaredNorm (); // PN0.0
179
179
Real r_ab = sqrt (r_ab2); // PN1.0
180
- Real r_ab3 = r_ab * r_ab2; // PN1.0
180
+ // Real r_ab3 = r_ab * r_ab2; // PN1.0
181
181
// Compute the normal vector, and mass to distance ratios
182
- Vec n_ab = x_ab / r_ab; // PN0.0
183
- Real m_a_over_r_ab = a.mass / r_ab; // PN1.0
184
- Real ma_am_b_over_r_ab = a.mass * b.mass / r_ab; // PN1.0
185
- Real ma_m_b_over_r_ab3 = a.mass * b.mass / r_ab3; // PN1.0
186
- Real ma_m_b_over_r_ab2 = a.mass * b.mass / r_ab2; // PN1.0
187
- Real m_b_over_r_ab = b.mass / r_ab; // PN1.0
182
+ Vec n_ab = x_ab / r_ab; // PN0.0
183
+ Real m_a_over_r_ab = a.mass / r_ab; // PN1.0
184
+ // Real ma_am_b_over_r_ab = a.mass * b.mass / r_ab; // PN1.0
185
+ // Real ma_m_b_over_r_ab3 = a.mass * b.mass / r_ab3; // PN1.0
186
+ // Real m_a_m_b_over_r_ab2 = a.mass * b.mass / r_ab2; // PN1.0
187
+ Real m_b_over_r_ab = b.mass / r_ab; // PN1.0
188
188
189
189
// Compute the inner products used in the equations
190
- Real v_b_dot_n_ab = b.vel .Dot (n_ab); // PN1.0
191
- Real v_a_dot_v_b = a.vel .Dot (b.vel ); // PN1.0
192
- Real v_b_dot_v_b = b.vel .SquaredNorm (); // PN1.0
193
- Real square_v_b_dot_v_b = b.vel .SquaredNorm () * b.vel .SquaredNorm (); // PN1.0
194
- Real v_a_dot_n_ab = a.vel .Dot (n_ab); // PN2.0
195
- Real square_v_b_dot_n_ab = v_b_dot_n_ab * v_b_dot_n_ab; // PN2.0
196
- Real square_v_a_dot_n_ab = v_a_dot_n_ab * v_a_dot_n_ab; // PN2.0
197
- Real cubed_v_a_dot_n_ab = square_v_a_dot_n_ab * v_a_dot_n_ab; // PN2.0
198
- Real quarted_v_b_dot_n_ab = square_v_b_dot_n_ab * square_v_b_dot_n_ab; // PN2.0
190
+ Real v_b_dot_n_ab = b.vel .Dot (n_ab); // PN1.0
191
+ Real v_a_dot_v_b = a.vel .Dot (b.vel ); // PN1.0
192
+ Real v_b_dot_v_b = b.vel .SquaredNorm (); // PN1.0
193
+ // Real square_v_b_dot_v_b = b.vel.SquaredNorm() * b.vel.SquaredNorm(); // PN1.0
194
+ Real v_a_dot_n_ab = a.vel .Dot (n_ab); // PN2.0
195
+ Real square_v_b_dot_n_ab = v_b_dot_n_ab * v_b_dot_n_ab; // PN2.0
196
+ Real square_v_a_dot_n_ab = v_a_dot_n_ab * v_a_dot_n_ab; // PN2.0
197
+ Real cubed_v_a_dot_n_ab = square_v_a_dot_n_ab * v_a_dot_n_ab; // PN2.0
198
+ // Real quarted_v_b_dot_n_ab = square_v_b_dot_n_ab * square_v_b_dot_n_ab; // PN2.0
199
199
200
200
// 1PN terms
201
- Real energypart = 1.5 * va_dot_va * m_b_over_r_ab;
201
+ Real energypart = 1.5 * v_a_dot_v_a * m_b_over_r_ab;
202
202
Real energytemp = 7.0 * a.vel .Dot (b.vel );
203
203
energytemp += a.vel .Dot (n_ab) * b.vel .Dot (n_ab);
204
204
energypart -= 0.25 * m_b_over_r_ab * energytemp;
@@ -215,23 +215,23 @@ Real EIHPerturbation::GetEnergy(Model *model)
215
215
// ADD 1PN terms
216
216
energy += c2_recipr * energypart * a.mass ;
217
217
// 2PN terms
218
- Real energypart = -0.5 * m_a_over_r_ab * m_a_over_r_ab;
219
- energypart -= 2.375 * m_a_over_r_ab * m_b_over_r_ab;
220
- energypart += 0.375 * cubed_v_a_dot_n_ab * v_b_dot_n_ab;
221
- energypart += 0.1875 * square_v_a_dot_n_ab * square_v_b_dot_n_ab;
222
- energypart -= 1.125 * v_a_dot_n_ab * v_b_dot_n_ab * v_a_dot_v_a;
223
- energypart -= 1.625 * square_v_b_dot_n_ab * v_a_dot_v_a;
224
- energypart += 2.625 * square_v_a_dot_v_a;
225
- energypart += 1.625 * square_v_a_dot_n_ab * v_a_dot_v_b;
226
- energypart += 0.75 * v_a_dot_n_ab * v_b_dot_n_ab * v_a_dot_v_b;
227
- energypart -= 6.875 * v_a_dot_v_a * v_a_dot_v_b;
228
- energypart += 2.125 * v_a_dot_v_b * v_a_dot_v_b;
229
- energypart += 1.9375 * v_a_dot_v_a * v_b_dot_v_b;
230
- Real energytemp = 7.25 * square_v_a_dot_n_ab;
231
- energytemp -= 3.25 * v_a_dot_n_ab * v_b_dot_n_ab;
232
- energytemp += 0.5 * square_v_b_dot_n_ab;
233
- energytemp + = 1.5 * v_a_dot_v_a;
234
- energytemp += 1.75 * v_b_dot_v_b;
218
+ energypart = -0.5 * m_a_over_r_ab * m_a_over_r_ab; // 1/2
219
+ energypart -= 2.375 * m_a_over_r_ab * m_b_over_r_ab; // 9/8
220
+ energypart += 0.375 * cubed_v_a_dot_n_ab * v_b_dot_n_ab; // 3/8
221
+ energypart += 0.1875 * square_v_a_dot_n_ab * square_v_b_dot_n_ab; // 3/16
222
+ energypart -= 1.125 * v_a_dot_n_ab * v_b_dot_n_ab * v_a_dot_v_a; // 9/8
223
+ energypart -= 1.625 * square_v_b_dot_n_ab * v_a_dot_v_a; // 13/8
224
+ energypart += 2.625 * square_v_a_dot_v_a; // 21/8
225
+ energypart += 1.625 * square_v_a_dot_n_ab * v_a_dot_v_b; // 13/8
226
+ energypart += 0.75 * v_a_dot_n_ab * v_b_dot_n_ab * v_a_dot_v_b; // 3/4
227
+ energypart -= 6.875 * v_a_dot_v_a * v_a_dot_v_b; // 55/8
228
+ energypart += 2.125 * v_a_dot_v_b * v_a_dot_v_b; // 17/8
229
+ energypart += 1.9375 * v_a_dot_v_a * v_b_dot_v_b; // 31/16
230
+ energytemp = 7.25 * square_v_a_dot_n_ab; // 29/4
231
+ energytemp -= 3.25 * v_a_dot_n_ab * v_b_dot_n_ab; // 13/4
232
+ energytemp += 0.5 * square_v_b_dot_n_ab; // 1/2
233
+ energytemp - = 1.5 * v_a_dot_v_a; // 3/2
234
+ energytemp += 1.75 * v_b_dot_v_b; // 7/4
235
235
energypart += energytemp * m_a_over_r_ab;
236
236
energy += a.mass * c4_recipr * energypart * m_b_over_r_ab;
237
237
}
@@ -251,7 +251,7 @@ Vec EIHPerturbation::GetLinearMomentum(Model *model)
251
251
// 0PN
252
252
momentum += a.mass * a.vel ;
253
253
// 1PN
254
- Real v_a_dot_v_a = v_a .SquaredNorm ();
254
+ Real v_a_dot_v_a = a. vel .SquaredNorm ();
255
255
Vec momentumpart = v_a_dot_v_a * a.vel ;
256
256
for (Particle &b : all)
257
257
{
@@ -261,9 +261,9 @@ Vec EIHPerturbation::GetLinearMomentum(Model *model)
261
261
Real r_ab = x_ab.Norm ();
262
262
Vec n_ab = x_ab / r_ab;
263
263
momentumpart -= b.mass / r_ab * a.vel ;
264
- momentumpart -= b.mass / r_ab * (a.vel .Dot (n_ab)) * n_ab
264
+ momentumpart -= b.mass / r_ab * (a.vel .Dot (n_ab)) * n_ab;
265
265
}
266
- momentum += 0.5 * c2_recipr * a.mass * momentumpart
266
+ momentum += 0.5 * c2_recipr * a.mass * momentumpart;
267
267
}
268
268
return momentum;
269
269
}
0 commit comments