// AHFinderDirect.hh -- misc global-within-this-thorn stuff
// $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 structure holds various information about the interpolator
// used to interpolate the "geometry" Cactus 3-D gridfns ($g_{ij}$ and
// $K_{ij}$) to the apparent horizon position.
//
struct	geometry_interpolator_info
	{
	int operator_handle;		// Cactus handle to interpolation op
	int param_table_handle;		// Cactus handle to parameter table
					// for the interpolator
	};

//
// This structure holds all the information we need about the Cactus grid
// and gridfns outside the top-level driver.  At present we use the
//  CCTK_InterpLocalUniform()  local interpolator to access the Cactus
// gridfns.  Much of the information in this structure is specific to
// that API, and (alas) will need changing when we eventually switch to
// a "global" multiprocessor/distributed interpolator.
//
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];

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

//
// This struct holds the various parameters used in computing the Jacobian
// matrix.
//
struct	Jacobian_info
	{
	fp perturbation_amplitude;
	};

//
// This struct holds the various parameters used in solving the H(h) = 0
// equations.
//
struct	solver_info
	{
	// stuff for Newton_solve()
	int max_Newton_iterations;
	fp H_norm_for_convergence;
	};

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

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

// driver.cc
extern "C"
  void AHFinderDirect_driver(CCTK_ARGUMENTS);

// horizon_function.cc
void horizon_function(patch_system& ps,
		      const struct cactus_grid_info& cgi,
		      const struct geometry_interpolator_info& gii,
		      bool Jacobian_flag,
		      jtutil::norm<fp>* H_norms_ptr,
		      bool msg_flag = true);

// horizon_Jacobian.cc
Jacobian& create_Jacobian(patch_system& ps,
			  const char Jacobian_type[]);
void horizon_Jacobian(patch_system& ps,
		      const struct cactus_grid_info& cgi,
		      const struct geometry_interpolator_info& gii,
		      const struct Jacobian_info& Jac_info,
		      Jacobian& Jac);

// Newton.cc
void Newton_solve(patch_system& ps,
		  const struct cactus_grid_info& cgi,
		  const struct geometry_interpolator_info& gii,
		  const char Jacobian_type[],
		  fp perturbation_amplitude,
		  int max_Newton_iterations,
		  fp H_norm_for_convergence);