// driver.cc -- top level driver for finding apparent horizons
// $Id$
//
// <<<prototypes for functions local to this file>>>
// AHFinderDirect_driver - top-level driver
///
/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
/// setup_ellipsoid - setup up a coordiante ellipsoid in h
///

#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 "../elliptic/Jacobian.hh"

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

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

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

namespace {
void setup_Kerr_horizon(patch_system& ps,
			fp x_center, fp y_center, fp z_center,
			fp m, fp a,
			bool Kerr_Schild_flag);
void setup_ellipsoid(patch_system& ps,
		     fp x_center, fp y_center, fp z_center,
		     fp x_radius, fp y_radius, fp z_radius);
	  };

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

//
// 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, "ADMBase::gxx"));
cgi.g_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxy"));
cgi.g_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gxz"));
cgi.g_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyy"));
cgi.g_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gyz"));
cgi.g_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::gzz"));
cgi.K_dd_11_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxx"));
cgi.K_dd_12_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxy"));
cgi.K_dd_13_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kxz"));
cgi.K_dd_22_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyy"));
cgi.K_dd_23_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kyz"));
cgi.K_dd_33_data = static_cast<fp*>(CCTK_VarDataPtr(cctkGH, 0, "ADMBase::kzz"));


//
// create the patch system and initialize the xyz derivative coefficients
//
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, "ellipsoid"))
   then setup_ellipsoid(ps,
			initial_guess__ellipsoid__x_center,
			initial_guess__ellipsoid__y_center,
			initial_guess__ellipsoid__z_center,
			initial_guess__ellipsoid__x_radius,
			initial_guess__ellipsoid__y_radius,
			initial_guess__ellipsoid__z_radius);
else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr"))
   then setup_Kerr_horizon(ps,
			   initial_guess__Kerr_KerrSchild__x_center,
			   initial_guess__Kerr_KerrSchild__y_center,
			   initial_guess__Kerr_KerrSchild__z_center,
			   initial_guess__Kerr_KerrSchild__mass,
			   initial_guess__Kerr_KerrSchild__spin,
			   false);	// use Kerr coords
else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
   then setup_Kerr_horizon(ps,
			   initial_guess__Kerr_KerrSchild__x_center,
			   initial_guess__Kerr_KerrSchild__y_center,
			   initial_guess__Kerr_KerrSchild__z_center,
			   initial_guess__Kerr_KerrSchild__mass,
			   initial_guess__Kerr_KerrSchild__spin,
			   true);	// use Kerr-Schild coords
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown initial_guess_method=\"%s\"!",
		   initial_guess_method);			/*NOTREACHED*/
ps.print_ghosted_gridfn_with_xyz(ghosted_gfns::gfn__h,
				 true, ghosted_gfns::gfn__h,
				 "h.dat",
				 false);	// no ghost zones


//
// find the apparent horizon
//
jtutil::norm<fp> H_norms;
if      (STRING_EQUAL(method, "horizon function"))
   then {
	horizon_function(ps, cgi, gii, false, H_norms);
	CCTK_VInfo(CCTK_THORNSTRING,
		   "   H(h) rms-norm %.2e, infinity-norm %.2e\n",
		   H_norms.rms_norm(), H_norms.infinity_norm());
	ps.print_gridfn_with_xyz(nominal_gfns::gfn__H,
				 true, ghosted_gfns::gfn__h,
				 "H.dat");
	}
else if (STRING_EQUAL(method, "Jacobian"))
   then {
	Jacobian& Jac = create_Jacobian(ps, Jacobian_type);
	horizon_function(ps, cgi, gii, true, H_norms);
	horizon_Jacobian_SD(ps, Jac);
	print_Jacobian(Jacobian_file_name, Jac);
	}
else if (STRING_EQUAL(method, "Jacobian test"))
   then {
	Jacobian& SD_Jac = create_Jacobian(ps, Jacobian_type);
	horizon_function(ps, cgi, gii, true, H_norms);
	horizon_Jacobian_SD(ps, SD_Jac);

	Jacobian& NP_Jac = create_Jacobian(ps, Jacobian_type);
	horizon_function(ps, cgi, gii, true, H_norms);
	horizon_Jacobian_NP(ps, cgi, gii,
			    NP_Jac,
			    NP_Jacobian__perturbation_amplitude);

	print_Jacobians(Jacobian_file_name, SD_Jac, NP_Jac);
	}
else if (STRING_EQUAL(method, "Newton (NP Jacobian)"))
   then {
	Newton_solve(ps, cgi, gii,
		     Jacobian_type,
		     NP_Jacobian__perturbation_amplitude,
		     max_Newton_iterations,
		     H_norm_for_convergence);
	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 method=\"%s\"!",
		   method);					/*NOTREACHED*/
}

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

//
// This function sets up the horizon of a Kerr black hole in Kerr or
// 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.2 of my AHFinderDirect notes.
//
// Arguments:
// [xyz]_center = The position of the Kerr black hole.
// (m,a) = Describe the Kerr black hole.  Note that my convention has
//	   a=J/m^2 dimensionless, while MTW take a=J/m=m*(my a).
// Kerr_Schild_flag = false to use Kerr coordinates,
//		      true to use Kerr-Schild coordinates
//
namespace {
void setup_Kerr_horizon(patch_system& ps,
			fp x_center, fp y_center, fp z_center,
			fp m, fp a,
			bool Kerr_Schild_flag)
{
const char* const name = Kerr_Schild_flag ? "Kerr-Schild" : "Kerr";

CCTK_VInfo(CCTK_THORNSTRING,
	   "setting up horizon for Kerr in %s coords",
	   name);
CCTK_VInfo(CCTK_THORNSTRING,
	   "   mass=%g, spin=J/m^2=%g, posn=(%g,%g,%g)",
	   double(m), double(a),
	   double(x_center), double(y_center), double(z_center));

// horizon in Kerr coordinates is coordinate sphere
const fp r = m * (1.0 + sqrt(1.0 - a*a));

// horizon in Kerr-Schild coordinates is coordinate ellipsoid
const fp  z_radius = r;
const fp xy_radius = Kerr_Schild_flag ? r * sqrt(1.0 + a*a*m*m/(r*r)) : r;

CCTK_VInfo(CCTK_THORNSTRING,
	   "   horizon is coordinate %s",
	   Kerr_Schild_flag ? "ellipsoid" : "sphere");
setup_ellipsoid(ps,
		x_center, y_center, z_center,
		xy_radius, xy_radius, z_radius);
}
	  }

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

//
// This function sets up an ellipsoid in the gridfn h, using the
// formulas in "ellipsoid.maple" and the Maple-generated C code in
// "ellipsoid.c":
//
// ellipsoid has center (A,B,C), radius (a,b,c)
// angular coordinate system has center (U,V,W)
//
// direction cosines wrt angular coordinate center are (xcos,ycos,zcos)
// i.e. a point has coordinates (U+xcos*r, V+ycos*r, W+zcos*r)
//
// then the equation of the ellipsoid is
//	(U+xcos*r - A)^2     (V+ycos*r - B)^2     (W+zcos*r - C)^2
//	-----------------  +  ----------------  +  -----------------  =  1
//	        a^2                  b^2                   c^2
//
// to solve this, we introduce intermediate variables
//	AU = A - U
//	BV = B - V
//	CW = C - W
//
namespace {
void setup_ellipsoid(patch_system& ps,
		     fp x_center, fp y_center, fp z_center,
		     fp x_radius, fp y_radius, fp z_radius)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "setting h = ellipsoid: center=(%g,%g,%g)",
	   double(x_center), double(y_center), double(z_center));
CCTK_VInfo(CCTK_THORNSTRING,
	   "                       radius=(%g,%g,%g)",
	   double(x_radius), double(y_radius), double(z_radius));

	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 xcos, ycos, zcos;
		p.xyzcos_of_rho_sigma(rho,sigma, xcos,ycos,zcos);

		// set up variables used by Maple-generated code
		const fp AU = x_center - ps.origin_x();
		const fp BV = y_center - ps.origin_y();
		const fp CW = z_center - ps.origin_z();
		const fp a = x_radius;
		const fp b = y_radius;
		const fp c = z_radius;

		// compute the solutions r_plus and r_minus
		fp r_plus, r_minus;
		#include "ellipsoid.c"

		// exactly one of the solutions (call it r) should be positive
		fp r;
		if      ((r_plus > 0.0) && (r_minus < 0.0))
		   then r = r_plus;
		else if ((r_plus < 0.0) && (r_minus > 0.0))
		   then r = r_minus;
		else    CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
				   "\n"
"   expected exactly one r>0 solution to quadratic, got 0 or 2!\n"
"   %s patch (irho,isigma)=(%d,%d) ==> (rho,sigma)=(%g,%g)\n"
"   direction cosines (xcos,ycos,zcos)=(%g,%g,%g)\n"
"   ==> r_plus=%g r_minus=%g\n"
				   ,
				   p.name(), irho, isigma,
				   double(rho), double(sigma),
				   double(xcos), double(ycos), double(zcos),
				   double(r_plus), double(r_minus));
		   						/*NOTREACHED*/

		// r = horizon radius at this grid point
		p.ghosted_gridfn(ghosted_gfns::gfn__h, irho,isigma) = r;
		}
		}
	}
}
	  }