aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorjthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5>2002-11-21 20:39:51 +0000
committerjthorn <jthorn@f88db872-0e4f-0410-b76b-b9085cfa78c5>2002-11-21 20:39:51 +0000
commit03bd2968bc2cafe03a4d76c3609526bd23896a0f (patch)
treee53cca7d80221f9157510fdacf982cd4fb504273
parent47f5db5fc54dfed7c03e6d5264a25ee02e382347 (diff)
change in terminology/notation:
LHS of apparent horizon equation was formerly called $H$, changed to $\Theta$ because it is in fact precisely the expansion of the surface r = h(angle). This implies renaming a lot of variables & functions & a few parameters (which aren't specified in most par files outside my own for testing this thorn) git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@901 f88db872-0e4f-0410-b76b-b9085cfa78c5
Diffstat
-rw-r--r--param.ccl42
-rw-r--r--src/driver/Newton.cc75
-rw-r--r--src/driver/driver.hh14
-rw-r--r--src/driver/find_horizons.cc98
-rw-r--r--src/driver/io.cc4
-rw-r--r--src/driver/setup.cc24
-rw-r--r--src/gr.cg/expansion.c182
-rw-r--r--src/gr.cg/expansion_Jacobian.c628
-rw-r--r--src/gr.cg/horizon_Jacobian.c625
-rw-r--r--src/gr.cg/horizon_function.c182
-rw-r--r--src/gr/cg.hh36
-rw-r--r--src/gr/expansion.cc (renamed from src/gr/horizon_function.cc)121
-rw-r--r--src/gr/expansion_Jacobian.cc (renamed from src/gr/horizon_Jacobian.cc)209
-rw-r--r--src/gr/gfns.hh14
-rw-r--r--src/gr/gr.hh30
-rw-r--r--src/gr/gr_gfas.minc42
-rw-r--r--src/gr/horizon.maple195
-rw-r--r--src/gr/make.code.defn4
-rw-r--r--src/gr/maple.log529
-rw-r--r--src/gr/setup_gr_gfas.maple55
20 files changed, 1575 insertions, 1534 deletions
diff --git a/param.ccl b/param.ccl
index eb6c87c..9b8ebe5 100644
--- a/param.ccl
+++ b/param.ccl
@@ -50,11 +50,13 @@ boolean find_AHs_at_poststep \
keyword method "what should this thorn do for each apparent horizon?"
{
# these options are mostly for testing/debugging
-"horizon function" :: "evaluate the LHS function H(h)"
-"test Jacobian" :: "compute/print the J[H(h)] Jacobian matrix"
+"evaluate expansion" :: "evaluate the LHS function Theta(h)"
+"test expansion Jacobian" :: "compute/print the J[Theta(h)] Jacobian matrix \
+ (possibly in several ways, depending on \
+ the test_all_Jacobian_methods parameter"
# this is for normal apparent horizon finding
-"find horizon" :: "find the apparent horizon"
+"find horizon" :: "find the apparent horizon"
} "find horizon"
#
@@ -92,7 +94,7 @@ keyword verbose_level \
# 1 line for each horizon giving position/mass/area, + a summary line or two
"physics details" :: "more detailed physics messages"
-# 1 line giving H(h) norms at each Newton iteration
+# 1 line giving Theta(h) norms at each Newton iteration
"algorithm highlights" :: \
"physics details + a few messages about the AH-finding algorithm"
@@ -169,7 +171,7 @@ boolean test_all_Jacobian_methods \
################################################################################
#
-# ***** parameters for the Newton's-method solution of H(h) = 0 *****
+# ***** parameters for the Newton's-method solution of Theta(h) = 0 *****
#
#
@@ -207,7 +209,7 @@ real max_Delta_h_over_h \
#
# we declare convergence if *either* of the following two criteria are met
#
-real H_norm_for_convergence "declare convergence if ||H||_inf <= this"
+real Theta_norm_for_convergence "declare convergence if ||Theta||_inf <= this"
{
(0.0:* :: "any positive real number"
} 1.0e-8
@@ -222,9 +224,9 @@ real Delta_h_norm_for_convergence \
# etc. On the other hand, setting it to true probably only slows down
# the apparent horizon finder by a few percent.
#
-boolean final_H_update_if_exit_x_H_small \
- "should we do a final H(h) update after a h += Delta_h update which is\
- so small it meets the Delta_h_norm_for_convergence convergence criterion?"
+boolean final_Theta_update_if_Delta_h_converged \
+ "should we do a final Theta(h) update if we terminate the \
+ Newton iteration by the small-||Delta_h|| convergence criterion?"
{
} "false"
@@ -254,15 +256,15 @@ int how_often_to_output_h \
# setting this > 0 is probably only of interest if the Newton iteration
# fails to converge, or if you're debugging AHFinderDirect internals
-int how_often_to_output_H_of_h \
- "how often (in Cactus time steps) should we output the H(h) functions?"
+int how_often_to_output_Theta \
+ "how often (in Cactus time steps) should we output the Theta(h) functions?"
{
-0 :: "don't output H(h) at all"
+0 :: "don't output Theta(h) at all"
1:* :: "any integer >= 1"
} 0
keyword horizon_file_format \
- "what file format should we use for h and H(h) data files?"
+ "what file format should we use for h and Theta(h) data files?"
{
"ASCII (gnuplot)" :: "simple ASCII format, directly readable by gnuplot"
"HDF5" :: "HDF5 surface format (alas not implemented yet)"
@@ -303,10 +305,10 @@ string h_base_file_name \
.+ :: "any nonempty string"
} "h"
-string H_of_h_base_file_name "base file name for H(h) output file(s)"
+string Theta_base_file_name "base file name for Theta(h) output file(s)"
{
.+ :: "any nonempty string"
-} "H"
+} "Theta"
string Delta_h_base_file_name \
"base file name for horizon-shape-update Delta_h output file(s)"
@@ -368,10 +370,10 @@ boolean output_initial_guess \
} "false"
# for debugging convergence failures, we can optionally output
-# h, H, and delta_h at each Newton iteration
+# h, Theta, and delta_h at each Newton iteration
# (the file names are the usual ones with ".it%d" appended)
boolean debugging_output_at_each_Newton_iteration \
- "should we output {h, H, delta_h} at each Newton iteration?"
+ "should we output {h, Theta, delta_h} at each Newton iteration?"
{
} "false"
@@ -419,7 +421,7 @@ real origin_z[5] "global z coordinate of patch system origin"
} 0.0
#
-# The "(rotating)" patch system types are ok for evaluating H(h),
+# The "(rotating)" patch system types are ok for evaluating Theta(h),
# but don't work yet for apparent horizon finding
# (the Jacobian computation doesn't yet grok the nonlocal rotation BCs).
#
@@ -567,7 +569,7 @@ keyword geometry_method "how do we compute the slice's geometry?"
# - It must support taking at least 1st derivatives as part of the
# interpolation.
# - It should give at least $C^1$ interpolants for smooth data, otherwise
-# the H(h) function will have "spikes" and the Newton iteration may
+# the Theta(h) function will have "spikes" and the Newton iteration may
# fail to converge all the way down to tight error tolerances. $C^2$
# would be even better, but in practice a ($C^1$) Hermite interpolant
# works well.
@@ -645,7 +647,7 @@ real geometry__Schwarzschild_EF__Delta_xyz \
#
# These tests control whether we check that various angular gridfns
# are finite (neither NaN nor infinity) at various points in evaluating
-# the H(h) function. These are pretty cheap tests, and they're quite
+# the Theta(h) function. These are pretty cheap tests, and they're quite
# useful in catching assorted wierdness, so it's probably worth leaving
# them enabled unless you're trying to squeeze every last nanosecond...
#
diff --git a/src/driver/Newton.cc b/src/driver/Newton.cc
index 8bb92d5..f13a40e 100644
--- a/src/driver/Newton.cc
+++ b/src/driver/Newton.cc
@@ -1,7 +1,7 @@
-// Newton.cc -- solve H(h) = 0 via Newton's method
+// Newton.cc -- solve Theta(h) = 0 via Newton's method
// $Header$
//
-// Newton_solve - driver to solve H(h) = 0 via Newton's method
+// Newton_solve - driver to solve Theta(h) = 0 via Newton's method
//
#include <stdio.h>
@@ -41,7 +41,7 @@ using jtutil::error_exit;
//******************************************************************************
//
-// This function solves H(h) = 0 for h via Newton's method.
+// This function solves Theta(h) = 0 for h via Newton's method.
//
// Arguments:
// initial_find_flag = true ==> This is the first attempt to find this
@@ -54,7 +54,7 @@ using jtutil::error_exit;
//
// Results:
// This function returns a success flag: this is true if (and only if)
-// it found an h satisfying H(h) = 0 to within the error tolerances,
+// it found an h satisfying Theta(h) = 0 to within the error tolerances,
// otherwise it's false.
//
bool Newton_solve(patch_system& ps,
@@ -77,7 +77,7 @@ const int max_Newton_iterations
{
if (verbose_info.print_algorithm_details)
then CCTK_VInfo(CCTK_THORNSTRING,
- "Newton iteration %d of %d",
+ "Newton iteration %d/%d",
iteration, max_Newton_iterations);
if (solver_info.debugging_output_at_each_Newton_iteration)
then output_gridfn(ps, gfns::gfn__h,
@@ -86,50 +86,52 @@ const int max_Newton_iterations
iteration);
//
- // evaluate H(h)
+ // evaluate the expansion Theta(h)
//
- jtutil::norm<fp> H_norms;
- if (! horizon_function(ps,
- cgi, gi,
- true, // yes, we want the Jacobian coeffs
- verbose_info.print_algorithm_details,
- &H_norms))
+ jtutil::norm<fp> Theta_norms;
+ if (! expansion(ps,
+ cgi, gi,
+ true, // yes, we want the Jacobian coeffs
+ verbose_info.print_algorithm_details,
+ &Theta_norms))
then return false; // *** ERROR RETURN ***
if (verbose_info.print_algorithm_highlights)
then CCTK_VInfo(CCTK_THORNSTRING,
- " iteration %d: H ||rms||=%.1e, ||infinity||=%.1e",
+ " iteration %d: Theta ||rms||=%.1e, ||infty||=%.1e",
iteration,
- H_norms.rms_norm(), H_norms.infinity_norm());
+ Theta_norms.rms_norm(), Theta_norms.infinity_norm());
if (solver_info.debugging_output_at_each_Newton_iteration)
- then output_gridfn(ps, gfns::gfn__H,
- IO_info, IO_info.H_base_file_name,
+ then output_gridfn(ps, gfns::gfn__Theta,
+ IO_info, IO_info.Theta_base_file_name,
hn, verbose_info.print_algorithm_details,
iteration);
//
- // convergence test on ||H||
+ // convergence test on ||Theta||
//
- if (H_norms.infinity_norm() <= solver_info.H_norm_for_convergence)
+ if (Theta_norms.infinity_norm() <= solver_info
+ .Theta_norm_for_convergence)
then {
if (verbose_info.print_algorithm_details)
then CCTK_VInfo(CCTK_THORNSTRING,
- "==> finished (||H|| < %g)",
- double(solver_info.H_norm_for_convergence));
+ "==> finished (||Theta|| < %g)",
+ double(solver_info
+ .Theta_norm_for_convergence));
return true; // *** NORMAL RETURN ***
}
//
// compute the Newton step
//
- if (! horizon_Jacobian(ps, Jac,
- cgi, gi, Jacobian_info,
- verbose_info.print_algorithm_details))
+ if (! expansion_Jacobian(ps, Jac,
+ cgi, gi, Jacobian_info,
+ verbose_info.print_algorithm_details))
then return false; // *** ERROR RETURN ***
if (verbose_info.print_algorithm_details)
then CCTK_VInfo(CCTK_THORNSTRING,
"solving linear system for Delta_h (%d unknowns)",
Jac.NN());
- const fp rcond = Jac.solve_linear_system(gfns::gfn__H,
+ const fp rcond = Jac.solve_linear_system(gfns::gfn__Theta,
gfns::gfn__Delta_h);
if (rcond == 0.0)
then {
@@ -211,31 +213,30 @@ const int max_Newton_iterations
"==> finished (||Delta_h|| < %g)",
double(solver_info
.Delta_h_norm_for_convergence));
- if (solver_info.final_H_update_if_exit_x_H_small)
+ if (solver_info.final_Theta_update_if_Delta_h_converged)
then {
if (verbose_info.print_algorithm_details)
then CCTK_VInfo(CCTK_THORNSTRING,
- " doing final H(h) evaluation");
- if (! horizon_function(ps,
- cgi, gi,
- false, // no Jacobian coeffs
- verbose_info
- .print_algorithm_details))
+ " doing final Theta(h) evaluation");
+ if (! expansion(ps,
+ cgi, gi,
+ false, // no Jacobian coeffs
+ verbose_info.print_algorithm_details))
then return false; // *** ERROR RETURN ***
}
return true; // *** NORMAL RETURN ***
}
}
-if (solver_info.final_H_update_if_exit_x_H_small)
+if (solver_info.final_Theta_update_if_Delta_h_converged)
then {
if (verbose_info.print_algorithm_details)
then CCTK_VInfo(CCTK_THORNSTRING,
- " doing final H(h) evaluation");
- if (! horizon_function(ps,
- cgi, gi,
- false, // no Jacobian coeffs
- verbose_info.print_algorithm_details))
+ " doing final Theta(h) evaluation");
+ if (! expansion(ps,
+ cgi, gi,
+ false, // no Jacobian coeffs
+ verbose_info.print_algorithm_details))
then return false; // *** ERROR RETURN ***
}
return false; // *** ERROR RETURN ***
diff --git a/src/driver/driver.hh b/src/driver/driver.hh
index d152ae3..87afdf5 100644
--- a/src/driver/driver.hh
+++ b/src/driver/driver.hh
@@ -14,8 +14,8 @@
//
enum method
{
- method__horizon_function,
- method__test_Jacobian,
+ method__evaluate_expansion,
+ method__test_expansion_Jacobian,
method__find_horizon // no comma
};
@@ -81,7 +81,7 @@ struct initial_guess_info
};
//
-// This struct holds parameters for solving the H(h) = 0 equations.
+// This struct holds parameters for solving the Theta(h) = 0 equations.
//
struct solver_info
{
@@ -89,9 +89,9 @@ struct solver_info
int max_Newton_iterations__initial,
max_Newton_iterations__subsequent;
fp max_Delta_h_over_h;
- fp H_norm_for_convergence;
+ fp Theta_norm_for_convergence;
fp Delta_h_norm_for_convergence;
- bool final_H_update_if_exit_x_H_small;
+ bool final_Theta_update_if_Delta_h_converged;
};
//
@@ -123,12 +123,12 @@ struct IO_info
enum horizon_file_format horizon_file_format;
bool output_initial_guess;
int how_often_to_output_h,
- how_often_to_output_H;
+ how_often_to_output_Theta;
bool output_ghost_zones_for_h;
const char* ASCII_gnuplot_file_name_extension;
const char* HDF5_file_name_extension;
const char* h_base_file_name;
- const char* H_base_file_name;
+ const char* Theta_base_file_name;
const char* Delta_h_base_file_name;
const char* Jacobian_base_file_name;
diff --git a/src/driver/find_horizons.cc b/src/driver/find_horizons.cc
index 35b5085..c9ef0b5 100644
--- a/src/driver/find_horizons.cc
+++ b/src/driver/find_horizons.cc
@@ -9,8 +9,8 @@
/// Cactus_gridfn_data_ptr - get a single data pointer from a variable index
///
/// find_horizon - find a horizon
-/// do_horizon_function
-/// do_test_Jacobian
+/// do_evaluate_expansion
+/// do_test_expansion_Jacobian
/// do_find_horizon
///
/// compute_BH_diagnostics - compute BH diagnostics for a horizon
@@ -70,13 +70,13 @@ const CCTK_REAL* Cactus_gridfn_data_ptr(const cGH *GH, int varindex,
const char gridfn_name[],
bool check_for_NULL = true);
-bool do_horizon_function
+bool do_evaluate_expansion
(const struct verbose_info& verbose_info, int timer_handle,
- struct IO_info& IO_info, bool output_h, bool output_H,
+ struct IO_info& IO_info, bool output_h, bool output_Theta,
struct cactus_grid_info& cgi, struct geometry_info& gi,
patch_system& ps,
int hn, int N_horizons);
-bool do_test_Jacobian
+bool do_test_expansion_Jacobian
(const struct verbose_info& verbose_info, int timer_handle,
struct IO_info& IO_info,
bool test_all_Jacobian_methods,
@@ -86,7 +86,7 @@ bool do_test_Jacobian
int hn, int N_horizons);
bool do_find_horizon
(const struct verbose_info& verbose_info, int timer_handle,
- struct IO_info& IO_info, bool output_h, bool output_H,
+ struct IO_info& IO_info, bool output_h, bool output_Theta,
const struct solver_info& solver_info, bool initial_find_flag,
const struct Jacobian_info& Jac_info,
struct cactus_grid_info& cgi, struct geometry_info& gi,
@@ -141,9 +141,9 @@ IO_info.time = cctk_time;
const bool output_h
= (IO_info.how_often_to_output_h > 0)
&& ((IO_info.time_iteration % IO_info.how_often_to_output_h) == 0);
-const bool output_H
- = (IO_info.how_often_to_output_H > 0)
- && ((IO_info.time_iteration % IO_info.how_often_to_output_H) == 0);
+const bool output_Theta
+ = (IO_info.how_often_to_output_Theta > 0)
+ && ((IO_info.time_iteration % IO_info.how_often_to_output_Theta) == 0);
//
// we need to re-fetch the Cactus data pointers at least each time step,
@@ -183,27 +183,27 @@ setup_Cactus_gridfn_data_ptrs(cctkGH, state.cgi);
AH_info.AH_found = false;
switch (state.method)
{
- case method__horizon_function:
- do_horizon_function(verbose_info, state.timer_handle,
- IO_info, output_h, output_H,
- state.cgi, state.gi,
- ps,
- hn, state.N_horizons);
+ case method__evaluate_expansion:
+ do_evaluate_expansion(verbose_info, state.timer_handle,
+ IO_info, output_h, output_Theta,
+ state.cgi, state.gi,
+ ps,
+ hn, state.N_horizons);
break;
- case method__test_Jacobian:
- do_test_Jacobian(verbose_info, state.timer_handle,
- IO_info,
- test_all_Jacobian_methods,
- state.Jac_info, state.cgi, state.gi,
- ps, AH_info.Jac_ptr,
- hn, state.N_horizons);
+ case method__test_expansion_Jacobian:
+ do_test_expansion_Jacobian(verbose_info, state.timer_handle,
+ IO_info,
+ test_all_Jacobian_methods,
+ state.Jac_info, state.cgi, state.gi,
+ ps, AH_info.Jac_ptr,
+ hn, state.N_horizons);
break;
case method__find_horizon:
AH_info.AH_found
= do_find_horizon(verbose_info, state.timer_handle,
- IO_info, output_h, output_H,
+ IO_info, output_h, output_Theta,
solver_info, initial_find_flag,
state.Jac_info, state.cgi, state.gi,
ps, AH_info.Jac_ptr,
@@ -356,9 +356,9 @@ return data_ptr;
//
// This function implements AHFinderDirect::method == "horizon function",
-// that is, it evaluates the H(h) function for a single apparent horizon.
+// that is, it evaluates the Theta(h) function for a single apparent horizon.
//
-// Note that if we decide to output h, we output it *after* any H(h)
+// Note that if we decide to output h, we output it *after* any Theta(h)
// evaluation or horizon finding has been done, to ensure that all the
// ghost zones are filled in in case we need to print them.
//
@@ -376,35 +376,35 @@ return data_ptr;
// otherwise.
//
namespace {
-bool do_horizon_function
+bool do_evaluate_expansion
(const struct verbose_info& verbose_info, int timer_handle,
- struct IO_info& IO_info, bool output_h, bool output_H,
+ struct IO_info& IO_info, bool output_h, bool output_Theta,
struct cactus_grid_info& cgi, struct geometry_info& gi,
patch_system& ps,
int hn, int N_horizons)
{
-jtutil::norm<fp> H_norms;
+jtutil::norm<fp> Theta_norms;
if (timer_handle >= 0)
then CCTK_TimerStartI(timer_handle);
-const bool status = horizon_function(ps, cgi, gi, false, true, &H_norms);
+const bool status = expansion(ps, cgi, gi, false, true, &Theta_norms);
if (timer_handle >= 0)
then CCTK_TimerStopI(timer_handle);
if (!status)
then return false; // *** ERROR RETURN ***
-if (H_norms.is_nonempty()) // might be empty if H(h) eval failed
+if (Theta_norms.is_nonempty()) // might be empty if Theta(h) eval failed
then CCTK_VInfo(CCTK_THORNSTRING,
- " H(h) rms-norm %.2e, infinity-norm %.2e",
- H_norms.rms_norm(), H_norms.infinity_norm());
+ " Theta(h) rms-norm %.2e, infinity-norm %.2e",
+ Theta_norms.rms_norm(), Theta_norms.infinity_norm());
if (output_h)
then output_gridfn(ps, gfns::gfn__h,
IO_info, IO_info.h_base_file_name,
hn, verbose_info.print_algorithm_details);
-if (output_H)
- then output_gridfn(ps, gfns::gfn__H,
- IO_info, IO_info.H_base_file_name,
+if (output_Theta)
+ then output_gridfn(ps, gfns::gfn__Theta,
+ IO_info, IO_info.Theta_base_file_name,
hn, verbose_info.print_algorithm_details);
return true; // *** NORMAL RETURN ***
@@ -415,7 +415,7 @@ return true; // *** NORMAL RETURN ***
//
// This function implements AHFinderDirect::method == "test Jacobian",
-// that is, it computes and prints the Jacobian matrix J[H(h)] for a
+// that is, it computes and prints the Jacobian matrix J[Theta(h)] for a
// single apparent horizon. The computation may optionally be done
// in several different ways, in which case all the resulting Jacobian
// matrices are printed, as are their differences. Alternatively, only
@@ -442,7 +442,7 @@ return true; // *** NORMAL RETURN ***
// false otherwise.
//
namespace {
-bool do_test_Jacobian
+bool do_test_expansion_Jacobian
(const struct verbose_info& verbose_info, int timer_handle,
struct IO_info& IO_info,
bool test_all_Jacobian_methods,
@@ -453,12 +453,12 @@ bool do_test_Jacobian
{
// numerical perturbation
Jacobian* Jac_NP_ptr = Jac_ptr;
-if (! horizon_function(ps, cgi, gi, true))
+if (! expansion(ps, cgi, gi, true))
then return false; // *** ERROR RETURN ***
Jac_info.Jacobian_method = Jacobian_method__numerical_perturb;
-if (! horizon_Jacobian(ps, *Jac_NP_ptr,
- cgi, gi, Jac_info,
- true))
+if (! expansion_Jacobian(ps, *Jac_NP_ptr,
+ cgi, gi, Jac_info,
+ true)) // print msgs
then return false; // *** ERROR RETURN ***
Jacobian* Jac_SD_FDdr_ptr = NULL;
@@ -466,12 +466,12 @@ if (test_all_Jacobian_methods)
then {
// symbolic differentiation with finite diff d/dr
Jac_SD_FDdr_ptr = & new_Jacobian(ps, Jac_info.Jacobian_storage_method);
- if (! horizon_function(ps, cgi, gi, true))
+ if (! expansion(ps, cgi, gi, true))
then return false; // *** ERROR RETURN ***
Jac_info.Jacobian_method = Jacobian_method__symbolic_diff_with_FD_dr;
- if (! horizon_Jacobian(ps, *Jac_SD_FDdr_ptr,
- cgi, gi, Jac_info,
- true))
+ if (! expansion_Jacobian(ps, *Jac_SD_FDdr_ptr,
+ cgi, gi, Jac_info,
+ true)) // print msgs
then return false; // *** ERROR RETURN ***
}
@@ -515,7 +515,7 @@ return true; // *** NORMAL RETURN ***
namespace {
bool do_find_horizon
(const struct verbose_info& verbose_info, int timer_handle,
- struct IO_info& IO_info, bool output_h, bool output_H,
+ struct IO_info& IO_info, bool output_h, bool output_Theta,
const struct solver_info& solver_info, bool initial_find_flag,
const struct Jacobian_info& Jac_info,
struct cactus_grid_info& cgi, struct geometry_info& gi,
@@ -545,9 +545,9 @@ if (output_h)
then output_gridfn(ps, gfns::gfn__h,
IO_info, IO_info.h_base_file_name,
hn, verbose_info.print_algorithm_details);
-if (output_H)
- then output_gridfn(ps, gfns::gfn__H,
- IO_info, IO_info.H_base_file_name,
+if (output_Theta)
+ then output_gridfn(ps, gfns::gfn__Theta,
+ IO_info, IO_info.Theta_base_file_name,
hn, verbose_info.print_algorithm_details);
return true; // *** NORMAL RETURN ***
diff --git a/src/driver/io.cc b/src/driver/io.cc
index e4c0ee7..a557257 100644
--- a/src/driver/io.cc
+++ b/src/driver/io.cc
@@ -139,10 +139,10 @@ case horizon_file_format__ASCII_gnuplot:
IO_info.output_ghost_zones_for_h); // should we include
// ghost zones?
break;
- case gfns::gfn__H:
+ case gfns::gfn__Theta:
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING,
- " writing H to \"%s\"", file_name);
+ " writing Theta to \"%s\"", file_name);
ps.print_gridfn(unknown_gfn, file_name);
break;
case gfns::gfn__Delta_h:
diff --git a/src/driver/setup.cc b/src/driver/setup.cc
index 033cbce..d48fbe7 100644
--- a/src/driver/setup.cc
+++ b/src/driver/setup.cc
@@ -122,13 +122,13 @@ state.timer_handle = (print_timing_stats != 0) ? CCTK_TimerCreateI() : -1;
struct IO_info& IO_info = state.IO_info;
IO_info.horizon_file_format = decode_horizon_file_format(horizon_file_format);
IO_info.output_initial_guess = (output_initial_guess != 0);
-IO_info.how_often_to_output_h = how_often_to_output_h;
-IO_info.how_often_to_output_H = how_often_to_output_H_of_h;
-IO_info.output_ghost_zones_for_h = (output_ghost_zones_for_h != 0);
+IO_info.how_often_to_output_h = how_often_to_output_h;
+IO_info.how_often_to_output_Theta = how_often_to_output_Theta;
+IO_info.output_ghost_zones_for_h = (output_ghost_zones_for_h != 0);
IO_info.ASCII_gnuplot_file_name_extension = ASCII_gnuplot_file_name_extension;
IO_info.HDF5_file_name_extension = HDF5_file_name_extension;
IO_info.h_base_file_name = h_base_file_name;
-IO_info.H_base_file_name = H_of_h_base_file_name;
+IO_info.Theta_base_file_name = Theta_base_file_name;
IO_info.Delta_h_base_file_name = Delta_h_base_file_name;
IO_info.Jacobian_base_file_name = Jacobian_base_file_name;
IO_info.output_BH_diagnostics = (output_BH_diagnostics != 0);
@@ -151,10 +151,10 @@ solver_info.max_Newton_iterations__initial
solver_info.max_Newton_iterations__subsequent
= max_Newton_iterations__subsequent;
solver_info.max_Delta_h_over_h = max_Delta_h_over_h;
-solver_info.H_norm_for_convergence = H_norm_for_convergence;
+solver_info.Theta_norm_for_convergence = Theta_norm_for_convergence;
solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
-solver_info.final_H_update_if_exit_x_H_small
- = (final_H_update_if_exit_x_H_small != 0);
+solver_info.final_Theta_update_if_Delta_h_converged
+ = (final_Theta_update_if_Delta_h_converged != 0);
//
@@ -297,7 +297,7 @@ state.AH_info_ptrs.push_back(NULL);
if (verbose_info.print_algorithm_details)
then CCTK_VInfo(CCTK_THORNSTRING,
" constructing Jacobian matrix");
- AH_ptr->Jac_ptr = (state.method == method__horizon_function)
+ AH_ptr->Jac_ptr = (state.method == method__evaluate_expansion)
? NULL // no Jacobian used for this case
: &new_Jacobian(ps, Jac_info.Jacobian_storage_method);
@@ -384,10 +384,10 @@ namespace {
enum method
decode_method(const char method_string[])
{
-if (STRING_EQUAL(method_string, "horizon function"))
- then return method__horizon_function;
-else if (STRING_EQUAL(method_string, "test Jacobian"))
- then return method__test_Jacobian;
+if (STRING_EQUAL(method_string, "evaluate expansion"))
+ then return method__evaluate_expansion;
+else if (STRING_EQUAL(method_string, "test expansion Jacobian"))
+ then return method__test_expansion_Jacobian;
else if (STRING_EQUAL(method_string, "find horizon"))
then return method__find_horizon;
else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
diff --git a/src/gr.cg/expansion.c b/src/gr.cg/expansion.c
new file mode 100644
index 0000000..65f3677
--- /dev/null
+++ b/src/gr.cg/expansion.c
@@ -0,0 +1,182 @@
+/*
+ * inputs = {r, partial_d_ln_sqrt_g, partial_d_g_uu, X_ud, X_udd, g_uu, K_uu, h}
+ * outputs = {Theta_A, Theta_B, Theta_C, Theta_D}
+ * cost = 134*assignments+401*multiplications+3*divisions+5*functions+173*additions
+ */
+fp t1, t2, t3, t5, t6, t8, t9, t11, t12, t14;
+fp t15, t17, t19, t25, t26, t27, t29, t31, t34, t35;
+fp t37, t39, t40, t42, t44, t46, t47, t49, t56, t61;
+fp t63, t65, t66, t67, t82, t93, t98, t100, t102, t106;
+fp t107, t110, t111, t112, t116, t119, t120, t121, t123, t124;
+fp t127, t128, t129, t130, t131, t133, t134, t135, t137, t138;
+fp t139, t141, t142, t143, t148, t149, t150, t153, t154, t155;
+fp t158, t159, t160, t163, t164, t167, t168, t171, t172, t177;
+fp t181, t182, t185, t186, t189, t191, t197, t198, t200, t205;
+fp t220, t224, t232, t239, t266, t273, t276, t280, t283, t289;
+fp t292, t302, t303, t306, t307, t310, t311, t314, t317, t326;
+fp t330, t334, t337, t340, t343, t353, t355, t356, t360, t362;
+fp t366, t382, t387, t394, t431, t440, t444, t447, t450, t465;
+ t1 = g_uu_13;
+ t2 = t1*t1;
+ t3 = 1/r;
+ t5 = X_ud_13;
+ t6 = PARTIAL_RHO(h);
+ t8 = X_ud_23;
+ t9 = PARTIAL_SIGMA(h);
+ t11 = zz*t3-t5*t6-t8*t9;
+ t12 = t11*t11;
+ t14 = yy*yy;
+ t15 = zz*zz;
+ t17 = r*r;
+ t19 = 1/t17/r;
+ t25 = X_ud_11;
+ t26 = t25*t25;
+ t27 = PARTIAL_RHO_RHO(h);
+ t29 = X_ud_21;
+ t31 = PARTIAL_RHO_SIGMA(h);
+ t34 = t29*t29;
+ t35 = PARTIAL_SIGMA_SIGMA(h);
+ t37 = (t14+t15)*t19-X_udd_111*t6-X_udd_211*t9-t26*t27-2.0*t29*t25*t31-t34
+*t35;
+ t39 = g_uu_23;
+ t40 = t39*t39;
+ t42 = X_ud_12;
+ t44 = X_ud_22;
+ t46 = yy*t3-t42*t6-t44*t9;
+ t47 = t46*t46;
+ t49 = xx*xx;
+ t56 = t5*t5;
+ t61 = t8*t8;
+ t63 = (t49+t14)*t19-X_udd_133*t6-X_udd_233*t9-t56*t27-2.0*t8*t5*t31-t61*
+t35;
+ t65 = t1*t11;
+ t66 = g_uu_22;
+ t67 = t66*t46;
+ t82 = -xx*yy*t19-X_udd_112*t6-X_udd_212*t9-t25*t42*t27-t29*t42*t31-t25*
+t44*t31-t29*t44*t35;
+ t93 = t42*t42;
+ t98 = t44*t44;
+ t100 = (t49+t15)*t19-X_udd_122*t6-X_udd_222*t9-t93*t27-2.0*t44*t42*t31-
+t98*t35;
+ t102 = t39*t11;
+ t106 = t1*t12;
+ t107 = partial_d_g_uu_123;
+ t110 = g_uu_12;
+ t111 = t110*t47;
+ t112 = partial_d_g_uu_112;
+ t116 = xx*t3-t25*t6-t29*t9;
+ t119 = t66*t47;
+ t120 = partial_d_g_uu_212;
+ t121 = t120*t116;
+ t123 = t39*t47;
+ t124 = partial_d_g_uu_312;
+ t127 = g_uu_11;
+ t128 = t116*t116;
+ t129 = t127*t128;
+ t130 = partial_d_g_uu_113;
+ t131 = t130*t11;
+ t133 = t1*t128;
+ t134 = partial_d_g_uu_313;
+ t135 = t134*t11;
+ t137 = g_uu_33;
+ t138 = t137*t12;
+ t139 = t134*t116;
+ t141 = -t2*t12*t37-t40*t47*t63-2.0*t65*t67*t82-t40*t12*t100-2.0*t102*t67*
+t100-t106*t107*t46-t111*t112*t116-t119*t121-t123*t124*t116-t129*t131-t133*t135-
+t138*t139;
+ t142 = t39*t12;
+ t143 = partial_d_g_uu_213;
+ t148 = t1*t116;
+ t149 = partial_d_g_uu_322;
+ t150 = t149*t47;
+ t153 = t110*t116;
+ t154 = partial_d_g_uu_222;
+ t155 = t154*t47;
+ t158 = t127*t116;
+ t159 = partial_d_g_uu_122;
+ t160 = t159*t47;
+ t163 = partial_d_g_uu_333;
+ t164 = t163*t12;
+ t167 = partial_d_g_uu_133;
+ t168 = t167*t12;
+ t171 = partial_d_g_uu_233;
+ t172 = t171*t12;
+ t177 = t110*t46;
+ t181 = partial_d_g_uu_323;
+ t182 = t181*t11;
+ t185 = t137*t11;
+ t186 = t124*t46;
+ t189 = -t142*t143*t116-t106*t130*t116+RATIONAL(-1.0,2.0)*t148*t150+
+RATIONAL(-1.0,2.0)*t153*t155+RATIONAL(-1.0,2.0)*t158*t160+RATIONAL(-1.0,2.0)*
+t148*t164+RATIONAL(-1.0,2.0)*t158*t168+RATIONAL(-1.0,2.0)*t153*t172+RATIONAL(
+-1.0,2.0)*t65*t160-2.0*t65*t177*t37-t148*t182*t46-t185*t186*t116;
+ t191 = t127*t127;
+ t197 = t110*t128;
+ t198 = t143*t11;
+ t200 = t137*t137;
+ t205 = t39*t46;
+ t220 = -xx*zz*t19-X_udd_113*t6-X_udd_213*t9-t25*t5*t27-t29*t5*t31-t25*t8*
+t31-t29*t8*t35;
+ t224 = t12*t11;
+ t232 = t1*t220;
+ t239 = -t191*t128*t37-2.0*t142*t1*t82-t197*t198-t200*t12*t63-t177*t131*
+t116-2.0*t65*t205*t220+RATIONAL(-1.0,2.0)*t39*t224*t171-t67*t198*t116-t205*t135
+*t116-2.0*t138*t232+RATIONAL(-1.0,2.0)*t205*t164+RATIONAL(-1.0,2.0)*t177*t168;
+ t266 = -yy*zz*t19-X_udd_123*t6-X_udd_223*t9-t42*t5*t27-t44*t5*t31-t42*t8*
+t31-t44*t8*t35;
+ t273 = t110*t110;
+ t276 = t47*t46;
+ t280 = t39*t266;
+ t283 = t158*t37;
+ t289 = t148*t266;
+ t292 = RATIONAL(-1.0,2.0)*t67*t172+RATIONAL(-1.0,2.0)*t185*t150+RATIONAL(
+-1.0,2.0)*t102*t155-2.0*t197*t127*t82-2.0*t133*t127*t220-2.0*t133*t110*t266+
+RATIONAL(-1.0,2.0)*t1*t224*t167-t273*t128*t100+RATIONAL(-1.0,2.0)*t39*t276*t149
+-2.0*t138*t280-2.0*t65*t283+RATIONAL(-1.0,2.0)*t110*t276*t159-2.0*t67*t289;
+ t302 = partial_d_g_uu_311;
+ t303 = t302*t128;
+ t306 = partial_d_g_uu_211;
+ t307 = t306*t128;
+ t310 = partial_d_g_uu_111;
+ t311 = t310*t128;
+ t314 = t148*t63;
+ t317 = t153*t266;
+ t326 = t107*t11;
+ t330 = RATIONAL(-1.0,2.0)*t66*t276*t154-2.0*t273*t46*t116*t82+RATIONAL(
+-1.0,2.0)*t205*t303+RATIONAL(-1.0,2.0)*t67*t307+RATIONAL(-1.0,2.0)*t177*t311
+-2.0*t205*t314-2.0*t205*t317+RATIONAL(-1.0,2.0)*t185*t303+RATIONAL(-1.0,2.0)*
+t102*t307+RATIONAL(-1.0,2.0)*t65*t311-t111*t326-t158*t326*t46;
+ t334 = t158*t82;
+ t337 = t110*t82;
+ t340 = t158*t220;
+ t343 = t153*t100;
+ t353 = t112*t46;
+ t355 = partial_d_g_uu_223;
+ t356 = t355*t11;
+ t360 = t120*t46;
+ t362 = -2.0*t177*t148*t220-2.0*t67*t334-2.0*t119*t337-2.0*t205*t340-2.0*
+t67*t343+RATIONAL(-1.0,2.0)*t137*t224*t163-t2*t128*t63-t273*t47*t37-t129*t353-
+t119*t356-t123*t182-t133*t186-t197*t360;
+ t366 = t181*t46;
+ t382 = t66*t66;
+ t387 = t128*t116;
+ t394 = -t142*t355*t46-t138*t366-2.0*t177*t283-2.0*t123*t110*t220-2.0*t123
+*t66*t266-t153*t356*t46-t65*t353*t116-t102*t360*t116-t382*t47*t100-2.0*t185*
+t317+RATIONAL(-1.0,2.0)*t127*t387*t310+RATIONAL(-1.0,2.0)*t110*t387*t306;
+ t431 = RATIONAL(-1.0,2.0)*t1*t387*t302-2.0*t2*t11*t116*t220-2.0*t185*t314
+-2.0*t102*t289-2.0*t65*t153*t82-2.0*t185*t205*t63-2.0*t40*t11*t46*t266-2.0*t102
+*t343-2.0*t102*t334-2.0*t185*t340-2.0*t102*t177*t82-2.0*t185*t67*t266-2.0*t185*
+t177*t220;
+ Theta_A = t141+t189+t239+t292+t330+t362+t394+t431;
+ t440 = t310*t116+t121+t139+t353+t154*t46+t366+t131+t356+t163*t11+t127*t37
++2.0*t337+2.0*t232;
+ t444 = partial_d_ln_sqrt_g_1;
+ t447 = partial_d_ln_sqrt_g_2;
+ t450 = partial_d_ln_sqrt_g_3;
+ t465 = t66*t100+2.0*t280+t137*t63+t127*t444*t116+t110*t447*t116+t1*t450*
+t116+t110*t444*t46+t66*t447*t46+t39*t450*t46+t1*t444*t11+t39*t447*t11+t137*t450
+*t11;
+ Theta_B = t440+t465;
+ Theta_C = K_uu_11*t128+2.0*K_uu_12*t46*t116+2.0*K_uu_13*t11*t116+K_uu_22*
+t47+2.0*K_uu_23*t11*t46+K_uu_33*t12;
+ Theta_D = t129+2.0*t177*t116+2.0*t65*t116+t119+2.0*t102*t46+t138;
diff --git a/src/gr.cg/expansion_Jacobian.c b/src/gr.cg/expansion_Jacobian.c
new file mode 100644
index 0000000..b5d2427
--- /dev/null
+++ b/src/gr.cg/expansion_Jacobian.c
@@ -0,0 +1,628 @@
+/*
+ * inputs = {r, partial_d_ln_sqrt_g, partial_d_g_uu, X_ud, X_udd, g_uu, K_uu, Theta_A, Theta_B, Theta_C, Theta_D, h}
+ * outputs = {partial_Theta_wrt_partial_d_h, partial_Theta_wrt_partial_dd_h}
+ * cost = 458*assignments+1692*multiplications+10*divisions+7*functions+729*additions
+ */
+fp t1, t2, t3, t4, t5, t7, t8, t10, t11, t13;
+fp t14, t16, t18, t20, t22, t24, t26, t28, t29, t31;
+fp t32, t35, t37, t38, t41, t42, t43, t46, t48, t52;
+fp t54, t55, t59, t60, t63, t67, t68, t69, t70, t71;
+fp t74, t76, t78, t80, t83, t85, t86, t92, t93, t94;
+fp t98, t99, t102, t103, t104, t107, t108, t112, t113, t114;
+fp t115, t116, t118, t119, t120, t122, t123, t126, t127, t128;
+fp t133, t136, t140, t141, t142, t143, t153, t156, t158, t160;
+fp t162, t165, t167, t168, t171, t172, t173, t174, t179, t183;
+fp t185, t189, t190, t193, t194, t195, t197, t198, t202, t205;
+fp t208, t209, t212, t216, t217, t218, t220, t222, t223, t224;
+fp t226, t227, t232, t235, t236, t237, t238, t240, t247, t248;
+fp t249, t254, t259, t263, t266, t267, t275, t278, t281, t284;
+fp t287, t288, t291, t296, t297, t298, t300, t307, t309, t311;
+fp t314, t316, t317, t322, t325, t326, t329, t334, t335, t336;
+fp t340, t346, t350, t351, t352, t354, t357, t358, t359, t361;
+fp t364, t365, t366, t368, t370, t373, t374, t376, t381, t385;
+fp t386, t392, t398, t401, t404, t405, t407, t408, t411, t414;
+fp t416, t417, t419, t421, t422, t424, t428, t431, t432, t434;
+fp t437, t440, t442, t449, t454, t458, t461, t467, t470, t471;
+fp t474, t475, t481, t485, t489, t494, t498, t503, t504, t505;
+fp t507, t514, t518, t534, t536, t542, t545, t548, t551, t552;
+fp t559, t561, t562, t565, t569, t571, t572, t573, t575, t576;
+fp t588, t589, t590, t593, t594, t599, t601, t605, t608, t609;
+fp t612, t613, t627, t632, t633, t640, t644, t652, t656, t664;
+fp t669, t672, t677, t678, t680, t694, t704, t707, t712, t716;
+fp t723, t738, t741, t746, t748, t750, t774, t776, t780, t785;
+fp t787, t792, t796, t797, t799, t800, t802, t803, t805, t807;
+fp t809, t811, t813, t815, t817, t819, t822, t824, t827, t829;
+fp t832, t835, t837, t840, t843, t847, t860, t869, t871, t876;
+fp t882, t886, t890, t891, t897, t899, t900, t902, t904, t905;
+fp t907, t913, t920, t929, t930, t933, t938, t944, t947, t949;
+fp t962, t970, t971, t976, t979, t983, t996, t997, t1000, t1001;
+fp t1004, t1010, t1012, t1015, t1033, t1036, t1039, t1047, t1048, t1050;
+fp t1062, t1065, t1070, t1074, t1075, t1078, t1080, t1082, t1087, t1093;
+fp t1095, t1097, t1103, t1107, t1112, t1114, t1138, t1139, t1141, t1145;
+fp t1150, t1163, t1166, t1169, t1174, t1186, t1189, t1192, t1200, t1214;
+fp t1234, t1266, t1281, t1289, t1300, t1301, t1308, t1335, t1342, t1345;
+fp t1364, t1370, t1405, t1414, t1427, t1457, t1460, t1463, t1465, t1469;
+fp t1475, t1476, t1477, t1483, t1486, t1487, t1491, t1492, t1493, t1497;
+fp t1505, t1508, t1510, t1513, t1516, t1517, t1520, t1526, t1536, t1547;
+fp t1552, t1555, t1558, t1561, t1572, t1580, t1594, t1600, t1606, t1610;
+fp t1622, t1629, t1639, t1641, t1643, t1645, t1648, t1655, t1659, t1660;
+fp t1666, t1667, t1684, t1697, t1704, t1718, t1721, t1739, t1748, t1751;
+fp t1757, t1760, t1761, t1768, t1771, t1783, t1785, t1788, t1791, t1803;
+fp t1809, t1812, t1825;
+ t1 = g_uu_13;
+ t2 = X_ud_13;
+ t3 = t1*t2;
+ t4 = g_uu_12;
+ t5 = 1/r;
+ t7 = X_ud_11;
+ t8 = PARTIAL_RHO(h);
+ t10 = X_ud_21;
+ t11 = PARTIAL_SIGMA(h);
+ t13 = xx*t5-t7*t8-t10*t11;
+ t14 = t4*t13;
+ t16 = r*r;
+ t18 = 1/t16/r;
+ t20 = X_udd_112;
+ t22 = X_udd_212;
+ t24 = X_ud_12;
+ t26 = PARTIAL_RHO_RHO(h);
+ t28 = t10*t24;
+ t29 = PARTIAL_RHO_SIGMA(h);
+ t31 = X_ud_22;
+ t32 = t7*t31;
+ t35 = PARTIAL_SIGMA_SIGMA(h);
+ t37 = -xx*yy*t18-t20*t8-t22*t11-t7*t24*t26-t28*t29-t32*t29-t10*t31*t35;
+ t38 = t14*t37;
+ t41 = g_uu_22;
+ t42 = t41*t24;
+ t43 = t1*t13;
+ t46 = X_udd_123;
+ t48 = X_udd_223;
+ t52 = t31*t2;
+ t54 = X_ud_23;
+ t55 = t24*t54;
+ t59 = -yy*zz*t18-t46*t8-t48*t11-t24*t2*t26-t52*t29-t55*t29-t31*t54*t35;
+ t60 = t43*t59;
+ t63 = g_uu_23;
+ t67 = yy*t5-t24*t8-t31*t11;
+ t68 = t63*t67;
+ t69 = t1*t7;
+ t70 = xx*xx;
+ t71 = yy*yy;
+ t74 = X_udd_133;
+ t76 = X_udd_233;
+ t78 = t2*t2;
+ t80 = t54*t2;
+ t83 = t54*t54;
+ t85 = (t70+t71)*t18-t74*t8-t76*t11-t78*t26-2.0*t80*t29-t83*t35;
+ t86 = t69*t85;
+ t92 = zz*t5-t2*t8-t54*t11;
+ t93 = t63*t92;
+ t94 = t4*t67;
+ t98 = t41*t67;
+ t99 = t69*t59;
+ t102 = g_uu_33;
+ t103 = t102*t92;
+ t104 = t43*t74;
+ t107 = t1*t92;
+ t108 = t4*t7;
+ t112 = g_uu_11;
+ t113 = t112*t13;
+ t114 = partial_d_g_uu_123;
+ t115 = t114*t2;
+ t116 = t115*t67;
+ t118 = partial_d_g_uu_211;
+ t119 = t118*t13;
+ t120 = t119*t7;
+ t122 = t63*t2;
+ t123 = t94*t37;
+ t126 = partial_d_g_uu_122;
+ t127 = t126*t67;
+ t128 = t127*t24;
+ t133 = t98*t37;
+ t136 = X_udd_113;
+ t140 = 2.0*t3*t38+2.0*t42*t60+2.0*t68*t86+2.0*t93*t94*t20+2.0*t98*t99+2.0
+*t103*t104+2.0*t107*t108*t37+t113*t116+t93*t120+2.0*t122*t123+t113*t128+2.0*
+t107*t14*t20+2.0*t3*t133+2.0*t107*t68*t136;
+ t141 = partial_d_g_uu_311;
+ t142 = t141*t13;
+ t143 = t142*t7;
+ t153 = zz*zz;
+ t156 = X_udd_122;
+ t158 = X_udd_222;
+ t160 = t24*t24;
+ t162 = t31*t24;
+ t165 = t31*t31;
+ t167 = (t70+t153)*t18-t156*t8-t158*t11-t160*t26-2.0*t162*t29-t165*t35;
+ t168 = t108*t167;
+ t171 = t13*t13;
+ t172 = t112*t171;
+ t173 = partial_d_g_uu_112;
+ t174 = t173*t24;
+ t179 = X_udd_213;
+ t183 = t10*t2;
+ t185 = t7*t54;
+ t189 = -xx*zz*t18-t136*t8-t179*t11-t7*t2*t26-t183*t29-t185*t29-t10*t54*
+t35;
+ t190 = t68*t189;
+ t193 = t112*t7;
+ t194 = t114*t92;
+ t195 = t194*t67;
+ t197 = t4*t4;
+ t198 = t197*t67;
+ t202 = t108*t59;
+ t205 = t193*t37;
+ t208 = t102*t2;
+ t209 = t14*t59;
+ t212 = t63*t24;
+ t216 = t63*t63;
+ t217 = t92*t92;
+ t218 = t216*t217;
+ t220 = t103*t143+2.0*t94*t43*t136+2.0*t107*t98*t20+2.0*t68*t104+2.0*t93*
+t168+t172*t174+2.0*t3*t190+t193*t195+2.0*t198*t7*t37+2.0*t103*t202+2.0*t93*t205
++2.0*t208*t209+2.0*t107*t212*t189+t218*t156;
+ t222 = t1*t1;
+ t223 = t222*t217;
+ t224 = X_udd_111;
+ t226 = t102*t102;
+ t227 = t226*t217;
+ t232 = t113*t189;
+ t235 = t67*t67;
+ t236 = t41*t235;
+ t237 = partial_d_g_uu_223;
+ t238 = t237*t2;
+ t240 = t194*t24;
+ t247 = partial_d_g_uu_333;
+ t248 = t247*t92;
+ t249 = t248*t2;
+ t254 = t113*t136;
+ t259 = t1*t171;
+ t263 = t193*t189;
+ t266 = t223*t224+t227*t74+2.0*t107*t42*t37+2.0*t208*t232+t236*t238+t113*
+t240+2.0*t93*t98*t156+2.0*t68*t202+t43*t249+2.0*t93*t42*t167+2.0*t103*t254+2.0*
+t212*t209+2.0*t259*t4*t46+2.0*t103*t263;
+ t267 = t98*t167;
+ t275 = t14*t46;
+ t278 = t43*t46;
+ t281 = t113*t224;
+ t284 = t113*t37;
+ t287 = t102*t217;
+ t288 = t63*t46;
+ t291 = t113*t20;
+ t296 = partial_d_g_uu_312;
+ t297 = t296*t67;
+ t298 = t297*t13;
+ t300 = t222*t92;
+ t307 = X_udd_211;
+ t309 = t7*t7;
+ t311 = t10*t7;
+ t314 = t10*t10;
+ t316 = (t71+t153)*t18-t224*t8-t307*t11-t309*t26-2.0*t311*t29-t314*t35;
+ t317 = t113*t316;
+ t322 = 2.0*t122*t267+2.0*t94*t69*t189+4.0*t43*t263+2.0*t103*t275+2.0*t98*
+t278+2.0*t107*t281+2.0*t122*t284+2.0*t287*t288+2.0*t93*t291+2.0*t68*t275+t208*
+t298+2.0*t300*t7*t189+2.0*t3*t317+2.0*t103*t86;
+ t325 = t4*t24;
+ t326 = t325*t189;
+ t329 = t43*t85;
+ t334 = partial_d_g_uu_313;
+ t335 = t334*t92;
+ t336 = t335*t13;
+ t340 = t335*t7;
+ t346 = t63*t59;
+ t350 = partial_d_g_uu_111;
+ t351 = t350*t13;
+ t352 = t351*t7;
+ t354 = t193*t316;
+ t357 = partial_d_g_uu_113;
+ t358 = t357*t2;
+ t359 = t358*t13;
+ t361 = t94*t189;
+ t364 = partial_d_g_uu_323;
+ t365 = t364*t2;
+ t366 = t365*t67;
+ t368 = 2.0*t103*t326+2.0*t208*t329+2.0*t212*t329+t212*t336+4.0*t68*t326+
+t68*t340+2.0*t93*t278+4.0*t43*t202+4.0*t103*t346*t2+t94*t352+2.0*t107*t354+t94*
+t359+2.0*t208*t361+t43*t366;
+ t370 = t41*t59*t24;
+ t373 = t357*t92;
+ t374 = t373*t13;
+ t376 = t1*t189;
+ t381 = t63*t235;
+ t385 = partial_d_g_uu_133;
+ t386 = t385*t217;
+ t392 = t4*t20;
+ t398 = t350*t171;
+ t401 = t118*t171;
+ t404 = t334*t2;
+ t405 = t404*t13;
+ t407 = t4*t37;
+ t408 = t407*t24;
+ t411 = t43*t189;
+ t414 = 4.0*t68*t370+t325*t374+4.0*t103*t376*t2+t98*t120+2.0*t381*t41*t46+
+RATIONAL(1.0,2.0)*t193*t386+2.0*t381*t4*t136+2.0*t236*t392+2.0*t259*t112*t136+
+RATIONAL(1.0,2.0)*t3*t398+RATIONAL(1.0,2.0)*t122*t401+t68*t405+4.0*t98*t408+2.0
+*t325*t411;
+ t416 = t364*t92;
+ t417 = t416*t67;
+ t419 = t297*t7;
+ t421 = t296*t24;
+ t422 = t421*t13;
+ t424 = t1*t37;
+ t428 = t94*t316;
+ t431 = t41*t41;
+ t432 = t431*t235;
+ t434 = t126*t235;
+ t437 = t247*t217;
+ t440 = t416*t24;
+ t442 = t373*t7;
+ t449 = t431*t67;
+ t454 = t69*t417+t103*t419+t103*t422+4.0*t93*t424*t2+2.0*t3*t428+t432*t156
++RATIONAL(1.0,2.0)*t193*t434+RATIONAL(1.0,2.0)*t69*t437+t43*t440+t94*t442+2.0*
+t300*t13*t136+t381*t296*t7+2.0*t449*t167*t24+t259*t421;
+ t458 = t350*t7;
+ t461 = t4*t235;
+ t467 = t13*t189;
+ t470 = t237*t92;
+ t471 = t470*t24;
+ t474 = t385*t92;
+ t475 = t474*t2;
+ t481 = t13*t37;
+ t485 = t67*t59;
+ t489 = t238*t67;
+ t494 = RATIONAL(3.0,2.0)*t259*t141*t7+RATIONAL(3.0,2.0)*t172*t458+t461*
+t115+2.0*t198*t13*t20+2.0*t222*t2*t467+2.0*t98*t471+t113*t475+2.0*t107*t94*t224
++2.0*t197*t24*t481+2.0*t216*t2*t485+t68*t249+t14*t489+t107*t128+2.0*t93*t99;
+ t498 = t470*t67;
+ t503 = partial_d_g_uu_233;
+ t504 = t503*t92;
+ t505 = t504*t2;
+ t507 = t4*t171;
+ t514 = t216*t92;
+ t518 = t334*t7;
+ t534 = t108*t498+2.0*t103*t94*t136+t14*t505+RATIONAL(3.0,2.0)*t507*t118*
+t7+2.0*t107*t325*t316+2.0*t514*t24*t59+t287*t518+t259*t404+RATIONAL(3.0,2.0)*
+t461*t126*t24+2.0*t514*t67*t46+RATIONAL(1.0,2.0)*t3*t434+2.0*t68*t440+t172*t358
++2.0*t68*t422;
+ t536 = partial_d_g_uu_213;
+ t542 = t98*t59;
+ t545 = t68*t85;
+ t548 = t216*t235;
+ t551 = t536*t13;
+ t552 = t551*t2;
+ t559 = t174*t13;
+ t561 = t536*t92;
+ t562 = t561*t7;
+ t565 = t226*t92;
+ t569 = t94*t475+t507*t536*t2+2.0*t43*t340+t14*t471+2.0*t208*t542+2.0*t208
+*t545+t548*t74+t98*t505+2.0*t93*t552+2.0*t94*t240+2.0*t113*t442+t107*t559+2.0*
+t14*t562+2.0*t565*t85*t2;
+ t571 = partial_d_g_uu_322;
+ t572 = t571*t67;
+ t573 = t572*t24;
+ t575 = t173*t67;
+ t576 = t575*t13;
+ t588 = partial_d_g_uu_212;
+ t589 = t588*t24;
+ t590 = t589*t13;
+ t593 = t588*t67;
+ t594 = t593*t13;
+ t599 = t575*t7;
+ t601 = t63*t217;
+ t605 = t141*t171;
+ t608 = t43*t573+t3*t576+2.0*t103*t405+2.0*t43*t419+t103*t573+2.0*t107*
+t359+2.0*t514*t167*t2+t93*t590+t381*t365+t122*t594+2.0*t103*t98*t46+t107*t599+
+2.0*t601*t1*t20+RATIONAL(1.0,2.0)*t208*t605;
+ t609 = t593*t7;
+ t612 = partial_d_g_uu_222;
+ t613 = t612*t24;
+ t627 = t588*t7;
+ t632 = t612*t67;
+ t633 = t632*t24;
+ t640 = t216*t67;
+ t644 = 2.0*t14*t609+RATIONAL(3.0,2.0)*t236*t613+t93*t609+2.0*t113*t599+
+RATIONAL(1.0,2.0)*t42*t401+2.0*t107*t116+RATIONAL(1.0,2.0)*t325*t398+2.0*t103*
+t366+t236*t627+2.0*t103*t212*t85+t14*t633+2.0*t93*t489+RATIONAL(3.0,2.0)*t381*
+t571*t24+2.0*t640*t85*t24;
+ t652 = t364*t24;
+ t656 = t1*t217;
+ t664 = t247*t2;
+ t669 = t1*t136;
+ t672 = t503*t217;
+ t677 = t112*t112;
+ t678 = t677*t171;
+ t680 = 4.0*t14*t205+2.0*t103*t68*t74+t287*t652+t461*t173*t7+t656*t114*t24
++t601*t237*t24+t507*t589+t601*t536*t7+RATIONAL(3.0,2.0)*t287*t664+RATIONAL(1.0,
+2.0)*t212*t437+2.0*t287*t669+RATIONAL(1.0,2.0)*t108*t672+RATIONAL(1.0,2.0)*t42*
+t672+t678*t224;
+ t694 = t677*t13;
+ t704 = t571*t235;
+ t707 = t612*t235;
+ t712 = t222*t13;
+ t716 = 2.0*t98*t590+2.0*t300*t316*t2+2.0*t94*t559+t98*t562+2.0*t122*t60+
+t93*t633+2.0*t103*t370+2.0*t694*t316*t7+RATIONAL(3.0,2.0)*t656*t385*t2+RATIONAL
+(3.0,2.0)*t601*t503*t2+RATIONAL(1.0,2.0)*t208*t704+RATIONAL(1.0,2.0)*t122*t707+
+RATIONAL(1.0,2.0)*t69*t704+2.0*t712*t85*t7;
+ t723 = t197*t13;
+ t738 = t14*t167;
+ t741 = t14*t156;
+ t746 = t561*t13;
+ t748 = t197*t235;
+ t750 = 2.0*t198*t316*t24+t656*t357*t7+2.0*t723*t167*t7+t68*t143+2.0*t507*
+t112*t20+2.0*t94*t354+t98*t552+RATIONAL(1.0,2.0)*t108*t707+RATIONAL(1.0,2.0)*
+t212*t605+2.0*t122*t738+2.0*t98*t741+2.0*t93*t408+t42*t746+t748*t224;
+ t774 = t197*t171;
+ t776 = t222*t171;
+ t780 = 2.0*t94*t281+2.0*t42*t284+2.0*t98*t168+t107*t352+2.0*t212*t232+2.0
+*t93*t741+RATIONAL(1.0,2.0)*t325*t386+2.0*t42*t738+2.0*t98*t205+2.0*t98*t291+
+2.0*t325*t317+2.0*t68*t254+t774*t156+t776*t74+2.0*t68*t263;
+ t785 = pow(Theta_D,1.0*RATIONAL(1.0,2.0));
+ t787 = 1/t785/Theta_D;
+ t792 = -t458-t627-t518-t174-t613-t652-t358-t238-t664-t112*t224-2.0*t392
+-2.0*t669;
+ t796 = partial_d_ln_sqrt_g_1;
+ t797 = t112*t796;
+ t799 = partial_d_ln_sqrt_g_2;
+ t800 = t4*t799;
+ t802 = partial_d_ln_sqrt_g_3;
+ t803 = t1*t802;
+ t805 = t4*t796;
+ t807 = t41*t799;
+ t809 = t63*t802;
+ t811 = t1*t796;
+ t813 = t63*t799;
+ t815 = t102*t802;
+ t817 = -t41*t156-2.0*t288-t102*t74-t797*t7-t800*t7-t803*t7-t805*t24-t807*
+t24-t809*t24-t811*t2-t813*t2-t815*t2;
+ t819 = 1/t785;
+ t822 = K_uu_11*t13;
+ t824 = K_uu_12;
+ t827 = t824*t67;
+ t829 = K_uu_13;
+ t832 = t829*t92;
+ t835 = K_uu_22*t67;
+ t837 = K_uu_23;
+ t840 = t837*t92;
+ t843 = K_uu_33*t92;
+ t847 = 1/Theta_D;
+ t860 = Theta_D*Theta_D;
+ t869 = RATIONAL(3.0,2.0)*Theta_A/t785/t860+RATIONAL(1.0,2.0)*Theta_B*t787
++Theta_C/t860;
+ partial_Theta_wrt_partial_d_h_1 = (t140+t220+t266+t322+t368+t414+t454+
+t494+t534+t569+t608+t644+t680+t716+t750+t780)*t787+(t792+t817)*t819+(-2.0*t822*
+t7-2.0*t824*t24*t13-2.0*t827*t7-2.0*t829*t2*t13-2.0*t832*t7-2.0*t835*t24-2.0*
+t837*t2*t67-2.0*t840*t24-2.0*t843*t2)*t847-(-2.0*t113*t7-2.0*t325*t13-2.0*t94*
+t7-2.0*t3*t13-2.0*t107*t7-2.0*t98*t24-2.0*t122*t67-2.0*t93*t24-2.0*t103*t2)*
+t869;
+ t871 = t113*t22;
+ t876 = t63*t54;
+ t882 = t551*t54;
+ t886 = t561*t10;
+ t890 = t112*t10;
+ t891 = t890*t316;
+ t897 = t334*t10;
+ t899 = 2.0*t93*t871+t381*t296*t10+t876*t594+2.0*t93*t94*t22+t432*t158+2.0
+*t93*t882+t218*t158+2.0*t14*t886+t748*t307+2.0*t94*t891+t890*t195+t548*t76+t223
+*t307+t287*t897;
+ t900 = t194*t31;
+ t902 = t334*t54;
+ t904 = t114*t54;
+ t905 = t904*t67;
+ t907 = t63*t31;
+ t913 = t102*t54;
+ t920 = t14*t48;
+ t929 = t4*t10;
+ t930 = t929*t59;
+ t933 = t335*t10;
+ t938 = t113*t900+t259*t902+t113*t905+2.0*t907*t209+2.0*t300*t13*t179+2.0*
+t913*t545+t507*t536*t54+t601*t536*t10+2.0*t68*t920+2.0*t712*t85*t10+2.0*t449*
+t167*t31+2.0*t68*t930+t68*t933+2.0*t197*t31*t481;
+ t944 = t1*t54;
+ t947 = t588*t31;
+ t949 = t113*t307;
+ t962 = t364*t54;
+ t970 = t4*t31;
+ t971 = t970*t189;
+ t976 = t913*t298+2.0*t103*t907*t85+2.0*t944*t317+t507*t947+2.0*t107*t949+
+RATIONAL(3.0,2.0)*t507*t118*t10+2.0*t259*t4*t48+t259*t296*t31+2.0*t107*t891+
+t381*t962+2.0*t198*t13*t22+2.0*t103*t68*t76+2.0*t103*t971+2.0*t876*t60;
+ t979 = t416*t31;
+ t983 = t351*t10;
+ t996 = t1*t10;
+ t997 = t996*t85;
+ t1000 = t41*t31;
+ t1001 = t1000*t59;
+ t1004 = t996*t59;
+ t1010 = t142*t10;
+ t1012 = 2.0*t907*t329+t43*t979+2.0*t913*t361+t107*t983+2.0*t944*t38+4.0*
+t93*t424*t54+2.0*t107*t929*t37+2.0*t103*t94*t179+2.0*t68*t997+2.0*t103*t1001+
+2.0*t98*t1004+t970*t374+2.0*t913*t542+t103*t1010;
+ t1015 = t119*t10;
+ t1033 = t43*t48;
+ t1036 = t297*t10;
+ t1039 = t373*t10;
+ t1047 = t357*t54;
+ t1048 = t1047*t13;
+ t1050 = t93*t1015+2.0*t107*t1000*t37+2.0*t259*t112*t179+2.0*t970*t411+2.0
+*t944*t133+2.0*t93*t1004+2.0*t103*t98*t48+t774*t158+2.0*t98*t1033+2.0*t43*t1036
++2.0*t113*t1039+2.0*t1000*t60+2.0*t94*t43*t179+t94*t1048;
+ t1062 = t43*t76;
+ t1065 = t113*t179;
+ t1070 = t470*t31;
+ t1074 = t237*t54;
+ t1075 = t1074*t67;
+ t1078 = t504*t54;
+ t1080 = t474*t54;
+ t1082 = 2.0*t107*t98*t22+2.0*t94*t996*t189+t98*t1015+2.0*t970*t317+2.0*
+t876*t123+2.0*t68*t1062+2.0*t68*t1065+2.0*t944*t428+t14*t1070+2.0*t94*t949+t14*
+t1075+t94*t1039+t14*t1078+t113*t1080;
+ t1087 = t112*t189*t10;
+ t1093 = t248*t54;
+ t1095 = t127*t31;
+ t1097 = t572*t31;
+ t1103 = t296*t13*t31;
+ t1107 = t407*t31;
+ t1112 = t632*t31;
+ t1114 = 4.0*t43*t930+4.0*t43*t1087+t227*t76+4.0*t68*t971+t68*t1093+t107*
+t1095+t103*t1097+4.0*t68*t1001+t98*t1078+2.0*t68*t1103+t678*t307+4.0*t98*t1107+
+RATIONAL(1.0,2.0)*t1000*t672+t93*t1112;
+ t1138 = t173*t31;
+ t1139 = t1138*t13;
+ t1141 = t4*t22;
+ t1145 = 2.0*t381*t41*t48+2.0*t216*t54*t485+t103*t1103+RATIONAL(1.0,2.0)*
+t996*t704+t103*t1036+t601*t237*t31+t172*t1047+RATIONAL(3.0,2.0)*t601*t503*t54+
+2.0*t514*t31*t59+2.0*t514*t67*t48+t776*t76+t107*t1139+2.0*t236*t1141+t996*t417;
+ t1150 = t593*t10;
+ t1163 = t947*t13;
+ t1166 = t962*t67;
+ t1169 = t575*t10;
+ t1174 = t944*t576+t93*t1150+2.0*t723*t167*t10+RATIONAL(1.0,2.0)*t890*t386
++RATIONAL(1.0,2.0)*t929*t672+RATIONAL(1.0,2.0)*t944*t434+RATIONAL(1.0,2.0)*t907
+*t437+t93*t1163+t98*t882+t43*t1166+t1000*t746+t107*t1169+t43*t1097+2.0*t107*
+t1048;
+ t1186 = t112*t37*t10;
+ t1189 = t929*t167;
+ t1192 = t14*t158;
+ t1200 = t43*t1093+t907*t336+2.0*t98*t1070+t113*t1095+2.0*t113*t1169+t68*
+t1010+t98*t886+t94*t1080+4.0*t14*t1186+2.0*t93*t1189+2.0*t93*t1192+2.0*t913*
+t232+2.0*t876*t738+t14*t1112;
+ t1214 = t902*t13;
+ t1234 = 2.0*t107*t14*t22+2.0*t1000*t738+2.0*t876*t267+2.0*t107*t68*t179+
+2.0*t94*t900+t68*t1214+2.0*t103*t920+2.0*t944*t190+2.0*t68*t1087+2.0*t93*t98*
+t158+2.0*t907*t232+2.0*t93*t1000*t167+2.0*t98*t1192+2.0*t98*t1189;
+ t1266 = 2.0*t103*t1065+t94*t983+2.0*t1000*t284+2.0*t198*t316*t31+2.0*t107
+*t907*t189+RATIONAL(1.0,2.0)*t890*t434+2.0*t103*t1166+2.0*t43*t933+2.0*t103*
+t930+2.0*t94*t1139+2.0*t98*t1163+2.0*t98*t1186+RATIONAL(3.0,2.0)*t259*t141*t10+
+2.0*t222*t54*t467;
+ t1281 = t588*t10;
+ t1289 = t247*t54;
+ t1300 = 2.0*t300*t10*t189+RATIONAL(3.0,2.0)*t656*t385*t54+t461*t904+t172*
+t1138+t236*t1074+2.0*t694*t316*t10+t236*t1281+RATIONAL(3.0,2.0)*t381*t571*t31+
+RATIONAL(3.0,2.0)*t461*t126*t31+RATIONAL(3.0,2.0)*t287*t1289+2.0*t103*t1087+2.0
+*t107*t905+2.0*t68*t979+2.0*t98*t871;
+ t1301 = t612*t31;
+ t1308 = t364*t31;
+ t1335 = RATIONAL(3.0,2.0)*t236*t1301+2.0*t93*t1075+2.0*t14*t1150+t287*
+t1308+t656*t114*t31+t461*t173*t10+RATIONAL(1.0,2.0)*t1000*t401+t656*t357*t10+
+2.0*t300*t316*t54+2.0*t507*t112*t22+2.0*t640*t85*t31+2.0*t601*t1*t22+RATIONAL(
+1.0,2.0)*t876*t707+2.0*t565*t85*t54;
+ t1342 = t63*t48;
+ t1345 = t1*t179;
+ t1364 = t350*t10;
+ t1370 = RATIONAL(1.0,2.0)*t913*t704+2.0*t514*t167*t54+2.0*t287*t1342+2.0*
+t287*t1345+RATIONAL(1.0,2.0)*t970*t398+RATIONAL(1.0,2.0)*t996*t437+RATIONAL(1.0
+,2.0)*t907*t605+RATIONAL(1.0,2.0)*t944*t398+RATIONAL(1.0,2.0)*t876*t401+
+RATIONAL(1.0,2.0)*t913*t605+RATIONAL(1.0,2.0)*t929*t707+RATIONAL(1.0,2.0)*t970*
+t386+RATIONAL(3.0,2.0)*t172*t1364+2.0*t198*t10*t37;
+ t1405 = 4.0*t103*t376*t54+2.0*t103*t1214+4.0*t103*t346*t54+2.0*t93*t1107+
+2.0*t107*t94*t307+2.0*t107*t970*t316+2.0*t913*t209+2.0*t381*t4*t179+t929*t498+
+2.0*t93*t1033+2.0*t103*t997+2.0*t103*t1062+2.0*t913*t329+2.0*t876*t284+2.0*t93*
+t1186;
+ t1414 = -t1364-t1281-t897-t1138-t1301-t1308-t1047-t1074-t1289-t112*t307
+-2.0*t1141-2.0*t1345;
+ t1427 = -t41*t158-2.0*t1342-t102*t76-t797*t10-t800*t10-t803*t10-t805*t31-
+t807*t31-t809*t31-t811*t54-t813*t54-t815*t54;
+ partial_Theta_wrt_partial_d_h_2 = (t899+t938+t976+t1012+t1050+t1082+t1114
++t1145+t1174+t1200+t1234+t1266+t1300+t1335+t1370+t1405)*t787+(t1414+t1427)*t819
++(-2.0*t822*t10-2.0*t824*t31*t13-2.0*t827*t10-2.0*t829*t54*t13-2.0*t832*t10-2.0
+*t835*t31-2.0*t837*t54*t67-2.0*t840*t31-2.0*t843*t54)*t847-(-2.0*t113*t10-2.0*
+t970*t13-2.0*t94*t10-2.0*t944*t13-2.0*t107*t10-2.0*t98*t31-2.0*t876*t67-2.0*t93
+*t31-2.0*t103*t54)*t869;
+ t1457 = t14*t160;
+ t1460 = t69*t2;
+ t1463 = t68*t4;
+ t1465 = t13*t24*t2;
+ t1469 = t43*t78;
+ t1475 = t103*t4;
+ t1476 = t67*t7;
+ t1477 = t1476*t2;
+ t1483 = t212*t2;
+ t1486 = t107*t41;
+ t1487 = t1476*t24;
+ t1491 = 2.0*t98*t1457+2.0*t287*t1460+2.0*t1463*t1465+t774*t160+2.0*t68*
+t1469+2.0*t107*t94*t309+2.0*t1475*t1477+2.0*t381*t108*t2+2.0*t287*t1483+2.0*
+t1486*t1487+t218*t160;
+ t1492 = t13*t7;
+ t1493 = t1492*t24;
+ t1497 = t113*t309;
+ t1505 = t98*t1;
+ t1508 = t103*t41;
+ t1510 = t67*t24*t2;
+ t1513 = t93*t4;
+ t1516 = t94*t1;
+ t1517 = t1492*t2;
+ t1520 = t107*t4;
+ t1526 = 2.0*t198*t1493+t223*t309+2.0*t107*t1497+2.0*t259*t193*t2+2.0*t93*
+t1457+2.0*t1505*t1465+2.0*t1508*t1510+2.0*t1513*t1487+2.0*t1516*t1517+2.0*t1520
+*t1493+2.0*t103*t68*t78;
+ t1536 = t93*t112;
+ t1547 = t93*t1;
+ t1552 = t432*t160+2.0*t93*t98*t160+t548*t78+2.0*t514*t1510+t748*t309+2.0*
+t1536*t1493+2.0*t300*t1517+2.0*t381*t42*t2+2.0*t507*t193*t24+2.0*t1547*t1465+
+2.0*t1475*t1465;
+ t1555 = t103*t112;
+ t1558 = t98*t112;
+ t1561 = t108*t24;
+ t1572 = t68*t112;
+ t1580 = 2.0*t1547*t1477+2.0*t1555*t1517+2.0*t1558*t1493+2.0*t236*t1561+
+2.0*t259*t325*t2+t776*t78+t678*t309+t227*t78+2.0*t94*t1497+2.0*t1572*t1517+2.0*
+t103*t1469+2.0*t601*t69*t24;
+ partial_Theta_wrt_partial_dd_h_11 = (t1491+t1526+t1552+t1580)*t787+(-t112
+*t309-2.0*t1561-2.0*t1460-t41*t160-2.0*t1483-t102*t78)*t819;
+ t1594 = -t183-t185;
+ t1600 = t67*t10;
+ t1606 = -t28-t32;
+ t1610 = t1*t1594;
+ t1622 = 2.0*t218*t162-2.0*t107*t68*t1594+2.0*t432*t162+4.0*t1520*t1600*t7
++2.0*t776*t80-2.0*t601*t1*t1606-2.0*t287*t1610+2.0*t223*t311+2.0*t748*t311-2.0*
+t93*t94*t1606+2.0*t774*t162;
+ t1629 = -t52-t55;
+ t1639 = t113*t1606;
+ t1641 = t113*t1594;
+ t1643 = t14*t1629;
+ t1645 = -t381*t4*t1594-t300*t13*t1594-t107*t98*t1606-t259*t4*t1629-t103*
+t94*t1594-t107*t14*t1606+t678*t311-t507*t112*t1606-t93*t1639-t103*t1641-t103*
+t1643;
+ t1648 = t43*t1629;
+ t1655 = t13*t54*t2;
+ t1659 = t13*t10;
+ t1660 = t1659*t7;
+ t1666 = t13*t31;
+ t1667 = t1666*t24;
+ t1684 = -2.0*t93*t1648+2.0*t227*t80+4.0*t103*t1*t1655+4.0*t94*t112*t1660
+-2.0*t198*t13*t1606+4.0*t1513*t1667-2.0*t514*t67*t1629-2.0*t94*t43*t1594+4.0*
+t68*t1*t1655+4.0*t98*t4*t1667-2.0*t98*t1639;
+ t1697 = t4*t1606;
+ t1704 = t67*t31;
+ t1718 = t63*t1629;
+ t1721 = -2.0*t259*t112*t1594-2.0*t68*t1643-2.0*t68*t1641-2.0*t98*t1648+
+4.0*t107*t112*t1660-2.0*t236*t1697-2.0*t381*t41*t1629+4.0*t93*t41*t1704*t24-2.0
+*t103*t98*t1629+2.0*t548*t80+4.0*t103*t63*t67*t54*t2-2.0*t287*t1718;
+ partial_Theta_wrt_partial_dd_h_12 = (t1622+2.0*t1645+t1684+t1721)*t787+(
+-2.0*t890*t7+2.0*t1697+2.0*t1610-2.0*t1000*t24+2.0*t1718-2.0*t913*t2)*t819;
+ t1739 = t996*t54;
+ t1748 = t1704*t54;
+ t1751 = t1600*t54;
+ t1757 = t1600*t31;
+ t1760 = 2.0*t507*t890*t31+2.0*t93*t98*t165+t227*t83+t548*t83+2.0*t287*
+t1739+2.0*t259*t970*t54+2.0*t601*t996*t31+2.0*t514*t1748+2.0*t1547*t1751+2.0*
+t107*t94*t314+2.0*t1486*t1757;
+ t1761 = t907*t54;
+ t1768 = t1659*t31;
+ t1771 = t113*t314;
+ t1783 = 2.0*t287*t1761+t748*t314+t774*t165+t678*t314+t223*t314+2.0*t198*
+t1768+2.0*t94*t1771+2.0*t1513*t1757+2.0*t1520*t1768+2.0*t1475*t1751+2.0*t103*
+t68*t83;
+ t1785 = t1666*t54;
+ t1788 = t1659*t54;
+ t1791 = t43*t83;
+ t1803 = t14*t165;
+ t1809 = 2.0*t1547*t1785+2.0*t1516*t1788+2.0*t68*t1791+2.0*t1558*t1768+2.0
+*t259*t890*t54+2.0*t107*t1771+2.0*t1463*t1785+2.0*t98*t1803+t218*t165+t776*t83+
+t432*t165;
+ t1812 = t929*t31;
+ t1825 = t1572*t1788+t1505*t1785+t236*t1812+t103*t1791+t300*t1788+t381*
+t1000*t54+t381*t929*t54+t93*t1803+t1555*t1788+t1536*t1768+t1475*t1785+t1508*
+t1748;
+ partial_Theta_wrt_partial_dd_h_22 = (t1760+t1783+t1809+2.0*t1825)*t787+(-
+t112*t314-2.0*t1812-2.0*t1739-t41*t165-2.0*t1761-t102*t83)*t819;
diff --git a/src/gr.cg/horizon_Jacobian.c b/src/gr.cg/horizon_Jacobian.c
deleted file mode 100644
index eff036b..0000000
--- a/src/gr.cg/horizon_Jacobian.c
+++ /dev/null
@@ -1,625 +0,0 @@
-/*
- * inputs = {r, partial_d_ln_sqrt_g, partial_d_g_uu, X_ud, X_udd, g_uu, K_uu, h, HA, HB, HC, HD}
- * outputs = {partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h}
- * cost = 457*assignments+10*divisions+7*functions+1685*multiplications+729*additions
- */
-fp t1, t2, t4, t5, t7, t8, t10, t11, t12, t14;
-fp t16, t18, t19, t20, t21, t23, t24, t25, t26, t29;
-fp t30, t33, t35, t36, t37, t38, t40, t42, t44, t46;
-fp t48, t49, t51, t52, t55, t56, t58, t59, t62, t63;
-fp t66, t68, t72, t74, t78, t82, t83, t86, t87, t88;
-fp t91, t92, t93, t94, t95, t97, t98, t101, t104, t107;
-fp t108, t109, t112, t113, t115, t119, t121, t125, t127, t131;
-fp t132, t135, t136, t137, t138, t142, t143, t146, t147, t150;
-fp t151, t152, t154, t157, t160, t162, t164, t166, t169, t171;
-fp t175, t182, t186, t188, t194, t196, t199, t206, t207, t210;
-fp t213, t215, t216, t219, t220, t221, t224, t226, t229, t235;
-fp t236, t239, t240, t242, t243, t244, t247, t248, t255, t259;
-fp t261, t265, t266, t270, t271, t273, t274, t275, t277, t280;
-fp t281, t282, t283, t286, t289, t290, t292, t295, t296, t299;
-fp t303, t309, t313, t319, t320, t323, t327, t333, t335, t337;
-fp t339, t342, t344, t345, t348, t351, t354, t359, t360, t361;
-fp t363, t366, t370, t371, t374, t375, t378, t381, t384, t385;
-fp t387, t390, t396, t398, t399, t403, t404, t406, t416, t417;
-fp t419, t423, t426, t429, t432, t440, t443, t445, t453, t454;
-fp t458, t460, t462, t464, t465, t467, t468, t470, t473, t479;
-fp t484, t485, t487, t488, t491, t494, t495, t497, t499, t503;
-fp t510, t514, t518, t521, t522, t537, t538, t540, t557, t560;
-fp t561, t566, t570, t573, t579, t581, t582, t589, t596, t606;
-fp t611, t619, t620, t622, t626, t630, t631, t636, t645, t647;
-fp t649, t650, t654, t659, t660, t662, t663, t672, t675, t676;
-fp t685, t686, t688, t697, t715, t719, t722, t725, t730, t747;
-fp t750, t751, t752, t758, t764, t777, t779, t786, t791, t793;
-fp t798, t802, t803, t805, t806, t808, t809, t811, t813, t815;
-fp t817, t819, t821, t823, t825, t828, t830, t833, t835, t838;
-fp t841, t843, t846, t849, t853, t866, t875, t877, t880, t883;
-fp t884, t886, t889, t895, t897, t900, t902, t905, t907, t909;
-fp t911, t912, t914, t924, t930, t936, t946, t948, t950, t954;
-fp t956, t957, t960, t963, t966, t971, t974, t978, t980, t985;
-fp t987, t988, t993, t994, t999, t1005, t1017, t1022, t1028, t1038;
-fp t1043, t1044, t1048, t1051, t1057, t1066, t1068, t1073, t1081, t1088;
-fp t1092, t1094, t1104, t1107, t1116, t1120, t1122, t1129, t1153, t1163;
-fp t1165, t1187, t1209, t1226, t1252, t1258, t1288, t1292, t1302, t1307;
-fp t1310, t1315, t1327, t1344, t1359, t1361, t1368, t1373, t1377, t1383;
-fp t1394, t1409, t1418, t1431, t1461, t1464, t1465, t1466, t1472, t1474;
-fp t1477, t1478, t1479, t1483, t1492, t1493, t1494, t1497, t1500, t1505;
-fp t1507, t1513, t1518, t1526, t1528, t1536, t1540, t1542, t1546, t1548;
-fp t1550, t1552, t1556, t1559, t1562, t1582, t1583, t1586, t1590, t1594;
-fp t1595, t1598, t1607, t1612, t1615, t1618, t1623, t1626, t1644, t1656;
-fp t1657, t1677, t1679, t1686, t1695, t1714, t1730, t1742, t1745, t1746;
-fp t1749, t1768, t1774, t1776, t1784, t1791, t1794, t1797, t1800, t1803;
-fp t1806, t1818;
- t1 = g_uu_23;
- t2 = 1/r;
- t4 = X_ud_12;
- t5 = PARTIAL_RHO(h);
- t7 = X_ud_22;
- t8 = PARTIAL_SIGMA(h);
- t10 = yy*t2-t4*t5-t7*t8;
- t11 = t1*t10;
- t12 = partial_d_g_uu_313;
- t14 = X_ud_13;
- t16 = X_ud_23;
- t18 = zz*t2-t14*t5-t16*t8;
- t19 = t12*t18;
- t20 = X_ud_11;
- t21 = t19*t20;
- t23 = g_uu_13;
- t24 = t18*t18;
- t25 = t23*t24;
- t26 = partial_d_g_uu_113;
- t29 = t1*t14;
- t30 = g_uu_12;
- t33 = X_ud_21;
- t35 = xx*t2-t20*t5-t33*t8;
- t36 = t30*t35;
- t37 = xx*xx;
- t38 = zz*zz;
- t40 = r*r;
- t42 = 1/t40/r;
- t44 = X_udd_122;
- t46 = X_udd_222;
- t48 = t4*t4;
- t49 = PARTIAL_RHO_RHO(h);
- t51 = t7*t4;
- t52 = PARTIAL_RHO_SIGMA(h);
- t55 = t7*t7;
- t56 = PARTIAL_SIGMA_SIGMA(h);
- t58 = (t37+t38)*t42-t44*t5-t46*t8-t48*t49-2.0*t51*t52-t55*t56;
- t59 = t36*t58;
- t62 = t23*t18;
- t63 = t1*t4;
- t66 = X_udd_113;
- t68 = X_udd_213;
- t72 = t33*t14;
- t74 = t20*t16;
- t78 = -xx*zz*t42-t66*t5-t68*t8-t20*t14*t49-t72*t52-t74*t52-t33*t16*t56;
- t82 = partial_d_g_uu_333;
- t83 = t82*t24;
- t86 = t1*t18;
- t87 = t30*t20;
- t88 = t87*t58;
- t91 = g_uu_11;
- t92 = t91*t35;
- t93 = partial_d_g_uu_122;
- t94 = t93*t10;
- t95 = t94*t4;
- t97 = t35*t35;
- t98 = t23*t97;
- t101 = t23*t35;
- t104 = t23*t20;
- t107 = t91*t20;
- t108 = partial_d_g_uu_133;
- t109 = t108*t24;
- t112 = t1*t1;
- t113 = t112*t24;
- t115 = t30*t30;
- t119 = X_udd_112;
- t121 = X_udd_212;
- t125 = t33*t4;
- t127 = t20*t7;
- t131 = -xx*yy*t42-t119*t5-t121*t8-t20*t4*t49-t125*t52-t127*t52-t33*t7*t56
-;
- t132 = t35*t131;
- t135 = t10*t10;
- t136 = t1*t135;
- t137 = g_uu_22;
- t138 = X_udd_123;
- t142 = t11*t21+t25*t26*t20+2.0*t29*t59+2.0*t62*t63*t78+RATIONAL(1.0,2.0)*
-t63*t83+2.0*t86*t88+t92*t95+t98*t12*t14+2.0*t101*t21+RATIONAL(1.0,2.0)*t104*t83
-+RATIONAL(1.0,2.0)*t107*t109+t113*t44+2.0*t115*t4*t132+2.0*t136*t137*t138;
- t143 = t30*t4;
- t146 = t137*t135;
- t147 = t30*t119;
- t150 = g_uu_33;
- t151 = t150*t18;
- t152 = partial_d_g_uu_323;
- t154 = t152*t10*t14;
- t157 = yy*yy;
- t160 = X_udd_133;
- t162 = X_udd_233;
- t164 = t14*t14;
- t166 = t16*t14;
- t169 = t16*t16;
- t171 = (t37+t157)*t42-t160*t5-t162*t8-t164*t49-2.0*t166*t52-t169*t56;
- t175 = t115*t10;
- t182 = t23*t131;
- t186 = partial_d_g_uu_223;
- t188 = t186*t10*t14;
- t194 = t23*t23;
- t196 = t35*t78;
- t199 = t194*t18;
- t206 = t150*t14;
- t207 = t101*t171;
- t210 = t101*t160;
- t213 = RATIONAL(1.0,2.0)*t143*t109+2.0*t146*t147+2.0*t151*t154+2.0*t151*
-t63*t171+2.0*t175*t20*t131+RATIONAL(3.0,2.0)*t25*t108*t14+4.0*t86*t182*t14+2.0*
-t86*t188+2.0*t175*t35*t119+2.0*t194*t14*t196+2.0*t199*t20*t78+2.0*t199*t35*t66+
-2.0*t206*t207+2.0*t151*t210;
- t215 = t26*t35;
- t216 = t215*t14;
- t219 = t91*t97;
- t220 = partial_d_g_uu_111;
- t221 = t220*t20;
- t224 = partial_d_g_uu_123;
- t226 = t224*t10*t14;
- t229 = t104*t171;
- t235 = t137*t10;
- t236 = t107*t131;
- t239 = partial_d_g_uu_312;
- t240 = t239*t4;
- t242 = partial_d_g_uu_213;
- t243 = t242*t35;
- t244 = t243*t14;
- t247 = t150*t150;
- t248 = t247*t24;
- t255 = X_udd_223;
- t259 = t7*t14;
- t261 = t4*t16;
- t265 = -yy*zz*t42-t138*t5-t255*t8-t4*t14*t49-t259*t52-t261*t52-t7*t16*t56
-;
- t266 = t1*t265;
- t270 = partial_d_g_uu_212;
- t271 = t270*t20;
- t273 = partial_d_g_uu_322;
- t274 = t273*t10;
- t275 = t274*t4;
- t277 = t1*t24;
- t280 = 2.0*t62*t216+RATIONAL(3.0,2.0)*t219*t221+2.0*t62*t226+2.0*t151*
-t229+2.0*t62*t11*t66+2.0*t235*t236+t98*t240+2.0*t86*t244+t248*t160+t136*t239*
-t20+4.0*t151*t266*t14+t146*t271+t101*t275+t277*t186*t4;
- t281 = t137*t4;
- t282 = partial_d_g_uu_211;
- t283 = t282*t97;
- t286 = t36*t138;
- t289 = t30*t10;
- t290 = partial_d_g_uu_112;
- t292 = t290*t35*t4;
- t295 = partial_d_g_uu_311;
- t296 = t295*t97;
- t299 = t23*t78;
- t303 = t92*t119;
- t309 = t247*t18;
- t313 = t87*t265;
- t319 = t150*t24;
- t320 = t82*t14;
- t323 = t30*t135;
- t327 = t112*t18;
- t333 = X_udd_111;
- t335 = X_udd_211;
- t337 = t20*t20;
- t339 = t33*t20;
- t342 = t33*t33;
- t344 = (t157+t38)*t42-t333*t5-t335*t8-t337*t49-2.0*t339*t52-t342*t56;
- t345 = t107*t344;
- t348 = RATIONAL(1.0,2.0)*t281*t283+2.0*t11*t286+2.0*t289*t292+RATIONAL(
-1.0,2.0)*t63*t296+4.0*t151*t299*t14+2.0*t235*t303+RATIONAL(3.0,2.0)*t136*t273*
-t4+2.0*t309*t171*t14+2.0*t11*t313+2.0*t98*t30*t138+RATIONAL(3.0,2.0)*t319*t320+
-RATIONAL(3.0,2.0)*t323*t93*t4+2.0*t327*t58*t14+2.0*t289*t345;
- t351 = t273*t135;
- t354 = t92*t333;
- t359 = partial_d_g_uu_222;
- t360 = t359*t10;
- t361 = t360*t4;
- t363 = t359*t4;
- t366 = partial_d_g_uu_233;
- t370 = t23*t14;
- t371 = t11*t78;
- t374 = t152*t18;
- t375 = t374*t4;
- t378 = t92*t344;
- t381 = t92*t78;
- t384 = t242*t18;
- t385 = t384*t35;
- t387 = t36*t265;
- t390 = t30*t97;
- t396 = RATIONAL(1.0,2.0)*t206*t351+2.0*t289*t354+RATIONAL(1.0,2.0)*t104*
-t351+t36*t361+RATIONAL(3.0,2.0)*t146*t363+RATIONAL(3.0,2.0)*t277*t366*t14+2.0*
-t370*t371+2.0*t11*t375+2.0*t143*t378+2.0*t63*t381+t281*t385+2.0*t63*t387+
-RATIONAL(3.0,2.0)*t390*t282*t20+t323*t290*t20;
- t398 = t186*t18;
- t399 = t398*t4;
- t403 = t270*t10;
- t404 = t403*t35;
- t406 = t107*t78;
- t416 = t290*t10;
- t417 = t416*t20;
- t419 = t186*t14;
- t423 = t1*t138;
- t426 = t23*t66;
- t429 = t220*t97;
- t432 = t62*t292+2.0*t235*t399+t235*t244+t29*t404+2.0*t11*t406+2.0*t390*
-t91*t119+t289*t216+2.0*t289*t104*t78+t62*t417+t146*t419+t136*t152*t14+2.0*t319*
-t423+2.0*t319*t426+RATIONAL(1.0,2.0)*t143*t429;
- t440 = t36*t131;
- t443 = t416*t35;
- t445 = t36*t44;
- t453 = t270*t4;
- t454 = t453*t35;
- t458 = t194*t97;
- t460 = t115*t135;
- t462 = t115*t97;
- t464 = t91*t91;
- t465 = t464*t97;
- t467 = 2.0*t277*t23*t119+RATIONAL(3.0,2.0)*t98*t295*t20+2.0*t370*t440+
-t370*t443+2.0*t86*t445+2.0*t281*t59+RATIONAL(1.0,2.0)*t206*t296+t92*t226+t86*
-t454+2.0*t235*t445+t458*t160+t460*t333+t462*t44+t465*t333;
- t468 = t384*t20;
- t470 = t92*t66;
- t473 = t101*t78;
- t479 = t403*t20;
- t484 = t282*t35;
- t485 = t484*t20;
- t487 = t220*t35;
- t488 = t487*t20;
- t491 = t12*t35*t14;
- t494 = t239*t10;
- t495 = t494*t35;
- t497 = t494*t20;
- t499 = t101*t265;
- t503 = t235*t468+2.0*t11*t470+2.0*t143*t473+2.0*t289*t101*t66+t86*t479+
-2.0*t235*t88+t86*t361+t235*t485+t289*t488+2.0*t151*t491+t206*t495+t151*t497+2.0
-*t281*t499+t62*t488;
- t510 = t101*t138;
- t514 = t112*t10;
- t518 = t104*t265;
- t521 = t295*t35;
- t522 = t521*t20;
- t537 = t62*t95+4.0*t36*t236+2.0*t235*t510+t86*t485+2.0*t514*t171*t4+2.0*
-t235*t518+t11*t522+2.0*t86*t303+2.0*t63*t207+2.0*t36*t479+2.0*t11*t210+2.0*t11*
-t229+2.0*t92*t417+t151*t522;
- t538 = t240*t35;
- t540 = t289*t131;
- t557 = t289*t78;
- t560 = t137*t265;
- t561 = t560*t4;
- t566 = t464*t35;
- t570 = t93*t135;
- t573 = t359*t135;
- t579 = t151*t538+2.0*t29*t540+4.0*t101*t406+2.0*t136*t30*t66+2.0*t98*t91*
-t66+2.0*t62*t36*t119+2.0*t62*t87*t131+2.0*t206*t557+4.0*t11*t561+t323*t224*t14+
-2.0*t566*t344*t20+RATIONAL(1.0,2.0)*t370*t570+RATIONAL(1.0,2.0)*t29*t573+2.0*
-t175*t344*t4;
- t581 = t30*t78;
- t582 = t581*t4;
- t589 = t143*t131;
- t596 = t374*t10;
- t606 = t115*t35;
- t611 = t10*t265;
- t619 = 4.0*t11*t582+2.0*t86*t289*t119+t151*t275+2.0*t86*t589+4.0*t101*
-t313+4.0*t235*t589+t104*t596+2.0*t62*t289*t333+2.0*t151*t289*t66+2.0*t151*t582+
-2.0*t606*t58*t20+2.0*t112*t14*t611+2.0*t327*t4*t265+2.0*t101*t497;
- t620 = t26*t14;
- t622 = t92*t131;
- t626 = t289*t344;
- t630 = t224*t18;
- t631 = t630*t4;
- t636 = t194*t35;
- t645 = t630*t10;
- t647 = t194*t24;
- t649 = t26*t18;
- t650 = t649*t20;
- t654 = t219*t620+2.0*t29*t622+t101*t375+2.0*t370*t626+t101*t154+t92*t631+
-2.0*t327*t10*t138+2.0*t636*t171*t20+2.0*t62*t143*t344+2.0*t235*t454+t107*t645+
-t647*t333+t289*t650+2.0*t86*t236;
- t659 = t108*t18;
- t660 = t659*t14;
- t662 = t366*t18;
- t663 = t662*t14;
- t672 = t398*t10;
- t675 = t82*t18;
- t676 = t675*t14;
- t685 = 2.0*t289*t631+t92*t660+t235*t663+2.0*t281*t622+2.0*t11*t538+t36*
-t188+RATIONAL(1.0,2.0)*t370*t429+t87*t672+t36*t399+t101*t676+t289*t660+t11*t676
-+2.0*t62*t354+2.0*t62*t281*t131;
- t686 = t290*t4;
- t688 = t19*t35;
- t697 = t235*t58;
- t715 = t219*t686+t63*t688+t36*t663+2.0*t370*t378+RATIONAL(1.0,2.0)*t29*
-t283+2.0*t206*t381+2.0*t29*t697+2.0*t62*t345+2.0*t86*t281*t58+2.0*t36*t468+2.0*
-t86*t235*t44+t390*t453+2.0*t151*t313+2.0*t151*t561;
- t719 = t366*t24;
- t722 = t235*t131;
- t725 = t12*t20;
- t730 = t235*t265;
- t747 = t152*t4;
- t750 = 2.0*t151*t470+RATIONAL(1.0,2.0)*t281*t719+2.0*t370*t722+t319*t725+
-2.0*t62*t235*t119+2.0*t206*t730+2.0*t206*t387+2.0*t151*t11*t160+2.0*t92*t650+
-t390*t242*t14+2.0*t199*t344*t14+2.0*t151*t286+t319*t747+t11*t491;
- t751 = t137*t137;
- t752 = t751*t10;
- t758 = t649*t35;
- t764 = t11*t171;
- t777 = t112*t135;
- t779 = t751*t135;
- t786 = 2.0*t752*t58*t4+RATIONAL(1.0,2.0)*t87*t719+t143*t758+2.0*t86*t518+
-t277*t242*t20+2.0*t206*t764+2.0*t151*t406+t25*t224*t4+2.0*t29*t499+2.0*t86*t510
-+RATIONAL(1.0,2.0)*t107*t570+t777*t160+t779*t44+2.0*t151*t235*t138+RATIONAL(1.0
-,2.0)*t87*t573;
- t791 = pow(HD,1.0*RATIONAL(1.0,2.0));
- t793 = 1/t791/HD;
- t798 = -t221-t271-t725-t686-t363-t747-t620-t419-t320-t91*t333-2.0*t147
--2.0*t426;
- t802 = partial_d_ln_sqrt_g_1;
- t803 = t91*t802;
- t805 = partial_d_ln_sqrt_g_2;
- t806 = t30*t805;
- t808 = partial_d_ln_sqrt_g_3;
- t809 = t23*t808;
- t811 = t30*t802;
- t813 = t137*t805;
- t815 = t1*t808;
- t817 = t23*t802;
- t819 = t1*t805;
- t821 = t150*t808;
- t823 = -t137*t44-2.0*t423-t150*t160-t803*t20-t806*t20-t809*t20-t811*t4-
-t813*t4-t815*t4-t817*t14-t819*t14-t821*t14;
- t825 = 1/t791;
- t828 = K_uu_11*t35;
- t830 = K_uu_12;
- t833 = t830*t10;
- t835 = K_uu_13;
- t838 = t835*t18;
- t841 = K_uu_22*t10;
- t843 = K_uu_23;
- t846 = t843*t18;
- t849 = K_uu_33*t18;
- t853 = 1/HD;
- t866 = HD*HD;
- t875 = RATIONAL(3.0,2.0)*HA/t791/t866+RATIONAL(1.0,2.0)*HB*t793+HC/t866;
- partial_H_wrt_partial_d_h_1 = (t142+t213+t280+t348+t396+t432+t467+t503+
-t537+t579+t619+t654+t685+t715+t750+t786)*t793+(t798+t823)*t825+(-2.0*t828*t20
--2.0*t830*t4*t35-2.0*t833*t20-2.0*t835*t14*t35-2.0*t838*t20-2.0*t841*t4-2.0*
-t843*t14*t10-2.0*t846*t4-2.0*t849*t14)*t853-(-2.0*t92*t20-2.0*t143*t35-2.0*t289
-*t20-2.0*t370*t35-2.0*t62*t20-2.0*t235*t4-2.0*t29*t10-2.0*t86*t4-2.0*t151*t14)*
-t875;
- t877 = t215*t16;
- t880 = t23*t16;
- t883 = t290*t7;
- t884 = t883*t35;
- t886 = t220*t33;
- t889 = t30*t33;
- t895 = t662*t16;
- t897 = t243*t16;
- t900 = t675*t16;
- t902 = t150*t16;
- t905 = t270*t33;
- t907 = t239*t7;
- t909 = t659*t16;
- t911 = t224*t16;
- t912 = t911*t10;
- t914 = 2.0*t62*t877+2.0*t880*t440+t62*t884+RATIONAL(3.0,2.0)*t219*t886+
-RATIONAL(1.0,2.0)*t889*t573+2.0*t327*t58*t16+t36*t895+2.0*t86*t897+t11*t900+2.0
-*t902*t207+t146*t905+t98*t907+t92*t909+t92*t912;
- t924 = t82*t16;
- t930 = t359*t7;
- t936 = t23*t33;
- t946 = t416*t33;
- t948 = t26*t16;
- t950 = 2.0*t289*t884+t101*t900+RATIONAL(3.0,2.0)*t25*t108*t16+RATIONAL(
-3.0,2.0)*t277*t366*t16+RATIONAL(3.0,2.0)*t319*t924+RATIONAL(3.0,2.0)*t323*t93*
-t7+RATIONAL(3.0,2.0)*t146*t930+RATIONAL(3.0,2.0)*t136*t273*t7+RATIONAL(1.0,2.0)
-*t936*t351+2.0*t880*t378+2.0*t309*t171*t16+t390*t242*t16+t62*t946+t219*t948;
- t954 = t630*t7;
- t956 = t91*t33;
- t957 = t956*t344;
- t960 = t936*t265;
- t963 = t36*t255;
- t966 = t889*t265;
- t971 = t92*t121;
- t974 = t1*t16;
- t978 = t487*t33;
- t980 = t403*t33;
- t985 = t1*t7;
- t987 = t136*t239*t33+t92*t954+2.0*t62*t957+2.0*t86*t960+2.0*t151*t963+2.0
-*t151*t966+2.0*t92*t946+2.0*t235*t971+2.0*t974*t499+t974*t404+t62*t978+2.0*t36*
-t980+2.0*t289*t957+t985*t688;
- t988 = t30*t7;
- t993 = t12*t16;
- t994 = t993*t35;
- t999 = t92*t335;
- t1005 = t956*t131;
- t1017 = t936*t171;
- t1022 = 2.0*t988*t378+2.0*t11*t966+t11*t994+2.0*t985*t387+t880*t443+2.0*
-t289*t999+2.0*t62*t889*t131+2.0*t235*t1005+2.0*t902*t387+2.0*t11*t963+2.0*t151*
-t985*t171+2.0*t86*t971+2.0*t151*t1017+2.0*t86*t1005;
- t1028 = t137*t7;
- t1038 = t956*t78;
- t1043 = t186*t16;
- t1044 = t1043*t10;
- t1048 = t36*t46;
- t1051 = t889*t58;
- t1057 = 2.0*t902*t381+t219*t883+2.0*t1028*t499+t86*t980+2.0*t151*t11*t162
-+2.0*t289*t101*t68+2.0*t151*t1038+2.0*t985*t381+t36*t1044+2.0*t1028*t622+2.0*
-t235*t1048+2.0*t235*t1051+2.0*t62*t912+t889*t672;
- t1066 = t398*t7;
- t1068 = t92*t68;
- t1073 = t101*t162;
- t1081 = t581*t7;
- t1088 = 2.0*t289*t954+t289*t909+2.0*t988*t473+2.0*t327*t10*t255+t36*t1066
-+2.0*t11*t1068+2.0*t11*t1038+2.0*t11*t1073+2.0*t11*t1017+2.0*t985*t207+t235*
-t895+4.0*t11*t1081+2.0*t151*t1068+2.0*t974*t540;
- t1092 = t494*t33;
- t1094 = t560*t7;
- t1104 = t988*t131;
- t1107 = t101*t255;
- t1116 = t384*t33;
- t1120 = t649*t33;
- t1122 = 2.0*t1028*t59+t151*t1092+4.0*t11*t1094+2.0*t235*t960+2.0*t289*
-t936*t78+2.0*t151*t1073+2.0*t86*t1104+2.0*t235*t1107+2.0*t86*t289*t121+4.0*t151
-*t299*t16+t235*t1116+t235*t897+t1028*t385+t289*t1120;
- t1129 = t484*t33;
- t1153 = 2.0*t151*t994+RATIONAL(1.0,2.0)*t956*t109+2.0*t880*t626+t86*t1129
-+t235*t1129+2.0*t902*t557+2.0*t752*t58*t7+RATIONAL(1.0,2.0)*t936*t83+RATIONAL(
-1.0,2.0)*t889*t719+2.0*t62*t289*t335+2.0*t62*t988*t344+t289*t978+2.0*t62*t999+
-2.0*t62*t11*t68;
- t1163 = t521*t33;
- t1165 = t907*t35;
- t1187 = 2.0*t136*t30*t68+4.0*t151*t266*t16+t11*t1163+t151*t1165+RATIONAL(
-1.0,2.0)*t974*t573+RATIONAL(1.0,2.0)*t880*t570+t151*t1163+2.0*t86*t1051+2.0*
-t606*t58*t33+2.0*t974*t59+t956*t645+2.0*t902*t730+2.0*t390*t91*t121+2.0*t86*
-t1048;
- t1209 = t30*t121;
- t1226 = 2.0*t62*t985*t78+2.0*t194*t16*t196+2.0*t199*t33*t78+2.0*t199*t35*
-t68+2.0*t112*t16*t611+2.0*t327*t7*t265+2.0*t98*t30*t255+2.0*t146*t1209+RATIONAL
-(1.0,2.0)*t902*t351+RATIONAL(1.0,2.0)*t988*t109+RATIONAL(1.0,2.0)*t1028*t719+
-2.0*t974*t697+2.0*t98*t91*t68+2.0*t151*t289*t68;
- t1252 = t152*t7;
- t1258 = 2.0*t62*t235*t121+2.0*t566*t344*t33+2.0*t880*t722+t777*t162+t146*
-t1043+2.0*t636*t171*t33+2.0*t86*t235*t46+2.0*t86*t1028*t58+t779*t46+RATIONAL(
-1.0,2.0)*t985*t83+RATIONAL(1.0,2.0)*t956*t570+t319*t1252+t277*t186*t7+t25*t224*
-t7;
- t1288 = t152*t16;
- t1292 = 2.0*t136*t137*t255+t323*t911+2.0*t151*t1081+2.0*t151*t235*t255+
-4.0*t86*t182*t16+2.0*t880*t371+t458*t162+RATIONAL(1.0,2.0)*t880*t429+2.0*t62*
-t1028*t131+2.0*t115*t7*t132+2.0*t175*t33*t131+2.0*t175*t35*t121+t136*t1288+
-RATIONAL(1.0,2.0)*t974*t283;
- t1302 = t12*t33;
- t1307 = t23*t68;
- t1310 = t1*t255;
- t1315 = t374*t7;
- t1327 = RATIONAL(1.0,2.0)*t902*t296+t390*t270*t7+2.0*t199*t344*t16+t319*
-t1302+2.0*t277*t23*t121+2.0*t319*t1307+2.0*t319*t1310+RATIONAL(1.0,2.0)*t988*
-t429+2.0*t11*t1315+2.0*t902*t764+RATIONAL(1.0,2.0)*t1028*t283+RATIONAL(1.0,2.0)
-*t985*t296+t98*t993+t25*t26*t33;
- t1344 = t19*t33;
- t1359 = t274*t7;
- t1361 = t277*t242*t33+t323*t290*t33+2.0*t514*t171*t7+2.0*t175*t344*t7+
-RATIONAL(3.0,2.0)*t98*t295*t33+RATIONAL(3.0,2.0)*t390*t282*t33+t11*t1344+2.0*
-t974*t622+2.0*t62*t36*t121+2.0*t86*t1107+2.0*t151*t1094+2.0*t235*t1066+4.0*t101
-*t1038+t151*t1359;
- t1368 = t1288*t10;
- t1373 = t360*t7;
- t1377 = t94*t7;
- t1383 = t936*t596+t289*t877+t460*t335+t462*t46+t465*t335+t101*t1368+t101*
-t1315+2.0*t101*t1344+t86*t1373+4.0*t36*t1005+t62*t1377+t902*t495+t988*t758+4.0*
-t101*t966;
- t1394 = t270*t35*t7;
- t1409 = 4.0*t235*t1104+t92*t1377+2.0*t151*t1368+2.0*t86*t1044+2.0*t101*
-t1092+2.0*t235*t1394+2.0*t11*t1165+2.0*t92*t1120+2.0*t36*t1116+t36*t1373+t101*
-t1359+t248*t162+t113*t46+t647*t335+t86*t1394;
- t1418 = -t886-t905-t1302-t883-t930-t1252-t948-t1043-t924-t91*t335-2.0*
-t1209-2.0*t1307;
- t1431 = -t137*t46-2.0*t1310-t150*t162-t803*t33-t806*t33-t809*t33-t811*t7-
-t813*t7-t815*t7-t817*t16-t819*t16-t821*t16;
- partial_H_wrt_partial_d_h_2 = (t914+t950+t987+t1022+t1057+t1088+t1122+
-t1153+t1187+t1226+t1258+t1292+t1327+t1361+t1383+t1409)*t793+(t1418+t1431)*t825+
-(-2.0*t828*t33-2.0*t830*t7*t35-2.0*t833*t33-2.0*t835*t16*t35-2.0*t838*t33-2.0*
-t841*t7-2.0*t843*t16*t10-2.0*t846*t7-2.0*t849*t16)*t853-(-2.0*t92*t33-2.0*t988*
-t35-2.0*t289*t33-2.0*t880*t35-2.0*t62*t33-2.0*t235*t7-2.0*t974*t10-2.0*t86*t7
--2.0*t151*t16)*t875;
- t1461 = t36*t48;
- t1464 = t11*t91;
- t1465 = t35*t20;
- t1466 = t1465*t14;
- t1472 = t11*t30;
- t1474 = t35*t4*t14;
- t1477 = t151*t30;
- t1478 = t10*t20;
- t1479 = t1478*t14;
- t1483 = t101*t164;
- t1492 = 2.0*t235*t1461+2.0*t1464*t1466+2.0*t62*t289*t337+2.0*t1472*t1474+
-2.0*t1477*t1479+t779*t48+2.0*t11*t1483+t462*t48+t460*t337+2.0*t277*t104*t4+t458
-*t164;
- t1493 = t86*t30;
- t1494 = t1478*t4;
- t1497 = t235*t23;
- t1500 = t63*t14;
- t1505 = t151*t137;
- t1507 = t10*t4*t14;
- t1513 = t104*t14;
- t1518 = 2.0*t1493*t1494+2.0*t1497*t1474+2.0*t319*t1500+t465*t337+t647*
-t337+2.0*t1505*t1507+t248*t164+t113*t48+t777*t164+2.0*t319*t1513+2.0*t199*t1466
-;
- t1526 = t289*t23;
- t1528 = t62*t137;
- t1536 = t62*t1;
- t1540 = t136*t87*t14+t86*t235*t48+t136*t281*t14+t1526*t1466+t1528*t1494+
-t151*t11*t164+t86*t1461+t327*t1507+t98*t107*t14+t1536*t1479+t390*t107*t4;
- t1542 = t87*t4;
- t1546 = t1465*t4;
- t1548 = t92*t337;
- t1550 = t62*t30;
- t1552 = t86*t91;
- t1556 = t151*t91;
- t1559 = t235*t91;
- t1562 = t146*t1542+t98*t143*t14+t175*t1546+t62*t1548+t1550*t1546+t1552*
-t1546+t151*t1483+t289*t1548+t1556*t1466+t1536*t1474+t1559*t1546+t1477*t1474;
- partial_H_wrt_partial_dd_h_11 = (t1492+t1518+2.0*t1540+2.0*t1562)*t793+(-
-t91*t337-2.0*t1542-2.0*t1513-t137*t48-2.0*t1500-t150*t164)*t825;
- t1582 = -t72-t74;
- t1583 = t92*t1582;
- t1586 = -t125-t127;
- t1590 = t92*t1586;
- t1594 = t35*t33;
- t1595 = t1594*t20;
- t1598 = -t259-t261;
- t1607 = t30*t1586;
- t1612 = 2.0*t777*t166+4.0*t151*t1*t10*t16*t14-2.0*t151*t1583-2.0*t175*t35
-*t1586-2.0*t86*t1590+4.0*t62*t91*t1595-2.0*t98*t30*t1598-2.0*t390*t91*t1586-2.0
-*t235*t1590-2.0*t146*t1607+2.0*t113*t51;
- t1615 = t35*t16*t14;
- t1618 = t36*t1598;
- t1623 = t23*t1582;
- t1626 = t101*t1598;
- t1644 = 4.0*t151*t23*t1615-2.0*t151*t1618+2.0*t460*t339-2.0*t319*t1623
--2.0*t86*t1626+2.0*t462*t51+2.0*t458*t166-2.0*t98*t91*t1582+4.0*t289*t91*t1595
--2.0*t11*t1583-2.0*t199*t35*t1582;
- t1656 = t35*t7;
- t1657 = t1656*t4;
- t1677 = -2.0*t11*t1618+2.0*t647*t339-2.0*t136*t137*t1598-2.0*t235*t1626+
-4.0*t235*t30*t1657+4.0*t11*t23*t1615-2.0*t136*t30*t1582-2.0*t62*t235*t1586-2.0*
-t277*t23*t1586-2.0*t151*t289*t1582+2.0*t465*t339;
- t1679 = t10*t7;
- t1686 = t10*t33;
- t1695 = t1*t1598;
- t1714 = 4.0*t86*t137*t1679*t4-2.0*t86*t289*t1586+4.0*t1550*t1686*t20+2.0*
-t248*t166-2.0*t151*t235*t1598-2.0*t319*t1695-2.0*t327*t10*t1598-2.0*t289*t101*
-t1582-2.0*t62*t11*t1582+4.0*t1493*t1657+2.0*t779*t51-2.0*t62*t36*t1586;
- partial_H_wrt_partial_dd_h_12 = (t1612+t1644+t1677+t1714)*t793+(-2.0*t956
-*t20+2.0*t1607+2.0*t1623-2.0*t1028*t4+2.0*t1695-2.0*t902*t14)*t825;
- t1730 = t1686*t7;
- t1742 = t1656*t16;
- t1745 = t458*t169+t248*t169+t465*t342+t113*t55+t462*t55+t460*t342+2.0*
-t1528*t1730+2.0*t62*t289*t342+2.0*t136*t1028*t16+2.0*t136*t889*t16+2.0*t1497*
-t1742;
- t1746 = t889*t7;
- t1749 = t1594*t16;
- t1768 = t101*t169;
- t1774 = 2.0*t146*t1746+2.0*t1526*t1749+2.0*t390*t956*t7+2.0*t199*t1749+
-2.0*t1493*t1730+t647*t342+2.0*t98*t988*t16+2.0*t98*t956*t16+2.0*t1464*t1749+2.0
-*t11*t1768+2.0*t277*t936*t7;
- t1776 = t36*t55;
- t1784 = t1594*t7;
- t1791 = t1686*t16;
- t1794 = t1679*t16;
- t1797 = t985*t16;
- t1800 = t92*t342;
- t1803 = 2.0*t86*t1776+2.0*t151*t11*t169+2.0*t1556*t1749+2.0*t175*t1784+
-2.0*t235*t1776+t779*t55+t777*t169+2.0*t1477*t1791+2.0*t327*t1794+2.0*t319*t1797
-+2.0*t62*t1800;
- t1806 = t936*t16;
- t1818 = t86*t235*t55+t319*t1806+t1505*t1794+t1550*t1784+t1536*t1791+t151*
-t1768+t1559*t1784+t289*t1800+t1472*t1742+t1477*t1742+t1552*t1784+t1536*t1742;
- partial_H_wrt_partial_dd_h_22 = (t1745+t1774+t1803+2.0*t1818)*t793+(-t91*
-t342-2.0*t1746-2.0*t1806-t137*t55-2.0*t1797-t150*t169)*t825;
diff --git a/src/gr.cg/horizon_function.c b/src/gr.cg/horizon_function.c
deleted file mode 100644
index 7e6e717..0000000
--- a/src/gr.cg/horizon_function.c
+++ /dev/null
@@ -1,182 +0,0 @@
-/*
- * inputs = {r, partial_d_ln_sqrt_g, partial_d_g_uu, X_ud, X_udd, g_uu, K_uu, h}
- * outputs = {HA, HB, HC, HD}
- * cost = 134*assignments+3*divisions+5*functions+401*multiplications+173*additions
- */
-fp t1, t2, t4, t5, t7, t8, t10, t11, t12, t14;
-fp t16, t18, t19, t20, t23, t24, t25, t26, t29, t30;
-fp t31, t33, t35, t37, t38, t40, t41, t42, t43, t45;
-fp t46, t47, t48, t50, t51, t52, t53, t55, t58, t59;
-fp t61, t64, t67, t68, t70, t72, t73, t75, t77, t83;
-fp t84, t87, t90, t91, t93, t95, t112, t130, t134, t137;
-fp t139, t143, t147, t150, t152, t159, t164, t166, t168, t176;
-fp t181, t183, t199, t203, t206, t209, t213, t216, t217, t219;
-fp t228, t229, t232, t235, t236, t237, t240, t243, t247, t251;
-fp t254, t259, t260, t262, t263, t269, t272, t273, t276, t279;
-fp t282, t283, t290, t293, t296, t299, t302, t305, t312, t315;
-fp t318, t324, t325, t329, t330, t338, t358, t381, t392, t395;
-fp t396, t397, t408, t411, t431, t440, t444, t447, t450, t465;
- t1 = g_uu_33;
- t2 = 1/r;
- t4 = X_ud_13;
- t5 = PARTIAL_RHO(h);
- t7 = X_ud_23;
- t8 = PARTIAL_SIGMA(h);
- t10 = zz*t2-t4*t5-t7*t8;
- t11 = t1*t10;
- t12 = partial_d_g_uu_311;
- t14 = X_ud_11;
- t16 = X_ud_21;
- t18 = xx*t2-t14*t5-t16*t8;
- t19 = t18*t18;
- t20 = t12*t19;
- t23 = g_uu_23;
- t24 = t23*t10;
- t25 = partial_d_g_uu_211;
- t26 = t25*t19;
- t29 = g_uu_13;
- t30 = t29*t19;
- t31 = partial_d_g_uu_312;
- t33 = X_ud_12;
- t35 = X_ud_22;
- t37 = yy*t2-t33*t5-t35*t8;
- t38 = t31*t37;
- t40 = g_uu_12;
- t41 = t40*t19;
- t42 = partial_d_g_uu_212;
- t43 = t42*t37;
- t45 = g_uu_11;
- t46 = t45*t19;
- t47 = partial_d_g_uu_112;
- t48 = t47*t37;
- t50 = g_uu_22;
- t51 = t37*t37;
- t52 = t50*t51;
- t53 = t42*t18;
- t55 = t40*t51;
- t58 = partial_d_g_uu_113;
- t59 = t58*t10;
- t61 = t23*t51;
- t64 = t29*t10;
- t67 = partial_d_g_uu_213;
- t68 = t67*t10;
- t70 = t45*t45;
- t72 = yy*yy;
- t73 = zz*zz;
- t75 = r*r;
- t77 = 1/t75/r;
- t83 = t14*t14;
- t84 = PARTIAL_RHO_RHO(h);
- t87 = PARTIAL_RHO_SIGMA(h);
- t90 = t16*t16;
- t91 = PARTIAL_SIGMA_SIGMA(h);
- t93 = (t72+t73)*t77-X_udd_111*t5-X_udd_211*t8-t83*t84-2.0*t16*t14*t87-t90
-*t91;
- t95 = RATIONAL(-1.0,2.0)*t11*t20+RATIONAL(-1.0,2.0)*t24*t26-t30*t38-t41*
-t43-t46*t48-t52*t53-t55*t47*t18-t46*t59-t61*t31*t18-t64*t48*t18-t41*t68-t70*t19
-*t93;
- t112 = -xx*yy*t77-X_udd_112*t5-X_udd_212*t8-t14*t33*t84-t16*t33*t87-t14*
-t35*t87-t16*t35*t91;
- t130 = -xx*zz*t77-X_udd_113*t5-X_udd_213*t8-t14*t4*t84-t16*t4*t87-t14*t7*
-t87-t16*t7*t91;
- t134 = t40*t37;
- t137 = t51*t37;
- t139 = partial_d_g_uu_322;
- t143 = partial_d_g_uu_122;
- t147 = partial_d_g_uu_222;
- t150 = t40*t40;
- t152 = xx*xx;
- t159 = t33*t33;
- t164 = t35*t35;
- t166 = (t152+t73)*t77-X_udd_122*t5-X_udd_222*t8-t159*t84-2.0*t35*t33*t87-
-t164*t91;
- t168 = t29*t29;
- t176 = t4*t4;
- t181 = t7*t7;
- t183 = (t152+t72)*t77-X_udd_133*t5-X_udd_233*t8-t176*t84-2.0*t7*t4*t87-
-t181*t91;
- t199 = -yy*zz*t77-X_udd_123*t5-X_udd_223*t8-t33*t4*t84-t35*t4*t87-t33*t7*
-t87-t35*t7*t91;
- t203 = t40*t112;
- t206 = t50*t37;
- t209 = -t24*t43*t18-2.0*t41*t45*t112-2.0*t30*t45*t130-t134*t59*t18+
-RATIONAL(-1.0,2.0)*t23*t137*t139+RATIONAL(-1.0,2.0)*t40*t137*t143+RATIONAL(-1.0
-,2.0)*t50*t137*t147-t150*t19*t166-t168*t19*t183-2.0*t30*t40*t199-2.0*t52*t203-
-t206*t68*t18;
- t213 = t50*t50;
- t216 = t10*t10;
- t217 = t216*t10;
- t219 = partial_d_g_uu_333;
- t228 = partial_d_g_uu_123;
- t229 = t228*t10;
- t232 = partial_d_g_uu_233;
- t235 = t23*t37;
- t236 = partial_d_g_uu_313;
- t237 = t236*t10;
- t240 = t23*t23;
- t243 = t29*t18;
- t247 = t243*t183;
- t251 = partial_d_g_uu_133;
- t254 = -t150*t51*t93-t213*t51*t166+RATIONAL(-1.0,2.0)*t1*t217*t219-2.0*
-t61*t40*t130-2.0*t61*t50*t199-t55*t229+RATIONAL(-1.0,2.0)*t23*t217*t232-t235*
-t237*t18-t240*t51*t183-2.0*t134*t243*t130-2.0*t235*t247+RATIONAL(-1.0,2.0)*t29*
-t217*t251;
- t259 = partial_d_g_uu_323;
- t260 = t259*t10;
- t262 = partial_d_g_uu_223;
- t263 = t262*t10;
- t269 = t243*t199;
- t272 = t40*t18;
- t273 = t272*t166;
- t276 = t272*t199;
- t279 = t219*t216;
- t282 = t45*t18;
- t283 = t251*t216;
- t290 = t232*t216;
- t293 = t23*t216;
- t296 = -2.0*t150*t37*t18*t112-t61*t260-t52*t263-2.0*t168*t10*t18*t130-2.0
-*t206*t269-2.0*t206*t273-2.0*t235*t276+RATIONAL(-1.0,2.0)*t243*t279+RATIONAL(
--1.0,2.0)*t282*t283+RATIONAL(-1.0,2.0)*t235*t279+RATIONAL(-1.0,2.0)*t134*t283+
-RATIONAL(-1.0,2.0)*t206*t290-t293*t262*t37;
- t299 = t29*t216;
- t302 = t147*t51;
- t305 = t139*t51;
- t312 = t282*t93;
- t315 = t282*t130;
- t318 = t282*t112;
- t324 = t1*t216;
- t325 = t236*t18;
- t329 = -t299*t228*t37+RATIONAL(-1.0,2.0)*t272*t302+RATIONAL(-1.0,2.0)*t11
-*t305+RATIONAL(-1.0,2.0)*t24*t302+RATIONAL(-1.0,2.0)*t272*t290-2.0*t134*t312
--2.0*t235*t315-2.0*t206*t318-t30*t237-2.0*t24*t269-t324*t325-2.0*t11*t247;
- t330 = t259*t37;
- t338 = t143*t51;
- t358 = -t324*t330-2.0*t11*t315-t272*t263*t37-2.0*t11*t276+RATIONAL(-1.0,
-2.0)*t282*t338-t243*t260*t37-2.0*t64*t312-t282*t229*t37-t299*t58*t18-2.0*t24*
-t273-2.0*t24*t318+RATIONAL(-1.0,2.0)*t64*t338-2.0*t64*t235*t130;
- t381 = t29*t130;
- t392 = t23*t199;
- t395 = -2.0*t240*t10*t37*t199-2.0*t11*t235*t183-2.0*t11*t134*t130-2.0*t64
-*t134*t93-2.0*t24*t206*t166-t11*t38*t18-2.0*t293*t29*t112-2.0*t324*t381-2.0*t11
-*t206*t199-2.0*t64*t206*t112-t168*t216*t93-2.0*t324*t392;
- t396 = partial_d_g_uu_111;
- t397 = t396*t19;
- t408 = t1*t1;
- t411 = t19*t18;
- t431 = RATIONAL(-1.0,2.0)*t134*t397+RATIONAL(-1.0,2.0)*t206*t26+RATIONAL(
--1.0,2.0)*t235*t20+RATIONAL(-1.0,2.0)*t64*t397-t240*t216*t166-t408*t216*t183+
-RATIONAL(-1.0,2.0)*t45*t411*t396+RATIONAL(-1.0,2.0)*t40*t411*t25+RATIONAL(-1.0,
-2.0)*t29*t411*t12-2.0*t24*t134*t112+RATIONAL(-1.0,2.0)*t243*t305-t293*t67*t18
--2.0*t64*t272*t112;
- HA = t95+t209+t254+t296+t329+t358+t395+t431;
- t440 = t396*t18+t53+t325+t48+t147*t37+t330+t59+t263+t219*t10+t45*t93+2.0*
-t203+2.0*t381;
- t444 = partial_d_ln_sqrt_g_1;
- t447 = partial_d_ln_sqrt_g_2;
- t450 = partial_d_ln_sqrt_g_3;
- t465 = t50*t166+2.0*t392+t1*t183+t45*t444*t18+t40*t447*t18+t29*t450*t18+
-t40*t444*t37+t50*t447*t37+t23*t450*t37+t29*t444*t10+t23*t447*t10+t1*t450*t10;
- HB = t440+t465;
- HC = K_uu_11*t19+2.0*K_uu_12*t37*t18+2.0*K_uu_13*t10*t18+K_uu_22*t51+2.0*
-K_uu_23*t10*t37+K_uu_33*t216;
- HD = t46+2.0*t134*t18+2.0*t64*t18+t52+2.0*t24*t37+t324;
diff --git a/src/gr/cg.hh b/src/gr/cg.hh
index 02bf370..affd54b 100644
--- a/src/gr/cg.hh
+++ b/src/gr/cg.hh
@@ -89,20 +89,20 @@
#define partial_d_g_dd_323 p.gridfn(gfns::gfn__partial_d_g_dd_323, irho,isigma)
#define partial_d_g_dd_333 p.gridfn(gfns::gfn__partial_d_g_dd_333, irho,isigma)
-#define H p.gridfn(gfns::gfn__H, irho,isigma)
-
-#define partial_H_wrt_partial_d_h_1 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1, irho,isigma)
-#define partial_H_wrt_partial_d_h_2 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2, irho,isigma)
-#define partial_H_wrt_partial_dd_h_11 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11, irho,isigma)
-#define partial_H_wrt_partial_dd_h_12 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12, irho,isigma)
-#define partial_H_wrt_partial_dd_h_22 \
- p.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22, irho,isigma)
-
-#define save_H p.gridfn(gfns::gfn__save_H, irho,isigma)
+#define Theta p.gridfn(gfns::gfn__Theta, irho,isigma)
+
+#define partial_Theta_wrt_partial_d_h_1 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_d_h_1, irho,isigma)
+#define partial_Theta_wrt_partial_d_h_2 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_d_h_2, irho,isigma)
+#define partial_Theta_wrt_partial_dd_h_11 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_11, irho,isigma)
+#define partial_Theta_wrt_partial_dd_h_12 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_12, irho,isigma)
+#define partial_Theta_wrt_partial_dd_h_22 \
+ p.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_22, irho,isigma)
+
+#define save_Theta p.gridfn(gfns::gfn__save_Theta, irho,isigma)
#define Delta_h p.gridfn(gfns::gfn__Delta_h, irho,isigma)
//******************************************************************************
@@ -148,7 +148,7 @@ fp partial_d_g_uu_322;
fp partial_d_g_uu_323;
fp partial_d_g_uu_333;
-fp HA;
-fp HB;
-fp HC;
-fp HD;
+fp Theta_A;
+fp Theta_B;
+fp Theta_C;
+fp Theta_D;
diff --git a/src/gr/horizon_function.cc b/src/gr/expansion.cc
index 0821eaf..9cafc2b 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/expansion.cc
@@ -1,8 +1,8 @@
-// horizon_function.cc -- evaluate LHS function H(h)
+// expansion.cc -- evaluate expansion Theta of trial horizon surface
// $Header$
//
// <<<prototypes for functions local to this file>>>
-// horizon_function - top-level driver
+// expansion - top-level driver
// decode_geometry_method - decode the geometry_method parameter
//
/// setup_xyz_posns - setup global xyz posns of grid points
@@ -12,7 +12,7 @@
/// h_is_finite - does function h contain any NaNs/infinities?
/// geometry_is_finite - do geometry vars contain NaN/infty?
///
-/// compute_H - compute H(h) given earlier setup
+/// compute_Theta - compute Theta(h) given earlier setup
///
#include <stdio.h>
@@ -73,10 +73,10 @@ bool h_is_finite(patch_system& ps, bool print_msg_flag);
bool geometry_is_finite(patch_system& ps,
bool print_msg_flag);
-bool compute_H(patch_system& ps,
- bool Jacobian_flag,
- jtutil::norm<fp>* H_norms_ptr,
- bool print_msg_flag);
+bool compute_Theta(patch_system& ps,
+ bool Jacobian_flag,
+ jtutil::norm<fp>* Theta_norms_ptr,
+ bool print_msg_flag);
}
//******************************************************************************
@@ -84,7 +84,7 @@ bool compute_H(patch_system& ps,
//******************************************************************************
//
-// This function computes the LHS function H(h), and optionally also
+// This function computes the LHS function Theta(h), and optionally also
// its Jacobian coefficients (from which the Jacobian matrix may be
// computed later).
//
@@ -113,17 +113,18 @@ bool compute_H(patch_system& ps,
// g_dd_{11,12,13,22,23,33} # $g_{ij}$
// K_dd_{11,12,13,22,23,33} # $K_{ij}$
// partial_d_g_dd_{1,2,3}{11,12,13,22,23,33} # $\partial_k g_{ij}$
-// H # $H = H(h)$
+// ## computed at the nominal grid points
+// Theta # $\Theta = \Theta(h)$
//
// Arguments:
// Jacobian_flag = true to compute the Jacobian coefficients,
// false to skip this.
// print_msg_flag = true to print status messages,
// false to skip this.
-// H_norms_ptr = (out) If this pointer is non-NULL, the norm object it
-// points to is updated with all the H values in the
-// grid. This norm object can then be used to compute
-// various (gridwise) norms of H.
+// Theta_norms_ptr = (out) If this pointer is non-NULL, the norm object it
+// points to is updated with all the Theta values
+// in the grid. This norm object can then be used
+// to compute various (gridwise) norms of Theta.
//
// Results:
// This function returns true for a successful computation, or false
@@ -132,17 +133,17 @@ bool compute_H(patch_system& ps,
// or from an excised region
// FIXME: excision isn't implemented yet :(
// - NaNs are found in h or in the interpolate geometry variables
-// - HD <= 0 (this means the interpolated g_ij aren't positive definite)
+// - Theta_D <= 0 (this means the interpolated g_ij aren't positive definite)
//
-bool horizon_function(patch_system& ps,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- bool Jacobian_flag /* = false */,
- bool print_msg_flag /* = false */,
- jtutil::norm<fp>* H_norms_ptr /* = NULL */)
+bool expansion(patch_system& ps,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ bool Jacobian_flag /* = false */,
+ bool print_msg_flag /* = false */,
+ jtutil::norm<fp>* Theta_norms_ptr /* = NULL */)
{
if (print_msg_flag)
- then CCTK_VInfo(CCTK_THORNSTRING, " horizon function");
+ then CCTK_VInfo(CCTK_THORNSTRING, " expansion");
// fill in values of all ghosted gridfns in ghost zones
ps.synchronize();
@@ -168,7 +169,7 @@ case geometry__Schwarzschild_EF:
break;
default:
error_exit(PANIC_EXIT,
-"***** horizon_function(): unknown gi.geometry_method=(int)%d!\n"
+"***** expansion(): unknown gi.geometry_method=(int)%d!\n"
" (this should never happen!)\n"
,
int(gi.geometry_method)); /*NOTREACHED*/
@@ -178,9 +179,10 @@ if (gi.check_that_geometry_is_finite
&& !geometry_is_finite(ps, print_msg_flag))
then return false; // *** ERROR RETURN ***
-// compute remaining gridfns --> $H$ and optionally Jacobian coefficients
+// compute remaining gridfns --> $\Theta$
+// and optionally also the Jacobian coefficients
// by algebraic ops and angular finite differencing
-if (!compute_H(ps, Jacobian_flag, H_norms_ptr, print_msg_flag))
+if (!compute_Theta(ps, Jacobian_flag, Theta_norms_ptr, print_msg_flag))
then return false; // *** ERROR RETURN ***
return true; // *** NORMAL RETURN ***
@@ -586,7 +588,7 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
if (print_msg_flag)
then {
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
" the trial-horizon-surface point"
@@ -743,7 +745,7 @@ if (print_msg_flag)
// or more NaNs or infinities.
//
// Results:
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
// This function returns true if all the h values are finite, false
// otherwise (i.e. if it contains any NaNs or infinities).
#else
@@ -757,7 +759,7 @@ bool h_is_finite(patch_system& ps, bool print_msg_flag)
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING, " checking that h is finite");
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
patch& p = ps.ith_patch(pn);
@@ -775,7 +777,7 @@ if (print_msg_flag)
const fp sigma = p.sigma_of_isigma(isigma);
const fp drho = jtutil::degrees_of_radians(rho);
const fp dsigma = jtutil::degrees_of_radians(sigma);
- CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+ CCTK_VWarn(expansion_warning_levels::nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
" h=%g isn't finite!\n"
@@ -791,7 +793,7 @@ if (print_msg_flag)
}
return true; // *** all values finite ***
#else
-CCTK_VWarn(horizon_function_warning_levels::skip_nonfinite,
+CCTK_VWarn(expansion_warning_levels::skip_nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
" no finite() fn ==> skipping is-h-finite check");
return true; // *** no check possible ***
@@ -810,7 +812,7 @@ return true; // *** no check possible ***
// or more NaNs or infinities.
//
// Results:
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
// This function returns true if all the geometry variables are finite,
// false otherwise (i.e. if they contain any NaNs or infinities).
#else
@@ -824,7 +826,7 @@ bool geometry_is_finite(patch_system& ps, bool print_msg_flag)
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING, " checking that geometry is finite");
-#ifdef HAVE_FINITE
+#ifdef Theta_AVE_FINITE
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
patch& p = ps.ith_patch(pn);
@@ -921,7 +923,7 @@ if (print_msg_flag)
const fp global_x = ps.global_x_of_local_x(local_x);
const fp global_y = ps.global_y_of_local_y(local_y);
const fp global_z = ps.global_z_of_local_z(local_z);
- CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+ CCTK_VWarn(expansion_warning_levels::nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
" geometry isn't finite at %s patch\n"
@@ -973,7 +975,7 @@ if (print_msg_flag)
}
return true; // *** no NaNs found ***
#else
-CCTK_VWarn(horizon_function_warning_levels::skip_nonfinite,
+CCTK_VWarn(expansion_warning_levels::skip_nonfinite,
__LINE__, __FILE__, CCTK_THORNSTRING,
" no finite() ==> skipping is-geometry-finite check");
return true; // *** no check possible ***
@@ -986,10 +988,11 @@ return true; // *** no check possible ***
//******************************************************************************
//
-// This function computes H(h), and optionally its Jacobian coefficients,
-// (from which the Jacobian matrix may be computed later). This function
-// uses a mixture of algebraic operations and (rho,sigma) finite differencing.
-// The computation is done (entirely) on the nominal angular grid.
+// This function computes the expansion Theta(h), and optionally also
+// its Jacobian coefficients, (from which the Jacobian matrix may be
+// computed later). This function uses a mixture of algebraic operations
+// and (rho,sigma) finite differencing. The computation is done entirely
+// on the nominal angular grid.
//
// N.b. This function #includes "cg.hh", which defines "dangerous" macros
// which will stay in effect for the rest of this compilation unit!
@@ -1000,17 +1003,17 @@ return true; // *** no check possible ***
//
// Results:
// This function returns true for a successful computation, or false
-// if the computation failed because HD <= 0 (this means the interpolated
+// if the computation failed because Theta_D <= 0 (this means the interpolated
// g_ij isn't positive definite).
//
namespace {
-bool compute_H(patch_system& ps,
- bool Jacobian_flag,
- jtutil::norm<fp>* H_norms_ptr,
- bool print_msg_flag)
+bool compute_Theta(patch_system& ps,
+ bool Jacobian_flag,
+ jtutil::norm<fp>* Theta_norms_ptr,
+ bool print_msg_flag)
{
if (print_msg_flag)
- then CCTK_VInfo(CCTK_THORNSTRING, " computing H(h)");
+ then CCTK_VInfo(CCTK_THORNSTRING, " computing Theta(h)");
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
{
@@ -1084,37 +1087,41 @@ if (print_msg_flag)
}
{
- // HA, HB, HC, HD
- #include "../gr.cg/horizon_function.c"
+ // Theta_A, Theta_B, Theta_C, Theta_D
+ #include "../gr.cg/expansion.c"
}
- if (HD <= 0)
+ if (Theta_D <= 0)
then {
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
- "HD <= 0 at %s patch rho=%g sigma=%g!",
+ "Theta_D <= 0 at %s patch rho=%g sigma=%g!",
p.name(), double(rho), double(sigma));
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
"(this probably means the interpolated g_ij");
- CCTK_VWarn(horizon_function_warning_levels::failure,
+ CCTK_VWarn(expansion_warning_levels::failure,
__LINE__, __FILE__, CCTK_THORNSTRING,
" isn't positive definite)");
return false; // *** ERROR RETURN ***
}
// compute H via equation (14) of my 1996 horizon finding paper
- const fp sqrt_HD = sqrt(HD);
- H = HA/(HD*sqrt_HD) + HB/sqrt_HD + HC/HD - K;
+ const fp sqrt_Theta_D = sqrt(Theta_D);
+ Theta = + Theta_A/(Theta_D*sqrt_Theta_D)
+ + Theta_B/sqrt_Theta_D
+ + Theta_C/Theta_D
+ - K;
- // update running norms of H(h) function
- if (H_norms_ptr != NULL)
- then H_norms_ptr->data(H);
+ // update running norms of Theta(h) function
+ if (Theta_norms_ptr != NULL)
+ then Theta_norms_ptr->data(Theta);
if (Jacobian_flag)
then {
- // partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h
- #include "../gr.cg/horizon_Jacobian.c"
+ // partial_Theta_wrt_partial_d_h,
+ // partial_Theta_wrt_partial_dd_h
+ #include "../gr.cg/expansion_Jacobian.c"
}
}
}
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/expansion_Jacobian.cc
index 2d4e8e3..c9cb8bc 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/expansion_Jacobian.cc
@@ -1,14 +1,14 @@
-// horizon_Jacobian.cc -- evaluate Jacobian matrix of LHS function H(h)
+// expansion_Jacobian.cc -- evaluate Jacobian matrix of LHS function Theta(h)
// $Header$
//
// <<<prototypes for functions local to this file>>>
//
-// horizon_Jacobian - top-level driver to compute the Jacobian
+// expansion_Jacobian - top-level driver to compute the Jacobian
// decode_geometry_method - decode the geometry_method parameter
///
-/// horizon_Jacobian_NP - compute the Jacobian by numerical perturbation
+/// expansion_Jacobian_NP - compute the Jacobian by numerical perturbation
///
-/// horizon_Jacobian_partial_SD - compute the partial-deriv terms: symbolic diff
+/// expansion_Jacobian_partial_SD - compute the partial-deriv terms: symbolic diff
/// add_ghost_zone_Jacobian - add ghost zone dependencies to Jacobian
///
@@ -51,58 +51,58 @@ using jtutil::error_exit;
//
namespace {
-bool horizon_Jacobian_NP(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& ii,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag);
-
-void horizon_Jacobian_partial_SD(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag);
+bool expansion_Jacobian_NP(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& ii,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag);
+
+void expansion_Jacobian_partial_SD(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag);
void add_ghost_zone_Jacobian(const patch_system& ps,
Jacobian& Jac,
fp mol,
const patch& xp, const ghost_zone& xmgz,
int x_II,
int xm_irho, int xm_isigma);
-bool horizon_Jacobian_dr_FD(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag);
+bool expansion_Jacobian_dr_FD(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag);
}
//******************************************************************************
//
// This function is a top-level driver to compute the Jacobian matrix
-// J[H(h)] of the LHS function H(h). It just decodes the Jacobian_method
-// parameter and calls the appropriate subfunction.
+// J[Theta(h)] of the expansion Theta(h). It just decodes the
+// Jacobian_method parameter and calls the appropriate subfunction.
//
-// We assume that H(h) has already been computed.
+// We assume that Theta(h) has already been computed.
//
// Results:
// This function returns true for a successful computation, or false
-// for failure. See horizon_function() (in "horizon_function.cc")
+// for failure. See expansion() (in "expansion.cc")
// for details of the various ways the computation may fail.
//
-bool horizon_Jacobian(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag)
+bool expansion_Jacobian(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag)
{
switch (Jacobian_info.Jacobian_method)
{
case Jacobian_method__numerical_perturb:
- if (! horizon_Jacobian_NP(ps, Jac,
+ if (! expansion_Jacobian_NP(ps, Jac,
cgi, gi, Jacobian_info,
print_msg_flag))
then return false; // *** ERROR RETURN ***
@@ -110,22 +110,22 @@ case Jacobian_method__numerical_perturb:
case Jacobian_method__symbolic_diff:
CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
-" horizon_Jacobian(): Jacobian_method == \"symbolic differentiation\"\n"
-" isn't implemented (yet); reverting to\n"
-" \"symbolic differentiation with finite diff d/dr\"");
+" expansion_Jacobian(): Jacobian_method == \"symbolic differentiation\"\n"
+" isn't implemented (yet); reverting to\n"
+" \"symbolic differentiation with finite diff d/dr\"");
// fall through
case Jacobian_method__symbolic_diff_with_FD_dr:
- horizon_Jacobian_partial_SD(ps, Jac,
- cgi, gi, Jacobian_info,
- print_msg_flag);
- if (! horizon_Jacobian_dr_FD(ps, Jac,
- cgi, gi, Jacobian_info,
- print_msg_flag))
+ expansion_Jacobian_partial_SD(ps, Jac,
+ cgi, gi, Jacobian_info,
+ print_msg_flag);
+ if (! expansion_Jacobian_dr_FD(ps, Jac,
+ cgi, gi, Jacobian_info,
+ print_msg_flag))
then return false; // *** ERROR RETURN ***
break;
default:
error_exit(PANIC_EXIT,
-"***** horizon_Jacobian(): unknown Jacobian_info.Jacobian_method=(int)%d!\n"
+"***** expansion_Jacobian(): unknown Jacobian_info.Jacobian_method=(int)%d!\n"
" (this should never happen!)\n"
,
int(Jacobian_info.Jacobian_method)); /*NOTREACHED*/
@@ -161,51 +161,52 @@ else CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
//******************************************************************************
//
-// This function computes the Jacobian matrix of the LHS function H(h)
-// by numerical perturbation of the H(h) function. The algorithm is as
+// This function computes the Jacobian matrix of the expansion Theta(h)
+// by numerical perturbation of the Theta(h) function. The algorithm is as
// follows:
//
-// we assume that H = H(h) has already been evaluated
-// save_H = H
+// we assume that Theta = Theta(h) has already been evaluated
+// save_Theta = Theta
// for each point (y,JJ)
// {
// const fp save_h_y = h at y;
// h at y += perturbation_amplitude;
-// evaluate H(h) (silently)
+// evaluate Theta(h) (silently)
// for each point (x,II)
// {
-// Jac(II,JJ) = (H(II) - save_H(II)) / perturbation_amplitude;
+// Jac(II,JJ) = (Theta(II) - save_Theta(II))
+// / perturbation_amplitude;
// }
// h at y = save_h_y;
// }
-// H = save_H
+// Theta = save_Theta
//
// Inputs (angular gridfns, on ghosted grid):
-// h # shape of trial surface
-// H # H(h) assumed to already be computed
+// h # shape of trial surface
+// Theta # Theta(h) assumed to already be computed
//
// Outputs:
// The Jacobian matrix is stored in the Jacobian object Jac.
//
// Results:
// This function returns true for a successful computation, or false
-// for failure. See horizon_function() (in "horizon_function.cc")
+// for failure. See expansion() (in "expansion.cc")
// for details of the various ways the computation may fail.
//
namespace {
-bool horizon_Jacobian_NP(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& ii,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag)
+bool expansion_Jacobian_NP(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& ii,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag)
{
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING,
" horizon Jacobian (numerical perturbation)");
const fp epsilon = Jacobian_info.perturbation_amplitude;
-ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
+ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
{
@@ -225,7 +226,7 @@ ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
const fp save_h_y = yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma);
yp.ghosted_gridfn(gfns::gfn__h, y_irho,y_isigma) += epsilon;
- if (! horizon_function(ps, cgi, ii))
+ if (! expansion(ps, cgi, ii))
then return false; // *** ERROR RETURN ***
for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
{
@@ -240,9 +241,11 @@ ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
++x_isigma)
{
const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);
- const fp old_H = xp.gridfn(gfns::gfn__save_H, x_irho,x_isigma);
- const fp new_H = xp.gridfn(gfns::gfn__H, x_irho,x_isigma);
- Jac(II,JJ) = (new_H - old_H) / epsilon;
+ const fp old_Theta = xp.gridfn(gfns::gfn__save_Theta,
+ x_irho,x_isigma);
+ const fp new_Theta = xp.gridfn(gfns::gfn__Theta,
+ x_irho,x_isigma);
+ Jac(II,JJ) = (new_Theta - old_Theta) / epsilon;
}
}
}
@@ -251,7 +254,7 @@ ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
}
}
-ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
+ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
return true; // *** NORMAL RETURN ***
}
}
@@ -262,27 +265,27 @@ return true; // *** NORMAL RETURN ***
//
// This function computes the partial derivative terms in the Jacobian
-// matrix of the LHS function H(h), by symbolic differentiation from the
-// Jacobian coefficient (angular) gridfns. The computation is done a
+// matrix of the expansion Theta(h), by symbolic differentiation from
+// the Jacobian coefficient (angular) gridfns. The computation is done a
// Jacobian row at a time, using equation (25) of my 1996 apparent horizon
// finding paper.
//
// Inputs (angular gridfns, on ghosted grid):
-// h # shape of trial surface
-// H # H(h) assumed to already be computed
-// partial_H_wrt_partial_d_h # Jacobian coefficients
-// partial_H_wrt_partial_dd_h # (also assumed to already be computed)
+// h # shape of trial surface
+// Theta # Theta(h) assumed to already be computed
+// partial_Theta_wrt_partial_d_h # Jacobian coefficients
+// partial_Theta_wrt_partial_dd_h # (also assumed to already be computed)
//
// Outputs:
// The Jacobian matrix is stored in the Jacobian object Jac.
//
namespace {
-void horizon_Jacobian_partial_SD(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag)
+void expansion_Jacobian_partial_SD(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag)
{
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING,
@@ -303,7 +306,7 @@ ps.compute_synchronize_Jacobian();
{
//
// compute the main Jacobian terms for this grid point, i.e.
- // partial H(this point x, Jacobian row II)
+ // partial Theta(this point x, Jacobian row II)
// ---------------------------------------------
// partial h(other points y, Jacobian column JJ)
//
@@ -313,19 +316,19 @@ ps.compute_synchronize_Jacobian();
// Jacobian coefficients for this point
const fp Jacobian_coeff_rho
- = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_1,
+ = xp.gridfn(gfns::gfn__partial_Theta_wrt_partial_d_h_1,
x_irho, x_isigma);
const fp Jacobian_coeff_sigma
- = xp.gridfn(gfns::gfn__partial_H_wrt_partial_d_h_2,
+ = xp.gridfn(gfns::gfn__partial_Theta_wrt_partial_d_h_2,
x_irho, x_isigma);
const fp Jacobian_coeff_rho_rho
- = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_11,
+ = xp.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_11,
x_irho, x_isigma);
const fp Jacobian_coeff_rho_sigma
- = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_12,
+ = xp.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_12,
x_irho, x_isigma);
const fp Jacobian_coeff_sigma_sigma
- = xp.gridfn(gfns::gfn__partial_H_wrt_partial_dd_h_22,
+ = xp.gridfn(gfns::gfn__partial_Theta_wrt_partial_dd_h_22,
x_irho, x_isigma);
// partial_rho, partial_rho_rho
@@ -471,14 +474,14 @@ const patch& yp = ye.my_patch();
//
// This function sums the d/dr terms into the Jacobian matrix of the
-// LHS function H(h), computing those terms by finite differencing.
+// expansion Theta(h), computing those terms by finite differencing.
//
// The basic algorithm is that
-// Jac += diag[ (H(h+epsilon) - H(h)) / epsilon ]
+// Jac += diag[ (Theta(h+epsilon) - Theta(h)) / epsilon ]
//
// Inputs (angular gridfns, on ghosted grid):
-// h # shape of trial surface
-// H # H(h) assumed to already be computed
+// h # shape of trial surface
+// Theta # Theta(h) assumed to already be computed
//
// Outputs:
// Jac += d/dr terms
@@ -490,12 +493,12 @@ const patch& yp = ye.my_patch();
// FIXME: excision isn't implemented yet :(
//
namespace {
-bool horizon_Jacobian_dr_FD(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag)
+bool expansion_Jacobian_dr_FD(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag)
{
if (print_msg_flag)
then CCTK_VInfo(CCTK_THORNSTRING,
@@ -503,10 +506,10 @@ if (print_msg_flag)
const fp epsilon = Jacobian_info.perturbation_amplitude;
-// compute H(h+epsilon)
-ps.gridfn_copy(gfns::gfn__H, gfns::gfn__save_H);
+// compute Theta(h+epsilon)
+ps.gridfn_copy(gfns::gfn__Theta, gfns::gfn__save_Theta);
ps.add_to_ghosted_gridfn(epsilon, gfns::gfn__h);
-if (! horizon_function(ps, cgi, gi))
+if (! expansion(ps, cgi, gi))
then return false; // *** ERROR RETURN ***
for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
@@ -519,16 +522,18 @@ if (! horizon_function(ps, cgi, gi))
++isigma)
{
const int II = ps.gpn_of_patch_irho_isigma(p, irho,isigma);
- const fp old_H = p.gridfn(gfns::gfn__save_H, irho,isigma);
- const fp new_H = p.gridfn(gfns::gfn__H, irho,isigma);
- Jac(II,II) += (new_H - old_H) / epsilon;
+ const fp old_Theta = p.gridfn(gfns::gfn__save_Theta,
+ irho,isigma);
+ const fp new_Theta = p.gridfn(gfns::gfn__Theta,
+ irho,isigma);
+ Jac(II,II) += (new_Theta - old_Theta) / epsilon;
}
}
}
-// restore h and H
+// restore h and Theta
ps.add_to_ghosted_gridfn(-epsilon, gfns::gfn__h);
-ps.gridfn_copy(gfns::gfn__save_H, gfns::gfn__H);
+ps.gridfn_copy(gfns::gfn__save_Theta, gfns::gfn__Theta);
return true; // *** NORMAL RETURN ***
}
diff --git a/src/gr/gfns.hh b/src/gr/gfns.hh
index 5c3f6d3..f64fa56 100644
--- a/src/gr/gfns.hh
+++ b/src/gr/gfns.hh
@@ -60,14 +60,14 @@ enum {
gfn__partial_d_psi_2, // no access macro
gfn__partial_d_psi_3, // no access macro
- gfn__H,
- gfn__partial_H_wrt_partial_d_h_1,
- gfn__partial_H_wrt_partial_d_h_2,
- gfn__partial_H_wrt_partial_dd_h_11,
- gfn__partial_H_wrt_partial_dd_h_12,
- gfn__partial_H_wrt_partial_dd_h_22,
+ gfn__Theta,
+ gfn__partial_Theta_wrt_partial_d_h_1,
+ gfn__partial_Theta_wrt_partial_d_h_2,
+ gfn__partial_Theta_wrt_partial_dd_h_11,
+ gfn__partial_Theta_wrt_partial_dd_h_12,
+ gfn__partial_Theta_wrt_partial_dd_h_22,
gfn__Delta_h,
- gfn__save_H,
+ gfn__save_Theta,
gfn__one,
nominal_max_gfn = gfn__one // no comma
};
diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index 45771ef..8ba54f8 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -129,12 +129,12 @@ struct Jacobian_info
// prototypes for functions visible outside their source files
//
-// horizon_function.cc
+// expansion.cc
// ... returns true for successful computation,
// ... returns false for failure (eg horizon outside grid);
// if (print_msg_flag) then also reports CCTK_VWarn() at the
// appropriate warning level
-namespace horizon_function_warning_levels
+namespace expansion_warning_levels
{
enum {
failure = 2, // point outside grid etc
@@ -143,23 +143,23 @@ namespace horizon_function_warning_levels
// ==> skipping nonfinite check
};
}
-bool horizon_function(patch_system& ps,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- bool Jacobian_flag = false,
- bool print_msg_flag = false,
- jtutil::norm<fp>* H_norms_ptr = NULL);
+bool expansion(patch_system& ps,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ bool Jacobian_flag = false,
+ bool print_msg_flag = false,
+ jtutil::norm<fp>* H_norms_ptr = NULL);
enum geometry_method
decode_geometry_method(const char geometry_method_string[]);
-// horizon_Jacobian.cc
+// expansion_Jacobian.cc
// ... returns true for successful computation, false for failure
-bool horizon_Jacobian(patch_system& ps,
- Jacobian& Jac,
- const struct cactus_grid_info& cgi,
- const struct geometry_info& gi,
- const struct Jacobian_info& Jacobian_info,
- bool print_msg_flag);
+bool expansion_Jacobian(patch_system& ps,
+ Jacobian& Jac,
+ const struct cactus_grid_info& cgi,
+ const struct geometry_info& gi,
+ const struct Jacobian_info& Jacobian_info,
+ bool print_msg_flag);
enum Jacobian_method
decode_Jacobian_method(const char Jacobian_method_string[]);
diff --git a/src/gr/gr_gfas.minc b/src/gr/gr_gfas.minc
index 7579e4d..34da085 100644
--- a/src/gr/gr_gfas.minc
+++ b/src/gr/gr_gfas.minc
@@ -21,24 +21,24 @@ partial_d_s_d, partial_d_s_d__fnd,
n_u, n_u__fnd,
-HA, HA__fnd,
-HB, HB__fnd,
-HC, HC__fnd,
-HD, HD__fnd,
-H, H__fnd,
-
-partial_d_HA, partial_d_HA__fnd,
-partial_d_HB, partial_d_HB__fnd,
-partial_d_HC, partial_d_HC__fnd,
-partial_d_HD, partial_d_HD__fnd,
-partial_d_H, partial_d_H__fnd,
-
-partial_HA_wrt_partial_d_h, partial_HA_wrt_partial_d_h__fnd,
-partial_HB_wrt_partial_d_h, partial_HB_wrt_partial_d_h__fnd,
-partial_HC_wrt_partial_d_h, partial_HC_wrt_partial_d_h__fnd,
-partial_HD_wrt_partial_d_h, partial_HD_wrt_partial_d_h__fnd,
-partial_HA_wrt_partial_dd_h, partial_HA_wrt_partial_dd_h__fnd,
-partial_HB_wrt_partial_dd_h, partial_HB_wrt_partial_dd_h__fnd,
-
-partial_H_wrt_partial_d_h, partial_H_wrt_partial_d_h__fnd,
-partial_H_wrt_partial_dd_h, partial_H_wrt_partial_dd_h__fnd # no comma
+Theta_A, Theta_A__fnd,
+Theta_B, Theta_B__fnd,
+Theta_C, Theta_C__fnd,
+Theta_D, Theta_D__fnd,
+Theta, Theta__fnd,
+
+partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd,
+partial_d_Theta_C, partial_d_Theta_C__fnd,
+partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd,
+
+partial_Theta_A_wrt_partial_d_h, partial_Theta_A_wrt_partial_d_h__fnd,
+partial_Theta_B_wrt_partial_d_h, partial_Theta_B_wrt_partial_d_h__fnd,
+partial_Theta_C_wrt_partial_d_h, partial_Theta_C_wrt_partial_d_h__fnd,
+partial_Theta_D_wrt_partial_d_h, partial_Theta_D_wrt_partial_d_h__fnd,
+partial_Theta_A_wrt_partial_dd_h, partial_Theta_A_wrt_partial_dd_h__fnd,
+partial_Theta_B_wrt_partial_dd_h, partial_Theta_B_wrt_partial_dd_h__fnd,
+
+partial_Theta_wrt_partial_d_h, partial_Theta_wrt_partial_d_h__fnd,
+partial_Theta_wrt_partial_dd_h, partial_Theta_wrt_partial_dd_h__fnd # no comma
diff --git a/src/gr/horizon.maple b/src/gr/horizon.maple
index 433c723..7284253 100644
--- a/src/gr/horizon.maple
+++ b/src/gr/horizon.maple
@@ -4,8 +4,8 @@
# horizon - overall driver
## non_unit_normal - compute s_d
## non_unit_normal_deriv - compute partial_d_s_d
-## horizon_function - compute H[ABCD],H = fn(s_d__fnd) = fn(h)
-## horizon_Jacobian - compute Jacobian coeffs for H[ABCD], H
+## expansion - compute Theta_[ABCD],Theta = fn(s_d__fnd) = fn(h)
+## expansion_Jacobian - compute Jacobian coeffs for Theta_[ABCD], Theta
#
################################################################################
@@ -20,15 +20,17 @@
# h
#
# Outputs:
-# H = H__fnd(h, ...)
+# Theta = Theta__fnd(h, ...)
+# partial_Theta_wrt_partial_d_h__fnd = partial_Theta_wrt_partial_d_h(h, ...)
+# partial_Theta_wrt_partial_dd_h__fnd = partial_Theta_wrt_partial_dd_h(h, ...)
#
horizon :=
proc(cg_flag::boolean)
non_unit_normal();
non_unit_normal_deriv();
-horizon_function(cg_flag);
-horizon_Jacobian(cg_flag);
+expansion(cg_flag);
+expansion_Jacobian(cg_flag);
NULL;
end proc;
@@ -124,22 +126,21 @@ 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$, using equations (14) and (15[abcd]) in my 1996
-# apparent horizon finding paper, but using partial_d_g_uu[k,i,j]
-# instead of Diff(g_uu[i,j], x_xyz[k]).
+# This function computes the expansion $\Theta$ of a trial horizon surface
+# as a function of 1st and 2nd angular partial derivatives of the surface
+# radius $h$, using equations (14) and (15[abcd]) in my 1996 AH-finding
+# paper, but using partial_d_g_uu[k,i,j] instead of Diff(g_uu[i,j], x_xyz[k]).
#
-# These equations give H[ABCD] and H as functions of s_d, but here we
-# use s_d__fnd, which we assume is already given in terms of angular
-# derivatives of h. The result is that we compute H[ABCD] and H directly
-# in terms of 1st and 2nd angular derivatives of h, without s_d ever
-# appearing in our final results.
+# These equations give Theta_[ABCD] and Theta as functions of s_d, but
+# here we use s_d__fnd, which we assume is already given in terms of angular
+# derivatives of h. The result is that we compute Theta_[ABCD] and Theta
+# directly in terms of 1st and 2nd angular derivatives of h, without s_d
+# ever appearing in our final results.
#
# This function also optionally generates C code for the computation
-# of H[ABCD]. It does *not* compute C code for the computation of H
-# itself, since we may want to check that HD > 0 in the C code before
-# we compute H itself.
+# of Theta_[ABCD]. It does *not* compute C code for the computation of
+# Theta itself, since we may want to check that Theta_D > 0 in the C
+# code before we compute Theta itself.
#
# Inputs:
# s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) )
@@ -147,17 +148,17 @@ end proc;
# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
#
# Outputs:
-# HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HB = HB__fnd( x_xyz, X_ud, X_udd,
-# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
-# HD = HD__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
+# Theta_A = Theta_A__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_B = Theta_B__fnd( x_xyz, X_ud, X_udd,
+# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_C = Theta_C__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
+# Theta_D = Theta_D__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
# ---------------------------------------------------------------
-# H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd )
+# Theta = Theta__fnd(Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd)
#
-horizon_function :=
+expansion :=
proc(cg_flag::boolean)
global
@include "../maple/coords.minc",
@@ -180,33 +181,38 @@ assert_fnd_exists(partial_d_s_d, fnd);
#
# (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]
- * partial_d_g_uu[i,k,l] * s_d__fnd[k]*s_d__fnd[l]',
- 'i'=1..N, 'j'=1..N, 'k'=1..N, 'l'=1..N);
+Theta_A__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]
+ * partial_d_g_uu[i,k,l] * 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('partial_d_g_uu[i,i,j]*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);
+Theta_B__fnd
+ := + msum('partial_d_g_uu[i,i,j]*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);
+Theta_C__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);
+Theta_D__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;
+# denominator, which would make the equation *much* more complicated
+Theta__fnd := + Theta_A__fnd/Theta_D__fnd^(3/2)
+ + Theta_B__fnd/Theta_D__fnd^(1/2)
+ + Theta_C__fnd/Theta_D__fnd
+ - K;
if (cg_flag)
- then codegen2([ HA__fnd, HB__fnd, HC__fnd, HD__fnd],
- ['HA', 'HB', 'HC', 'HD' ],
- "../gr.cg/horizon_function.c");
+ then codegen2([ Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd],
+ ['Theta_A', 'Theta_B', 'Theta_C', 'Theta_D' ],
+ "../gr.cg/expansion.c");
fi;
NULL;
@@ -215,12 +221,12 @@ end proc;
################################################################################
#
-# This function computes the Jacobian coefficients for the LHS of the
-# apparent horizon equation with respect to 1st and 2nd angular partial
-# derivatives of the apparent horizon radius $h$. These coefficients
+# This function computes the Jacobian coefficients for the expansion
+# $\Theta$ of a trial horizon surface with respect to 1st and 2nd angular
+# partial derivatives of the surface radius $h$. These coefficients
# are (or at least should be :) the same as those in equation (A1) in
# my 1996 apparent horizon finding paper, bue here we compute them in
-# a different manner: we have Maple directly differentiate H__fnd with
+# a different manner: we have Maple directly differentiate Theta__fnd with
# respect to Diff(h,y_rs[u]) and Diff(h,y_rs[u],y_rs[v]). We use the
# Maple frontend() function to do this.
#
@@ -228,24 +234,24 @@ end proc;
# coefficients.
#
# Inputs:
-# HA = HA__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HB = HB__fnd( x_xyz, X_ud, X_udd,
-# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
-# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
-# HC = HC__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
-# HD = HD__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
-# H = H__fnd( HA__fnd, HB__fnd, HC__fnd, HD__fnd )
+# Theta_A = Theta_A__fnd( x_xyz, X_ud, X_udd, g_uu, partial_d_g_uu,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_B = Theta_B__fnd( x_xyz, X_ud, X_udd,
+# g_uu, partial_d_g_uu, partial_d_ln_sqrt_g,
+# Diff(h,y_rs), Diff(h,y_rs,y_rs) )
+# Theta_C = Theta_C__fnd( x_xyz, X_ud, K_uu, Diff(h,y_rs) )
+# Theta_D = Theta_D__fnd( x_xyz, X_ud, g_uu, Diff(h,y_rs) )
+# Theta = Theta__fnd(Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd)
#
# Outputs:
-# partial_HA_wrt_partial_d_h partial_HA_wrt_partial_dd_h
-# partial_HB_wrt_partial_d_h partial_HB_wrt_partial_dd_h
-# partial_HC_wrt_partial_d_h
-# partial_HD_wrt_partial_d_h
+# partial_Theta_A_wrt_partial_d_h partial_Theta_A_wrt_partial_dd_h
+# partial_Theta_B_wrt_partial_d_h partial_Theta_B_wrt_partial_dd_h
+# partial_Theta_C_wrt_partial_d_h
+# partial_Theta_D_wrt_partial_d_h
# ---------------------------------------------------------------
-# partial_H_wrt_partial_d_h partial_H_wrt_partial_dd_h
+# partial_Theta_wrt_partial_d_h partial_Theta_wrt_partial_dd_h
#
-horizon_Jacobian :=
+expansion_Jacobian :=
proc(cg_flag::boolean)
global
@include "../maple/coords.minc",
@@ -260,54 +266,57 @@ assert_fnd_exists(g_uu);
assert_fnd_exists(K_uu);
assert_fnd_exists(partial_d_ln_sqrt_g);
assert_fnd_exists(partial_d_g_uu);
-assert_fnd_exists(HA);
-assert_fnd_exists(HB);
-assert_fnd_exists(HC);
-assert_fnd_exists(HD);
+assert_fnd_exists(Theta_A);
+assert_fnd_exists(Theta_B);
+assert_fnd_exists(Theta_C);
+assert_fnd_exists(Theta_D);
-# Jacobian coefficients of H[ABCD] and H wrt Diff(h,y_rs[u])
+# Jacobian coefficients of Theta_[ABCD] and Theta wrt Diff(h,y_rs[u])
for u from 1 to N_ang
do
- partial_HA_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HA__fnd, Diff(h,y_rs[u])]);
- partial_HB_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HB__fnd, Diff(h,y_rs[u])]);
- partial_HC_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HC__fnd, Diff(h,y_rs[u])]);
- partial_HD_wrt_partial_d_h__fnd[u]
- := frontend('diff', [HD__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_A_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_A__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_B_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_B__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_C_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_C__fnd, Diff(h,y_rs[u])]);
+ partial_Theta_D_wrt_partial_d_h__fnd[u]
+ := frontend('diff', [Theta_D__fnd, Diff(h,y_rs[u])]);
# equation (A1a) in my 1996 apparent horizon finding paper
- temp := (3/2)*HA/HD^(5/2) + (1/2)*HB/HD^(3/2) + HC/HD^2;
- partial_H_wrt_partial_d_h__fnd[u]
- := + partial_HA_wrt_partial_d_h__fnd[u] / HD^(3/2)
- + partial_HB_wrt_partial_d_h__fnd[u] / HD^(1/2)
- + partial_HC_wrt_partial_d_h__fnd[u] / HD
- - partial_HD_wrt_partial_d_h__fnd[u] * temp;
+ temp := + (3/2)*Theta_A/Theta_D^(5/2)
+ + (1/2)*Theta_B/Theta_D^(3/2)
+ + Theta_C/Theta_D^2;
+ partial_Theta_wrt_partial_d_h__fnd[u]
+ := + partial_Theta_A_wrt_partial_d_h__fnd[u] / Theta_D^(3/2)
+ + partial_Theta_B_wrt_partial_d_h__fnd[u] / Theta_D^(1/2)
+ + partial_Theta_C_wrt_partial_d_h__fnd[u] / Theta_D
+ - partial_Theta_D_wrt_partial_d_h__fnd[u] * temp;
end do;
-# Jacobian coefficients of H[AB] and H wrt Diff(h,y_rs[u],y_rs[v])
+# Jacobian coefficients of Theta_[AB] and Theta wrt Diff(h,y_rs[u],y_rs[v])
for u from 1 to N_ang
do
for v from u to N_ang
do
- partial_HA_wrt_partial_dd_h__fnd[u,v]
- := frontend('diff', [HA__fnd, Diff(h,y_rs[u],y_rs[v])]);
- partial_HB_wrt_partial_dd_h__fnd[u,v]
- := frontend('diff', [HB__fnd, Diff(h,y_rs[u],y_rs[v])]);
+ partial_Theta_A_wrt_partial_dd_h__fnd[u,v]
+ := frontend('diff', [Theta_A__fnd, Diff(h,y_rs[u],y_rs[v])]);
+ partial_Theta_B_wrt_partial_dd_h__fnd[u,v]
+ := frontend('diff', [Theta_B__fnd, Diff(h,y_rs[u],y_rs[v])]);
# equation (A1b) in my 1996 apparent horizon finding paper
- partial_H_wrt_partial_dd_h__fnd[u,v]
- := + partial_HA_wrt_partial_dd_h__fnd[u,v] / HD^(3/2)
- + partial_HB_wrt_partial_dd_h__fnd[u,v] / HD^(1/2)
+ partial_Theta_wrt_partial_dd_h__fnd[u,v]
+ := + partial_Theta_A_wrt_partial_dd_h__fnd[u,v] / Theta_D^(3/2)
+ + partial_Theta_B_wrt_partial_dd_h__fnd[u,v] / Theta_D^(1/2)
end do;
end do;
if (cg_flag)
- then codegen2([partial_H_wrt_partial_d_h__fnd,
- partial_H_wrt_partial_dd_h__fnd],
- ['partial_H_wrt_partial_d_h', 'partial_H_wrt_partial_dd_h'],
- "../gr.cg/horizon_Jacobian.c");
+ then codegen2([partial_Theta_wrt_partial_d_h__fnd,
+ partial_Theta_wrt_partial_dd_h__fnd],
+ ['partial_Theta_wrt_partial_d_h',
+ 'partial_Theta_wrt_partial_dd_h'],
+ "../gr.cg/expansion_Jacobian.c");
fi;
NULL;
diff --git a/src/gr/make.code.defn b/src/gr/make.code.defn
index a371c50..e93414f 100644
--- a/src/gr/make.code.defn
+++ b/src/gr/make.code.defn
@@ -2,8 +2,8 @@
# $Header$
# Source files in this directory
-SRCS = horizon_function.cc \
- horizon_Jacobian.cc \
+SRCS = expansion.cc \
+ expansion_Jacobian.cc \
Schwarzschild_EF.cc
# Subdirectories containing source files
diff --git a/src/gr/maple.log b/src/gr/maple.log
index 0a28bb5..54c762f 100644
--- a/src/gr/maple.log
+++ b/src/gr/maple.log
@@ -586,19 +586,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
make_gfa('g_dd', {inert}, [1 .. N, 1 .. N], symmetric);
make_gfa('K_dd', {inert}, [1 .. N, 1 .. N], symmetric);
make_gfa('g_uu', {inert, fnd}, [1 .. N, 1 .. N], symmetric);
@@ -612,31 +613,31 @@ partial_H_wrt_partial_dd_h__fnd;
make_gfa('h', {inert, fnd}, [], none);
make_gfa('s_d', {inert, fnd}, [1 .. N], none);
make_gfa('partial_d_s_d', {inert, fnd}, [1 .. N, 1 .. N], none);
- make_gfa('HA', {inert, fnd}, [], none);
- make_gfa('HB', {inert, fnd}, [], none);
- make_gfa('HC', {inert, fnd}, [], none);
- make_gfa('HD', {inert, fnd}, [], none);
- make_gfa('H', {inert, fnd}, [], none);
- make_gfa('partial_d_HA', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_HB', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_HC', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_HD', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_d_H', {inert, fnd}, [1 .. N], none);
- make_gfa('partial_HA_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HB_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HC_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HD_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_HA_wrt_partial_dd_h', {inert, fnd},
+ make_gfa('Theta_A', {inert, fnd}, [], none);
+ make_gfa('Theta_B', {inert, fnd}, [], none);
+ make_gfa('Theta_C', {inert, fnd}, [], none);
+ make_gfa('Theta_D', {inert, fnd}, [], none);
+ make_gfa('Theta', {inert, fnd}, [], none);
+ make_gfa('partial_d_Theta_A', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta_B', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta_C', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta_D', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_d_Theta', {inert, fnd}, [1 .. N], none);
+ make_gfa('partial_Theta_A_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_B_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_C_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_D_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_A_wrt_partial_dd_h', {inert, fnd},
[1 .. N_ang, 1 .. N_ang], symmetric);
- make_gfa('partial_HB_wrt_partial_dd_h', {inert, fnd},
+ make_gfa('partial_Theta_B_wrt_partial_dd_h', {inert, fnd},
[1 .. N_ang, 1 .. N_ang], symmetric);
- make_gfa('partial_H_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang], none)
- ;
- make_gfa('partial_H_wrt_partial_dd_h', {inert, fnd},
+ make_gfa('partial_Theta_wrt_partial_d_h', {inert, fnd}, [1 .. N_ang],
+ none);
+ make_gfa('partial_Theta_wrt_partial_dd_h', {inert, fnd},
[1 .. N_ang, 1 .. N_ang], symmetric);
NULL
end proc
@@ -648,19 +649,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
option remember;
var_list := [args[2 .. nargs]];
if type(operand, indexed) and op(0, operand) = 'g_dd' and
@@ -683,19 +685,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_dd);
assert_fnd_exists(g_uu, fnd);
@@ -713,19 +716,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_uu);
assert_fnd_exists(K_dd);
@@ -756,19 +760,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_dd);
assert_fnd_exists(g_uu);
@@ -793,19 +798,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_dd);
assert_fnd_exists(g_uu);
@@ -823,8 +829,8 @@ end proc
horizon := proc(cg_flag::boolean)
non_unit_normal();
non_unit_normal_deriv();
- horizon_function(cg_flag);
- horizon_Jacobian(cg_flag);
+ expansion(cg_flag);
+ expansion_Jacobian(cg_flag);
NULL
end proc
@@ -835,19 +841,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(h);
assert_fnd_exists(X_ud);
@@ -865,19 +872,20 @@ rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(h);
assert_fnd_exists(X_ud);
@@ -897,26 +905,27 @@ partial_H_wrt_partial_dd_h__fnd;
NULL
end proc
-horizon_function := proc(cg_flag::boolean)
+expansion := proc(cg_flag::boolean)
local i, j, k, l;
global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_uu);
assert_fnd_exists(K_uu);
@@ -924,88 +933,94 @@ partial_H_wrt_partial_dd_h__fnd;
assert_fnd_exists(partial_d_g_uu);
assert_fnd_exists(s_d, fnd);
assert_fnd_exists(partial_d_s_d, fnd);
- HA__fnd := -msum('g_uu[i, k]*s_d__fnd[k]*g_uu[j, l]*s_d__fnd[l]*
+ Theta_A__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]*
partial_d_g_uu[i, k, l]*s_d__fnd[k]*s_d__fnd[l]', 'i' = 1 .. N,
'j' = 1 .. N, 'k' = 1 .. N, 'l' = 1 .. N);
- HB__fnd := msum('partial_d_g_uu[i, i, j]*s_d__fnd[j]', 'i' = 1 .. N,
- 'j' = 1 .. N) + msum('g_uu[i, j]*partial_d_s_d__fnd[i, j]',
+ Theta_B__fnd := msum('partial_d_g_uu[i, i, j]*s_d__fnd[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,
+ '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);
+ Theta_C__fnd := msum('K_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
'j' = 1 .. N);
- HC__fnd := msum('K_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
+ Theta_D__fnd := msum('g_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
'j' = 1 .. N);
- HD__fnd := msum('g_uu[i, j]*s_d__fnd[i]*s_d__fnd[j]', 'i' = 1 .. N,
- 'j' = 1 .. N);
- H__fnd :=
- HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K
- ;
- if cg_flag then codegen2([HA__fnd, HB__fnd, HC__fnd, HD__fnd],
- ['HA', 'HB', 'HC', 'HD'], "../gr.cg/horizon_function.c")
+ Theta__fnd := Theta_A__fnd/Theta_D__fnd^(3/2)
+ + Theta_B__fnd/Theta_D__fnd^(1/2) + Theta_C__fnd/Theta_D__fnd - K;
+ if cg_flag then codegen2(
+ [Theta_A__fnd, Theta_B__fnd, Theta_C__fnd, Theta_D__fnd],
+ ['Theta_A', 'Theta_B', 'Theta_C', 'Theta_D'],
+ "../gr.cg/expansion.c")
end if;
NULL
end proc
-horizon_Jacobian := proc(cg_flag::boolean)
+expansion_Jacobian := proc(cg_flag::boolean)
local u, v, temp;
global delta, N, N_ang, xx, yy, zz, x_xyz, x_xyz_list, x_xyz_set, r, r__fnd,
rho, sigma, y_rs, y_rs_list, y_rs_set, xy_all_list, xy_all_set, inert, none,
fnd, symmetric3_23, X_ud, X_ud__fnd, X_udd, X_udd__fnd, g_dd, K_dd, g_uu,
g_uu__fnd, K_uu, K_uu__fnd, K, K__fnd, partial_d_g_dd, partial_d_ln_sqrt_g,
partial_d_ln_sqrt_g__fnd, partial_d_g_uu, partial_d_g_uu__fnd, h, h__fnd,
-s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, HA,
-HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd, H, H__fnd, partial_d_HA,
-partial_d_HA__fnd, partial_d_HB, partial_d_HB__fnd, partial_d_HC,
-partial_d_HC__fnd, partial_d_HD, partial_d_HD__fnd, partial_d_H,
-partial_d_H__fnd, partial_HA_wrt_partial_d_h,
-partial_HA_wrt_partial_d_h__fnd, partial_HB_wrt_partial_d_h,
-partial_HB_wrt_partial_d_h__fnd, partial_HC_wrt_partial_d_h,
-partial_HC_wrt_partial_d_h__fnd, partial_HD_wrt_partial_d_h,
-partial_HD_wrt_partial_d_h__fnd, partial_HA_wrt_partial_dd_h,
-partial_HA_wrt_partial_dd_h__fnd, partial_HB_wrt_partial_dd_h,
-partial_HB_wrt_partial_dd_h__fnd, partial_H_wrt_partial_d_h,
-partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h,
-partial_H_wrt_partial_dd_h__fnd;
+s_d, s_d__fnd, partial_d_s_d, partial_d_s_d__fnd, n_u, n_u__fnd, Theta_A,
+Theta_A__fnd, Theta_B, Theta_B__fnd, Theta_C, Theta_C__fnd, Theta_D,
+Theta_D__fnd, Theta, Theta__fnd, partial_d_Theta_A, partial_d_Theta_A__fnd,
+partial_d_Theta_B, partial_d_Theta_B__fnd, partial_d_Theta_C,
+partial_d_Theta_C__fnd, partial_d_Theta_D, partial_d_Theta_D__fnd,
+partial_d_Theta_, partial_d_Theta__fnd, partial_Theta_A_wrt_partial_d_h,
+partial_Theta_A_wrt_partial_d_h__fnd, partial_Theta_B_wrt_partial_d_h,
+partial_Theta_B_wrt_partial_d_h__fnd, partial_Theta_C_wrt_partial_d_h,
+partial_Theta_C_wrt_partial_d_h__fnd, partial_Theta_D_wrt_partial_d_h,
+partial_Theta_D_wrt_partial_d_h__fnd, partial_Theta_A_wrt_partial_dd_h,
+partial_Theta_A_wrt_partial_dd_h__fnd, partial_Theta_B_wrt_partial_dd_h,
+partial_Theta_B_wrt_partial_dd_h__fnd, partial_Theta_wrt_partial_d_h,
+partial_Theta_wrt_partial_d_h__fnd, partial_Theta_wrt_partial_dd_h,
+partial_Theta_wrt_partial_dd_h__fnd;
printf("%a...\n", procname);
assert_fnd_exists(g_uu);
assert_fnd_exists(K_uu);
assert_fnd_exists(partial_d_ln_sqrt_g);
assert_fnd_exists(partial_d_g_uu);
- assert_fnd_exists(HA);
- assert_fnd_exists(HB);
- assert_fnd_exists(HC);
- assert_fnd_exists(HD);
+ assert_fnd_exists(Theta_A);
+ assert_fnd_exists(Theta_B);
+ assert_fnd_exists(Theta_C);
+ assert_fnd_exists(Theta_D);
for u to N_ang do
- partial_HA_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HA__fnd, Diff(h, y_rs[u])]);
- partial_HB_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HB__fnd, Diff(h, y_rs[u])]);
- partial_HC_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HC__fnd, Diff(h, y_rs[u])]);
- partial_HD_wrt_partial_d_h__fnd[u] :=
- frontend('diff', [HD__fnd, Diff(h, y_rs[u])]);
- temp := 3/2*HA/HD^(5/2) + 1/2*HB/HD^(3/2) + HC/HD^2;
- partial_H_wrt_partial_d_h__fnd[u] :=
- partial_HA_wrt_partial_d_h__fnd[u]/HD^(3/2)
- + partial_HB_wrt_partial_d_h__fnd[u]/HD^(1/2)
- + partial_HC_wrt_partial_d_h__fnd[u]/HD
- - partial_HD_wrt_partial_d_h__fnd[u]*temp
+ partial_Theta_A_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_A__fnd, Diff(h, y_rs[u])]);
+ partial_Theta_B_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_B__fnd, Diff(h, y_rs[u])]);
+ partial_Theta_C_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_C__fnd, Diff(h, y_rs[u])]);
+ partial_Theta_D_wrt_partial_d_h__fnd[u] :=
+ frontend('diff', [Theta_D__fnd, Diff(h, y_rs[u])]);
+ temp := 3/2*Theta_A/Theta_D^(5/2) + 1/2*Theta_B/Theta_D^(3/2)
+ + Theta_C/Theta_D^2;
+ partial_Theta_wrt_partial_d_h__fnd[u] :=
+ partial_Theta_A_wrt_partial_d_h__fnd[u]/Theta_D^(3/2)
+ + partial_Theta_B_wrt_partial_d_h__fnd[u]/Theta_D^(1/2)
+ + partial_Theta_C_wrt_partial_d_h__fnd[u]/Theta_D
+ - partial_Theta_D_wrt_partial_d_h__fnd[u]*temp
end do;
for u to N_ang do for v from u to N_ang do
- partial_HA_wrt_partial_dd_h__fnd[u, v] :=
- frontend('diff', [HA__fnd, Diff(h, y_rs[u], y_rs[v])]);
- partial_HB_wrt_partial_dd_h__fnd[u, v] :=
- frontend('diff', [HB__fnd, Diff(h, y_rs[u], y_rs[v])]);
- partial_H_wrt_partial_dd_h__fnd[u, v] :=
- partial_HA_wrt_partial_dd_h__fnd[u, v]/HD^(3/2)
- + partial_HB_wrt_partial_dd_h__fnd[u, v]/HD^(1/2)
+ partial_Theta_A_wrt_partial_dd_h__fnd[u, v] :=
+ frontend('diff', [Theta_A__fnd, Diff(h, y_rs[u], y_rs[v])])
+ ;
+ partial_Theta_B_wrt_partial_dd_h__fnd[u, v] :=
+ frontend('diff', [Theta_B__fnd, Diff(h, y_rs[u], y_rs[v])])
+ ;
+ partial_Theta_wrt_partial_dd_h__fnd[u, v] :=
+ partial_Theta_A_wrt_partial_dd_h__fnd[u, v]/Theta_D^(3/2)
+ +
+ partial_Theta_B_wrt_partial_dd_h__fnd[u, v]/Theta_D^(1/2)
end do
end do;
- if cg_flag then codegen2(
- [partial_H_wrt_partial_d_h__fnd, partial_H_wrt_partial_dd_h__fnd],
- ['partial_H_wrt_partial_d_h', 'partial_H_wrt_partial_dd_h'],
- "../gr.cg/horizon_Jacobian.c")
+ if cg_flag then codegen2([partial_Theta_wrt_partial_d_h__fnd,
+ partial_Theta_wrt_partial_dd_h__fnd],
+ ['partial_Theta_wrt_partial_d_h', 'partial_Theta_wrt_partial_dd_h']
+ , "../gr.cg/expansion_Jacobian.c")
end if;
NULL
end proc
@@ -1025,7 +1040,7 @@ codegen2(g_uu) --> "../gr.cg/inverse_metric.c"
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=1006876, alloc=917336, time=0.09
+bytes used=1000088, alloc=917336, time=0.12
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
--> `codegen2/unindex`
@@ -1038,7 +1053,7 @@ codegen2([K, K_uu]) --> "../gr.cg/extrinsic_curvature_trace_raise.c"
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=2007252, alloc=1441528, time=0.17
+bytes used=2000352, alloc=1441528, time=0.17
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
@@ -1046,39 +1061,39 @@ bytes used=2007252, alloc=1441528, time=0.17
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
--> `codegen2/unindex`
-bytes used=3007472, alloc=1638100, time=0.25
+bytes used=3000584, alloc=1638100, time=0.22
convert p/q --> RATIONAL(p/q)
--> `codegen2/fix_rationals`
writing C code
> curvature(true);
inverse_metric_gradient...
-bytes used=4007696, alloc=1638100, time=0.33
+bytes used=4000860, alloc=1703624, time=0.28
codegen2(partial_d_g_uu) --> "../gr.cg/inverse_metric_gradient.c"
--> `codegen2/input`
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=5007880, alloc=1703624, time=0.43
+bytes used=5001020, alloc=1703624, time=0.38
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=6008088, alloc=1769148, time=0.49
+bytes used=6001324, alloc=1769148, time=0.44
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
--> `codegen2/unindex`
-bytes used=7008488, alloc=1769148, time=0.56
+bytes used=7001528, alloc=1769148, time=0.51
convert p/q --> RATIONAL(p/q)
--> `codegen2/fix_rationals`
writing C code
-bytes used=8008672, alloc=1769148, time=0.62
+bytes used=8001896, alloc=1769148, time=0.57
metric_det_gradient...
codegen2(partial_d_ln_sqrt_g) --> "../gr.cg/metric_det_gradient.c"
--> `codegen2/input`
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=9008904, alloc=1769148, time=0.75
+bytes used=9002080, alloc=1769148, time=0.64
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
@@ -1093,90 +1108,90 @@ codegen/C/expression: Unknown function: RATIONAL will be left as is.
> horizon(true);
non_unit_normal...
non_unit_normal_deriv...
-bytes used=10009180, alloc=1769148, time=0.83
-horizon_function...
-bytes used=11009748, alloc=1834672, time=0.91
-codegen2([HA, HB, HC, HD]) --> "../gr.cg/horizon_function.c"
+bytes used=10002316, alloc=1769148, time=0.69
+expansion...
+bytes used=11002748, alloc=1834672, time=0.77
+codegen2([Theta_A, Theta_B, Theta_C, Theta_D]) --> "../gr.cg/expansion.c"
--> `codegen2/input`
convert --> equation list
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=12010468, alloc=1834672, time=0.98
-bytes used=13010628, alloc=2031244, time=1.08
-bytes used=14010940, alloc=2031244, time=1.16
+bytes used=12003104, alloc=1834672, time=0.83
+bytes used=13003356, alloc=2031244, time=0.95
+bytes used=14003520, alloc=2031244, time=1.05
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=15011240, alloc=2096768, time=1.23
+bytes used=15003896, alloc=2096768, time=1.11
--> `codegen2/fix_Diff`
convert R_dd[2,3] --> R_dd_23 etc
-bytes used=16011692, alloc=2096768, time=1.31
+bytes used=16004136, alloc=2096768, time=1.18
--> `codegen2/unindex`
-bytes used=17012040, alloc=2096768, time=1.37
+bytes used=17004388, alloc=2096768, time=1.25
convert p/q --> RATIONAL(p/q)
-bytes used=18012780, alloc=2096768, time=1.44
+bytes used=18004580, alloc=2096768, time=1.31
--> `codegen2/fix_rationals`
writing C code
codegen/C/expression: Unknown function: PARTIAL_RHO will be left as is.
codegen/C/expression: Unknown function: PARTIAL_SIGMA will be left as is.
-bytes used=19012960, alloc=2096768, time=1.52
codegen/C/expression: Unknown function: PARTIAL_RHO_RHO
will be left as is.
codegen/C/expression: Unknown function: PARTIAL_RHO_SIGMA
will be left as is.
codegen/C/expression: Unknown function: PARTIAL_SIGMA_SIGMA
will be left as is.
-bytes used=20013320, alloc=2096768, time=1.66
-horizon_Jacobian...
-bytes used=21013548, alloc=2096768, time=1.73
-codegen2([partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h]) --> "../gr.cg/horizon_Jacobian.c"
+bytes used=19004732, alloc=2096768, time=1.37
+bytes used=20004972, alloc=2096768, time=1.47
+expansion_Jacobian...
+bytes used=21005184, alloc=2096768, time=1.56
+codegen2([partial_Theta_wrt_partial_d_h, partial_Theta_wrt_partial_dd_h]) --> "../gr.cg/expansion_Jacobian.c"
--> `codegen2/input`
convert --> equation list
-bytes used=22013808, alloc=2096768, time=1.79
+bytes used=22005440, alloc=2096768, time=1.63
--> `codegen2/eqnlist`
optimizing computation sequence
-bytes used=23014336, alloc=2096768, time=1.87
-bytes used=24014548, alloc=2293340, time=1.92
-bytes used=25014936, alloc=2293340, time=1.98
-bytes used=26015132, alloc=2293340, time=2.08
-bytes used=27015308, alloc=2293340, time=2.20
-bytes used=28015956, alloc=2293340, time=2.36
-bytes used=29016184, alloc=2293340, time=2.48
-bytes used=30016340, alloc=2424388, time=2.62
-bytes used=31016540, alloc=2424388, time=2.72
-bytes used=32016792, alloc=2555436, time=2.82
+bytes used=23005788, alloc=2162292, time=1.71
+bytes used=24006060, alloc=2293340, time=1.78
+bytes used=25006304, alloc=2293340, time=1.84
+bytes used=26006456, alloc=2293340, time=1.94
+bytes used=27007216, alloc=2293340, time=2.10
+bytes used=28007500, alloc=2293340, time=2.26
+bytes used=29007664, alloc=2293340, time=2.39
+bytes used=30007868, alloc=2424388, time=2.50
+bytes used=31008068, alloc=2424388, time=2.60
+bytes used=32009692, alloc=2555436, time=2.71
--> `codegen2/optimize`
find temporary variables
--> `codegen2/temps`
convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=33016952, alloc=2555436, time=2.89
-bytes used=34017300, alloc=2555436, time=2.96
-bytes used=35017536, alloc=2555436, time=3.03
+bytes used=33010028, alloc=2555436, time=2.77
+bytes used=34010432, alloc=2555436, time=2.84
+bytes used=35010740, alloc=2555436, time=2.91
--> `codegen2/fix_Diff`
+bytes used=36021172, alloc=2555436, time=2.99
convert R_dd[2,3] --> R_dd_23 etc
-bytes used=36018544, alloc=2555436, time=3.10
-bytes used=37018944, alloc=2555436, time=3.17
-bytes used=38019240, alloc=2555436, time=3.24
-bytes used=39019452, alloc=2555436, time=3.31
+bytes used=37021516, alloc=2555436, time=3.05
+bytes used=38021792, alloc=2555436, time=3.12
+bytes used=39022208, alloc=2555436, time=3.20
--> `codegen2/unindex`
-bytes used=40019836, alloc=2555436, time=3.38
-bytes used=41019988, alloc=2555436, time=3.45
-bytes used=42020192, alloc=2555436, time=3.51
-bytes used=43020452, alloc=2555436, time=3.58
+bytes used=40022676, alloc=2555436, time=3.27
+bytes used=41022924, alloc=2555436, time=3.33
+bytes used=42023096, alloc=2555436, time=3.40
+bytes used=43023256, alloc=2555436, time=3.47
+bytes used=44023688, alloc=2555436, time=3.53
convert p/q --> RATIONAL(p/q)
-bytes used=44021232, alloc=2555436, time=3.65
-bytes used=45021924, alloc=2555436, time=3.71
-bytes used=46022240, alloc=2555436, time=3.77
-bytes used=47023068, alloc=2555436, time=3.84
-bytes used=48023224, alloc=2555436, time=3.92
+bytes used=45023996, alloc=2555436, time=3.60
+bytes used=46024268, alloc=2555436, time=3.67
+bytes used=47024688, alloc=2555436, time=3.73
+bytes used=48025108, alloc=2555436, time=3.80
--> `codegen2/fix_rationals`
writing C code
-bytes used=49023460, alloc=2555436, time=3.99
-bytes used=50023664, alloc=2555436, time=4.04
-bytes used=51023844, alloc=2555436, time=4.20
-bytes used=52024240, alloc=2555436, time=4.32
-bytes used=53024400, alloc=2555436, time=4.48
-bytes used=54024568, alloc=2555436, time=4.58
+bytes used=49025260, alloc=2555436, time=3.86
+bytes used=50025744, alloc=2555436, time=3.91
+bytes used=51026008, alloc=2555436, time=4.03
+bytes used=52026168, alloc=2555436, time=4.18
+bytes used=53026328, alloc=2555436, time=4.32
+bytes used=54026584, alloc=2555436, time=4.44
> quit
-bytes used=54791456, alloc=2555436, time=4.67
+bytes used=54895808, alloc=2555436, time=4.56
diff --git a/src/gr/setup_gr_gfas.maple b/src/gr/setup_gr_gfas.maple
index ff39681..6af4aed 100644
--- a/src/gr/setup_gr_gfas.maple
+++ b/src/gr/setup_gr_gfas.maple
@@ -44,39 +44,38 @@ make_gfa('h', {inert,fnd}, [], none);
make_gfa('s_d', {inert,fnd}, [1..N], none);
make_gfa('partial_d_s_d', {inert,fnd}, [1..N, 1..N], none);
-# subexpressions for computing LHS of apparent horizon equation
-# ... these are defined by (15) in my 1996 apparent horizon finding paper
-make_gfa('HA', {inert,fnd}, [], none);
-make_gfa('HB', {inert,fnd}, [], none);
-make_gfa('HC', {inert,fnd}, [], none);
-make_gfa('HD', {inert,fnd}, [], none);
-
-# LHS of apparent horizon equation
-make_gfa('H', {inert,fnd}, [], none);
-
-# 1st derivatives of H[ABCD] and of H
-make_gfa('partial_d_HA', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_HB', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_HC', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_HD', {inert,fnd}, [1..N], none);
-make_gfa('partial_d_H', {inert,fnd}, [1..N], none);
-
-# Jacobian coefficients for H[ABCD]
-# these are the partial derivatives of H[ABCD]
+# expansion of horizon and subexpressions for computing it
+# ... these are defined by (14) and (15) in my 1996 AH-finding paper
+# except I now use Theta instead of H for the LHS
+make_gfa('Theta_A', {inert,fnd}, [], none);
+make_gfa('Theta_B', {inert,fnd}, [], none);
+make_gfa('Theta_C', {inert,fnd}, [], none);
+make_gfa('Theta_D', {inert,fnd}, [], none);
+make_gfa('Theta', {inert,fnd}, [], none);
+
+# 1st derivatives of Theta[ABCD] and of Theta
+make_gfa('partial_d_Theta_A', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta_B', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta_C', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta_D', {inert,fnd}, [1..N], none);
+make_gfa('partial_d_Theta', {inert,fnd}, [1..N], none);
+
+# Jacobian coefficients for Theta[ABCD]
+# these are the partial derivatives of Theta[ABCD]
# ... wrt Diff(h,x_rs[x])
-make_gfa('partial_HA_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_HB_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_HC_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_HD_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_A_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_B_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_C_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_D_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
# ... wrt Diff(h,x_rs[x],x_rs[y])
-make_gfa('partial_HA_wrt_partial_dd_h', {inert,fnd},
+make_gfa('partial_Theta_A_wrt_partial_dd_h', {inert,fnd},
[1..N_ang, 1..N_ang], symmetric);
-make_gfa('partial_HB_wrt_partial_dd_h', {inert,fnd},
+make_gfa('partial_Theta_B_wrt_partial_dd_h', {inert,fnd},
[1..N_ang, 1..N_ang], symmetric);
-# Jacobian coefficients for H itself
-make_gfa('partial_H_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
-make_gfa('partial_H_wrt_partial_dd_h', {inert,fnd},
+# Jacobian coefficients for Theta itself
+make_gfa('partial_Theta_wrt_partial_d_h', {inert,fnd}, [1..N_ang], none);
+make_gfa('partial_Theta_wrt_partial_dd_h', {inert,fnd},
[1..N_ang, 1..N_ang], symmetric);
NULL;
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-20 12:12:54 +0000