diff options
Diffstat (limited to 'src/gr/LHS.maple')
| -rw-r--r-- | src/gr/LHS.maple | 131 |
1 files changed, 131 insertions, 0 deletions
diff --git a/src/gr/LHS.maple b/src/gr/LHS.maple new file mode 100644 index 0000000..2bd34be --- /dev/null +++ b/src/gr/LHS.maple @@ -0,0 +1,131 @@ +# H.maple - evaluate $H$ = LHS of horizon function +# $Id$ +# +# 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 "../maple/coeffs.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 "../maple/coeffs.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 "../maple/coeffs.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; |