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
180
181
182
183
184
185
|
// TwoPunctures: File "Equations.c"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <ctype.h>
#include "cctk_Parameters.h"
#include "TP_utilities.h"
#include "TwoPunctures.h"
// U.d0[ivar] = U[ivar]; (ivar = 0..nvar-1)
// U.d1[ivar] = U[ivar]_x;
// U.d2[ivar] = U[ivar]_y;
// U.d3[ivar] = U[ivar]_z;
// U.d11[ivar] = U[ivar]_xx;
// U.d12[ivar] = U[ivar]_xy;
// U.d13[ivar] = U[ivar]_xz;
// U.d22[ivar] = U[ivar]_yy;
// U.d23[ivar] = U[ivar]_yz;
// U.d33[ivar] = U[ivar]_zz;
double
BY_KKofxyz (double x, double y, double z)
{
DECLARE_CCTK_PARAMETERS;
int i, j;
double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
Aij, AijAij, n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
r_plus = sqrt (r2_plus);
r_minus = sqrt (r2_minus);
r3_plus = r_plus * r2_plus;
r3_minus = r_minus * r2_minus;
n_plus[0] = (x - par_b) / r_plus;
n_minus[0] = (x + par_b) / r_minus;
n_plus[1] = y / r_plus;
n_minus[1] = y / r_minus;
n_plus[2] = z / r_plus;
n_minus[2] = z / r_minus;
// dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-)
np_Pp = 0;
nm_Pm = 0;
for (i = 0; i < 3; i++)
{
np_Pp += n_plus[i] * par_P_plus[i];
nm_Pm += n_minus[i] * par_P_minus[i];
}
// cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i
np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
AijAij = 0;
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{ // Bowen-York-Curvature :
Aij =
+ 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
+ np_Pp * n_plus[i] * n_plus[j]) / r2_plus
+ 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
+ nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
if (i == j)
Aij -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
AijAij += Aij * Aij;
}
}
return AijAij;
}
void
BY_Aijofxyz (double x, double y, double z, double Aij[3][3])
{
DECLARE_CCTK_PARAMETERS;
int i, j;
double r_plus, r2_plus, r3_plus, r_minus, r2_minus, r3_minus, np_Pp, nm_Pm,
n_plus[3], n_minus[3], np_Sp[3], nm_Sm[3];
r2_plus = (x - par_b) * (x - par_b) + y * y + z * z;
r2_minus = (x + par_b) * (x + par_b) + y * y + z * z;
r_plus = sqrt (r2_plus);
r_minus = sqrt (r2_minus);
r3_plus = r_plus * r2_plus;
r3_minus = r_minus * r2_minus;
n_plus[0] = (x - par_b) / r_plus;
n_minus[0] = (x + par_b) / r_minus;
n_plus[1] = y / r_plus;
n_minus[1] = y / r_minus;
n_plus[2] = z / r_plus;
n_minus[2] = z / r_minus;
// dot product: np_Pp = (n_+).(P_+); nm_Pm = (n_-).(P_-)
np_Pp = 0;
nm_Pm = 0;
for (i = 0; i < 3; i++)
{
np_Pp += n_plus[i] * par_P_plus[i];
nm_Pm += n_minus[i] * par_P_minus[i];
}
// cross product: np_Sp[i] = [(n_+) x (S_+)]_i; nm_Sm[i] = [(n_-) x (S_-)]_i
np_Sp[0] = n_plus[1] * par_S_plus[2] - n_plus[2] * par_S_plus[1];
np_Sp[1] = n_plus[2] * par_S_plus[0] - n_plus[0] * par_S_plus[2];
np_Sp[2] = n_plus[0] * par_S_plus[1] - n_plus[1] * par_S_plus[0];
nm_Sm[0] = n_minus[1] * par_S_minus[2] - n_minus[2] * par_S_minus[1];
nm_Sm[1] = n_minus[2] * par_S_minus[0] - n_minus[0] * par_S_minus[2];
nm_Sm[2] = n_minus[0] * par_S_minus[1] - n_minus[1] * par_S_minus[0];
for (i = 0; i < 3; i++)
{
for (j = 0; j < 3; j++)
{ // Bowen-York-Curvature :
Aij[i][j] =
+ 1.5 * (par_P_plus[i] * n_plus[j] + par_P_plus[j] * n_plus[i]
+ np_Pp * n_plus[i] * n_plus[j]) / r2_plus
+ 1.5 * (par_P_minus[i] * n_minus[j] + par_P_minus[j] * n_minus[i]
+ nm_Pm * n_minus[i] * n_minus[j]) / r2_minus
- 3.0 * (np_Sp[i] * n_plus[j] + np_Sp[j] * n_plus[i]) / r3_plus
- 3.0 * (nm_Sm[i] * n_minus[j] + nm_Sm[j] * n_minus[i]) / r3_minus;
if (i == j)
Aij[i][j] -= +1.5 * (np_Pp / r2_plus + nm_Pm / r2_minus);
}
}
}
//---------------------------------------------------------------
//******* Nonlinear Equations **********
//---------------------------------------------------------------
void
NonLinEquations (double A, double B, double X, double R,
double x, double r, double phi,
double y, double z, derivs U, double *values)
{
DECLARE_CCTK_PARAMETERS;
double r_plus, r_minus, psi, psi2, psi4, psi7;
double mu;
r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
psi =
1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
psi2 = psi * psi;
psi4 = psi2 * psi2;
psi7 = psi * psi2 * psi4;
values[0] =
U.d11[0] + U.d22[0] + U.d33[0] + 0.125 * BY_KKofxyz (x, y, z) / psi7;
}
//---------------------------------------------------------------
//******* Linear Equations **********
//---------------------------------------------------------------
void
LinEquations (double A, double B, double X, double R,
double x, double r, double phi,
double y, double z, derivs dU, derivs U, double *values)
{
DECLARE_CCTK_PARAMETERS;
double r_plus, r_minus, psi, psi2, psi4, psi8;
r_plus = sqrt ((x - par_b) * (x - par_b) + y * y + z * z);
r_minus = sqrt ((x + par_b) * (x + par_b) + y * y + z * z);
psi =
1. + 0.5 * par_m_plus / r_plus + 0.5 * par_m_minus / r_minus + U.d0[0];
psi2 = psi * psi;
psi4 = psi2 * psi2;
psi8 = psi4 * psi4;
values[0] = dU.d11[0] + dU.d22[0] + dU.d33[0]
- 0.875 * BY_KKofxyz (x, y, z) / psi8 * dU.d0[0];
}
//---------------------------------------------------------------
|
