// BH_diagnostics.cc -- compute/print BH diagnostics
// $Header$
//
// BH_diagnostics::BH_diagnostics - initialize a  struct BH_diagnostics
//
// BH_diagnostics::compute - compute BH diagnostics after an AH has been found
/// BH_diagnostics::surface_integral - integrate gridfn over the 2-sphere
//
// print - print a line or two summarizing the diagnostics
// setup_output_file - create/open output file, write header describing fields
// output - write a (long) line of all the diagnostics
//

#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 "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"

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

//
// This function initializes a  struct BH_diagnostics  to all zeros.
//
BH_diagnostics::BH_diagnostics()
	: centroid_x(0.0), centroid_y(0.0), centroid_z(0.0),
	  min_radius(0.0), max_radius(0.0), mean_radius(0.0),
	  circumference_xy(0.0), circumference_xz(0.0), circumference_yz(0.0),
	  area(0.0),
	  m_irreducible(0.0)
{ }

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

//
// Given that an apparent horizon has been found, this function computes
// various black hole diagnostics.
//
// Inputs (gridfns)
// h		# ghosted
// one		# nominal
// global_[xyz]	# nominal
//
void BH_diagnostics::compute
	(const patch_system& ps,
	 const struct BH_diagnostics_info& BH_diagnostics_info)
{
//
// min/max radius of horizon
//
jtutil::norm<fp> h_norms;
ps.ghosted_gridfn_norms(gfns::gfn__h, h_norms);
min_radius = h_norms.min_abs_value();
max_radius = h_norms.max_abs_value();

//
// surface integrals
//
const fp integral_one = surface_integral(ps,
					 gfns::gfn__one, true, true, true,
					 BH_diagnostics_info.integral_method);
const fp integral_h = surface_integral(ps,
				       gfns::gfn__h, true, true, true,
				       BH_diagnostics_info.integral_method);
const fp integral_x = surface_integral(ps,
				       gfns::gfn__global_x, true, true, false,
				       BH_diagnostics_info.integral_method);
const fp integral_y = surface_integral(ps,
				       gfns::gfn__global_y, true, false, true,
				       BH_diagnostics_info.integral_method);
const fp integral_z = surface_integral(ps,
				       gfns::gfn__global_z, false, true, true,
				       BH_diagnostics_info.integral_method);

//
// centroids
//
centroid_x = integral_x / integral_one;
centroid_y = integral_y / integral_one;
centroid_z = integral_z / integral_one;

//
// area, mean radius, and mass
//
area = integral_one;
mean_radius = integral_h / integral_one;
m_irreducible = sqrt(area / (16.0*PI));

//
// circumferences
//
circumference_xy
  = ps.circumference("xy", gfns::gfn__h,
		     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
					 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
							     gfns::gfn__g_dd_33,
		     BH_diagnostics_info.integral_method);
circumference_xz
  = ps.circumference("xz", gfns::gfn__h,
		     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
					 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
							     gfns::gfn__g_dd_33,
		     BH_diagnostics_info.integral_method);
circumference_yz
  = ps.circumference("yz", gfns::gfn__h,
		     gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
					 gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
							     gfns::gfn__g_dd_33,
		     BH_diagnostics_info.integral_method);
}

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

//
// This function computes the surface integral of a gridfn over the
// horizon.
//
//static
  fp BH_diagnostics::surface_integral
	(const patch_system& ps,
	 int src_gfn, bool src_gfn_is_even_across_xy_plane,
		      bool src_gfn_is_even_across_xz_plane,
		      bool src_gfn_is_even_across_yz_plane,
	 enum patch::integration_method method)
{
return ps.integrate_gridfn
	   (src_gfn, src_gfn_is_even_across_xy_plane,
		     src_gfn_is_even_across_xz_plane,
		     src_gfn_is_even_across_yz_plane,
	    gfns::gfn__h,
	    gfns::gfn__g_dd_11, gfns::gfn__g_dd_12, gfns::gfn__g_dd_13,
				gfns::gfn__g_dd_22, gfns::gfn__g_dd_23,
						    gfns::gfn__g_dd_33,
	    method);
}

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

//
// This function prints a line or two summarizing the diagnostics,
// using CCTK_VInfo().
//
void BH_diagnostics::print(int N_horizons, int hn)
	const
{
CCTK_VInfo(CCTK_THORNSTRING,
	   "AH %d/%d: r=%g at (%f,%f,%f)",
	   hn, N_horizons,
	   double(mean_radius),
	   double(centroid_x), double(centroid_y), double(centroid_z));
CCTK_VInfo(CCTK_THORNSTRING,
	   "AH %d/%d: area=%.10g m_irreducible=%.10g",
	   hn, N_horizons,
	   double(area), double(m_irreducible));
}

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

//
// This function creates a BH-diagnostics output file, writes a suitable
// header comment identifying the fields to be written by  output() ,
// flushes the stream (to help in examining the output while Cactus is
// still running), and finally returns a stdio file pointer which can be
// used by  output()  to output data to the file.
//
FILE* BH_diagnostics::setup_output_file(const struct IO_info& IO_info,
					int N_horizons, int hn)
	const
{
const int N_file_name_buffer = 200;
char file_name_buffer[N_file_name_buffer];

snprintf(file_name_buffer, N_file_name_buffer,
	 "%s.ah%d.%s",
	 IO_info.BH_diagnostics_base_file_name,
	 hn, IO_info.BH_diagnostics_file_name_extension);
FILE *fileptr = fopen(file_name_buffer, "w");
if (fileptr == NULL)
   then CCTK_VWarn(FATAL_ERROR, __LINE__, __FILE__, CCTK_THORNSTRING,
"\n"
"   BH_diagnostics::setup_output_file():\n"
"        can't open BH-diagnostics output file\n"
"        \"%s\"!"
		   ,
		   file_name_buffer);				/*NOTREACHED*/

fprintf(fileptr, "# apparent horizon %d/%d\n", hn, N_horizons);
fprintf(fileptr, "#\n");
fprintf(fileptr, "# column  1 = cctk_iteration\n");
fprintf(fileptr, "# column  2 = cctk_time\n");
fprintf(fileptr, "# column  3 = centroid_x\n");
fprintf(fileptr, "# column  4 = centroid_y\n");
fprintf(fileptr, "# column  5 = centroid_z\n");
fprintf(fileptr, "# column  6 = min radius\n");
fprintf(fileptr, "# column  7 = max radius\n");
fprintf(fileptr, "# column  8 = mean radius\n");
fprintf(fileptr, "# column  9 = xy-plane circumference\n");
fprintf(fileptr, "# column 10 = xz-plane circumference\n");
fprintf(fileptr, "# column 11 = yz-plane circumference\n");
fprintf(fileptr, "# column 12 = area\n");
fprintf(fileptr, "# column 13 = m_irreducible\n");
fflush(fileptr);

return fileptr;
}

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

//
// This function outputs a BH-diagnostics line to a stdio stream, then
// flushes the stream (to help in examining the output while Cactus is
// still running).
//
// Arguments:
// BH_diagnostics = The BH diagnostics to be written
// fileptr = The stdio file pointer to append to
//
void BH_diagnostics::output(FILE*fileptr, const struct IO_info& IO_info)
	const
{
assert(fileptr != NULL);

fprintf(fileptr,
//  cctk_iteration        mean radius     max radius                      area    m_irreducible
//  ==  cctk_time         ======  min radius      {xy,xz,yz}-plane        ======  ======
//  ==  ====  centroid_[xyz]      ======  ======  circumferences          ======  ======
//  ==  ====  ==========  ======  ======  ======  ======================  ======  ======
   "%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
	IO_info.time_iteration, double(IO_info.time),
	centroid_x, centroid_y, centroid_z,
	min_radius, max_radius, mean_radius,
	circumference_xy, circumference_xz, circumference_yz,
	area, m_irreducible);

fflush(fileptr);
}