// horizon_radius.cc -- provide horizon radius & other info for other thorns
// $Header$
//
// <<<access to persistent data>>>
// AHFinderDirect_local_coordinate_origin - provide our local coordinate origin
// AHFinderDirect_horizon_was_found - query if a given horizon was found
// AHFinderDirect_radius_in_direction - provide r(angle) function
//

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

#include "util_Table.h"
#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"
#include "cctk_Functions.h"

#include "config.h"
#include "stdc.h"
#include "../jtutil/util.hh"
#include "../jtutil/array.hh"
#include "../jtutil/cpm_map.hh"
#include "../jtutil/linear_map.hh"
using jtutil::error_exit;

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

#include "../elliptic/Jacobian.hh"

#include "../gr/gfns.hh"
#include "../gr/gr.hh"

#include "horizon_sequence.hh"
#include "BH_diagnostics.hh"
#include "driver.hh"

// all the code in this file is inside this namespace
namespace AHFinderDirect
	  {

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

//
// ***** access to persistent data *****
//
extern struct state state;

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

//
// This function is called (via the Cactus flesh function-aliasing mechanism)
// by other thorns to find out our local coordinate origin for a given AH.
//
// Results:
// This function returns 0 for ok, or -1 if the horizon number is invalid.
//
extern "C"
  CCTK_INT AHFinderDirect_local_coordinate_origin
    (CCTK_INT horizon_number,
     CCTK_REAL* origin_x_ptr, CCTK_REAL* origin_y_ptr, CCTK_REAL* origin_z_ptr)
{
const struct verbose_info& verbose_info = state.verbose_info;

if (!  ((horizon_number >= 1) && (horizon_number <= state.N_horizons))  )
   then {
	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_local_coordinate_origin():\n"
"        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
		   ,
		   int(horizon_number), state.N_horizons);
	return -1;					// *** ERROR RETURN ***
	}

assert(state.AH_data_array[horizon_number] != NULL);
const struct AH_data& AH_data = *state.AH_data_array[horizon_number];

assert(AH_data.ps_ptr != NULL);
const patch_system& ps = *AH_data.ps_ptr;

*origin_x_ptr = ps.origin_x();
*origin_y_ptr = ps.origin_y();
*origin_z_ptr = ps.origin_z();

return 0;						// *** NORMAL RETURN ***
}

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

//
// This function is called (via the Cactus flesh function-aliasing mechanism)
// by other thorns to query whether or not the specified horizon was found
// the last time we searched for it.
//
// Results:
// This function returns
//  1 if the horizon was found
//  0 if the horizon was not found
//  -1 if the horizon number is invalid.
//
extern "C"
  CCTK_INT AHFinderDirect_horizon_was_found(CCTK_INT horizon_number)
{
const struct verbose_info& verbose_info = state.verbose_info;

if (!  ((horizon_number >= 1) && (horizon_number <= state.N_horizons))  )
   then {
	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_horizon_was_found():\n"
"        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
		   ,
		   int(horizon_number), state.N_horizons);
	return -1;					// *** ERROR RETURN ***
	}

assert(state.AH_data_array[horizon_number] != NULL);
const struct AH_data& AH_data = *state.AH_data_array[horizon_number];

return AH_data.found_flag ? 1 : 0;
}

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

//
// This function is called (via the Cactus flesh function-aliasing mechanism)
// by other thorns to find out a given AH's radius in the direction from
// its local coordinate origin to a given (x,y,z) coordinate or coordinates.
//
// FIXME:
// The current implementation is quite inefficient: it does a full 2-D
// interpolation (to find the horizon radius) for each of the xyz points.
// It would be more efficient to batch the interpolations.
//
// Arguments:
// horizon_number = must be in the range 1 to N_horizons
// N_points = should be >= 0
// x[], y[], z[] = these give the (x,y,z) coordinates
// radius[] = this is set to the horizon radius values (Euclidean distance
//	      from the local coordinate origin), or to all -1.0 if we didn't
//	      find this horizon the last time we looked for it
//
// Results:
// This function returns 0 for ok, or -1 if the horizon number is invalid.
//
extern "C"
  CCTK_INT AHFinderDirect_radius_in_direction
    (CCTK_INT horizon_number,
     CCTK_INT N_points,
     const CCTK_REAL* const x, const CCTK_REAL* const y, const CCTK_REAL* const z,
     CCTK_REAL* const radius)
{
const struct verbose_info& verbose_info = state.verbose_info;

if (! ((horizon_number >= 1) && (horizon_number <= state.N_horizons)) )
   then {
	CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
"AHFinderDirect_distance_outside_thorn():\n"
"        horizon_number=%d must be in the range [1,N_horizons=%d]!\n"
		   ,
		   int(horizon_number), state.N_horizons);
	return -1;					// *** ERROR RETURN ***
	}

assert(state.AH_data_array[horizon_number] != NULL);
const struct AH_data& AH_data = *state.AH_data_array[horizon_number];

	for (int point = 0 ; point < N_points ; ++point)
	{
	if (AH_data.found_flag)
	   then {
		assert(AH_data.ps_ptr != NULL);
		const patch_system& ps = *AH_data.ps_ptr;

		const fp local_x = x[point] - ps.origin_x();
		const fp local_y = y[point] - ps.origin_y();
		const fp local_z = z[point] - ps.origin_z();
		radius[point] = ps.radius_in_local_xyz_direction
				      (gfns::gfn__h,
				       local_x, local_y, local_z);
		}
	   else radius[point] = -1.0;
	}

return 0;						// *** NORMAL RETURN ***
}

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

	  }	// namespace AHFinderDirect