* expand NaN checks in gridfns to "is finite" checks

  (also catchs +/-infinity)
* move final stage of calculation of LHS function H(h)
  out of Maple-generated code (in src/gr.cg/) into hand-written code
  in src/gr/horizon_function.c (comput_H()) so we can check that HD > 0
  before taking its sqrt.
	[HD <= 0 probably means the inerpolated g_ij
	isn't positive definite]
* also systematize CCTK_VWarn() printing for various error conditions
  in horizon_function.c


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

4 files changed, 367 insertions, 252 deletions

diff --git a/src/gr/gr.hh b/src/gr/gr.hh
index 9109e7b..3f86254 100644
--- a/src/gr/gr.hh
+++ b/src/gr/gr.hh
@@ -101,8 +101,8 @@ struct	geometry_info
 						// computation of partial_k g_ij
 
 	// should we check various gridfns for NaNs?
-	bool check_h_for_NaNs;
-	bool check_geometry_for_NaNs;
+	bool check_that_h_is_finite;
+	bool check_that_geometry_is_finite;
 	};
 
 //
@@ -122,7 +122,19 @@ struct	Jacobian_info
 //
 
 // horizon_function.cc
-// ... returns true for successful computation, false for failure
+// ... 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
+	  {
+	  enum {
+	       failure = 1,		// point outside grid etc
+	       nonfinite = 0,		// nonfinite gridfn detected
+	       skip_nonfinite = 3	// no finite() function available
+					// ==> skipping nonfinite check
+	       };
+	  }
 bool horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
 		      const struct geometry_info& gi,
diff --git a/src/gr/horizon.maple b/src/gr/horizon.maple
index e0213a7..433c723 100644
--- a/src/gr/horizon.maple
+++ b/src/gr/horizon.maple
@@ -124,19 +124,22 @@ end proc;
 ################################################################################
 
 #
-# This function computes the LHS of the apparent horizon equation
-# as a function of 1st and 2nd angular partial derivatives of the
-# apparent horizon radius $h$, in the Maple gridfn array H__fnd,
-# using equation (15) in my 1996 apparent horizon finding paper,
-# but using partial_d_g_uu[k,i,j] instead of Diff(g_uu[i,j], x_xyz[k]).
+# 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]).
 #
-# These equation gives H[ABCD] and H as functions of s_d, but here we
+# 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.
 #
-# This function also optionally generates C code for this computation.
+# 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.
 #
 # Inputs:
 #   s_d = s_d__fnd( x_xyz, X_ud, Diff(h,y_rs) )
@@ -150,7 +153,7 @@ end proc;
 #		  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) )
-#   HC = HC__fnd( x_xyz, X_ud, g_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 )
 #
@@ -201,8 +204,8 @@ 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, H__fnd],
-		 ['HA',    'HB',    'HC',    'HD',    'H'    ],
+   then codegen2([ HA__fnd, HB__fnd, HC__fnd, HD__fnd],
+		 ['HA',    'HB',    'HC',    'HD'    ],
 		 "../gr.cg/horizon_function.c");
 fi;
 
@@ -231,7 +234,7 @@ end proc;
 #		  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) )
-#   HC = HC__fnd( x_xyz, X_ud, g_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 )
 #
 # Outputs:
diff --git a/src/gr/horizon_function.cc b/src/gr/horizon_function.cc
index 92bc397..88ebfc7 100644
--- a/src/gr/horizon_function.cc
+++ b/src/gr/horizon_function.cc
@@ -9,8 +9,8 @@
 /// interpolate_geometry - interpolate g_ij and K_ij from Cactus 3-D grid
 /// convert_conformal_to_physical - convert conformal gij to physical
 ///
-/// h_contains_NaNs - does the horizon shape function h contain any NaNs?
-/// geometry_contains_NaNs - do the geometry variables contain any NaNs?
+/// 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
 ///
@@ -69,10 +69,11 @@ void Schwarzschild_EF_geometry(patch_system& ps,
 			       const struct geometry_info& gi,
 			       bool print_msg_flag);
 
-bool h_contains_NaNs(patch_system& ps, bool print_msg_flag);
-bool geometry_contains_NaNs(patch_system& ps, bool print_msg_flag);
+bool h_is_finite(patch_system& ps, bool print_msg_flag);
+bool geometry_is_finite(patch_system& ps,
+					  bool print_msg_flag);
 
-void compute_H(patch_system& ps,
+bool compute_H(patch_system& ps,
 	       bool Jacobian_flag,
 	       jtutil::norm<fp>* H_norms_ptr,
 	       bool print_msg_flag);
@@ -131,6 +132,7 @@ void 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)
 //
 bool horizon_function(patch_system& ps,
 		      const struct cactus_grid_info& cgi,
@@ -145,8 +147,8 @@ if (print_msg_flag)
 // fill in values of all ghosted gridfns in ghost zones
 ps.synchronize();
 
-if (gi.check_h_for_NaNs
-    && h_contains_NaNs(ps, print_msg_flag))
+if (gi.check_that_h_is_finite
+    && !h_is_finite(ps, print_msg_flag))
    then return false;					// *** ERROR RETURN ***
 
 // set up xyz positions of grid points
@@ -172,13 +174,14 @@ default:
 		   int(gi.geometry_method));			/*NOTREACHED*/
 	}
 
-if (gi.check_geometry_for_NaNs
-    && geometry_contains_NaNs(ps, print_msg_flag))
+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
 // by algebraic ops and angular finite differencing
-compute_H(ps, Jacobian_flag, H_norms_ptr, print_msg_flag);
+if (!compute_H(ps, Jacobian_flag, H_norms_ptr, print_msg_flag))
+   then return false;					// *** ERROR RETURN ***
 
 return true;						// *** NORMAL RETURN ***
 }
@@ -583,15 +586,19 @@ if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
 
 	if (print_msg_flag)
 	   then {
-		CCTK_VInfo(CCTK_THORNSTRING,
-		   "*** the trial-horizon-surface point");
-		CCTK_VInfo(CCTK_THORNSTRING,
-		   "***    (%g,%g,%g)",
+		CCTK_VWarn(horizon_function_warning_levels::failure,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "the trial-horizon-surface point");
+		CCTK_VWarn(horizon_function_warning_levels::failure,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "   (%g,%g,%g)",
 			   global_x, global_y, global_z);
-		CCTK_VInfo(CCTK_THORNSTRING,
-		   "*** is outside the grid (or too close to the boundary)");
-		CCTK_VInfo(CCTK_THORNSTRING,
-		   "*** in the %c%c direction!",
+		CCTK_VWarn(horizon_function_warning_levels::failure,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "is outside the grid, or too close to the boundary");
+		CCTK_VWarn(horizon_function_warning_levels::failure,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "in the %c%c direction!",
 			   end, axis);
 		}
 	return false;					// *** ERROR RETURN ***
@@ -737,18 +744,25 @@ if (print_msg_flag)
 
 //
 // This function whether the horizon shape function  h  (nominal grid only)
-// contain any NaNs.
+// contains all finite floating-point numbers, or whether it contain one
+// or more NaNs or infinities.
 //
 // Results:
-// This function returns true if any NaNs are found, false otherwise.
+#ifdef HAVE_FINITE
+// This function returns true if all the h values are finite, false
+// otherwise (i.e. if it contains any NaNs or infinities).
+#else
+// This function is a no-op, and just prints a warning message and
+// returns true.
+#endif
 //
 namespace {
-bool h_contains_NaNs(patch_system& ps, bool print_msg_flag)
+bool h_is_finite(patch_system& ps, bool print_msg_flag)
 {
 if (print_msg_flag)
-   then CCTK_VInfo(CCTK_THORNSTRING, "      checking h for NaNs");
+   then CCTK_VInfo(CCTK_THORNSTRING, "      checking that h is finite");
 
-#ifdef HAVE_ISNAN
+#ifdef HAVE_FINITE
 	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
 	{
 	patch& p = ps.ith_patch(pn);
@@ -760,23 +774,27 @@ if (print_msg_flag)
 		     ++isigma)
 		{
 		const fp h = p.ghosted_gridfn(gfns::gfn__h, irho,isigma);
-		if (isnan(h))
+		if (!finite(h))
 		   then {
 			const fp rho = p.rho_of_irho(irho);
 			const fp sigma = p.sigma_of_isigma(isigma);
-			CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-				   "h = NaN at %s patch rho=%g sigma=%g!",
-				   p.name(), double(rho), double(sigma));
-			return true;			// *** found a NaN ***
+			CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+				   __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "h=%g isn't finite at %s patch rho=%g sigma=%g!",
+				   double(h),
+				   p.name(),
+				   double(rho), double(sigma));
+			return false;			// *** found a NaN ***
 			}
 		}
 		}
 	}
-return false;						// *** no NaNs found ***
+return true;					// *** all values finite ***
 #else
-CCTK_VWarn(3, __LINE__, __FILE__, CCTK_THORNSTRING,
-           "      no isnan() ==>  skipping NaN check of h");
-return false;					// *** no check possible ***
+CCTK_VWarn(horizon_function_warning_levels::skip_nonfinite,
+	   __LINE__, __FILE__, CCTK_THORNSTRING,
+           "      no finite() fn ==>  skipping is-h-finite check");
+return true;					// *** no check possible ***
 #endif
 }
 	  }
@@ -788,18 +806,25 @@ return false;					// *** no check possible ***
 //	g_dd_ij
 //	partial_d_g_dd_kij
 //	K_dd_ij
-// contain any NaNs.
+// are all finite floating-point numbers, or whether they contain one
+// or more NaNs or infinities.
 //
 // Results:
-// This function returns true if any NaNs are found, false otherwise.
+#ifdef HAVE_FINITE
+// This function returns true if all the geometry variables are finite,
+// false otherwise (i.e. if they contain any NaNs or infinities).
+#else
+// This function is a no-op, and just prints a warning message and
+// returns true.
+#endif
 //
 namespace {
-bool geometry_contains_NaNs(patch_system& ps, bool print_msg_flag)
+bool geometry_is_finite(patch_system& ps, bool print_msg_flag)
 {
 if (print_msg_flag)
-   then CCTK_VInfo(CCTK_THORNSTRING, "      checking geometry for NaNs");
+   then CCTK_VInfo(CCTK_THORNSTRING, "      checking that geometry is finite");
 
-#ifdef HAVE_ISNAN
+#ifdef HAVE_FINITE
     for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
     {
     patch& p = ps.ith_patch(pn);
@@ -816,16 +841,17 @@ if (print_msg_flag)
 	const fp g_dd_22 = p.gridfn(gfns::gfn__g_dd_22, irho,isigma);
 	const fp g_dd_23 = p.gridfn(gfns::gfn__g_dd_23, irho,isigma);
 	const fp g_dd_33 = p.gridfn(gfns::gfn__g_dd_33, irho,isigma);
-	if (isnan(g_dd_11) || isnan(g_dd_12) || isnan(g_dd_13)
-			   || isnan(g_dd_22) || isnan(g_dd_23)
-					     || isnan(g_dd_33))
+	if (!finite(g_dd_11) || !finite(g_dd_12) || !finite(g_dd_13)
+			     || !finite(g_dd_22) || !finite(g_dd_23)
+						 || !finite(g_dd_33))
 	   then {
 		const fp rho = p.rho_of_irho(irho);
 		const fp sigma = p.sigma_of_isigma(isigma);
-		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "NaN in g_ij at %s patch rho=%g sigma=%g!",
+		CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "g_ij isn't finite at %s patch rho=%g sigma=%g!",
 			   p.name(), double(rho), double(sigma));
-		return true;			// *** found a NaN ***
+		return false;			// *** found a NaN ***
 		}
 
 	const fp partial_d_g_dd_111 = p.gridfn(gfns::gfn__partial_d_g_dd_111,
@@ -864,31 +890,32 @@ if (print_msg_flag)
 					       irho,isigma);
 	const fp partial_d_g_dd_333 = p.gridfn(gfns::gfn__partial_d_g_dd_333,
 					       irho,isigma);
-	if (    isnan(partial_d_g_dd_111)
-	     || isnan(partial_d_g_dd_112)
-	     || isnan(partial_d_g_dd_113)
-	     || isnan(partial_d_g_dd_122)
-	     || isnan(partial_d_g_dd_123)
-	     || isnan(partial_d_g_dd_133)
-	     || isnan(partial_d_g_dd_211)
-	     || isnan(partial_d_g_dd_212)
-	     || isnan(partial_d_g_dd_213)
-	     || isnan(partial_d_g_dd_222)
-	     || isnan(partial_d_g_dd_223)
-	     || isnan(partial_d_g_dd_233)
-	     || isnan(partial_d_g_dd_311)
-	     || isnan(partial_d_g_dd_312)
-	     || isnan(partial_d_g_dd_313)
-	     || isnan(partial_d_g_dd_322)
-	     || isnan(partial_d_g_dd_323)
-	     || isnan(partial_d_g_dd_333)    )
+	if (    !finite(partial_d_g_dd_111)
+	     || !finite(partial_d_g_dd_112)
+	     || !finite(partial_d_g_dd_113)
+	     || !finite(partial_d_g_dd_122)
+	     || !finite(partial_d_g_dd_123)
+	     || !finite(partial_d_g_dd_133)
+	     || !finite(partial_d_g_dd_211)
+	     || !finite(partial_d_g_dd_212)
+	     || !finite(partial_d_g_dd_213)
+	     || !finite(partial_d_g_dd_222)
+	     || !finite(partial_d_g_dd_223)
+	     || !finite(partial_d_g_dd_233)
+	     || !finite(partial_d_g_dd_311)
+	     || !finite(partial_d_g_dd_312)
+	     || !finite(partial_d_g_dd_313)
+	     || !finite(partial_d_g_dd_322)
+	     || !finite(partial_d_g_dd_323)
+	     || !finite(partial_d_g_dd_333)    )
 	   then {
 		const fp rho = p.rho_of_irho(irho);
 		const fp sigma = p.sigma_of_isigma(isigma);
-		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "NaN in partial_k g_ij at %s patch rho=%g sigma=%g!",
+		CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+		   "partial_k g_ij isn't finite at %s patch rho=%g sigma=%g!",
 			   p.name(), double(rho), double(sigma));
-		return true;			// *** found a NaN ***
+		return false;			// *** found a NaN ***
 		}
 
 	const fp K_dd_11 = p.gridfn(gfns::gfn__K_dd_11, irho,isigma);
@@ -897,25 +924,27 @@ if (print_msg_flag)
 	const fp K_dd_22 = p.gridfn(gfns::gfn__K_dd_22, irho,isigma);
 	const fp K_dd_23 = p.gridfn(gfns::gfn__K_dd_23, irho,isigma);
 	const fp K_dd_33 = p.gridfn(gfns::gfn__K_dd_33, irho,isigma);
-	if (isnan(K_dd_11) || isnan(K_dd_12) || isnan(K_dd_13)
-			   || isnan(K_dd_22) || isnan(K_dd_23)
-					     || isnan(K_dd_33))
+	if (!finite(K_dd_11) || !finite(K_dd_12) || !finite(K_dd_13)
+			     || !finite(K_dd_22) || !finite(K_dd_23)
+						 || !finite(K_dd_33))
 	   then {
 		const fp rho = p.rho_of_irho(irho);
 		const fp sigma = p.sigma_of_isigma(isigma);
-		CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
-			   "NaN in K_ij at %s patch rho=%g sigma=%g!",
+		CCTK_VWarn(horizon_function_warning_levels::nonfinite,
+			   __LINE__, __FILE__, CCTK_THORNSTRING,
+			   "K_ij isn't finite at %s patch rho=%g sigma=%g!",
 			   p.name(), double(rho), double(sigma));
-		return true;			// *** found a NaN ***
+		return false;			// *** found a NaN ***
 		}
 	}
 	}
     }
-return false;						// *** no NaNs found ***
+return true;						// *** no NaNs found ***
 #else
-CCTK_VWarn(3, __LINE__, __FILE__, CCTK_THORNSTRING,
-           "      no isnan() ==>  skipping NaN check of geometry");
-return false;					// *** no check possible ***
+CCTK_VWarn(horizon_function_warning_levels::skip_nonfinite,
+	   __LINE__, __FILE__, CCTK_THORNSTRING,
+           "      no finite() ==>  skipping is-geometry-finite check");
+return true;					// *** no check possible ***
 #endif
 }
 	  }
@@ -937,8 +966,13 @@ return false;					// *** no check possible ***
 // Jacobian_flag = true to compute the Jacobian coefficients,
 //		   false to skip this.
 //
+// Results:
+// This function returns true for a successful computation, or false
+// if the computation failed because HD <= 0 (this means the interpolated
+// g_ij isn't positive definite).
+//
 namespace {
-void compute_H(patch_system& ps,
+bool compute_H(patch_system& ps,
 	       bool Jacobian_flag,
 	       jtutil::norm<fp>* H_norms_ptr,
 	       bool print_msg_flag)
@@ -1018,10 +1052,29 @@ if (print_msg_flag)
 		  }
 
 		  {
-		// HA, HB, HC, HD, H
+		// HA, HB, HC, HD
 		#include "../gr.cg/horizon_function.c"
 		  }
 
+		if (HD <= 0)
+		   then {
+			CCTK_VWarn(horizon_function_warning_levels::failure,
+				   __LINE__, __FILE__, CCTK_THORNSTRING,
+				   "HD <= 0 at %s patch rho=%g sigma=%g!",
+				   p.name(), double(rho), double(sigma));
+			CCTK_VWarn(horizon_function_warning_levels::failure,
+				   __LINE__, __FILE__, CCTK_THORNSTRING,
+				   "(this probably means the interpolated g_ij");
+			CCTK_VWarn(horizon_function_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;
+
 		// update running norms of H(h) function
 		if (H_norms_ptr != NULL)
 		   then H_norms_ptr->data(H);
@@ -1034,5 +1087,7 @@ if (print_msg_flag)
 		}
 		}
 	}
+
+return true;						// *** NORMAL RETURN ***
 }
 	  }
diff --git a/src/gr/maple.log b/src/gr/maple.log
index 568cbe6..0a28bb5 100644
--- a/src/gr/maple.log
+++ b/src/gr/maple.log
@@ -4,7 +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 $
+# $Header: /numrelcvs/AEIDevelopment/AHFinderDirect/src/gr/doit.maple,v 1.5 2002/09/13 14:12:18 jthorn Exp $
 > 
 > read "../maple/setup.mm";
 msum := proc(fn::algebraic)
@@ -205,8 +205,19 @@ 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;
+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;
 option remember;
     var_list := [args[2 .. nargs]];
     if type(operand, indexed) and op(0, operand) = 'X_ud' and
@@ -575,16 +586,19 @@ 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;
+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;
     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);
@@ -603,6 +617,11 @@ partial_H_wrt_partial_dd_h, partial_H_wrt_partial_dd_h__fnd;
     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)
@@ -629,16 +648,19 @@ 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;
+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;
 option remember;
     var_list := [args[2 .. nargs]];
     if type(operand, indexed) and op(0, operand) = 'g_dd' and
@@ -661,16 +683,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_dd);
     assert_fnd_exists(g_uu, fnd);
@@ -688,16 +713,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_uu);
     assert_fnd_exists(K_dd);
@@ -728,16 +756,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_dd);
     assert_fnd_exists(g_uu);
@@ -762,16 +793,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_dd);
     assert_fnd_exists(g_uu);
@@ -801,16 +835,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(h);
     assert_fnd_exists(X_ud);
@@ -828,16 +865,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(h);
     assert_fnd_exists(X_ud);
@@ -864,16 +904,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_uu);
     assert_fnd_exists(K_uu);
@@ -898,8 +941,8 @@ partial_H_wrt_partial_dd_h, partial_H_wrt_partial_dd_h__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([HA__fnd, HB__fnd, HC__fnd, HD__fnd, H__fnd],
-        ['HA', 'HB', 'HC', 'HD', 'H'], "../gr.cg/horizon_function.c")
+    if cg_flag then codegen2([HA__fnd, HB__fnd, HC__fnd, HD__fnd],
+        ['HA', 'HB', 'HC', 'HD'], "../gr.cg/horizon_function.c")
     end if;
     NULL
 end proc
@@ -911,16 +954,19 @@ 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;
+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;
     printf("%a...\n", procname);
     assert_fnd_exists(g_uu);
     assert_fnd_exists(K_uu);
@@ -979,7 +1025,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=1009876, alloc=917336, time=0.09
+bytes used=1006876, alloc=917336, time=0.09
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
@@ -992,7 +1038,7 @@ codegen2([K, K_uu]) --> "../gr.cg/extrinsic_curvature_trace_raise.c"
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=2010048, alloc=1441528, time=0.16
+bytes used=2007252, alloc=1441528, time=0.17
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -1000,39 +1046,39 @@ bytes used=2010048, alloc=1441528, time=0.16
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=3010268, alloc=1638100, time=0.23
+bytes used=3007472, alloc=1638100, time=0.25
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
 > curvature(true);
 inverse_metric_gradient...
-bytes used=4010500, alloc=1638100, time=0.30
+bytes used=4007696, alloc=1638100, time=0.33
 codegen2(partial_d_g_uu) --> "../gr.cg/inverse_metric_gradient.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-bytes used=5010984, alloc=1638100, time=0.37
+bytes used=5007880, alloc=1703624, time=0.43
       --> `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
+bytes used=6008088, alloc=1769148, time=0.49
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
       --> `codegen2/unindex`
-bytes used=7011816, alloc=1769148, time=0.50
+bytes used=7008488, alloc=1769148, time=0.56
    convert p/q --> RATIONAL(p/q)
       --> `codegen2/fix_rationals`
    writing C code
-bytes used=8012448, alloc=1769148, time=0.56
+bytes used=8008672, alloc=1769148, time=0.62
 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=9012608, alloc=1769148, time=0.64
+bytes used=9008904, alloc=1769148, time=0.75
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
@@ -1047,91 +1093,90 @@ codegen/C/expression:   Unknown function:   RATIONAL   will be left as is.
 > horizon(true);
 non_unit_normal...
 non_unit_normal_deriv...
-bytes used=10012896, alloc=1769148, time=0.71
+bytes used=10009180, alloc=1769148, time=0.83
 horizon_function...
-bytes used=11013596, alloc=1834672, time=0.79
-codegen2([HA, HB, HC, HD, H]) --> "../gr.cg/horizon_function.c"
+bytes used=11009748, alloc=1834672, time=0.91
+codegen2([HA, HB, HC, HD]) --> "../gr.cg/horizon_function.c"
       --> `codegen2/input`
    convert --> equation list
       --> `codegen2/eqnlist`
    optimizing computation sequence
-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
+bytes used=12010468, alloc=1834672, time=0.98
+bytes used=13010628, alloc=2031244, time=1.08
+bytes used=14010940, alloc=2031244, time=1.16
       --> `codegen2/optimize`
    find temporary variables
       --> `codegen2/temps`
    convert Diff(expr,rho,sigma) --> PARTIAL_RHO_SIGMA(expr) etc
-bytes used=16015736, alloc=2162292, time=1.21
+bytes used=15011240, alloc=2096768, time=1.23
       --> `codegen2/fix_Diff`
    convert R_dd[2,3] --> R_dd_23 etc
-bytes used=17016120, alloc=2162292, time=1.27
+bytes used=16011692, alloc=2096768, time=1.31
       --> `codegen2/unindex`
-bytes used=18016512, alloc=2162292, time=1.34
+bytes used=17012040, alloc=2096768, time=1.37
    convert p/q --> RATIONAL(p/q)
-bytes used=19017260, alloc=2162292, time=1.40
+bytes used=18012780, alloc=2096768, time=1.44
       --> `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=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=21017636, alloc=2162292, time=1.62
+bytes used=20013320, alloc=2096768, time=1.66
 horizon_Jacobian...
-bytes used=22135056, alloc=2227816, time=1.70
+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"
       --> `codegen2/input`
    convert --> equation list
-bytes used=23135472, alloc=2227816, time=1.77
+bytes used=22013808, alloc=2096768, time=1.79
       --> `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
+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
       --> `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
+bytes used=33016952, alloc=2555436, time=2.89
+bytes used=34017300, alloc=2555436, time=2.96
+bytes used=35017536, alloc=2555436, time=3.03
       --> `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
+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
       --> `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
+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
    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
+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
       --> `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
+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
 > quit
-bytes used=56126644, alloc=2620960, time=4.42
+bytes used=54791456, alloc=2555436, time=4.67