PeTar
N-body code for collisional gravitational systems
pn_BH.h
Go to the documentation of this file.
1 /***************************************************************************/
2 /*
3  Coded by : Peter Berczik (on the base of Gabor Kupi original PN code)
4  Version number : 1.0
5  Last redaction : 2011.V.6. 11:00
6 */
7 
8 inline double SQR(const double x) {return x*x;}
9 extern bool is_root();
10 
11 int calc_force_pn_BH(double m1, double xx1[], double vv1[], double spin1[],
12  double m2, double xx2[], double vv2[], double spin2[],
13  double CCC_NB, int usedOrNot[], double dt_bh,
14  double a_pn1[][3], double adot_pn1[][3],
15  double a_pn2[][3], double adot_pn2[][3])
16 {
17  /*
18  INPUT
19 
20  m1 - mass of the 1 BH
21  xx1[0,1,2] - coordinate of the 1 BH
22  vv1[0,1,2] - velocity of the 1 BH
23  spin1[0,1,2] - normalized spin of the 1 BH
24 
25  m2 - mass of the 2 BH
26  xx2[0,1,2] - coordinate of the 2 BH
27  vv2[0,1,2] - velocity of the 2 BH
28  spin2[0,1,2] - normalized spin of the 2 BH
29 
30  CCC_NB - Speed of light "c" in internal units
31  usedOrNot[PN1, PN2, PN2.5, PN3, PN3.5, SPIN] - different PN term usage: PN1, PN2, PN2.5, PN3, PN3.5, SPIN
32  0 1 2 3 4 5
33 
34  OUTPUT
35 
36  a_pn1 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5] [3] for the 1 BH
37  adot_pn1[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5] [3] for the 1 BH
38 
39  a_pn2 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5] [3] for the 2 BH
40  adot_pn2[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5] [3] for the 2 BH
41 
42  return - 0 if everything OK
43  - 505 if BH's separation < 4 x (RSwarch1 + RSwarch2)
44  */
45 
46 
47  int j, k;
48  double PI2 = 9.86960440108935;
49 
50  double c_1, c_2, c_4, c_5, c_6, c_7, RS_DIST;
51 
52  double M, eta, r, r2, r3, MOR;
53  double V1_V22,VWHOLE, RP, RPP, VA;
54  double N[3], x[3], v[3], A[3];
55  double A1, B1, A2, B2, A2_5, B2_5, AK2, BK2, AK4, BK4, AK5, BK5;
56  double A1D, A2D, A2_5D, B1D, B2D, B2_5D, ADK2, BDK2, ADK4, BDK4, ADK5, BDK5;
57 
58  double A3, B3, A3_5, B3_5, AK6, BK6, AK7, BK7;
59  double A3D, A3_5D, B3D, B3_5D, ADK6, BDK6, ADK7, BDK7;
60 
61  int Van_Spin=0;
62  double DM, S1[3], SPIN[3][2], S2[3], KSS[3], KSSIG[3], XS[3], XA[3], NCV[3], NCS[3], NCSIG[3],
63  VCS[3], VCSIG[3], SDNCV, SIGDNCV, NDV, XS2, XA2, NXA, NXS, VDS, VDSIG, NDS, NDSIG,
64  C1_5[3], C2[3], C2_5[3];
65  double LABS, LU[3], S1DLU, S2DLU, SU1[3], SV1[3], SS1[3], SS2[3], SU2[3], SV2[3];
66  double AT[3], NDOT[3], NVDOT, NDOTCV[3], NCA[3];
67  double SS1aux[3],SS2aux[3],SU[3],SV[3],XAD[3],XSD[3];
68  double NDOTCS[3], NCSU[3], NDOTCSIG[3], NCSV[3], ACS[3], VCSU[3], ACSIG[3], VCSV[3], SNVDOT,
69  SIGNVDOT, NSDOT, NSIGDOT, VSDOT, VSIGDOT, NXSDOT, NXADOT;
70  double C1_5D[3],C2D[3], C2_5D[3];
71  double ADK, BDK, AD[3], KSAK, KSBK;
72 
73  Van_Spin = usedOrNot[5]; // Van vagy nincs SPIN szamolas...
74 
75  for(k=0;k<3;k++)
76  {
77  SPIN[k][0] = spin1[k];
78  SPIN[k][1] = spin2[k];
79  }
80 
81 
82  for(j=0;j<4;j++)
83  {
84  for(k=0;k<3;k++)
85  {
86  a_pn1[j][k] = 0.0; adot_pn1[j][k] = 0.0;
87  a_pn2[j][k] = 0.0; adot_pn2[j][k] = 0.0;
88  }
89  }
90 
91 
92  // Speed of light "c" and its powers
93 
94  c_1 = CCC_NB;
95 
96  c_2 = SQR(c_1);
97  c_4 = SQR(c_2);
98  c_5 = c_4*c_1;
99  c_6 = c_5*c_1;
100  c_7 = c_6*c_1;
101 
102 
103  // Mass parameters
104 
105  M = m1+m2;
106  eta = m1*m2/(M*M);
107 
108  for(k=0;k<3;k++)
109  {
110  x[k] = xx1[k] - xx2[k];
111  v[k] = vv1[k] - vv2[k];
112  }
113 
114 
115  r2 = SQR(x[0]) + SQR(x[1]) + SQR(x[2]);
116  r = sqrt(r2);
117  r3 = r2*r;
118 
119  MOR = M/r;
120  V1_V22 = v[0]*v[0]+v[1]*v[1]+v[2]*v[2];
121  VWHOLE = sqrt(V1_V22);
122  RP = (x[0]*v[0]+x[1]*v[1]+x[2]*v[2])/r;
123 
124 
125 
126  // Newton accelerations
127 
128  for(k=0;k<3;k++) N[k] = x[k]/r;
129 
130  // PN accelerations
131 
132  AK2 = 0.0;
133  BK2 = 0.0;
134  AK4 = 0.0;
135  BK4 = 0.0;
136  AK5 = 0.0;
137  BK5 = 0.0;
138  AK6 = 0.0;
139  BK6 = 0.0;
140  AK7 = 0.0;
141  BK7 = 0.0;
142  C1_5[0] = 0.0;
143  C1_5[1] = 0.0;
144  C1_5[2] = 0.0;
145  C2[0] = 0.0;
146  C2[1] = 0.0;
147  C2[2] = 0.0;
148  C2_5[0] = 0.0;
149  C2_5[1] = 0.0;
150  C2_5[2] = 0.0;
151 
152  if(usedOrNot[0] == 1) // PN1 ~1/c^2
153  {
154  A1 = 2.0*(2.0+eta)*MOR-(1.0+3.0*eta)*V1_V22 +1.5*eta*RP*RP;
155  B1 = 2.0*(2.0-eta)*RP;
156 
157  AK2 = A1/c_2;
158  BK2 = B1/c_2;
159  }
160 
161  if(usedOrNot[1] == 1) // PN2 ~1/c^4
162  {
163  A2 = -0.75*(12.0+29.0*eta)*MOR*MOR-eta*(3.0-4.0*eta)*V1_V22*V1_V22-1.875*eta*(1.0-3.0*eta)*RP*RP*RP*RP+0.5*eta*(13.0-4.0*eta)*MOR*V1_V22+(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RP+1.5*eta*(3.0-4.0*eta)*V1_V22*RP*RP;
164  B2 = -0.5*RP*((4.0+41.0*eta+8.0*eta*eta)*MOR-eta*(15.0+4.0*eta)*V1_V22+3.0*eta*(3.0+2.0*eta)*RP*RP);
165 
166  AK4 = A2/c_4;
167  BK4 = B2/c_4;
168  }
169 
170  if(usedOrNot[2] == 1) // PN2.5 ~1/c^5
171  {
172  A2_5 = 1.6*eta*MOR*RP*(17.0*MOR/3.0+3.0*V1_V22);
173  B2_5 = -1.6*eta*MOR*(3.0*MOR+V1_V22);
174 
175  AK5 = A2_5/c_5;
176  BK5 = B2_5/c_5;
177  }
178 
179  if(usedOrNot[3] == 1) // PN3 ~1/c^6
180  {
181  A3 = MOR*MOR*MOR*(16.0+(1399.0/12.0-41.0*PI2/16.0)*eta+
182  71.0*eta*eta/2.0)+eta*(20827.0/840.0+123.0*PI2/64.0-eta*eta)
183  *MOR*MOR*V1_V22-(1.0+(22717.0/168.0+615.0*PI2/64.0)*eta+
184  11.0*eta*eta/8.0-7.0*eta*eta*eta)*MOR*MOR*RP*RP-
185  0.25*eta*(11.0-49.0*eta+52.0*eta*eta)*V1_V22*V1_V22*V1_V22+
186  35.0*eta*(1.0-5.0*eta+5.0*eta*eta)*RP*RP*RP*RP*RP*RP/16.0-
187  0.25*eta*(75.0+32.0*eta-40.0*eta*eta)*MOR*V1_V22*V1_V22-
188  0.5*eta*(158.0-69.0*eta-60.0*eta*eta)*MOR*RP*RP*RP*RP+
189  eta*(121.0-16.0*eta-20.0*eta*eta)*MOR*V1_V22*RP*RP+
190  3.0*eta*(20.0-79.0*eta+60.0*eta*eta)*V1_V22*V1_V22*RP*RP/8.0-
191  15.0*eta*(4.0-18.0*eta+17.0*eta*eta)*V1_V22*RP*RP*RP*RP/8.0;
192 
193  B3 = RP*((4.0+(5849.0/840.0+123.0*PI2/32.0)*eta-25.0*eta*eta-
194  8.0*eta*eta*eta)*MOR*MOR+eta*(65.0-152.0*eta-48.0*eta*eta)*
195  V1_V22*V1_V22/8.0+15.0*eta*(3.0-8.0*eta-2.0*eta*eta)*RP*RP*RP*RP/8.0+
196  eta*(15.0+27.0*eta+10.0*eta*eta)*MOR*V1_V22-eta*(329.0+177.0*eta+
197  108.0*eta*eta)*MOR*RP*RP/6.0-
198  3.0*eta*(16.0-37.0*eta-16.0*eta*eta)*V1_V22*RP*RP/4.0);
199 
200  AK6 = A3/c_6;
201  BK6 = B3/c_6;
202  }
203 
204  if(usedOrNot[4] == 1) // PN3.5 ~1/c^7
205  {
206  A3_5 = MOR*eta*(V1_V22*V1_V22*(-366.0/35.0-12.0*eta)+V1_V22*RP*RP*(114.0+12.0*eta)-112.0*RP*RP*RP*RP+MOR*(V1_V22*(-692.0/35.0+724.0*eta/15.0)+RP*RP*(-294.0/5.0-376.0*eta/5.0)+MOR*(-3956.0/35.0-184.0*eta/5.0)));
207  B3_5 = 8.0*eta*MOR*((1325.0+546.0*eta)*MOR*MOR/42.0+(313.0+42.0*eta)*V1_V22*V1_V22/28.0+75.0*RP*RP*RP*RP-(205.0+777.0*eta)*MOR*V1_V22/42.0+(205.0+424.0*eta)*MOR*RP*RP/12.0-3.0*(113.0+2.0*eta)*V1_V22*RP*RP/4.0)/5.0;
208 
209  AK7 = A3_5/c_7;
210  BK7 = B3_5/c_7;
211  }
212 
213 
214  // Spin accelerations
215 
216  if(Van_Spin==1)
217  {
218  DM = m1 - m2;
219 
220 
221  for(k=0;k<3;k++)
222  {
223  S1[k] = SPIN[k][0]*m1*m1/c_1; // fizikai spin
224  S2[k] = SPIN[k][1]*m2*m2/c_1;
225  KSS[k] = S1[k]+S2[k];
226  KSSIG[k] = M*(S2[k]/m2-S1[k]/m1);
227  XS[k] = 0.5*(SPIN[k][0]+SPIN[k][1]);
228  XA[k] = 0.5*(SPIN[k][0]-SPIN[k][1]);
229  }
230 
231 
232  //NCV crossproduct of N[k] and relative v = N[k]Xv[j]
233  NCV[0] = N[1]*v[2] - N[2]*v[1];
234  NCV[1] = N[2]*v[0] - N[0]*v[2];
235  NCV[2] = N[0]*v[1] - N[1]*v[0];
236 
237  //NCS crossproduct of N[k] and KSS = N[k]XKSS
238  NCS[0] = N[1]*KSS[2] - N[2]*KSS[1];
239  NCS[1] = N[2]*KSS[0] - N[0]*KSS[2];
240  NCS[2] = N[0]*KSS[1] - N[1]*KSS[0];
241 
242  //NCSIG crossproduct of N[k] and KSSIG = N[k]XKSSIG
243  NCSIG[0] = N[1]*KSSIG[2] - N[2]*KSSIG[1];
244  NCSIG[1] = N[2]*KSSIG[0] - N[0]*KSSIG[2];
245  NCSIG[2] = N[0]*KSSIG[1] - N[1]*KSSIG[0];
246 
247  //VCS crossproduct of v[k] and KSS = v[k]XKSS
248  VCS[0] = v[1]*KSS[2] - v[2]*KSS[1];
249  VCS[1] = v[2]*KSS[0] - v[0]*KSS[2];
250  VCS[2] = v[0]*KSS[1] - v[1]*KSS[0];
251 
252  //VCSIG crossproduct of v[k] and KSSIG = v[k]XKSSIG
253  VCSIG[0] = v[1]*KSSIG[2] - v[2]*KSSIG[1];
254  VCSIG[1] = v[2]*KSSIG[0] - v[0]*KSSIG[2];
255  VCSIG[2] = v[0]*KSSIG[1] - v[1]*KSSIG[0];
256 
257  SDNCV = KSS[0]*NCV[0]+KSS[1]*NCV[1]+KSS[2]*NCV[2];
258 
259  SIGDNCV = KSSIG[0]*NCV[0]+KSSIG[1]*NCV[1]+KSSIG[2]*NCV[2];
260 
261  NDV = N[0]*v[0] + N[1]*v[1] + N[2]*v[2];
262 
263  XS2 = XS[0]*XS[0]+XS[1]*XS[1]+XS[2]*XS[2];
264  XA2 = XA[0]*XA[0]+XA[1]*XA[1]+XA[2]*XA[2];
265 
266  NXA = N[0]*XA[0]+N[1]*XA[1]+N[2]*XA[2];
267  NXS = N[0]*XS[0]+N[1]*XS[1]+N[2]*XS[2];
268 
269  VDS = v[0]*KSS[0]+v[1]*KSS[1]+v[2]*KSS[2];
270  VDSIG = v[0]*KSSIG[0]+v[1]*KSSIG[1]+v[2]*KSSIG[2];
271 
272  NDS = N[0]*KSS[0]+N[1]*KSS[1]+N[2]*KSS[2];
273  NDSIG = N[0]*KSSIG[0]+N[1]*KSSIG[1]+N[2]*KSSIG[2];
274 
275  for(k=0;k<3;k++)
276  {
277  C1_5[k] = (N[k]*(12.0*SDNCV+6.0*DM*SIGDNCV/M)+9.0*NDV*NCS[k]+3.0*DM*NDV*NCSIG[k]/M -7.0*VCS[k]-3.0*DM*VCSIG[k]/M)/r3;
278  C2[k] = -MOR*MOR*MOR/r*3.0*eta*(N[k]*(XS2-XA2-5.0*NXS*NXS+5.0*NXA*NXA)+2.0*(XS[k]*NXS-XA[k]*NXA));
279  C2_5[k] = (N[k]*(SDNCV*(-30.0*eta*NDV*NDV+24.0*eta*V1_V22-MOR*(38.0+25.0*eta))+DM/M*SIGDNCV*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22
280  -MOR*(18.0+14.5*eta)))+NDV*v[k]*(SDNCV*(-9.0+9.0*eta)+DM/M* SIGDNCV*(-3.0+6.0*eta))+NCV[k]*(NDV*VDS*(-3.0+3.0*eta)
281  -8.0*MOR*eta*NDS-DM/M*(4.0*MOR*eta*NDSIG+3.0*NDV*VDSIG))+NDV*NCS[k]*(-22.5*eta*NDV*NDV+21.0*eta*V1_V22-MOR*(25.0+15.0*eta))
282  +DM/M*NDV*NCSIG[k]*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22-MOR*(9.0+8.5*eta))+VCS[k]*(16.5*eta*NDV*NDV+MOR*(21.0+9.0*eta)
283  -14.0*eta*V1_V22)+DM/M*VCSIG[k]*(9.0*eta*NDV*NDV-7.0*eta*V1_V22+MOR*(9.0+4.5*eta)))/r3;
284  }
285 
286  } /* if(Van_Spin==1) */
287 
288 
289 
290 
291 
292  for(k=0;k<3;k++)
293  {
294  // Newton ~1/c^0
295 
296  a_pn1[0][k] = -m2*x[k]/r3;
297  a_pn2[0][k] = m1*x[k]/r3;
298 
299  if(usedOrNot[0] == 1) // PN1 ~1/c^2
300  {
301  a_pn1[1][k] = ((AK2*N[k] + BK2*v[k])/r2)*m2;
302  a_pn2[1][k] = -((AK2*N[k] + BK2*v[k])/r2)*m1;
303  }
304 
305  if(usedOrNot[1] == 1) // PN2 ~1/c^4
306  {
307  a_pn1[2][k] = ((AK4*N[k] + BK4*v[k])/r2)*m2;
308  a_pn2[2][k] = -((AK4*N[k] + BK4*v[k])/r2)*m1;
309  }
310 
311  if(usedOrNot[2] == 1) // PN2.5 ~1/c^5
312  {
313  a_pn1[3][k] = ((AK5*N[k] + BK5*v[k])/r2)*m2;
314  a_pn2[3][k] = -((AK5*N[k] + BK5*v[k])/r2)*m1;
315  }
316 
317  if(usedOrNot[3] == 1) // PN3 ~1/c^6
318  {
319  a_pn1[4][k] = ((AK6*N[k] + BK6*v[k])/r2)*m2;
320  a_pn2[4][k] = -((AK6*N[k] + BK6*v[k])/r2)*m1;
321  }
322 
323  if(usedOrNot[4] == 1) // PN3.5 ~1/c^7
324  {
325  a_pn1[5][k] = ((AK7*N[k] + BK7*v[k])/r2)*m2;
326  a_pn2[5][k] = -((AK7*N[k] + BK7*v[k])/r2)*m1;
327  }
328 
329  A[k] = MOR*((AK2+AK4+AK5+AK6+AK7)*N[k] + (BK2+BK4+BK5+BK6+BK7)*v[k])/r + C1_5[k]/c_2 + C2[k]/c_4 + C2_5[k]/c_4;
330  }
331 
332  // PN accelerations
333 
334 
335 
336 
337  // PN jerks
338 
339 
340  AT[0] = A[0] - MOR*N[0]/r; // miert van AT - ?
341  AT[1] = A[1] - MOR*N[1]/r;
342  AT[2] = A[2] - MOR*N[2]/r;
343 
344  /*
345  AT[0] = A[0];
346  AT[1] = A[1];
347  AT[2] = A[2];
348  */
349 
350  RPP = V1_V22/r + AT[0]*N[0]+AT[1]*N[1] + AT[2]*N[2] - RP*RP/r;
351  VA = AT[0]*v[0] + AT[1]*v[1] + AT[2]*v[2];
352 
353  for(k=0;k<3;k++) NDOT[k] = (v[k]-N[k]*RP)/r;
354 
355  NVDOT = NDOT[0]*v[0]+NDOT[1]*v[1]+NDOT[2]*v[2]+N[0]*AT[0]+N[1]*AT[1]+N[2]*AT[2];
356 
357  //NDOTCV crossproduct of NDOT[k] and relative v = NDOT[k]Xv[j]
358  NDOTCV[0] = NDOT[1]*v[2] - NDOT[2]*v[1];
359  NDOTCV[1] = NDOT[2]*v[0] - NDOT[0]*v[2];
360  NDOTCV[2] = NDOT[0]*v[1] - NDOT[1]*v[0];
361 
362  //NCA crossproduct of N and AT = N[k]XAT[j]
363  NCA[0] = N[1]*AT[2] - N[2]*AT[1];
364  NCA[1] = N[2]*AT[0] - N[0]*AT[2];
365  NCA[2] = N[0]*AT[1] - N[1]*AT[0];
366 
367  ADK2 = 0.0;
368  BDK2 = 0.0;
369  ADK4 = 0.0;
370  BDK4 = 0.0;
371  ADK5 = 0.0;
372  BDK5 = 0.0;
373  ADK6 = 0.0;
374  BDK6 = 0.0;
375  ADK7 = 0.0;
376  BDK7 = 0.0;
377  C1_5D[0] = 0.0;
378  C1_5D[1] = 0.0;
379  C1_5D[2] = 0.0;
380  C2D[0] = 0.0;
381  C2D[1] = 0.0;
382  C2D[2] = 0.0;
383  C2_5D[0] = 0.0;
384  C2_5D[1] = 0.0;
385  C2_5D[2] = 0.0;
386 
387  if(usedOrNot[0] == 1) // PN1 ~1/c^2
388  {
389  A1D = -2.0*(2.0+eta)*MOR*RP/r - 2.0*(1.0+3.0*eta)*VA + 3.0*eta*RP*RPP;
390  B1D = 2.0*(2.0-eta)*RPP;
391 
392  ADK2 = A1D/c_2;
393  BDK2 = B1D/c_2;
394  }
395 
396  if(usedOrNot[1] == 1) // PN2 ~1/c^4
397  {
398  A2D = 1.5*(12.0+29.0*eta)*MOR*MOR*RP/r -eta*(3.0-4.0*eta)*4.0*V1_V22*VA - 7.5*eta*(1.0-3.0*eta)*RPP -0.5*eta*(13.0-4.0*eta)*MOR*RP*V1_V22/r+eta*(13.0-4.0*eta)*MOR*VA -(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RP*RP/r+2.0*(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RPP + 3.0*eta*(3.0-4.0*eta)*VA*RP*RP + 3.0*eta*(3.0-4.0*eta)*V1_V22*RP*RPP;
399  B2D = -0.5*RPP*((4.0+41.0*eta+8.0*eta*eta)*MOR - eta*(15.0+4.0*eta)*V1_V22+3.0*eta*(3.0+2.0*eta)*RP*RP) - 0.5*RP*(-(4.0+41.0*eta+8.0*eta*eta)*MOR*RP/r - 2.0*eta*(15.0+4.0*eta)*VA + 6.0*eta*(3.0+2.0*eta)*RP*RPP);
400 
401  ADK4 = A2D/c_4;
402  BDK4 = B2D/c_4;
403  }
404 
405  if(usedOrNot[2] == 1) // PN2.5 ~1/c^5
406  {
407  A2_5D = -1.6*eta*MOR*RP*RP*(17.0/3.0*MOR+3.0*V1_V22)/r +1.6*eta*MOR*RPP*(17.0/3.0*MOR+3.0*V1_V22)+1.6*eta*MOR*RP*(-17.0*MOR*RP/3.0/r+6.0*VA);
408  B2_5D = 1.6*eta*MOR*RP*(3.0*MOR+V1_V22)/r - 1.6*eta*MOR*(-3.0*MOR*RP/r+2.0*VA);
409 
410  ADK5 = A2_5D/c_5;
411  BDK5 = B2_5D/c_5;
412  }
413 
414  if(usedOrNot[3] == 1) // PN3 ~1/c^6
415  {
416  A3D = 6.0*eta*RP*RP*RP*RP*RP*RPP*(35.0-175.0*eta+175.0*eta*eta)/16.0 + eta*(4.0*RP*RP*RP*RPP*V1_V22 + 2.0*RP*RP*RP*RP*VA)*(-15.0+135.0*eta/2.0-255.0*eta*eta/4.0)/2.0 + eta*(2.0*RP*RPP*V1_V22*V1_V22+4.0*RP*RP*V1_V22*VA)/2.0*(15.0-237.0*eta/2.0+45.0*eta*eta) + 6.0*V1_V22*V1_V22*VA*eta*(-11.0/4.0-49.0*eta/4.0-13.0*eta*eta) + MOR*(4.0*RP*RP*RP*RPP*eta*(-79.0+69.0/2.0*eta+30.0*eta*eta) + eta*(2.0*RP*RPP*V1_V22+2.0*RP*RP*VA)*(121.0-16.0*eta-20.0*eta*eta)+4.0*V1_V22*VA*eta*(-75.0/4.0-8.0*eta+10.0*eta*eta)) - MOR*RP*((-79.0+69.0*eta/2.0+30.0*eta*eta)*RP*RP*RP*RP*eta+eta*RP*RP*V1_V22*(121.0-16.0*eta-20.0*eta*eta)+eta*V1_V22*V1_V22*(-75.0/4.0-8.0*eta+10.0*eta*eta))/r - 2.0*MOR*MOR*RP*(RP*RP*((-1.0-615.0*PI2*eta/64.0)-22717.0*eta/168.0-11.0*eta*eta/8.0+7.0*eta*eta*eta)+eta*V1_V22*((20827.0/840.0+123.0*PI2/64.0)-eta*eta))/r + MOR*MOR*(2.0*RP*RPP*((-1.0-615*PI2*eta/64.0)-22717.0*eta/168.0-11.0*eta*eta/8.0+7*eta*eta*eta)+2.0*eta*VA*((20827.0/840.0 +123.0*PI2/64.0)-eta*eta)) - 3.0*MOR*MOR*MOR*RP*(16.0+(1399.0/12.0-41.0*PI2/16.0)*eta+71.0*eta*eta/2.0)/r;
417  B3D = 75.0*RP*RP*RP*RP*RPP*eta*(3.0/8.0-eta-.25*eta*eta)+eta*(3.0*RP*RP*RPP*V1_V22+2.0*RP*RP*RP*VA)*(-12.0+111.0*eta/4.0+12.0*eta*eta)+eta*(RPP*V1_V22*V1_V22+4.0*RP*V1_V22*VA)*(65.0/8.0-19.0*eta-6.0*eta*eta)-MOR*RP*(RP*RP*RP*eta*(-329.0/6.0-59.0*eta/2.0-18.0*eta*eta)+RP*V1_V22*eta*(15.0+27.0*eta+10.0*eta*eta))/r+MOR*(3.0*RP*RP*RPP*eta*(-329.0/6.0-59.0*eta/2.0-18.0*eta*eta)+eta*(RPP*V1_V22+2.0*RP*VA)*(15.0+27.0*eta+10.0*eta*eta))-2.0*MOR*MOR*RP*(RP*((4.0+123.0*PI2*eta/32.0)+5849.0*eta/840.0-25.0*eta*eta-8.0*eta*eta*eta))/r+MOR*MOR*(RPP*((4.0+123.0*PI2*eta/32.0)+5849.0/840.0*eta-25.0*eta*eta-8.0*eta*eta*eta));
418 
419  ADK6 = A3D/c_6;
420  BDK6 = B3D/c_6;
421  }
422 
423  if(usedOrNot[4] == 1) // PN3.5 ~1/c^7
424  {
425  A3_5D = MOR*eta*(-RP*(V1_V22*V1_V22*(-366.0/35.0-12.0*eta)+V1_V22*RP*RP*(114.0+12.0*eta)+RP*RP*RP*RP*(-112.0))/r+4.0*V1_V22*VA*(-366.0/35.0-12.0*eta)+2.0*(VA*RP*RP+RP*RPP*V1_V22)*(114.0+12.0*eta)+4.0*RP*RP*RP*RPP*(-112.0)+MOR*(2.0*VA*(-692.0/35.0+724.0*eta/15.0)+2.0*RP*RPP*(-294.0/5.0-376.0*eta/5.0)-2.0*RP*(V1_V22*(-692.0/35.0+724.0*eta/15.0)+RP*RP*(-294.0/5.0-376.0*eta/5.0))/r-3.0*MOR*RP*(-3956.0/35.0-184.0*eta/5.0)/r));
426  B3_5D = MOR*eta*(4.0*V1_V22*VA*(626.0/35.0+12.0*eta/5.0)+2.0*(VA*RP*RP+V1_V22*RP*RPP)*(-678.0/5.0-12.0*eta/5.0)+4.0*RP*RP*RP*RPP*120.0-RP*(V1_V22*V1_V22*(626.0/35.0+12.0*eta/5.0)+V1_V22*RP*RP*(-678.0/5.0-12.0*eta/5.0)+120.0*RP*RP*RP*RP)/r+MOR*(2.0*VA*(-164.0/21.0-148.0*eta/5.0)+2*RP*RPP*(82.0/3.0+848.0*eta/15.0)-2.0*RP*(V1_V22*(-164.0/21-148.0*eta/5.0)+RP*RP*(82.0/3.0+848.0*eta/15.0))/r-3.0*MOR*RP*(1060.0/21.0+104.0*eta/5.0)/r));
427 
428  ADK7 = A3_5D/c_7;
429  BDK7 = B3_5D/c_7;
430  }
431 
432 
433  for(k=0;k<3;k++)
434  {
435  // Newton ~1/c^0
436 
437  adot_pn1[0][k] = -m2*(v[k]/r3 - 3.0*RP*x[k]/r2/r2);
438  adot_pn2[0][k] = m1*(v[k]/r3 - 3.0*RP*x[k]/r2/r2);
439 
440  if(usedOrNot[0] == 1) // PN1 ~1/c^2
441  {
442  adot_pn1[1][k] = (-2.0*MOR*RP*(AK2*N[k]+BK2*v[k])/r2 + MOR*(ADK2*N[k]+BDK2*v[k])/r + MOR*(AK2*(v[k]-N[k]*RP)/r+BK2*AT[k])/r)*m2/M;
443  adot_pn2[1][k] = -(-2.0*MOR*RP*(AK2*N[k]+BK2*v[k])/r2 + MOR*(ADK2*N[k]+BDK2*v[k])/r + MOR*(AK2*(v[k]-N[k]*RP)/r+BK2*AT[k])/r)*m1/M;
444  }
445 
446  if(usedOrNot[1] == 1) // PN2 ~1/c^4
447  {
448  adot_pn1[2][k] = (-2.0*MOR*RP*(AK4*N[k]+BK4*v[k])/r2 + MOR*(ADK4*N[k]+BDK4*v[k])/r + MOR*(AK4*(v[k]-N[k]*RP)/r+BK4*AT[k])/r)*m2/M;
449  adot_pn2[2][k] = -(-2.0*MOR*RP*(AK4*N[k]+BK4*v[k])/r2 + MOR*(ADK4*N[k]+BDK4*v[k])/r + MOR*(AK4*(v[k]-N[k]*RP)/r+BK4*AT[k])/r)*m1/M;
450  }
451 
452  if(usedOrNot[2] == 1) // PN2.5 ~1/c^5
453  {
454  adot_pn1[3][k] = (-2.0*MOR*RP*(AK5*N[k]+BK5*v[k])/r2 + MOR*(ADK5*N[k]+BDK5*v[k])/r + MOR*(AK5*(v[k]-N[k]*RP)/r+BK5*AT[k])/r)*m2/M;
455  adot_pn2[3][k] = -(-2.0*MOR*RP*(AK5*N[k]+BK5*v[k])/r2 + MOR*(ADK5*N[k]+BDK5*v[k])/r + MOR*(AK5*(v[k]-N[k]*RP)/r+BK5*AT[k])/r)*m1/M;
456  }
457 
458  if(usedOrNot[3] == 1) // PN3 ~1/c^6
459  {
460  adot_pn1[4][k] = (-2.0*MOR*RP*(AK6*N[k]+BK6*v[k])/r2 + MOR*(ADK6*N[k]+BDK6*v[k])/r + MOR*(AK6*(v[k]-N[k]*RP)/r+BK6*AT[k])/r)*m2/M;
461  adot_pn2[4][k] = -(-2.0*MOR*RP*(AK6*N[k]+BK6*v[k])/r2 + MOR*(ADK6*N[k]+BDK6*v[k])/r + MOR*(AK6*(v[k]-N[k]*RP)/r+BK6*AT[k])/r)*m1/M;
462  }
463 
464  if(usedOrNot[4] == 1) // PN3.5 ~1/c^7
465  {
466  adot_pn1[5][k] = (-2.0*MOR*RP*(AK7*N[k]+BK7*v[k])/r2 + MOR*(ADK7*N[k]+BDK7*v[k])/r + MOR*(AK7*(v[k]-N[k]*RP)/r+BK7*AT[k])/r)*m2/M;
467  adot_pn2[5][k] = -(-2.0*MOR*RP*(AK7*N[k]+BK7*v[k])/r2 + MOR*(ADK7*N[k]+BDK7*v[k])/r + MOR*(AK7*(v[k]-N[k]*RP)/r+BK7*AT[k])/r)*m1/M;
468  }
469 
470 
471  }
472 
473  // PN jerks
474 
475  if(Van_Spin==1)
476  {
477  double L[3];
478  //L crossproduct of x[k] and relative v = x[k]Xv[j]
479  L[0] = x[1]*v[2] - x[2]*v[1];
480  L[1] = x[2]*v[0] - x[0]*v[2];
481  L[2] = x[0]*v[1] - x[1]*v[0];
482 
483  LABS = sqrt(L[0]*L[0]+L[1]*L[1]+L[2]*L[2]);
484 
485  LU[0] = L[0]/LABS;
486  LU[1] = L[1]/LABS;
487  LU[2] = L[2]/LABS;
488 
489  S1DLU = S1[0]*LU[0]+S1[1]*LU[1]+S1[2]*LU[2];
490  S2DLU = S2[0]*LU[0]+S2[1]*LU[1]+S2[2]*LU[2];
491 
492  for(k=0;k<3;k++)
493  {
494  SU1[k] = MOR*eta*(N[k]*(-4.0*VDS-2.0*DM/M*VDSIG)+ v[k]*(3.0*NDS+DM/M*NDSIG)+NDV*(2.0*KSS[k]+DM/M*KSSIG[k])) /r;
495  SV1[k] = MOR*(N[k]*(VDSIG*(-2.0+4.0*eta)-2.0*DM/M*VDS)+ v[k]*(NDSIG*(1.0-eta)+DM/M*NDS)+NDV*(KSSIG[k]*(1.0- 2.0*eta)+ DM/M*KSS[k]))/r;
496 
497  SS1[k] = 0.5*(L[k]*(4.0+3.0*(m2/m1))+ (S2[k]-3.0*S2DLU*LU[k]))/r3;
498  SS2[k] = 0.5*(L[k]*(4.0+3.0*(m1/m2))+ (S1[k]-3.0*S1DLU*LU[k]))/r3;
499 
500  SU2[k] = MOR*eta/r*(N[k]*(VDS*(-2.0*V1_V22+3.0*NDV*NDV- 6.0*eta*NDV*NDV+7.0*MOR-8.0*eta*MOR)-14.0*MOR*NDS*NDV+ DM/M*VDSIG*eta*(-3.0*NDV*NDV-4.0*MOR)+DM/M*MOR*NDSIG*NDV* (2.0-eta/2.))+v[k]*(NDS*(2.0*V1_V22-4.0*eta*V1_V22-3.0*NDV* NDV+7.5*eta*NDV*NDV+4.0*MOR-6.0*eta*MOR)+VDS*NDV*(2.0- 6.0*eta)+ DM/M*NDSIG*(-1.5*eta*V1_V22+3.0*eta*NDV*NDV-MOR-3.5*eta* MOR)-3.0*DM/M*VDSIG*NDV*eta)+KSS[k]*NDV*(V1_V22-2.0*eta* V1_V22-1.5*NDV*NDV+3.0*eta*NDV*NDV-MOR+2.0*eta*MOR)+ DM/M*KSSIG[k]*NDV*(-eta*V1_V22+1.5*eta*NDV*NDV+ (eta-1.)*MOR));
501  SV2[k] = MOR/r*(N[k]*(VDSIG*eta*(-2.0*V1_V22+6.0*eta*NDV* NDV+(3.0+8.0*eta)*MOR)+MOR*NDSIG*NDV*(2.0-22.5*eta+2.0* eta*eta)+ DM/M*VDS*eta*(-3.0*NDV*NDV-4.0*MOR)+DM/M*MOR*NDS*NDV*(2.0- 0.5*eta))+v[k]*(NDSIG*(0.5*eta*V1_V22+2.0*eta*eta*V1_V22- 4.5*eta*eta*NDV*NDV+(4.5*eta-1.0+8.0*eta*eta)*MOR)+VDSIG*NDV* eta*(6.0*eta-1.)-3.0*DM/M*VDS*NDV*eta+DM/M*NDS*(-1.5* eta*V1_V22+ 3.0*eta*NDV*NDV-(1.0+3.5*eta)*MOR))+KSSIG[k]*NDV*(2.0*eta*eta* V1_V22-3.0*eta*eta*NDV*NDV+(-1.0+4.0*eta-2.0*eta*eta)*MOR)+ DM/M*KSS[k]*NDV*(-eta*V1_V22+1.5*eta*NDV*NDV+(-1.0+eta)* MOR));
502  }
503 
504  //SS1 crossproduct of SS1 and S1 = SS1[k]XS1[j]
505  SS1aux[0] = SS1[1]*S1[2] - SS1[2]*S1[1];
506  SS1aux[1] = SS1[2]*S1[0] - SS1[0]*S1[2];
507  SS1aux[2] = SS1[0]*S1[1] - SS1[1]*S1[0];
508 
509  SS1[0] = SS1aux[0];
510  SS1[1] = SS1aux[1];
511  SS1[2] = SS1aux[2];
512 
513  //SS2 crossproduct of SS2 and S2 = SS2[k]XS2[j]
514  SS2aux[0] = SS2[1]*S2[2] - SS2[2]*S2[1];
515  SS2aux[1] = SS2[2]*S2[0] - SS2[0]*S2[2];
516  SS2aux[2] = SS2[0]*S2[1] - SS2[1]*S2[0];
517 
518  SS2[0] = SS2aux[0];
519  SS2[1] = SS2aux[1];
520  SS2[2] = SS2aux[2];
521 
522 
523  for(k=0;k<3;k++)
524  {
525  SU[k] = SU1[k]/c_2 + SU2[k]/c_4 + (SS1[k] + SS2[k])/c_2;
526  SV[k] = SV1[k]/c_2 + SV2[k]/c_4+M*(SS2[k]/m2-SS1[k]/ m1)/c_2;
527 
528  KSS[k] = KSS[k] + SU[k]*dt_bh; // integrate for dt_bh timestep
529  KSSIG[k] = KSSIG[k] + SV[k]*dt_bh;
530 
531  SPIN[k][0] = m1*(M*KSS[k]-m2*KSSIG[k])/M/M/m1/m1*c_1 ;
532  SPIN[k][1] = m2*(M*KSS[k]+m1*KSSIG[k])/M/M/m2/m2*c_1;
533  XAD[k] = 0.5/(M*M*m1*m2)*(-SU[k]*M*DM-SV[k]*(m1*m1+m2*m2));
534  XSD[k] = 0.5/(M*M*m1*m2)*(SU[k]*M*M+SV[k]*(m1*m1-m2*m2));
535 
536  }
537 
538 
539  //NDOTCS crossproduct of NDOT and KSS = NDOT[k]XKSS[j]
540  NDOTCS[0] = NDOT[1]*KSS[2] - NDOT[2]*KSS[1];
541  NDOTCS[1] = NDOT[2]*KSS[0] - NDOT[0]*KSS[2];
542  NDOTCS[2] = NDOT[0]*KSS[1] - NDOT[1]*KSS[0];
543  //NCSU crossproduct of N and SU = N[k]XSU[j]
544  NCSU[0] = N[1]*SU[2] - N[2]*SU[1];
545  NCSU[1] = N[2]*SU[0] - N[0]*SU[2];
546  NCSU[2] = N[0]*SU[1] - N[1]*SU[0];
547  //NDOTCSIG crossproduct of NDOT and KSSIG = NDOT[k]XKSSIG[j]
548  NDOTCSIG[0] = NDOT[1]*KSSIG[2] - NDOT[2]*KSSIG[1];
549  NDOTCSIG[1] = NDOT[2]*KSSIG[0] - NDOT[0]*KSSIG[2];
550  NDOTCSIG[2] = NDOT[0]*KSSIG[1] - NDOT[1]*KSSIG[0];
551  //NCSV crossproduct of N and SV = N[k]XSV[j]
552  NCSV[0] = N[1]*SV[2] - N[2]*SV[1];
553  NCSV[1] = N[2]*SV[0] - N[0]*SV[2];
554  NCSV[2] = N[0]*SV[1] - N[1]*SV[0];
555  //ACS crossproduct of AT and KSS = AT[k]XKSS[j]
556  ACS[0] = AT[1]*KSS[2] - AT[2]*KSS[1];
557  ACS[1] = AT[2]*KSS[0] - AT[0]*KSS[2];
558  ACS[2] = AT[0]*KSS[1] - AT[1]*KSS[0];
559  //VCSU crossproduct of relative v and SU = v[k]XSU[j]
560  VCSU[0] = v[1]*SU[2] - v[2]*SU[1];
561  VCSU[1] = v[2]*SU[0] - v[0]*SU[2];
562  VCSU[2] = v[0]*SU[1] - v[1]*SU[0];
563  //ACSIG crossproduct of AT and KSSIG = AT[k]XKSSIG[j]
564  ACSIG[0] = AT[1]*KSSIG[2] - AT[2]*KSSIG[1];
565  ACSIG[1] = AT[2]*KSSIG[0] - AT[0]*KSSIG[2];
566  ACSIG[2] = AT[0]*KSSIG[1] - AT[1]*KSSIG[0];
567  //VCSV crossproduct of relative v and SV = v[k]XSV[j]
568  VCSV[0] = v[1]*SV[2] - v[2]*SV[1];
569  VCSV[1] = v[2]*SV[0] - v[0]*SV[2];
570  VCSV[2] = v[0]*SV[1] - v[1]*SV[0];
571 
572  SNVDOT = SU[0]*NCV[0]+SU[1]*NCV[1]+SU[2]*NCV[2]+ KSS[0]*NDOTCV[0]+KSS[1]*NDOTCV[1]+KSS[2]*NDOTCV[2]+ KSS[0]*NCA[0]+KSS[1]*NCA[1]+KSS[2]*NCA[2];
573 
574  SIGNVDOT = SV[0]*NCV[0]+SV[1]*NCV[1]+SV[2]*NCV[2]+ KSSIG[0]*NDOTCV[0]+KSSIG[1]*NDOTCV[1]+KSSIG[2]*NDOTCV[2]+ KSSIG[0]*NCA[0]+KSSIG[1]*NCA[1]+KSSIG[2]*NCA[2];
575 
576  NSDOT = NDOT[0]*KSS[0]+NDOT[1]*KSS[1]+NDOT[2]*KSS[2]+ N[0]*SU[0]+N[1]*SU[1]+N[2]*SU[2];
577  NSIGDOT = NDOT[0]*KSSIG[0]+NDOT[1]*KSSIG[1]+NDOT[2]*KSSIG[2]+ N[0]*SV[0]+N[1]*SV[1]+N[2]*SV[2];
578  VSDOT = AT[0]*KSS[0]+AT[1]*KSS[1]+AT[2]*KSS[2]+ v[0]*SU[0]+v[1]*SU[1]+v[2]*SU[2];
579  VSIGDOT = AT[0]*KSSIG[0]+AT[1]*KSSIG[1]+AT[2]*KSSIG[2]+ v[0]*SV[0]+v[1]*SV[1]+v[2]*SV[2];
580 
581  NXSDOT = NDOT[0]*XS[0]+NDOT[1]*XS[1]+NDOT[2]*XS[2]+ N[0]*XSD[0]+N[1]*XSD[1]+N[2]*XSD[2];
582  NXADOT = NDOT[0]*XA[0]+NDOT[1]*XA[1]+NDOT[2]*XA[2]+ N[0]*XAD[0]+N[1]*XAD[1]+N[2]*XAD[2];
583 
584  for(k=0;k<3;k++)
585  {
586  C1_5D[k] = -3.0*RP/r*C1_5[k]+(NDOT[k]*(12.0*SDNCV+6.0*DM/M* SIGDNCV)+N[k]*(12.0*SNVDOT+6.0*DM/M*SIGNVDOT)+9.0*NVDOT* NCS[k]+9.0*NDV*(NDOTCS[k]+NCSU[k])+3.0*DM/M*(NVDOT*NCSIG[k]+ NDV*(NDOTCSIG[k]+NCSV[k]))-7.0*(ACS[k]+VCSU[k])-3.0*DM/M* (ACSIG[k]+VCSV[k]))/(r3);
587  C2D[k] = -4.0*RP/r*C2[k]-MOR*MOR*MOR*3.0*eta/r*(NDOT[k]* (XS2-XA2-5.0*NXS*NXS+5.0*NXA*NXA)+N[k]*(2.0*(XS[0]*XSD[0]+ XS[1]*XSD[1]+XS[2]*XSD[2]-XA[0]*XAD[0]-XA[1]*XAD[1]- XA[2]*XAD[2])-10.0*NXS*NXSDOT+10.0*NXA*NXADOT)+2.0*(XSD[k]* NXS+XS[k]*NXSDOT-XAD[k]*NXA-XA[k]*NXADOT));
588  C2_5D[k] = -3.0*RP/r*C2_5[k]+(NDOT[k]*(SDNCV*(-30.0*eta* NDV*NDV+24.0*eta*V1_V22-MOR*(38.0+25.0*eta))+DM/M*SIGDNCV* (-15.0*eta*NDV*NDV+12.0*eta*V1_V22-MOR*(18.0+14.5*eta)))+ N[k]*(SNVDOT*(-30.0*eta*NDV*NDV+24.0*eta*V1_V22-MOR* (38.0+25.0*eta))+SDNCV*(-60.0*eta*NDV*NVDOT+48.0*eta*VA+ MOR*RP/r*(38.0+25.0*eta))+DM/M*SIGNVDOT*(-15.0*eta*NDV* NDV+12.0*eta*V1_V22-MOR*(18.0+14.5*eta))+DM/M*SIGDNCV* (-30.0*eta*NDV*NVDOT+24.0*eta*VA+MOR*RP/r*(18.0+14.5*eta)))+ (NVDOT*v[k]+NDV*AT[k])*(SDNCV*(-9.0+9.0*eta)+DM/M*SIGDNCV* (-3.0+6.0*eta))+NDV*v[k]*(SNVDOT*(-9.0+9.0*eta)+DM/M* SIGNVDOT*(-3.0+6.0*eta))+(NDOTCV[k]+NCA[k])*(NDV*VDS*(-3.0+ 3.0*eta)-8.0*MOR*eta*NDS-DM/M*(4.0*MOR*eta*NDSIG+3.0*NDV*VDSIG) )+NCV[k]*((NVDOT*VDS+NDV*VSDOT)*(-3.0+3.0*eta)-8.0*eta*MOR* (NSDOT-RP/r*NDS)-DM/M*(4.0*eta*MOR*(NSIGDOT-RP/r*NDSIG)+ 3.0*(NVDOT*VDSIG+NDV*VSIGDOT)))+(NVDOT*NCS[k]+NDV* (NDOTCS[k]+NCSU[k]))*(-22.5*eta*NDV*NDV+21.0*eta*V1_V22- MOR*(25.0+15.0*eta))+NDV*NCS[k]*(-45.0*eta*NDV*NVDOT+42.0*eta* VA+MOR*RP/r*(25.0+15.0*eta))+DM/M*(NVDOT*NCSIG[k]+NDV* (NDOTCSIG[k]+NCSV[k]))*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22- MOR*(9.0+8.5*eta))+DM/M*NDV*NCSIG[k]*(-30.0*eta*NDV*NVDOT+ 24.0*eta*VA+MOR*RP/r*(9.0+8.5*eta))+(ACS[k]+VCSU[k])* (16.5*eta*NDV*NDV+MOR*(21.0+9.0*eta)-14.0*eta*V1_V22)+ VCS[k]*(33.0*eta*NDV*NVDOT-MOR*RP/r*(21.0+9.0*eta)- 28.0*eta*VA)+DM/M*(ACSIG[k]+VCSV[k])*(9.0*eta*NDV*NDV- 7.0*eta*V1_V22+MOR*(9.0+4.5*eta))+DM/M*VCSIG[k]*(18.0* eta*NDV*NVDOT-14.0*eta*VA-MOR*RP/r*(9.0+4.5*eta)))/ (r3);
589  }
590 
591 
592  } /* if(Van_Spin==1) */
593 
594  ADK = ADK2+ADK4+ADK5+ADK6+ADK7;
595  BDK = BDK2+BDK4+BDK5+BDK6+BDK7;
596 
597  KSAK = AK2+AK4+AK5+AK6+AK7;
598  KSBK = BK2+BK4+BK5+BK6+BK7;
599 
600  for(k=0;k<3;k++) AD[k] = -2.0*MOR*RP*(KSAK*N[k]+KSBK*v[k])/r2 + MOR*(ADK*N[k]+BDK*v[k])/r + MOR*(KSAK*(v[k]-N[k]*RP)/r+KSBK*AT[k])/r + C1_5D[k]/c_2 + C2D[k]/c_4 +C2_5D[k]/c_4;
601 
602 
603  for(k=0;k<3;k++) // new values of the BH's spins, returned back to the main program...
604  {
605  spin1[k] = SPIN[k][0];
606  spin2[k] = SPIN[k][1];
607  }
608 
609 
610 
611  // Check RS_DIST conditions !!!
612 
613  RS_DIST = 4.0*(2.0*m1/c_2 + 2.0*m2/c_2);
614 
615  if(r < RS_DIST)
616  {
617  //if(myRank == rootRank)
618  //if(is_root())
619  //{
620  // fprintf(stdout,"PN RSDIST: r = %.8E \t RS = %.8E \n", r, RS_DIST);
621  // fflush(stdout);
622  //}
623  return(505);
624  }
625  else
626  {
627  return(0);
628  }
629 
630 
631 }
632 /***************************************************************************/
calc_force_pn_BH
int calc_force_pn_BH(double m1, double xx1[], double vv1[], double spin1[], double m2, double xx2[], double vv2[], double spin2[], double CCC_NB, int usedOrNot[], double dt_bh, double a_pn1[][3], double adot_pn1[][3], double a_pn2[][3], double adot_pn2[][3])
Definition: pn_BH.h:11
SQR
double SQR(const double x)
Definition: pn_BH.h:8
is_root
bool is_root()