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
|
# horizon.maple - evaluate LHS of horizon function
#
# non_unit_normal - compute s_d
# non_unit_normal_gradient - compute partial_d_s_d
# horizon_function - compute HA,HB,HC,HD,H = fn(s_d)
#
################################################################################
################################################################################
################################################################################
#
# This function computes the non-unit outward-pointing normal vector
# (field) to the horizon, using the equation on p.7.1 of my apparent
# horizon finding notes.
#
non_unit_normal :=
proc()
global
@include "../maple/coords.minc",
@include "../gr/gr_gfas.minc";
local i, u;
assert_fnd_exists(h);
assert_fnd_exists(X_ud);
assert_fnd_exists(s_d, fnd);
for i from 1 to N
do
s_d__fnd[i] := x_xyz[i]/r
- sum('X_ud[u,i]*Diff(h,x_rrs[u])', 'u'=A_min..N);
end do;
NULL;
end proc;
################################################################################
#
# This function computes the gradient of the non-unit outward-pointing
# normal vector (field) to the horizon, using the equation on p.7.2 of
# my apparent horizon finding notes.
#
non_unit_normal_gradient :=
proc()
global
@include "../maple/coords.minc",
@include "../gr/gr_gfas.minc";
local temp, i, j, k, u, v;
assert_fnd_exists(h);
assert_fnd_exists(X_ud);
assert_fnd_exists(X_udd);
assert_fnd_exists(partial_d_s_d, fnd);
for i from 1 to N
do
for j from 1 to N
do
temp := `if`(i=j, sum('x_xyz[k]^2','k'=1..3), 0);
# simplify likes to put things over a common denominator,
# so we only apply it to the individual terms, not to
# the whole experssion
partial_d_s_d__fnd[i,j]
:= + simplify( (temp-x_xyz[i]*x_xyz[j])/r^3 )
- simplify( sum('X_udd[u,i,j]*Diff(h,x_rrs[u])', 'u'=A_min..N) )
- simplify( msum('X_ud[u,i]*X_ud[v,j]*Diff(h,x_rrs[u],x_rrs[v])',
'u'=A_min..N, 'v'=A_min..N) );
end do;
end do;
NULL;
end proc;
################################################################################
#
# This function computes the LHS of the apparent horizon equation
# as a function of 1st and 2nd angular partial derivatives of the
# apparent horizon radius $h$, in the Maple gridfn array H.fnd ,
# using equation (15) in my 1996 apparent horizon finding paper.
#
horizon_function :=
proc()
global
@include "../maple/coords.minc",
@include "../gr/gr_gfas.minc";
local i,j,k,l;
assert_fnd_exists(g_uu);
assert_fnd_exists(K_uu);
assert_fnd_exists(partial_d_ln_sqrt_g);
assert_fnd_exists(s_d, fnd);
assert_fnd_exists(partial_d_s_d, fnd);
#
# n.b. no simplify() in these subexpressions; empirically it makes
# them much *more* complicated (measured by codegen(...,optimized))
#
# (15a) in my paper
HA__fnd := - msum('g_uu[i,k]*s_d__fnd[k] * g_uu[j,l]*s_d__fnd[l]
* partial_d_s_d__fnd[i,j]',
'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N)
- (1/2)*msum('g_uu[i,j]*s_d__fnd[j]
* Diff(g_uu[k,l],x_xyz[i]) * s_d__fnd[k]*s_d__fnd[l]',
'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N);
# (15b) in my paper
HB__fnd := + msum('Diff(g_uu[i,j],x_xyz[i])*s_d__fnd[j]', 'i'=1..N, 'j'=1..N)
+ msum('g_uu[i,j]*partial_d_s_d__fnd[i,j]', 'i'=1..N, 'j'=1..N)
+ msum('g_uu[i,j]*partial_d_ln_sqrt_g[i]*s_d__fnd[j]',
'i'=1..N, 'j'=1..N);
# (15c) in my paper
HC__fnd := msum('K_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
# (15d) in my paper
HD__fnd := msum('g_uu[i,j]*s_d__fnd[i]*s_d__fnd[j]', 'i'=1..N, 'j'=1..N);
# n.b. no simplify() here, it would try to put things over a common
# denominator, which would make the equation much more complicated
H__fnd := HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K;
NULL;
end proc;
|
