// 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 enum holds the (a) decoded  Jacobian_method  parameter,
// i.e. it specifies how we compute the (a) Jacobian matrix
//
enum	Jacobian_method
	{
	numerical_perturbation,
	symbolic_differentiation_with_FD_dr,
	symbolic_differentiation // no comma
	};


//
// 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
	{
	//
	// parameters for hard-wiring Schwarzschild/EF geometry
	//
	bool hardwire_Schwarzschild_EF;		// should we hard-wire the
						// Schwarzschild/EF geometry?
	fp hardwire_Schwarzschild_EF__x_posn;	// x posn of Schwarzschild BH
	fp hardwire_Schwarzschild_EF__y_posn;	// y posn of Schwarzschild BH
	fp hardwire_Schwarzschild_EF__z_posn;	// z posn of Schwarzschild BH
	fp hardwire_Schwarzschild_EF__mass;	// mass of Schwarzschild BH
	fp hardwire_Schwarzschild_EF__epsilon;	// threshold for sin^2 theta
						// = (x^2+y^2)/r^2 below which
						// we use z axis limits
	fp Delta_xyz;		// pseudo-grid spacing for finite differencing
				// computation of $\partial_k g_{ij}$

	//
	// parameters for normal interpolation from Cactus grid
	//
	int operator_handle;		// Cactus handle to interpolation op
	int param_table_handle;		// Cactus handle to parameter table
					// for the interpolator

	// this doesn't really belong in this structure (it doesn't
	// have any logical connection to the geometry calculations),
	// but it's convenient to put it here anyway...
	};

//
// 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 parameters for computing the Jacobian matrix.
//
struct	Jacobian_info
	{
	enum Jacobian_method Jacobian_method;
	fp perturbation_amplitude;
	};

//
// This struct holds parameters for solving the H(h) = 0 equations.
//
struct	solver_info
	{
	// stuff for Newton_solve()
	int max_Newton_iterations;
	fp H_norm_for_convergence;
	fp Delta_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_info& gi,
		      bool Jacobian_flag = false,
		      bool msg_flag = false,
		      jtutil::norm<fp>* H_norms_ptr = NULL);

// horizon_Jacobian.cc
void horizon_Jacobian(patch_system& ps,
		      const struct cactus_grid_info& cgi,
		      const struct geometry_info& gi,
		      const struct Jacobian_info& Jacobian_info,
		      Jacobian& Jac);

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

// Newton.cc
// return true for success, false for failure to converge
bool Newton_solve(patch_system& ps,
		  const struct cactus_grid_info& cgi,
		  const struct geometry_info& gi,
		  const struct Jacobian_info& Jacobian_info,
		  const struct solver_info& solver_info,
		  Jacobian& Jac);