// driver.cc -- top level driver for finding apparent horizons
// $Id$
//
// <<<prototypes for functions local to this file>>>
// AHFinderDirect_driver - top-level driver
/// setup_Kerr_KerrSchild_initial_guess - set up Kerr/Kerr-Schild initial guess
//

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

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.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;

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

#include "gfn.hh"
#include "AHFinderDirect.hh"

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

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

namespace {
void setup_Kerr_KerrSchild_horizon(patch_system& ps,
				   fp mass, fp spin,
				   fp xposn, fp yposn, fp zposn);
	  };

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

//
// This function is the Cactus interface for the test driver.
//
extern "C"
  void AHFinderDirect_driver(CCTK_ARGUMENTS)
{
DECLARE_CCTK_ARGUMENTS
DECLARE_CCTK_PARAMETERS

CCTK_VInfo(CCTK_THORNSTRING, "initializing AHFinderDirect data structures");


//
// set up the geometry interpolator
//
struct geometry_interpolator_info gii;
CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
gii.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
if (gii.operator_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "couldn't find interpolator \"%s\"!",
		   geometry_interpolator_name);		/*NOTREACHED*/

gii.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
if (gii.param_table_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "bad geometry-interpolator parameter(s) \"%s\"!",
		   geometry_interpolator_pars);		/*NOTREACHED*/


//
// set up the interpatch interpolator
//
CCTK_VInfo(CCTK_THORNSTRING, "   setting up interpatch interpolator");
const int interp_handle = CCTK_InterpHandle(interpatch_interpolator_name);
if (interp_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "couldn't find interpolator \"%s\"!",
		   interpatch_interpolator_name);		/*NOTREACHED*/
const int interp_param_table_handle
	= Util_TableCreateFromString(interpatch_interpolator_pars);
if (interp_param_table_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "bad interpatch-interpolator parameter(s) \"%s\"!",
		   interpatch_interpolator_pars);		/*NOTREACHED*/


//
// set up the Cactus grid info
//
CCTK_VInfo(CCTK_THORNSTRING, "   setting up Cactus grid info");
struct cactus_grid_info cgi;
cgi.GH = cctkGH;
cgi.coord_origin[0] = cctk_origin_space[0];
cgi.coord_origin[1] = cctk_origin_space[1];
cgi.coord_origin[2] = cctk_origin_space[2];
cgi.coord_delta[0] = cctk_delta_space[0];
cgi.coord_delta[1] = cctk_delta_space[1];
cgi.coord_delta[2] = cctk_delta_space[2];
cgi.gridfn_dims[0] = cctk_lsh[0];
cgi.gridfn_dims[1] = cctk_lsh[1];
cgi.gridfn_dims[2] = cctk_lsh[2];
// n.b. The  cgi.[gK]_dd_??_data  are actually  const fp *  pointers,
//	since we won't modify the 3-D gridfn data!  But  static_cast<...>
//      won't change const modifiers, so we just cast to  fp*  and let
//	the assignment take care of the const part...
cgi.g_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxx"));
cgi.g_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxy"));
cgi.g_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gxz"));
cgi.g_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyy"));
cgi.g_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gyz"));
cgi.g_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::gzz"));
cgi.K_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxx"));
cgi.K_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxy"));
cgi.K_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kxz"));
cgi.K_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyy"));
cgi.K_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kyz"));
cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH,0, "einstein::kzz"));


//
// create the patch system and initialize the xyz derivative coefficients
//
CCTK_VInfo(CCTK_THORNSTRING, "   creating patch system");
patch_system ps(origin_x, origin_y, origin_z,
		patch_system::type_of_name(patch_system_type),
		N_ghost_points, N_overlap_points, delta_drho_dsigma,
		nominal_gfns::min_gfn, nominal_gfns::max_gfn,
		ghosted_gfns::min_gfn, ghosted_gfns::max_gfn,
		interp_handle, interp_param_table_handle);


//
// set up the initial guess for the apparent horizon shape
//
if	(STRING_EQUAL(initial_guess_method, "read from file"))
   then {
	CCTK_VInfo(CCTK_THORNSTRING,
		   "   reading initial guess from \"%s\"",
		   initial_guess__read_from_file__file_name);
	ps.read_ghosted_gridfn(ghosted_gfns::gfn__h,
			       initial_guess__read_from_file__file_name,
			       false);		// no ghost zones
	}
else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
   then {
	setup_Kerr_KerrSchild_horizon(ps,
				      initial_guess__Kerr_KerrSchild__mass,
				      initial_guess__Kerr_KerrSchild__spin,
				      initial_guess__Kerr_KerrSchild__xposn,
				      initial_guess__Kerr_KerrSchild__yposn,
				      initial_guess__Kerr_KerrSchild__zposn);
	ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
					 true, ghosted_gfns::gfn__h,
					 "h.dat",
					 false);	// no ghost zones
	}
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown initial_guess_method=\"%s\"!",
		   initial_guess_method);			/*NOTREACHED*/


//
// find the apparent horizon
//
if      (STRING_EQUAL(method, "horizon"))
   then {
	horizon_function(ps, cgi, gii, false);
	ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
				 true, ghosted_gfns::gfn__h,
				 "H.dat");
	}
if      (STRING_EQUAL(method, "horizon Jacobian"))
   then {
	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
	horizon_function(ps, cgi, gii, true);
	horizon_Jacobian(ps, Jac);
	}
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown method=\"%s\"!",
		   method);					/*NOTREACHED*/
}

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

//
// This function sets up the horizon of a Kerr black hole in Kerr-Schild
// coordinates, on the nominal grid, in the  h  gridfn.
//
// Kerr-Schild coordinates are described in MTW Exercise 33.8, page 903,
// and the horizon is worked out on page 13 of my AHFinderDirect notes.
//
// Arguments:
// (mass,spin) = Describe the Kerr black hole.
// [xyz]_posn = The position of the Kerr black hole.
//
namespace {
void setup_Kerr_KerrSchild_horizon(patch_system& ps,
				   fp mass, fp spin,
				   fp xposn, fp yposn, fp zposn)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "   setting up Kerr/Kerr-Schild horizon");
CCTK_VInfo(CCTK_THORNSTRING,
	   "      mass=%g, spin=%g, posn=(%g,%g,%g)",
	   double(mass), double(spin),
	   double(xposn), double(yposn), double(zposn));

// horizon in Kerr-Schild is coordinate sphere $r = (1 + \sqrt{1-a^2}) m$
const fp r = (1.0 + sqrt(1.0 - spin*spin)) * mass;

	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 rho   = p.rho_of_irho(irho);
		const fp sigma = p.sigma_of_isigma(isigma);

		fp theta, phi;
		p.theta_phi_of_rho_sigma(rho,sigma, theta,phi);
		fp Kerr_x, Kerr_y, Kerr_z;
		p.xyz_of_r_rho_sigma(r,rho,sigma, Kerr_x,Kerr_y,Kerr_z);

		const fp KS_x = Kerr_x - spin*sin(theta)*sin(phi);
		const fp KS_y = Kerr_y + spin*sin(theta)*cos(phi);
		const fp KS_z = Kerr_z;
		p.ghosted_gridfn(ghosted_gfns::gfn__h,irho,isigma)
			= jtutil::hypot3(KS_x, KS_y, KS_z);
		}
		}
	}
}
	  }