// BH_diagnostics.cc -- compute/print BH diagnostics
// $Header$
//
// BH_diagnostics::BH_diagnostics - initialize a  struct BH_diagnostics
//
// BH_diagnostics::copy_to_buffer - copy diagnostics to buffer
// BH_diagnostics::copy_from_buffer - copy buffer to 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),
	  min_x(0.0), max_x(0.0),
	  min_y(0.0), max_y(0.0),
	  min_z(0.0), max_z(0.0),
	  circumference_xy(0.0), circumference_xz(0.0), circumference_yz(0.0),
	  area(0.0),
	  m_irreducible(0.0)
{ }

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

//
// This function copies the diagnostics to a user-supplied buffer.
//
void BH_diagnostics::copy_to_buffer(CCTK_REAL buffer[N_buffer])
	const
{
buffer[posn__centroid_x] = centroid_x;
buffer[posn__centroid_y] = centroid_y;
buffer[posn__centroid_z] = centroid_z;

buffer[posn__min_radius]  = min_radius;
buffer[posn__max_radius]  = max_radius;
buffer[posn__mean_radius] = mean_radius;

buffer[posn__min_x] = min_x;
buffer[posn__max_x] = max_x;
buffer[posn__min_y] = min_y;
buffer[posn__max_y] = max_y;
buffer[posn__min_z] = min_z;
buffer[posn__max_z] = max_z;

buffer[posn__circumference_xy] = circumference_xy;
buffer[posn__circumference_xz] = circumference_xz;
buffer[posn__circumference_yz] = circumference_yz;
buffer[posn__area]             = area;
buffer[posn__m_irreducible]    = m_irreducible;
}

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

//
// This function copies a user-supplied buffer to the diagnostics.
//
void BH_diagnostics::copy_from_buffer(const CCTK_REAL buffer[N_buffer])
{
centroid_x = buffer[posn__centroid_x];
centroid_y = buffer[posn__centroid_y];
centroid_z = buffer[posn__centroid_z];

 min_radius = buffer[posn__min_radius];
 max_radius = buffer[posn__max_radius];
mean_radius = buffer[posn__mean_radius];

circumference_xy = buffer[posn__circumference_xy];
circumference_xz = buffer[posn__circumference_xz];
circumference_yz = buffer[posn__circumference_yz];
            area = buffer[posn__area];
   m_irreducible = buffer[posn__m_irreducible];
}

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

//
// Given that an apparent horizon has been found, this function computes
// various black hole diagnostics.
//
// Inputs (gridfns)
// h		# ghosted
// one		# nominal
// global_[xyz]	# nominal
//
// Bugs:
// The computation is rather inefficient -- we make many passes over the
// angular grid, instead of doing everything in one pass.
//
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();


//
// xyz bounding box of horizon
//

jtutil::norm<fp> x_norms;
ps.gridfn_norms(gfns::gfn__global_x, x_norms);
min_x = x_norms.min_value();
max_x = x_norms.max_value();

jtutil::norm<fp> y_norms;
ps.gridfn_norms(gfns::gfn__global_y, y_norms);
min_y = y_norms.min_value();
max_y = y_norms.max_value();

jtutil::norm<fp> z_norms;
ps.gridfn_norms(gfns::gfn__global_z, z_norms);
min_z = z_norms.min_value();
max_z = z_norms.max_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 = min x\n");
fprintf(fileptr, "# column 10 = max x\n");
fprintf(fileptr, "# column 11 = min y\n");
fprintf(fileptr, "# column 12 = max y\n");
fprintf(fileptr, "# column 13 = min z\n");
fprintf(fileptr, "# column 14 = max z\n");
fprintf(fileptr, "# column 15 = xy-plane circumference\n");
fprintf(fileptr, "# column 16 = xz-plane circumference\n");
fprintf(fileptr, "# column 17 = yz-plane circumference\n");
fprintf(fileptr, "# column 18 = area\n");
fprintf(fileptr, "# column 19 = 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        min radius      mean radius
     //  ==  cctk_time         ======  max radius
     //  ==  ====  centroid_[xyz]      ======  ======
     //  ==  ====  ==========  ======  ======  ======
	"%d\t%.3f\t%f\t%f\t%f\t%#.10g\t%#.10g\t%#.10g\t",
	IO_info.time_iteration, double(IO_info.time),
	double(centroid_x), double(centroid_y), double(centroid_z),
	double(min_radius), double(max_radius), double(mean_radius));

fprintf(fileptr,
     //  {min,max}_x     {min,max}_y     {min,max}_z
     //  ==============  ==============  ==============
	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\t",
	double(min_x), double(max_x),
	double(min_y), double(max_y),
	double(min_z), double(max_z));

fprintf(fileptr,
     //  {xy,xz,yz}-plane        area    m_irreducible
     //  circumferences          ======  ======
     //  ======================  ======  ======
	"%#.10g\t%#.10g\t%#.10g\t%#.10g\t%#.10g\n",
	double(circumference_xy),
	double(circumference_xz),
	double(circumference_yz),
	double(area), double(m_irreducible));

fflush(fileptr);
}