// gr.hh -- header file for general relativity code
// $Header$

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

//
// number of spatial dimensions in the main Cactus grid
// and in our trial-horizon-surface grid
//
enum	{ N_GRID_DIMS = 3, N_HORIZON_DIMS = 2 };

//
// this enum holds the (a) decoded  geometry_method  parameter,
// i.e. it specifies how we compute the (a) geometry
//
enum	geometry_method
	{
	geometry__interpolate_from_Cactus_grid,
	geometry__Schwarzschild_EF // no comma
	};

//
// this enum holds the (a) decoded  Jacobian_method  parameter,
// i.e. it specifies how we compute the (a) Jacobian matrix
//
enum	Jacobian_method
	{
	Jacobian_method__numerical_perturb,
	Jacobian_method__symbolic_diff_with_FD_dr,
	Jacobian_method__symbolic_diff // no comma
	};

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

//
// This structure holds all the information we need about the Cactus grid
// and gridfns.  At present we use the  CCTK_InterpLocalUniform()  local
// interpolator to access the Cactus gridfns.  Alas, much of the information
// in this structure is specific to that API, and (alas) will need changing
// when we eventually switch to the  CCTK_InterpGridArrays()  API.
//
struct	cactus_grid_info
	{
	cGH *GH;			// --> Cactus grid hierarchy

	// Cactus coordinate system
	fp coord_origin[N_GRID_DIMS];	// (x,y,z) of grid posn (0,0,0)
	fp coord_delta[N_GRID_DIMS];	// (x,y,z) grid spacing

	// dimensions of gridfn data, viewed as a 3-D array
	// n.b. storage ordering is Fortran,
	//	i.e. i is contiguous, j has stride Ni, k has stride Ni*Nj
	CCTK_INT gridfn_dims[N_GRID_DIMS];

	// this describes the semantics of the Cactus gij gridfns:
	// true ==> the Cactus gij are conformal values as per
	//          CactusEinstein/StaticConformal and the  psi  gridfn
	// false ==> the Cactus gij are the physical metric
	bool Cactus_conformal_metric;

	// --> Cactus gridfn data for grid posn (0,0,0)
	const fp* g_dd_11_data;
	const fp* g_dd_12_data;
	const fp* g_dd_13_data;
	const fp* g_dd_22_data;
	const fp* g_dd_23_data;
	const fp* g_dd_33_data;
	const fp* K_dd_11_data;
	const fp* K_dd_12_data;
	const fp* K_dd_13_data;
	const fp* K_dd_22_data;
	const fp* K_dd_23_data;
	const fp* K_dd_33_data;

	const fp* psi_data;		// NULL if !Cactus_conformal_metric
	};

//
// This structure holds information for computing the spacetime geometry.
// This is normally done by interpolating $g_{ij}$ and $K_{ij}$ from the
// usual Cactus grid, but can optionally instead by done by hard-wiring
// the Schwarzschild geometry in Eddington-Finkelstein coordinates.
//
struct	geometry_info
	{
	enum geometry_method geometry_method;

	// parameters for geometry_method = interpolate_from_Cactus_grid
	int operator_handle;		// Cactus handle to interpolation op
	int param_table_handle;		// Cactus handle to parameter table
					// for the interpolator

	// parameters for geometry_method = Schwarzschild_EF
	fp geometry__Schwarzschild_EF__x_posn;	// x posn of Schwarzschild BH
	fp geometry__Schwarzschild_EF__y_posn;	// y posn of Schwarzschild BH
	fp geometry__Schwarzschild_EF__z_posn;	// z posn of Schwarzschild BH
	fp geometry__Schwarzschild_EF__mass;	// mass of Schwarzschild BH
	fp geometry__Schwarzschild_EF__epsilon;	// use z axis limits if
						// (x^2+y^2)/r^2 <= this
	fp geometry__Schwarzschild_EF__Delta_xyz;// "grid spacing" for FD
						// computation of partial_k g_ij

	// should we check various gridfns for NaNs?
	bool check_h_for_NaNs;
	bool check_geometry_for_NaNs;
	};

//
// This struct holds parameters for computing the Jacobian matrix.
//
struct	Jacobian_info
	{
	enum Jacobian_method Jacobian_method;
	enum Jacobian_type   Jacobian_storage_method;
	fp perturbation_amplitude;
	};

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

//
// prototypes for functions visible outside their source files
//

// horizon_function.cc
// ... returns true for successful computation, false for failure
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);
enum geometry_method
  decode_geometry_method(const char geometry_method_string[]);

// horizon_Jacobian.cc
// ... returns true for successful computation, false for failure
bool horizon_Jacobian(patch_system& ps,
		      Jacobian& Jac,
		      const struct cactus_grid_info& cgi,
		      const struct geometry_info& gi,
		      const struct Jacobian_info& Jacobian_info,
		      bool print_msg_flag);
enum Jacobian_method
  decode_Jacobian_method(const char Jacobian_method_string[]);

// Schwarzschild_EF.cc
void Schwarzschild_EF_geometry(patch_system& ps,
			       const struct geometry_info& gi,
			       bool msg_flag);