// driver.cc -- top level driver for finding apparent horizons
// $Id$
//
// <<<prototypes for functions local to this file>>>
// AHFinderDirect_driver - top-level driver
/// setup_ellipsoid_initial_guess - set up ellipsoid in h gridfn
//

#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 "jt/stdc.h"
#include "jt/util.hh"
#include "jt/array.hh"
#include "jt/cpm_map.hh"
#include "jt/linear_map.hh"
using jtutil::error_exit;

#include "../config.hh"

#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 "gfn.hh"
#include "AHFinderDirect.hh"

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

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

namespace {
void setup_ellipsoid_initial_guess
	(patch_system& ps,
	 fp global_center_x, fp global_center_y, fp global_center_z,
	 fp radius_x, fp radius_y, fp radius_z);
	  };

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

//
// 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];
cgi.g_dd_11_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::gxx")
					  );
cgi.g_dd_12_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::gxy")
					  );
cgi.g_dd_13_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::gxz")
					  );
cgi.g_dd_22_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::gyy")
					  );
cgi.g_dd_23_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::gyz")
					  );
cgi.g_dd_33_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::gzz")
					  );
cgi.K_dd_11_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::kxx")
					  );
cgi.K_dd_12_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::kxy")
					  );
cgi.K_dd_13_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::kxz")
					  );
cgi.K_dd_22_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::kyy")
					  );
cgi.K_dd_23_data = static_cast<const fp *>(
			CCTK_VarDataPtr(cctkGH, 0, "einstein::kyz")
					  );
cgi.K_dd_33_data = static_cast<const 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, "ellipsoid"))
   then {
	setup_ellipsoid_initial_guess(ps,
				      initial_guess__ellipsoid__center_global_x,
				      initial_guess__ellipsoid__center_global_y,
				      initial_guess__ellipsoid__center_global_z,
				      initial_guess__ellipsoid__radius_x,
				      initial_guess__ellipsoid__radius_y,
				      initial_guess__ellipsoid__radius_z);
	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");
	}
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown method=\"%s\"!",
		   method);					/*NOTREACHED*/
}

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

//
// 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_initial_guess
	(patch_system& ps,
	 fp global_center_x, fp global_center_y, fp global_center_z,
	 fp radius_x, fp radius_y, fp radius_z)
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "   h = ellipsoid: global_center=(%g,%g,%g)",
	   global_center_x, global_center_y, global_center_z);
CCTK_VInfo(CCTK_THORNSTRING,
	   "                  radius=(%g,%g,%g)",
	   radius_x, radius_y, radius_z);

	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 = global_center_x - ps.origin_x();
		const fp BV = global_center_y - ps.origin_y();
		const fp CW = global_center_z - ps.origin_z();
		const fp a = radius_x;
		const fp b = radius_y;
		const fp c = radius_z;

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