path: root/src/gr/horizon_function.cc

// horizon_function.cc -- evaluate LHS function H(h)
// $Header$
//
// <<<prototypes for functions local to this file>>>
// horizon_function - top-level driver
/// setup_xyz_posns - setup global xyz posns of grid points
/// interpolate_geometry - interpolate g_ij and K_ij from Cactus 3-D grid
/// compute_H - compute H(h) given earlier setup
///

#include <stdio.h>
#include <assert.h>
#include <math.h>
#include <vector>

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"

#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;
using jtutil::pow2;
using jtutil::pow4;

#include "../patch/coords.hh"
#include "../patch/grid.hh"
#include "../patch/fd_grid.hh"
#include "../patch/patch.hh"
#include "../patch/patch_edge.hh"
#include "../patch/patch_interp.hh"
#include "../patch/ghost_zone.hh"
#include "../patch/patch_system.hh"

#include "../elliptic/Jacobian.hh"

#include "gfns.hh"
#include "gr.hh"

//******************************************************************************

//
// ***** prototypes for functions local to this file *****
//

namespace {
void setup_xyz_posns(patch_system& ps, bool print_msg_flag);
bool interpolate_geometry(patch_system& ps,
			  const struct cactus_grid_info& cgi,
			  const struct geometry_info& gi,
			  bool print_msg_flag);
void Schwarzschild_EF_geometry(patch_system& ps,
			       const struct cactus_grid_info& cgi,
			       const struct geometry_info& gi,
			       bool print_msg_flag);
void compute_H(patch_system& ps,
	       bool Jacobian_flag,
	       jtutil::norm<fp>* H_norms_ptr,
	       bool print_msg_flag);
	  }

//******************************************************************************

//
// This function computes the LHS function H(h), and optionally also
// its Jacobian coefficients (from which the Jacobian matrix may be
// computed later).
//
// Inputs (angular gridfns, on ghosted grid):
// ... defined on ghosted grid
// ... only values on nominal grid are actually used as input
//	h				# shape of trial surface
//
// Inputs (Cactus 3-D gridfns):
//	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
//	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
//
// Outputs (temporaries computed at each grid point)
//	## computed by hand-written code
//	global_[xyz]			# xyz positions of grid points
//	X_ud_*, X_udd_*			# xyz derivative coefficients
//	## computed by Maple-generated code
//	g_uu_{11,12,13,22,23,33}	# $g^{ij}$
//	K				# $K$
//	K_dd_{11,12,13,22,23,33}	# $K^{ij}$
//	partial_d_ln_sqrt_g_d		# $\partial_i \ln \sqrt{g}$
//	partial_d_g_uu_{1,2,3}{11,12,13,22,23,33}	# $\partial_k g^{ij}$
//
// Outputs (angular gridfns, all on the nominal grid):
//	## interpolated from 3-D Cactus grid
//	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)$
//
// 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.
//
// Results:
// This function returns true for a successful computation, or false
// if the computation failed because the geometry interpolation would
// need data outside the Cactus grid, or data from an excised region.
// FIXME: excision isn't implemented yet :(
//
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)
{
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING, "   horizon function");

// fill in values of all ghosted gridfns in ghost zones
ps.synchronize();

// set up xyz positions of grid points
setup_xyz_posns(ps, print_msg_flag);

// compute the "geometry" g_ij, K_ij, and partial_k g_ij
switch	(gi.geometry_method)
	{
case geometry__interpolate_from_Cactus_grid:
	if (! interpolate_geometry(ps, cgi, gi, print_msg_flag))
	   then return false;				// *** ERROR RETURN ***
	break;
case geometry__Schwarzschild_EF:
	Schwarzschild_EF_geometry(ps, gi, print_msg_flag);
	break;
default:
	error_exit(PANIC_EXIT,
"***** horizon_function(): unknown gi.geometry_method=(int)%d!\n"
"                          (this should never happen!)\n"
,
		   int(gi.geometry_method));			/*NOTREACHED*/
	}

// 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);

return true;						// *** NORMAL RETURN ***
}

//******************************************************************************

//
// This function sets up the global xyz positions of the grid points
// in the gridfns global_[xyz].  These will be used by interplate_geometry().
//
namespace {
void setup_xyz_posns(patch_system& ps, bool print_msg_flag)
{
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "      xyz positions and derivative coefficients");

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
		{
		for (int isigma = p.min_isigma() ;
		     isigma <= p.max_isigma() ;
		     ++isigma)
		{
		const fp r = p.ghosted_gridfn(gfns::gfn__h, irho,isigma);
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);

		fp local_x, local_y, local_z;
		p.xyz_of_r_rho_sigma(r,rho,sigma, local_x,local_y,local_z);

		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);

		p.gridfn(gfns::gfn__global_x, irho,isigma) = global_x;
		p.gridfn(gfns::gfn__global_y, irho,isigma) = global_y;
		p.gridfn(gfns::gfn__global_z, irho,isigma) = global_z;
		}
		}
	}
}
	  }

//******************************************************************************

//
// This function interpolates the Cactus gridfns $g_{ij}$ and $K_{ij}$
// to determine $g_{ij}$, $K_{ij}$, and the spatial derivatives
// $\partial_k g_{ij}$, on the trial horizon surface position given by h.
//
// Inputs (angular gridfns, on ghosted grid):
// ... defined on ghosted grid
// ... only values on nominal grid are actually used as input
//	h				# shape of trial surface
//
// Inputs (angular gridfns, all on the nominal grid):
//	global_[xyz]			# xyz positions of grid points
//
// Inputs (Cactus 3-D gridfns):
//	gxx,gxy,gxz,gyy,gyz,gzz		# 3-metric $g_{ij}$
//	kxx,kxy,kxz,kyy,kyz,kzz		# extrinsic curvature $K_{ij}$
//
// Outputs (angular gridfns, all on the nominal grid):
//	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}$
//
// On the first call this function also modifies the interpolator
// parameter table.
//
// Results:
// This function returns true for a successful computation, or false
// if the computation failed because the geometry interpolation would
// need data outside the Cactus grid, or data from an excised region.
// FIXME: excision isn't implemented yet :(
//
namespace {
bool interpolate_geometry(patch_system& ps,
			  const struct cactus_grid_info& cgi,
			  const struct geometry_info& gi,
			  bool print_msg_flag)
{
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "      interpolating g_ij and Kij from Cactus 3-D grid");

int status;

const int N_interp_points = ps.N_grid_points();
const int interp_coords_type_code = CCTK_VARIABLE_REAL;
const void* const interp_coords[N_GRID_DIMS]
  = {
    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_x)),
    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_y)),
    static_cast<const void*>(ps.gridfn_data(gfns::gfn__global_z)),
    };

const int N_INPUT_ARRAYS = 12;
const CCTK_INT input_array_type_codes[N_INPUT_ARRAYS]
	= {
	  // g_ij
	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
			      CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
						  CCTK_VARIABLE_REAL,
	  // K_ij
	  CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
			      CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
						  CCTK_VARIABLE_REAL,
	  };
const void* const input_arrays[N_INPUT_ARRAYS]
	= {
	  static_cast<const void*>(cgi.g_dd_11_data),
	  static_cast<const void*>(cgi.g_dd_12_data),
	  static_cast<const void*>(cgi.g_dd_13_data),
	  static_cast<const void*>(cgi.g_dd_22_data),
	  static_cast<const void*>(cgi.g_dd_23_data),
	  static_cast<const void*>(cgi.g_dd_33_data),
	  static_cast<const void*>(cgi.K_dd_11_data),
	  static_cast<const void*>(cgi.K_dd_12_data),
	  static_cast<const void*>(cgi.K_dd_13_data),
	  static_cast<const void*>(cgi.K_dd_22_data),
	  static_cast<const void*>(cgi.K_dd_23_data),
	  static_cast<const void*>(cgi.K_dd_33_data),
	  };

const int N_OUTPUT_ARRAYS = 30;

const CCTK_INT output_array_type_codes[N_OUTPUT_ARRAYS]
	= {
 // g_ij             partial_x g_ij      partial_y g_ij      partial_z g_ij
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
 // $K_{ij}$
 CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
		     CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
					 CCTK_VARIABLE_REAL,
	  };
const CCTK_INT operand_indices[N_OUTPUT_ARRAYS]
	= {
	  0, 0, 0, 0,		// g_dd_11
	  1, 1, 1, 1,		// g_dd_12
	  2, 2, 2, 2,		// g_dd_13
	  3, 3, 3, 3,		// g_dd_22
	  4, 4, 4, 4,		// g_dd_23
	  5, 5, 5, 5,		// g_dd_33
	   6,  7,  8,		// K_dd_{11,12,13}
	       9, 10,		// K_dd_{22,23}
		  11,		// K_dd_33
	  };
#define DERIV(x)	x
const CCTK_INT operation_codes[N_OUTPUT_ARRAYS]
	= {
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_11
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_12
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_13
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_22
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_23
	  DERIV(0), DERIV(1), DERIV(2), DERIV(3),	// g_dd_33
	  DERIV(0), DERIV(0), DERIV(0),			// K_dd_{11,12,13}
		    DERIV(0), DERIV(0),			// K_dd_{22,23}
			      DERIV(0),			// K_dd_{33}
	  };
void* const output_arrays[N_OUTPUT_ARRAYS]
  = {
    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_11)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_111)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_211)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_311)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_12)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_112)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_212)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_312)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_13)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_113)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_213)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_313)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_22)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_122)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_222)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_322)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_23)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_123)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_223)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_323)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__g_dd_33)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_133)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_233)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__partial_d_g_dd_333)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_11)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_12)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_13)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_22)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_23)),
    static_cast<void*>(ps.gridfn_data(gfns::gfn__K_dd_33)),
    };

bool first_time = true;
if (first_time)
   then {
	first_time = false;

	// store derivative info in interpolator parameter table
	if (print_msg_flag)
	   then CCTK_VInfo(CCTK_THORNSTRING,
			   "         setting up interpolator derivative info");

	status = Util_TableSetIntArray(gi.param_table_handle,
				       N_OUTPUT_ARRAYS, operand_indices,
				       "operand_indices");
	if (status < 0)
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        unable to set operand_indices in interpolator parameter table!\n"
"        Util_TableSetIntArray() status=%d\n"
,
			   status);				/*NOTREACHED*/

	status = Util_TableSetIntArray(gi.param_table_handle,
				       N_OUTPUT_ARRAYS, operation_codes,
				       "operation_codes");
	if (status < 0)
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        unable to set operation_codes in interpolator parameter table!\n"
"        Util_TableSetIntArray() status=%d\n"
,
			   status);				/*NOTREACHED*/
	}

if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING,
		   "         calling interpolator (%d points)",
		   N_interp_points);
status = CCTK_InterpLocalUniform(N_GRID_DIMS,
				 gi.operator_handle, gi.param_table_handle,
				 cgi.coord_origin, cgi.coord_delta,
				 N_interp_points,
				    interp_coords_type_code,
				    interp_coords,
				 N_INPUT_ARRAYS,
				    cgi.gridfn_dims,
				    input_array_type_codes,
				    input_arrays,
				 N_OUTPUT_ARRAYS,
				    output_array_type_codes,
				    output_arrays);
if (status == CCTK_ERROR_INTERP_POINT_X_RANGE)
   then {
	// look in interpolator output table entries
	// to see *which* point is out-of-range
	CCTK_INT out_of_range_pt, out_of_range_axis, out_of_range_end;
	if (    (Util_TableGetInt(gi.param_table_handle,
				  &out_of_range_pt,
				  "out_of_range_pt") < 0)
	     || (Util_TableGetInt(gi.param_table_handle,
				  &out_of_range_axis,
				  "out_of_range_axis") < 0)
	     || (Util_TableGetInt(gi.param_table_handle,
				  &out_of_range_end,
				  "out_of_range_end") < 0)    )
	   then error_exit(ERROR_EXIT,
"***** interpolate_geometry():\n"
"        point out of range when interpolating geometry info from 3-D grid!\n"
"        ==> the trial horizon surface is (at least partially)\n"
"            outside the grid and/or in an excised region!\n"
"        (unable to get info about which point is out of range:\n"
"         maybe an interpolator problem?)\n");			/*NOTREACHED*/

	assert(out_of_range_pt >= 0);
	assert(out_of_range_pt < ps.N_grid_points());
	const double global_x
	   = ps.gridfn_data(gfns::gfn__global_x)[out_of_range_pt];
	const double global_y
	   = ps.gridfn_data(gfns::gfn__global_y)[out_of_range_pt];
	const double global_z
	   = ps.gridfn_data(gfns::gfn__global_z)[out_of_range_pt];

	assert(out_of_range_axis >= 0);
	assert(out_of_range_axis < N_GRID_DIMS);
	const char axis = "xyz"[out_of_range_axis];

	assert((out_of_range_end == -1) || (out_of_range_end == +1));
	const char end = (out_of_range_end == -1) ? '-' : '+';

	CCTK_VInfo(CCTK_THORNSTRING,
		   "*** the trial-horizon-surface point");
	CCTK_VInfo(CCTK_THORNSTRING,
		   "***    (%g,%g,%g)",
		   global_x, global_y, global_z);
	CCTK_VInfo(CCTK_THORNSTRING,
		   "*** is outside the grid (or too to the grid boundary)");
	CCTK_VInfo(CCTK_THORNSTRING,
		   "*** in the %c%c direction!",
		   end, axis);
	return false;					// *** ERROR RETURN ***
	}
if (status < 0)
   then error_exit(ERROR_EXIT,
"***** interpolate_geometry(): error return from interpolator!\n"
"        CCTK_InterpLocalUniform() status=%d\n"
,
		   status);					/*NOTREACHED*/

return true;						// *** NORMAL RETURN ***
}
	  }

//******************************************************************************

//
// 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.
//
// Arguments:
// Jacobian_flag = true to compute the Jacobian coefficients,
//		   false to skip this.
//
namespace {
void compute_H(patch_system& ps,
	       bool Jacobian_flag,
	       jtutil::norm<fp>* H_norms_ptr,
	       bool print_msg_flag)
{
if (print_msg_flag)
   then CCTK_VInfo(CCTK_THORNSTRING, "      computing H(h)");

	for (int pn = 0 ; pn < ps.N_patches() ; ++pn)
	{
	patch& p = ps.ith_patch(pn);

		for (int irho = p.min_irho() ; irho <= p.max_irho() ; ++irho)
		{
		for (int isigma = p.min_isigma() ;
		     isigma <= p.max_isigma() ;
		     ++isigma)
		{
		//
		// compute the X_ud and X_udd derivative coefficients
		// ... n.b. this uses the *local* (x,y,z) coordinates
		//
		const fp r = p.ghosted_gridfn(gfns::gfn__h, irho,isigma);
		const fp rho = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);
		fp xx, yy, zz;
		p.xyz_of_r_rho_sigma(r,rho,sigma, xx, yy, zz);

		// 1st derivative coefficients X_ud
		const fp X_ud_11 = p.partial_rho_wrt_x(xx, yy, zz);
		const fp X_ud_12 = p.partial_rho_wrt_y(xx, yy, zz);
		const fp X_ud_13 = p.partial_rho_wrt_z(xx, yy, zz);
		const fp X_ud_21 = p.partial_sigma_wrt_x(xx, yy, zz);
		const fp X_ud_22 = p.partial_sigma_wrt_y(xx, yy, zz);
		const fp X_ud_23 = p.partial_sigma_wrt_z(xx, yy, zz);

		// 2nd derivative coefficient gridfns X_udd
		const fp X_udd_111 = p.partial2_rho_wrt_xx(xx, yy, zz);
		const fp X_udd_112 = p.partial2_rho_wrt_xy(xx, yy, zz);
		const fp X_udd_113 = p.partial2_rho_wrt_xz(xx, yy, zz);
		const fp X_udd_122 = p.partial2_rho_wrt_yy(xx, yy, zz);
		const fp X_udd_123 = p.partial2_rho_wrt_yz(xx, yy, zz);
		const fp X_udd_133 = p.partial2_rho_wrt_zz(xx, yy, zz);
		const fp X_udd_211 = p.partial2_sigma_wrt_xx(xx, yy, zz);
		const fp X_udd_212 = p.partial2_sigma_wrt_xy(xx, yy, zz);
		const fp X_udd_213 = p.partial2_sigma_wrt_xz(xx, yy, zz);
		const fp X_udd_222 = p.partial2_sigma_wrt_yy(xx, yy, zz);
		const fp X_udd_223 = p.partial2_sigma_wrt_yz(xx, yy, zz);
		const fp X_udd_233 = p.partial2_sigma_wrt_zz(xx, yy, zz);

		//
		// "call" the Maple-generated code
		// ... each cg/*.c file has a separate set of temp variables,
		//     and so must be inside its own set of { } braces
		//

		// gridfn #defins
		#include "cg.hh"

		  {
		// g_uu
		#include "../gr.cg/inverse_metric.c"
		  }

		  {
		// K, K_uu
		#include "../gr.cg/extrinsic_curvature_trace_raise.c"
		  }

		  {
		// partial_d_g_uu
		#include "../gr.cg/inverse_metric_gradient.c"
		  }

		  {
		// partial_d_ln_sqrt_g
		#include "../gr.cg/metric_det_gradient.c"
		  }

		  {
		// HA, HB, HC, HD, H
		#include "../gr.cg/horizon_function.c"
		  }

		// update running norms of H(h) function
		if (H_norms_ptr != NULL)
		   then H_norms_ptr->data(H);

		if (Jacobian_flag)
		   then {
			// partial_H_wrt_partial_d_h, partial_H_wrt_partial_dd_h
			#include "../gr.cg/horizon_Jacobian.c"
			}
		}
		}
	}
}
	  }