path: root/src/gr/driver.cc

// driver.cc -- top level driver for finding apparent horizons
// $Id$
//
// <<<prototypes for functions local to this file>>>
// AHFinderDirect_driver - top-level driver
/// decode_Jacobian_method - decode the Jacobian_method parameter
/// decode_geometry_method - decode the geometry_method parameter
///
/// setup_Kerr_horizon - set up Kerr horizon in h (Kerr or Kerr-Schild coords)
/// setup_ellipsoid - setup up a coordiante ellipsoid in h
///
/// print_Jacobians - print a pair of Jacobians
///

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

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

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

namespace {
enum Jacobian_method
  decode_Jacobian_method(const char Jacobian_method_string[]);
enum geometry_method
  decode_geometry_method(const char geometry_method_string[]);

void setup_Kerr_horizon(patch_system& ps,
			fp x_posn, fp y_posn, fp z_posn,
			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);

void print_Jacobians(const patch_system& ps,
		     const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr,
		     const char file_name[]);
	  };

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

//
// 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 parameters and the geometry interpolator
//
struct geometry_info gi;
gi.geometry_method = decode_geometry_method(geometry_method);

// parameters for geometry_method = "interpolate from Cactus grid"
CCTK_VInfo(CCTK_THORNSTRING, "   setting up geometry interpolator");
gi.operator_handle = CCTK_InterpHandle(geometry_interpolator_name);
if (gi.operator_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "couldn't find interpolator \"%s\"!",
		   geometry_interpolator_name);		/*NOTREACHED*/
gi.param_table_handle = Util_TableCreateFromString(geometry_interpolator_pars);
if (gi.param_table_handle < 0)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "bad geometry-interpolator parameter(s) \"%s\"!",
		   geometry_interpolator_pars);		/*NOTREACHED*/

// parameters for geometry_method = "Schwarzschild/EF"
gi.geometry__Schwarzschild_EF__mass     = geometry__Schwarzschild_EF__mass;
gi.geometry__Schwarzschild_EF__x_posn   = geometry__Schwarzschild_EF__x_posn;
gi.geometry__Schwarzschild_EF__y_posn   = geometry__Schwarzschild_EF__y_posn;
gi.geometry__Schwarzschild_EF__z_posn   = geometry__Schwarzschild_EF__z_posn;
gi.geometry__Schwarzschild_EF__epsilon  = geometry__Schwarzschild_EF__epsilon;
gi.geometry__Schwarzschild_EF__Delta_xyz= geometry__Schwarzschild_EF__Delta_xyz;


//
// 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,
		gfns::nominal_min_gfn, gfns::nominal_max_gfn,
		gfns::ghosted_min_gfn, gfns::ghosted_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 h from \"%s\"",
		   h_file_name);
	ps.read_ghosted_gridfn(gfns::gfn__h,
			       h_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__x_posn,
				   initial_guess__Kerr__y_posn,
				   initial_guess__Kerr__z_posn,
				   initial_guess__Kerr__mass,
				   initial_guess__Kerr__spin,
				   false);	// use Kerr coords
	else if (STRING_EQUAL(initial_guess_method, "Kerr/Kerr-Schild"))
	   then setup_Kerr_horizon(ps,
				   initial_guess__Kerr__x_posn,
				   initial_guess__Kerr__y_posn,
				   initial_guess__Kerr__z_posn,
				   initial_guess__Kerr__mass,
				   initial_guess__Kerr__spin,
				   true);	// use Kerr-Schild coords
	else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
			   "unknown initial_guess_method=\"%s\"!",
			   initial_guess_method);		/*NOTREACHED*/

	// write initial guess back to this file
	CCTK_VInfo(CCTK_THORNSTRING,
		   "writing initial guess h to \"%s\"",
		   h_file_name);
	ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
					 true, gfns::gfn__h,
					 h_file_name,
					 false);	// no ghost zones
	}


//
// decode/copy parameters into structures
//

struct Jacobian_info Jacobian_info;
Jacobian_info.Jacobian_method = decode_Jacobian_method(Jacobian_method);
Jacobian_info.perturbation_amplitude = Jacobian_perturbation_amplitude;

struct solver_info solver_info;
solver_info.max_Newton_iterations = max_Newton_iterations;
solver_info.output_h_and_H_at_each_Newton_iteration
	= (output_h_and_H_at_each_Newton_iteration != 0);
solver_info.h_file_name = h_file_name;
solver_info.H_of_h_file_name = H_of_h_file_name;
solver_info.H_norm_for_convergence = H_norm_for_convergence;
solver_info.Delta_h_norm_for_convergence = Delta_h_norm_for_convergence;
solver_info.final_H_update_if_exit_x_H_small
	= (final_H_update_if_exit_x_H_small != 0);

//
// find the apparent horizon
//
if      (STRING_EQUAL(method, "horizon function"))
   then {
	jtutil::norm<fp> H_norms;

	const int timer_handle = CCTK_TimerCreateI();
	CCTK_TimerStartI(timer_handle);
	  horizon_function(ps, cgi, gi, false, true, &H_norms);
	CCTK_TimerStopI(timer_handle);
	printf("timer stats for evaluating H(h):\n");
	CCTK_TimerPrintDataI(timer_handle, -1);

	CCTK_VInfo(CCTK_THORNSTRING,
		   "   H(h) rms-norm %.2e, infinity-norm %.2e",
		   H_norms.rms_norm(), H_norms.infinity_norm());
	CCTK_VInfo(CCTK_THORNSTRING,
		   "writing H(h) to \"%s\"",
		   H_of_h_file_name);
	ps.print_gridfn_with_xyz(gfns::gfn__H,
				 true, gfns::gfn__h,
				 H_of_h_file_name);
	}
else if (STRING_EQUAL(method, "Jacobian test"))
   then {
	// symbolic differentiation with finite diff d/dr
	Jacobian& Jac_SD_FDdr = new_Jacobian(ps, Jacobian_storage_method);
	horizon_function(ps, cgi, gi);
	Jacobian_info.Jacobian_method = symbolic_differentiation_with_FD_dr;
	horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_SD_FDdr);

	// numerical perturbation
	Jacobian& Jac_NP = new_Jacobian(ps, Jacobian_storage_method);
	horizon_function(ps, cgi, gi);
	Jacobian_info.Jacobian_method = numerical_perturbation;
	horizon_Jacobian(ps, cgi, gi, Jacobian_info, Jac_NP);

	CCTK_VInfo(CCTK_THORNSTRING,
		   "writing Jacobians to \"%s\"",
		   Jacobian_file_name);
	print_Jacobians(ps, Jac_NP, Jac_SD_FDdr, Jacobian_file_name);
	}
else if (STRING_EQUAL(method, "Newton solve"))
   then {
	Jacobian& Jac = new_Jacobian(ps, Jacobian_storage_method);

	const int timer_handle = CCTK_TimerCreateI();
	CCTK_TimerStartI(timer_handle);
	  Newton_solve(ps, cgi, gi, Jacobian_info, solver_info, Jac);
	CCTK_TimerStopI(timer_handle);
	printf("timer stats for finding the apparent horizon:\n");
	CCTK_TimerPrintDataI(timer_handle, -1);

	CCTK_VInfo(CCTK_THORNSTRING,
		   "writing final horizon shape h to \"%s\"",
		   h_file_name);
	ps.print_ghosted_gridfn_with_xyz(gfns::gfn__h,
					 true, gfns::gfn__h,
					 h_file_name,
					 false);	// no ghost zones
	CCTK_VInfo(CCTK_THORNSTRING,
		   "writing H(h) to \"%s\"",
		   H_of_h_file_name);
	ps.print_gridfn_with_xyz(gfns::gfn__H,
				 true, gfns::gfn__h,
				 H_of_h_file_name);
	}
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown method=\"%s\"!",
		   method);					/*NOTREACHED*/
}

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

//
// This function decodes the  Jacobian_method  parameter (string) into
// an internal enum for future use.
//
namespace {
enum Jacobian_method
  decode_Jacobian_method(const char Jacobian_method_string[])
{
if	(STRING_EQUAL(Jacobian_method_string,
		      "numerical perturbation"))
   then return numerical_perturbation;
else if (STRING_EQUAL(Jacobian_method_string,
		      "symbolic differentiation with finite diff d/dr"))
   then return symbolic_differentiation_with_FD_dr;
else if (STRING_EQUAL(Jacobian_method_string,
		      "symbolic differentiation"))
   then return symbolic_differentiation;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown Jacobian_method_string=\"%s\"!",
		   Jacobian_method_string);			/*NOTREACHED*/
}
	  }

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

//
// This function decodes the  geometry_method  parameter (string) into
// an internal enum for future use.
//
namespace {
enum geometry_method
  decode_geometry_method(const char geometry_method_string[])
{
if	(STRING_EQUAL(geometry_method_string, "interpolate from Cactus grid"))
   then return geometry__interpolate_from_Cactus_grid;
else if (STRING_EQUAL(geometry_method_string, "Schwarzschild/EF"))
   then return geometry__Schwarzschild_EF;
else	CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "unknown geometry_method_string=\"%s\"!",
		   geometry_method_string);			/*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]_posn = 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_posn, fp y_posn, fp z_posn,
			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,
	   "   posn=(%g,%g,%g) mass=%g spin=J/m^2=%g",
	   double(x_posn), double(y_posn), double(z_posn),
	   double(m), double(a));

// 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_posn, y_posn, z_posn,
		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(gfns::gfn__h, irho,isigma) = r;
		}
		}
	}
}
	  }

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

//
// This function prints a pair of Jacobian matrices (and their difference)
// to a named output file.
//
namespace {
void print_Jacobians(const patch_system& ps,
		     const Jacobian& Jac_NP, const Jacobian& Jac_SD_FDdr,
		     const char file_name[])
{
FILE *fileptr = fopen(file_name, "w");
if (fileptr == NULL)
   then CCTK_VWarn(-1, __LINE__, __FILE__, CCTK_THORNSTRING,
		   "print_Jacobians(): can't open output file \"%s\"!",
		   file_name);					/*NOTREACHED*/

fprintf(fileptr, "# column 1 = x II\n");
fprintf(fileptr, "# column 2 = x patch number\n");
fprintf(fileptr, "# column 3 = x irho\n");
fprintf(fileptr, "# column 4 = x isigma\n");
fprintf(fileptr, "# column 5 = y JJ\n");
fprintf(fileptr, "# column 6 = y patch number\n");
fprintf(fileptr, "# column 7 = y irho\n");
fprintf(fileptr, "# column 8 = y isigma\n");
fprintf(fileptr, "# column 9 = Jac_NP(II,JJ)\n");
fprintf(fileptr, "# column 10 = Jac_SD_FDdr(II,JJ)\n");
fprintf(fileptr, "# column 11 = abs error in Jac_SD_FDdr\n");
fprintf(fileptr, "# column 12 = rel error in Jac_SD_FDdr\n");

	for (int xpn = 0 ; xpn < ps.N_patches() ; ++xpn)
	{
	patch& xp = ps.ith_patch(xpn);

	for (int x_irho = xp.min_irho() ; x_irho <= xp.max_irho() ; ++x_irho)
	{
	for (int x_isigma = xp.min_isigma() ;
	     x_isigma <= xp.max_isigma() ;
	     ++x_isigma)
	{
	const int II = ps.gpn_of_patch_irho_isigma(xp, x_irho,x_isigma);

		for (int ypn = 0 ; ypn < ps.N_patches() ; ++ypn)
		{
		patch& yp = ps.ith_patch(ypn);

		for (int y_irho = yp.min_irho() ;
		     y_irho <= yp.max_irho() ;
		     ++y_irho)
		{
		for (int y_isigma = yp.min_isigma() ;
		     y_isigma <= yp.max_isigma() ;
		     ++y_isigma)
		{
		const int JJ = ps.gpn_of_patch_irho_isigma(yp, y_irho,y_isigma);

		if (! Jac_SD_FDdr.is_explicitly_stored(II,JJ))
		   then continue;			// skip sparse points

		const fp NP      = Jac_NP     (II,JJ);
		const fp SD_FDdr = Jac_SD_FDdr(II,JJ);

		if ((NP == 0.0) && (SD_FDdr == 0.0))
		   then continue;		// skip zero values (if == )

		const fp abs_NP      = jtutil::abs(NP     );
		const fp abs_SD_FDdr = jtutil::abs(SD_FDdr);
		const fp max_abs = jtutil::max(abs_SD_FDdr, abs_NP);
		const fp SD_FDdr_abs_error = SD_FDdr - NP;
		const fp SD_FDdr_rel_error = SD_FDdr_abs_error / max_abs;

		fprintf(fileptr,
			"%d %d %d %d\t%d %d %d %d\t%.10g\t%.10g\t%e\t%e\n",
			II, xpn, x_irho, x_isigma,
			JJ, ypn, y_irho, y_isigma,
			double(NP), double(SD_FDdr),
			double(SD_FDdr_abs_error), double(SD_FDdr_rel_error));
		}
		}
		}
	}
	}
	}

fclose(fileptr);
}
	  }