forked from Shallyn/pyWaveformGenerator
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpPrecRRForce.c
179 lines (171 loc) · 11 KB
/
pPrecRRForce.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
/**
* Writer: Xiaolin.liu
*
* This module contains basic functions for calculation.
* Functions list:
* Kernel:
* 20xx.xx.xx, LOC
**/
#include "pPrecRRForce.h"
void prec_CalculateRRForce(SEOBPrecVariables *vars,
REAL8 eta, REAL8 prTDot,
REAL8 fRRVec[], REAL8 fRRs1Vec[], REAL8 fRRs2Vec[])
{
REAL8 dm = sqrt(1. - 3.*eta);
REAL8 eta2, eta3, eta4;
eta2 = eta*eta;
eta3 = eta2*eta;
eta4 = eta2*eta2;
REAL8 xpS1, xpS2, pS1, pS2, xS1, xS2, pS12, pS22, xS12, xS22;
xpS1 = vars->xpS1;
xpS2 = vars->xpS2;
pS1 = vars->pTS1;
pS2 = vars->pTS2;
xS1 = vars->xS1;
xS2 = vars->xS2;
pS12 = pS1*pS1;
pS22 = pS2*pS2;
xS12 = xS1*xS1;
xS22 = xS2*xS2;
REAL8 S1S1, S1S2, S2S2;
S1S1 = vars->S1S1;
S1S2 = vars->S1S2;
S2S2 = vars->S2S2;
REAL8 r, sqr, x, x2, x3, x4, x5, x6, x7, x8, x9, x11, x12;
REAL8 y, y2, y3, y4, y6;
REAL8 z, z2, z3, z4, z5, z6, z7, z8, z9, z10, z11, z12, z13, z15;
INT i;
r = vars->r;
sqr = sqrt(r);
x = 1./sqr;
x2 = 1./r;
x3 = x2*x;
x4 = x2*x2;
x5 = x3*x2;
x6 = x4*x2;
x7 = x6*x;
x8 = x7*x;
x9 = x8*x2;
x11 = x9*x2;
x12 = x11*x;
y = sqr*vars->prT;
y2 = y*y;
y3 = y2*y;
y4 = y3*y;
y6 = y4*y2;
z3 = 1. + r*r*prTDot;
if (z3 > 0.)
z = cbrt(z3);
else
z = -cbrt(fabs(z3));
z2 = z*z;
z4 = z2*z2;
z5 = z4*z;
z6 = z3*z3;
z7 = z6*z;
z8 = z7*z;
z9 = z6*z3;
z10 = z5*z5;
z11 = z10*z;
z12 = z9*z3;
z13 = z12*z;
z15 = z13*z2;
REAL8 fxVec0, fxVec2, fxVec3, fxVec4;
REAL8 fpVec0, fpVec2, fpVec3, fpVec4;
REAL8 fLVec3, fLVec4;
fxVec0 = (-5*x3*y3)/(4*z2);
fxVec2 = (x5*y*(448*(11 - 5*z3)*z3 +
80*y4*(-7 - 40*z3 + 7*eta*(-1 + 2*z3)) +
y2*(-1120 + 3148*z3 - 6235*z4 - 3200*z6 +
4*eta*(280 + 4068*z3 - 1225*z4 + 280*z6))))/(1344.*z5);
fxVec3 = (-(x7*y3*(xpS1*(8 + 30*(3 + 2*eta)*y2 - 540*z3 -
360*eta*z3 + 75*z5 + 180*eta*z5 - dm*(8 + 90*y2 - 540*z3 + 75*z5)) +
xpS2*(8 + 30*(3 + 2*eta)*y2 - 540*z3 - 360*eta*z3 + 75*z5 + 180*eta*z5 +
dm*(8 + 90*y2 - 540*z3 + 75*z5)))))/(96.*eta*z5);
fxVec4 = (x5*y*(-105840*eta*y2*z3*(48*S2S2*x2 - 144*x6*xS12 + 288*x6*xS1*xS2 -
144*x6*xS22 - 168*S2S2*x2*z3 + 73*pS12*z4 - 142*pS1*pS2*z4 + 73*pS22*z4 +
99*S2S2*x2*z5 - 6*S1S2*x2*(16 - 56*z3 + 31*z5) + 3*S1S1*x2*(16 - 56*z3 + 33*z5) -
360*x6*xS12*z7 + 720*x6*xS1*xS2*z7 - 360*x6*xS22*z7) - 52920*y2*z3*(144*x6*xS12 -
144*dm*x6*xS12 + 144*x6*xS22 + 144*dm*x6*xS22 - 73*pS12*z4 + 73*dm*pS12*z4 -
73*pS22*z4 - 73*dm*pS22*z4 + 3*(-1 + dm)*S1S1*x2*(16 - 56*z3 + 33*z5) -
3*(1 + dm)*S2S2*x2*(16 - 56*z3 + 33*z5) + 360*x6*xS12*z7 - 360*dm*x6*xS12*z7 +
360*x6*xS22*z7 + 360*dm*x6*xS22*z7) - 336*eta4*x2*y2*(8400 - 36750*z10 -
19440*z12 + 82512*z3 - 14700*z4 + 240764*z6 + 60*y4*(35 - 119*z3 + 292*z6) -
419790*z7 + 74025*z8 + 167088*z9 - 6*y2*(1400 + 4496*z3 - 1225*z4 - 13666*z6 +
6125*z7 + 320*z9)) - 24*eta3*x2*(840*y6*(70 + 41*z3 + 32*z6) -
105*y4*z3*(10552 - 2227*z + 180488*z3 - 3985*z4 + 6920*z6) +
672*z3*(616 + 4021*z3 - 2695*z4 - 2953*z6 + 1225*z7 + 542*z9) -
y2*(235200 - 418425*z10 + 753480*z12 + 1238832*z3 + 56070*z4 +
11823112*z6 + 8326449*z7 - 688905*z8 + 21453264*z9)) +
eta2*(-51480576*(pS1 + pS2)*x3*(xS1 + xS2)*y3*z3 -
338688*pow(pS1 + pS2,2)*y2*z3*(-76 + 15*z4) +
338688*x6*pow(xS1 + xS2,2)*y2*z3*(15 + 86*y2 + 61*z3 + 15*z7) +
x2*(20160*y6*(-35 - 97*z3 + 204*z6) -
4032*z3*(-2464 + 10687*z3 - 13717*z4 - 5339*z6 + 6235*z7 + 566*z9) -
24*y4*(117600 + 504*(411 + 280*S1S1 + 560*S1S2 + 280*S2S2)*z3 +
130935*z4 + 468664*z6 + 1365465*z7 + 317520*z9) -
y2*(2822400 + 32771160*z10 + 11733120*z12 + 8064*(-1193 + 630*S1S1 + 1260*S1S2 + 630*S2S2)*z3 +
6284880*z4 + 36288*(1387 + 70*S1S1 + 140*S1S2 + 70*S2S2)*z6 - 25902684*z7 + 95013565*z8 + 21646464*z9)))))/(4064256.*eta2*z8);
fpVec0 = (22 + 29*y2 - 10*z3)/(12.*z2);
fpVec2 = (x2*(8*y4*(406 + 853*z3) + y2*(8960 - 46494*z3 + 36163*z4 + 292*z6) +
2*(2464 - 18079*z3 + 13717*z4 + 8699*z6 - 6235*z7 - 566*z9) + 4*eta*(y4*(812 + 68*z3) +
y2*(-1008 - 15512*z3 + 7105*z4 + 874*z6) - 2*(616 + 2635*z3 - 2695*z4 - 2323*z6 + 1225*z7 + 542*z9))))/(4032.*z5);
fpVec3 = (-(x4*(xpS1*(-4*eta*(-22 + 45*y4 - 119*z3 + 198*z5 + y2*(-14 - 267*z3 + 261*z5) + 33*z6 - 90*z8) +
3*dm*(90*y4 + y2*(-20 - 254*z3 + 145*z5) - 2*(22 + 6*z3 - 55*z5 + 2*z6 + 25*z8)) +
3*(-90*y4 + y2*(20 + 254*z3 - 145*z5) + 2*(22 + 6*z3 - 55*z5 + 2*z6 + 25*z8))) + xpS2*(-4*eta*(-22 + 45*y4 - 119*z3 + 198*z5 +
y2*(-14 - 267*z3 + 261*z5) + 33*z6 - 90*z8) + 3*(-90*y4 + y2*(20 + 254*z3 - 145*z5) + 2*(22 + 6*z3 - 55*z5 + 2*z6 + 25*z8)) +
3*dm*(-90*y4 + y2*(20 + 254*z3 - 145*z5) + 2*(22 + 6*z3 - 55*z5 + 2*z6 + 25*z8))))))/(288.*eta*z5);
fpVec4 = (x2*(10584*z3*(3168*x6*xS12 - 3168*dm*x6*xS12 + 3168*x6*xS22 + 3168*dm*x6*xS22 +
4176*x6*xS12*y2 - 4176*dm*x6*xS12*y2 + 4176*x6*xS22*y2 + 4176*dm*x6*xS22*y2 - 3600*x6*xS12*z10 +
3600*dm*x6*xS12*z10 - 3600*x6*xS22*z10 - 3600*dm*x6*xS22*z10 + 876*pS12*z3 - 876*dm*pS12*z3 + 876*pS22*z3 +
876*dm*pS22*z3 - 7488*x6*xS12*z3 + 7488*dm*x6*xS12*z3 - 7488*x6*xS22*z3 - 7488*dm*x6*xS22*z3 - 1606*pS12*z4 +
1606*dm*pS12*z4 - 1606*pS22*z4 - 1606*dm*pS22*z4 - 2117*pS12*y2*z4 + 2117*dm*pS12*y2*z4 - 2117*pS22*y2*z4 -
2117*dm*pS22*y2*z4 + 730*pS12*z7 - 730*dm*pS12*z7 + 730*pS22*z7 + 730*dm*pS22*z7 + 7920*x6*xS12*z7 - 7920*dm*x6*xS12*z7 +
7920*x6*xS22*z7 + 7920*dm*x6*xS22*z7 + 10440*x6*xS12*y2*z7 - 10440*dm*x6*xS12*y2*z7 + 10440*x6*xS22*y2*z7 + 10440*dm*x6*xS22*y2*z7 +
3*(-1 + dm)*S1S1*x2*(352 - 1090*z3 + 726*z5 + y2*(464 - 1728*z3 + 957*z5) + 342*z6 - 330*z8) - 3*(1 + dm)*S2S2*x2*(352 - 1090*z3 +
726*z5 + y2*(464 - 1728*z3 + 957*z5) + 342*z6 - 330*z8)) - 21168*eta*z3*(-1056*S2S2*x2 + 3168*x6*xS12 - 6336*x6*xS1*xS2 + 3168*x6*xS22 -
1392*S2S2*x2*y2 + 4176*x6*xS12*y2 - 8352*x6*xS1*xS2*y2 + 4176*x6*xS22*y2 - 3600*x6*xS12*z10 + 7200*x6*xS1*xS2*z10 - 3600*x6*xS22*z10 + 876*pS12*z3 -
1704*pS1*pS2*z3 + 876*pS22*z3 + 3270*S2S2*x2*z3 - 7488*x6*xS12*z3 + 14976*x6*xS1*xS2*z3 - 7488*x6*xS22*z3 + 5184*S2S2*x2*y2*z3 - 1606*pS12*z4 + 3124*pS1*pS2*z4 -
1606*pS22*z4 - 2117*pS12*y2*z4 + 4118*pS1*pS2*y2*z4 - 2117*pS22*y2*z4 - 2178*S2S2*x2*z5 - 2871*S2S2*x2*y2*z5 - 1026*S2S2*x2*z6 + 730*pS12*z7 - 1420*pS1*pS2*z7 +
730*pS22*z7 + 7920*x6*xS12*z7 - 15840*x6*xS1*xS2*z7 + 7920*x6*xS22*z7 + 10440*x6*xS12*y2*z7 - 20880*x6*xS1*xS2*y2*z7 + 10440*x6*xS22*y2*z7 - 3*S1S1*x2*(352 - 1090*z3 + 726*z5 +
y2*(464 - 1728*z3 + 957*z5) + 342*z6 - 330*z8) + 6*S1S2*x2*(352 - 1102*z3 + 682*z5 + y2*(464 - 1872*z3 + 899*z5) + 378*z6 - 310*z8) + 990*S2S2*x2*z8) + 1008*eta4*x2*(y6*(4060 - 7036*z3 + 79132*z6) +
2*y4*(-6580 - 19592*z3 + 7105*z4 + 118769*z6 - 5915*z7 + 43101*z9) - 2*(-6160 - 73465*z10 + 24675*z11 - 54084*z12 + 16520*z13 + 19912*z15 - 9432*z3 + 10780*z4 - 24320*z6 + 91735*z7 - 54285*z8 + 58124*z9) +
y2*(3920 + 21280*z10 - 834*z12 + 97080*z3 - 17640*z4 + 184884*z6 - 534590*z7 + 143115*z8 + 266762*z9)) + 72*eta3*x2*(-344960 + 4391988*z10 - 459270*z11 + 5221006*z12 - 658654*z13 - 4726694*z15 + 557032*z3 - 82236*z4 +
1794772*z6 + 28*y6*(4060 - 106*z3 + 324487*z6) - 7107782*z7 + 1010394*z8 + 404404*z9 + y4*(86240 - 1652084*z3 + 452081*z4 + 18918872*z6 + 426257*z7 + 12918402*z9) + y2*(-454720 + 284518*z10 + 17472*z12 -
996940*z3 + 234556*z4 + 20356*z6 - 14704354*z7 + 1331883*z8 + 28166278*z9)) - eta2*(-4064256*(pS1 + pS2)*x3*(xS1 + xS2)*y*(11 + 45*y2 - 15*z3)*z3 + 1016064*pow(pS1 + pS2,2)*z3*(y2*(90 - 29*z4) + 2*(11 - 5*z3 - 11*z4 + 5*z7)) +
338688*x6*pow(xS1 + xS2,2)*z3*(66 + 328*y4 - 30*z10 - 102*z3 + 66*z7 + y2*(197 + 10*z3 + 87*z7)) + x2*(2016*y6*(-2030 + 242*z3 + 60855*z6) - y2*(28788480 - 8664156*z10 + 4003776*z12 + 6048*(-18855 + 7336*S1S1 + 14672*S1S2 + 7336*S2S2)*z3 +
50279040*z4 + 6048*(-58049 + 2996*S1S1 + 5992*S1S2 + 2996*S2S2)*z6 - 474665562*z7 + 551078677*z8 - 418972176*z9) + 72*y4*(-270480 - 84*(-4821 + 3248*S1S1 + 6496*S1S2 + 3248*S2S2)*z3 - 253141*z4 + 8566796*z6 - 810550*z7 + 1964298*z9) +
2*(-6209280 - 81288189*z10 + 95013565*z11 - 11256840*z12 + 3209778*z13 - 27885816*z15 - 22176*(-1501 + 504*S1S1 + 1008*S1S2 + 504*S2S2)*z3 - 13826736*z4 + 336*(-152807 + 25704*S1S1 + 51408*S1S2 + 25704*S2S2)*z6 +
175875633*z7 - 209029843*z8 + 1008*(92713 + 2520*S1S1 + 5040*S1S2 + 2520*S2S2)*z9)))))/(12192768.*eta2*z8);
fLVec3 = (x4*y2*(45*(pS1 - dm*pS1 + pS2 + dm*pS2)*z3 - 10*eta2*(pS1*(1 + 3*y2 - 3*z3) + pS2*(1 + 3*y2 - 3*z3) -
x3*(xS1 + xS2)*y*(1 + 3*y2 + 3*z3)) + eta*(pS1*(-19 - 45*y2 + dm*(19 + 45*y2 - 75*z3) - 75*z3) -
x3*((-1 + dm)*xS1 - (1 + dm)*xS2)*y*(19 + 45*y2 + 45*z3) - pS2*(19 + 45*y2 + dm*(19 + 45*y2 - 75*z3) + 75*z3))))/(48.*eta2*z5);
fLVec4 =(-61*x6*(pS2*xpS2 + pS1*(xpS1 + 2*xpS2) - x3*(xpS1*xS1 + 2*xpS2*xS1 + xpS2*xS2)*y)*y2)/(12.*z5);
REAL8 fLVecS1, fxVecS1, fpVecS1, fsVecS1;
REAL8 fLVecS2, fxVecS2, fpVecS2, fsVecS2;
REAL8 w2;
w2 = r*vars->pT2;
fLVecS1 = (2*x11*y*(16 - 48*w2 + 2*eta*(-38 + 84*w2 - 135*y2) + dm*(-16 + 48*w2 - 75*y2) + 75*y2))/(15);
fxVecS1 = (x12*(4*eta*pS2*(-8 + 9*w2 - 45*y2) - 30*eta*x3*(xS1 + xS2)*y*(2 - 3*w2 + 7*y2) - 15*(-1 + dm)*x3*xS1*y*(2 - 21*w2 + 49*y2) - 2*pS1*(4 - 27*w2 + dm*(-4 + 27*w2 - 135*y2) + 135*y2 + eta*(16 - 18*w2 + 90*y2))))/(15.);
fpVecS1 = (x9*(-60*(-1 + dm)*pS1*y + x3*(2*eta*xS2*(8 - 9*w2 + 45*y2) + xS1*(-16 + 33*w2 - 165*y2 + eta*(16 - 18*w2 + 90*y2) + dm*(16 - 33*w2 + 165*y2)))))/(5.);
fsVecS1 = (-4*eta*x11*y*(11 - 18*w2 + 30*y2))/(5.);
fLVecS2 = (2*x11*y*(16 - 48*w2 + 2*eta*(-38 + 84*w2 - 135*y2) + 75*y2 + dm*(16 - 48*w2 + 75*y2)))/(15.);
fxVecS2 = (x12*(2*eta*(2*pS1*(-8 + 9*w2 - 45*y2) + 2*pS2*(-8 + 9*w2 - 45*y2) - 15*x3*(xS1 + xS2)*y*(2 - 3*w2 + 7*y2)) + (1 + dm)*(pS2*(-8 + 54*w2 - 270*y2) + 15*x3*xS2*y*(2 - 21*w2 + 49*y2))))/(15.);
fpVecS2 = (x9*((1 + dm)*(60*pS2*y + x3*xS2*(-16 + 33*w2 - 165*y2)) - 2*eta*x3*(xS1 + xS2)*(-8 + 9*w2 - 45*y2)))/(5.);
fsVecS2 = (-4*eta*x11*y*(11 - 18*w2 + 30*y2))/(5.);
// results
REAL8 tmpfRRs1[3], tmpfRRs2[3];
for (i=0; i<3; i++)
{
fRRVec[i] = (fxVec0 + fxVec2 + fxVec3 + fxVec4) * vars->xVec[i] + (fpVec0 + fpVec2 + fpVec3 + fpVec4) * vars->pTVec[i] + (fLVec3 + fLVec4) * vars->LVec[i];
tmpfRRs1[i] = fLVecS1 * vars->LVec[i] + fxVecS1 * vars->xVec[i] + fpVecS1 * vars->pTVec[i] + fsVecS1 * vars->s2Vec[i];
tmpfRRs2[i] = fLVecS2 * vars->LVec[i] + fxVecS2 * vars->xVec[i] + fpVecS2 * vars->pTVec[i] + fsVecS2 * vars->s1Vec[i];
}
cross_product3d(tmpfRRs1, vars->s1Vec, fRRs1Vec);
cross_product3d(tmpfRRs2, vars->s2Vec, fRRs2Vec);
// print_debug("fRRs1Vec = (%.5e, %.5e, %.5e), fRRs2Vec = (%.5e, %.5e, %.5e)\n",
// fRRs1Vec[0], fRRs1Vec[1], fRRs1Vec[2],
// fRRs2Vec[0], fRRs2Vec[1], fRRs2Vec[2]);
return;
}