changes so I can compile the full horizon test

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@629 f88db872-0e4f-0410-b76b-b9085cfa78c5

4 files changed, 258 insertions, 59 deletions

diff --git a/src/gr/driver.cc b/src/gr/driver.cc
index e95c783..3982be2 100644
--- a/src/gr/driver.cc
+++ b/src/gr/driver.cc
@@ -16,15 +16,14 @@
 #include "cctk_Arguments.h"
 #include "cctk_Parameters.h"
 
-#include "jt/stdc.h"
-#include "jt/util.hh"
-#include "jt/array.hh"
-#include "jt/cpm_map.hh"
-#include "jt/linear_map.hh"
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../config.hh"
-
 #include "../util/coords.hh"
 #include "../util/grid.hh"
 #include "../util/fd_grid.hh"
diff --git a/src/gr/horizon_Jacobian.cc b/src/gr/horizon_Jacobian.cc
index 39db258..23390ab 100644
--- a/src/gr/horizon_Jacobian.cc
+++ b/src/gr/horizon_Jacobian.cc
@@ -14,15 +14,14 @@
 #include "cctk.h"
 #include "cctk_Arguments.h"
 
-#include "jt/stdc.h"
-#include "jt/util.hh"
-#include "jt/array.hh"
-#include "jt/cpm_map.hh"
-#include "jt/linear_map.hh"
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../config.hh"
-
 #include "../util/coords.hh"
 #include "../util/grid.hh"
 #include "../util/fd_grid.hh"
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 400516b..50e9254 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -17,15 +17,14 @@
 #include "cctk.h"
 #include "cctk_Arguments.h"
 
-#include "jt/stdc.h"
-#include "jt/util.hh"
-#include "jt/array.hh"
-#include "jt/cpm_map.hh"
-#include "jt/linear_map.hh"
+#include "stdc.h"
+#include "config.hh"
+#include "../jtutil/util.hh"
+#include "../jtutil/array.hh"
+#include "../jtutil/cpm_map.hh"
+#include "../jtutil/linear_map.hh"
 using jtutil::error_exit;
 
-#include "../config.hh"
-
 #include "../util/coords.hh"
 #include "../util/grid.hh"
 #include "../util/fd_grid.hh"
@@ -49,7 +48,7 @@ void setup_xyz_posns(patch_system& ps);
 void interpolate_geometry(patch_system& ps,
 			  const struct cactus_grid_info& cgi,
 			  const struct geometry_interpolator_info& ii);
-void compute_H(patch_system& ps);
+void compute_H(patch_system& ps, bool Jacobian_flag);
 	  }
 
 //******************************************************************************
diff --git a/src/gr/maple.log b/src/gr/maple.log
index 2e5ddd0..568cbe6 100644
--- a/src/gr/maple.log
+++ b/src/gr/maple.log
@@ -4,6 +4,7 @@
  <____ ____>  Waterloo Maple Inc.
       |       Type ? for help.
 # top-level Maple file to read/run all code in this directory
+# $Id: doit.maple,v 1.4 2002/07/16 10:26:13 jthorn Exp $
 > 
 > read "../maple/setup.mm";
 msum := proc(fn::algebraic)
@@ -575,7 +576,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     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);
@@ -589,15 +598,31 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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('H', {inert, fnd}, [], 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_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},
+        [1 .. N_ang, 1 .. N_ang], symmetric);
+    make_gfa('partial_HB_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},
+        [1 .. N_ang, 1 .. N_ang], symmetric);
     NULL
 end proc
 
-Diff/gridfn := proc(operand)
+Diff/gridfn2 := proc(operand)
 local var_list, posn;
 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,
@@ -605,7 +630,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
 option remember;
     var_list := [args[2 .. nargs]];
     if type(operand, indexed) and op(0, operand) = 'g_dd' and
@@ -629,7 +662,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_dd);
     assert_fnd_exists(g_uu, fnd);
@@ -648,7 +689,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_uu);
     assert_fnd_exists(K_dd);
@@ -667,8 +716,8 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
     NULL
 end proc
 
-> read "metric_derivs.mm";
-metric_derivs := proc(cg_flag::boolean)
+> read "curvature.mm";
+curvature := proc(cg_flag::boolean)
     inverse_metric_gradient(cg_flag); metric_det_gradient(cg_flag); NULL
 end proc
 
@@ -680,7 +729,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_dd);
     assert_fnd_exists(g_uu);
@@ -706,7 +763,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_dd);
     assert_fnd_exists(g_uu);
@@ -725,6 +790,7 @@ horizon := proc(cg_flag::boolean)
     non_unit_normal();
     non_unit_normal_deriv();
     horizon_function(cg_flag);
+    horizon_Jacobian(cg_flag);
     NULL
 end proc
 
@@ -736,7 +802,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     printf("%a...\n", procname);
     assert_fnd_exists(h);
     assert_fnd_exists(X_ud);
@@ -755,7 +829,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     printf("%a...\n", procname);
     assert_fnd_exists(h);
     assert_fnd_exists(X_ud);
@@ -783,7 +865,15 @@ 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, H, H__fnd,
-HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_uu);
     assert_fnd_exists(K_uu);
@@ -808,7 +898,68 @@ HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__fnd;
     H__fnd :=
         HA__fnd/HD__fnd^(3/2) + HB__fnd/HD__fnd^(1/2) + HC__fnd/HD__fnd - K
         ;
-    if cg_flag then codegen2(H__fnd, 'H', "../gr.cg/horizon_function.c")
+    if cg_flag then codegen2([HA__fnd, HB__fnd, HC__fnd, HD__fnd, H__fnd],
+        ['HA', 'HB', 'HC', 'HD', 'H'], "../gr.cg/horizon_function.c")
+    end if;
+    NULL
+end proc
+
+horizon_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, H, H__fnd,
+HA, HA__fnd, HB, HB__fnd, HC, HC__fnd, HD, HD__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;
+    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);
+    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
+    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)
+        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")
     end if;
     NULL
 end proc
@@ -828,7 +979,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=1000148, alloc=917336, time=0.12
+bytes used=1009876, alloc=917336, time=0.09
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
@@ -841,7 +992,7 @@ codegen2([K, K_uu]) --> "../gr.cg/extrinsic_curvature_trace_raise.c"
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=2000796, alloc=1441528, time=0.22
+bytes used=2010048, alloc=1441528, time=0.16
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -849,39 +1000,39 @@ bytes used=2000796, alloc=1441528, time=0.22
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=3001776, alloc=1572576, time=0.28
+bytes used=3010268, alloc=1638100, time=0.23
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
-> metric_derivs(true);
+> curvature(true);
 inverse_metric_gradient...
-bytes used=4002196, alloc=1638100, time=0.36
+bytes used=4010500, alloc=1638100, time=0.30
 codegen2(partial_d_g_uu) --> "../gr.cg/inverse_metric_gradient.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=5002380, alloc=1638100, time=0.43
+bytes used=5010984, alloc=1638100, time=0.37
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
+bytes used=6011336, alloc=1769148, time=0.43
       --> `codegen2/fix_Diff`
-bytes used=6004052, alloc=1703624, time=0.50
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=7004228, alloc=1703624, time=0.57
+bytes used=7011816, alloc=1769148, time=0.50
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
-bytes used=8004548, alloc=1703624, time=0.63
+bytes used=8012448, alloc=1769148, time=0.56
 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=9004804, alloc=1703624, time=0.74
+bytes used=9012608, alloc=1769148, time=0.64
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -896,40 +1047,91 @@ codegen/C/expression:   Unknown function:   RATIONAL   will be left as is.
 > horizon(true);
 non_unit_normal...
 non_unit_normal_deriv...
-bytes used=10004956, alloc=1769148, time=0.85
+bytes used=10012896, alloc=1769148, time=0.71
 horizon_function...
-bytes used=11005112, alloc=1834672, time=0.93
-codegen2(H) --> "../gr.cg/horizon_function.c"
+bytes used=11013596, alloc=1834672, time=0.79
+codegen2([HA, HB, HC, HD, H]) --> "../gr.cg/horizon_function.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=12005688, alloc=1834672, time=1.01
-bytes used=13005864, alloc=2031244, time=1.12
-bytes used=14006148, alloc=2031244, time=1.21
+bytes used=12014936, alloc=1834672, time=0.86
+bytes used=13015120, alloc=2162292, time=0.92
+bytes used=14015272, alloc=2162292, time=1.04
+bytes used=15015576, alloc=2162292, time=1.14
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=15006312, alloc=2096768, time=1.28
+bytes used=16015736, alloc=2162292, time=1.21
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
-bytes used=16006476, alloc=2096768, time=1.35
+bytes used=17016120, alloc=2162292, time=1.27
       --> `codegen2/unindex`
-bytes used=17006632, alloc=2096768, time=1.41
+bytes used=18016512, alloc=2162292, time=1.34
    convert p/q --> RATIONAL(p/q)
-bytes used=18006900, alloc=2096768, time=1.48
+bytes used=19017260, alloc=2162292, time=1.40
       --> `codegen2/fix_rationals`
    writing C code
+bytes used=20017484, alloc=2162292, time=1.46
 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=19007116, alloc=2096768, time=1.54
 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=20007280, alloc=2096768, time=1.66
+bytes used=21017636, alloc=2162292, time=1.62
+horizon_Jacobian...
+bytes used=22135056, alloc=2227816, time=1.70
+codegen2([partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h]) --> "../gr.cg/horizon_Jacobian.c"
+      --> `codegen2/input`
+   convert --> equation list
+bytes used=23135472, alloc=2227816, time=1.77
+      --> `codegen2/eqnlist`
+   optimizing computation sequence
+bytes used=24135760, alloc=2227816, time=1.84
+bytes used=25136784, alloc=2293340, time=1.91
+bytes used=26136972, alloc=2293340, time=1.97
+bytes used=27137268, alloc=2293340, time=2.05
+bytes used=28137440, alloc=2293340, time=2.15
+bytes used=29137604, alloc=2293340, time=2.27
+bytes used=30138036, alloc=2293340, time=2.41
+bytes used=31138284, alloc=2424388, time=2.53
+bytes used=32138448, alloc=2424388, time=2.62
+bytes used=33138692, alloc=2489912, time=2.73
+      --> `codegen2/optimize`
+   find temporary variables
+      --> `codegen2/temps`
+   convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
+bytes used=34138856, alloc=2620960, time=2.80
+bytes used=35139568, alloc=2620960, time=2.86
+bytes used=36139824, alloc=2620960, time=2.93
+bytes used=37140068, alloc=2620960, time=2.99
+      --> `codegen2/fix_Diff`
+   convert R_dd[2,3] --> R_dd_23 etc
+bytes used=38140236, alloc=2620960, time=3.06
+bytes used=39140408, alloc=2620960, time=3.13
+bytes used=40141148, alloc=2620960, time=3.19
+      --> `codegen2/unindex`
+bytes used=41141672, alloc=2620960, time=3.26
+bytes used=42141908, alloc=2620960, time=3.32
+bytes used=43142740, alloc=2620960, time=3.39
+bytes used=44142896, alloc=2620960, time=3.45
+bytes used=45143236, alloc=2620960, time=3.52
+   convert p/q --> RATIONAL(p/q)
+bytes used=46143764, alloc=2620960, time=3.58
+bytes used=47143924, alloc=2620960, time=3.64
+bytes used=48144500, alloc=2620960, time=3.70
+bytes used=49144756, alloc=2620960, time=3.77
+      --> `codegen2/fix_rationals`
+   writing C code
+bytes used=50145012, alloc=2620960, time=3.82
+bytes used=51145184, alloc=2620960, time=3.87
+bytes used=52145484, alloc=2620960, time=3.99
+bytes used=53145800, alloc=2620960, time=4.10
+bytes used=54146020, alloc=2620960, time=4.21
+bytes used=55146332, alloc=2620960, time=4.31
 > quit
-bytes used=20238456, alloc=2096768, time=1.70
+bytes used=56126644, alloc=2620960, time=4.42