add support for explicitly getting the induced 2-D (rho,sigma) metric

add support for computing arc lengths of patches


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@890 f88db872-0e4f-0410-b76b-b9085cfa78c5

2 files changed, 407 insertions, 56 deletions

diff --git a/src/patch/patch.cc b/src/patch/patch.cc
index 931adf6..73bd0d3 100644
--- a/src/patch/patch.cc
+++ b/src/patch/patch.cc
@@ -8,9 +8,14 @@
 // x_patch::x_patch
 // y_patch::y_patch
 //
-// patch::Jacobian_xyz_wrt_rho_sigma
+// patch::rho_sigma_metric
 //
 // patch::decode_integration_method
+// patch::rho_arc_length
+// patch::sigma_arc_length
+// z_patch::plane_arc_length
+// x_patch::plane_arc_length
+// y_patch::plane_arc_length
 // patch::integrate_gridfn
 /// patch::integration_coeff
 //
@@ -161,23 +166,29 @@ y_patch::y_patch(patch_system &my_patch_system_in, int patch_number_in,
 //******************************************************************************
 
 //
-// This function computes the Jacobian of (x,y,z) with respect to (rho,sigma),
-//	det [g_{ij} (partial x^i/partial y^U) (partial y^j/partial y^V)]
-// as per p.33 of my apparent horizon finding notes.
+// This function computes the (rho,sigma) induced 2-D metric from the
+// 3-D (x,y,z) metric of the space containing the patch, as per p.33 of
+// my apparent horizon finding notes.
 //
 // Arguments:
 // (r,rho,sigma) = The coordinates where the Jacobian is wanted.
 // partial_surface_r_wrt_(rho,sigma)
 //	= The partial derivatives of the surface radius with respect to
 //	  the (rho,sigma) coordinates.
-// g_{xx,xy,xz,yy,yz,zz} = The xyz 3-metric components.
-//
-fp patch::Jacobian_xyz_wrt_rho_sigma(fp r, fp rho, fp sigma,
-				     fp partial_surface_r_wrt_rho,
-				     fp partial_surface_r_wrt_sigma,
-				     fp g_xx, fp g_xy, fp g_xz,
-					      fp g_yy, fp g_yz,
-						       fp g_zz)
+// g_{xx,xy,xz,yy,yz,zz} = The xyz 3-metric components $g_{ij}$.
+// g_{rho_rho,rho_sigma,sigma_sigma} = The (rho,sigma) induced 2-D metric.
+//
+// Results:
+// This function returns the Jacobian of the (rho,sigma) induced 2-D metric.
+//
+fp patch::rho_sigma_metric(fp r, fp rho, fp sigma,
+			   fp partial_surface_r_wrt_rho,
+			   fp partial_surface_r_wrt_sigma,
+			   fp g_xx, fp g_xy, fp g_xz,
+				    fp g_yy, fp g_yz,
+					     fp g_zz,
+			   fp& g_rho_rho, fp& g_rho_sigma,
+					  fp& g_sigma_sigma)
 	const
 {
 fp partial_x_wrt_r, partial_x_wrt_rho, partial_x_wrt_sigma;
@@ -202,27 +213,24 @@ const fp dz_wrt_rho   = partial_z_wrt_rho
 const fp dz_wrt_sigma = partial_z_wrt_sigma
 			+ partial_z_wrt_r*partial_surface_r_wrt_sigma;
 
-const fp g_rho_rho = +     dx_wrt_rho*dx_wrt_rho*g_xx
-		     + 2.0*dx_wrt_rho*dy_wrt_rho*g_xy
-		     + 2.0*dx_wrt_rho*dz_wrt_rho*g_xz
-		     +     dy_wrt_rho*dy_wrt_rho*g_yy
-		     + 2.0*dy_wrt_rho*dz_wrt_rho*g_yz
-		     +     dz_wrt_rho*dz_wrt_rho*g_zz;
-const fp g_rho_sigma = +     dx_wrt_rho*dx_wrt_sigma    *g_xx
-		       + (   dx_wrt_rho*dy_wrt_sigma
-			   + dy_wrt_rho*dx_wrt_sigma   )*g_xy
-		       + (   dx_wrt_rho*dz_wrt_sigma
-			   + dz_wrt_rho*dx_wrt_sigma   )*g_xz
-		       +     dy_wrt_rho*dy_wrt_sigma    *g_yy
-		       + (   dy_wrt_rho*dz_wrt_sigma
-			   + dz_wrt_rho*dy_wrt_sigma   )*g_yz
-		       +     dz_wrt_rho*dz_wrt_sigma    *g_zz;
-const fp g_sigma_sigma = +     dx_wrt_sigma*dx_wrt_sigma*g_xx
-			 + 2.0*dx_wrt_sigma*dy_wrt_sigma*g_xy
-			 + 2.0*dx_wrt_sigma*dz_wrt_sigma*g_xz
-			 +     dy_wrt_sigma*dy_wrt_sigma*g_yy
-			 + 2.0*dy_wrt_sigma*dz_wrt_sigma*g_yz
-			 +     dz_wrt_sigma*dz_wrt_sigma*g_zz;
+g_rho_rho     = +     dx_wrt_rho*dx_wrt_rho*g_xx
+		+ 2.0*dx_wrt_rho*dy_wrt_rho*g_xy
+		+ 2.0*dx_wrt_rho*dz_wrt_rho*g_xz
+		+     dy_wrt_rho*dy_wrt_rho*g_yy
+		+ 2.0*dy_wrt_rho*dz_wrt_rho*g_yz
+		+     dz_wrt_rho*dz_wrt_rho*g_zz;
+g_rho_sigma =   +  dx_wrt_rho*dx_wrt_sigma                           *g_xx
+		+ (dx_wrt_rho*dy_wrt_sigma + dy_wrt_rho*dx_wrt_sigma)*g_xy
+		+ (dx_wrt_rho*dz_wrt_sigma + dz_wrt_rho*dx_wrt_sigma)*g_xz
+		+  dy_wrt_rho*dy_wrt_sigma                           *g_yy
+		+ (dy_wrt_rho*dz_wrt_sigma + dz_wrt_rho*dy_wrt_sigma)*g_yz
+		+  dz_wrt_rho*dz_wrt_sigma                           *g_zz;
+g_sigma_sigma = + dx_wrt_sigma*dx_wrt_sigma*g_xx
+		+ 2.0*dx_wrt_sigma*dy_wrt_sigma*g_xy
+		+ 2.0*dx_wrt_sigma*dz_wrt_sigma*g_xz
+		+     dy_wrt_sigma*dy_wrt_sigma*g_yy
+		+ 2.0*dy_wrt_sigma*dz_wrt_sigma*g_yz
+		+     dz_wrt_sigma*dz_wrt_sigma*g_zz;
 
 return g_rho_rho*g_sigma_sigma - jtutil::pow2(g_rho_sigma);
 }
@@ -262,6 +270,289 @@ else	error_exit(ERROR_EXIT,
 //******************************************************************************
 
 //
+// This function computes an approximation to the arc length of a surface
+// over the patch's nominal bounds along the rho direction (i.e. in a
+// dsigma=constant plane where dsigma is a multiple of 90 degrees)
+//
+// Arguments:
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
+//
+fp patch::rho_arc_length(int ghosted_radius_gfn,
+			 int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+				       int g_yy_gfn, int g_yz_gfn,
+						     int g_zz_gfn,
+			 enum integration_method method)
+	const
+{
+fp dsigma;
+if	(is_valid_dsigma(  0.0)) then dsigma =   0.0;
+else if (is_valid_dsigma( 90.0)) then dsigma =  90.0;
+else if (is_valid_dsigma(180.0)) then dsigma = 180.0;
+else if (is_valid_dsigma(-90.0)) then dsigma = -90.0;
+else	error_exit(PANIC_EXIT,
+"***** patch::rho_arc_length(): can't find valid dsigma\n"
+"                               which is a multiple of 90 degrees!\n"
+"                               %s patch: [min,max]_dsigma()=[%g,%g]\n"
+		   ,
+		   name(), min_dsigma(), max_dsigma());
+const fp sigma = sigma_of_dsigma(dsigma);
+const int isigma = isigma_of_sigma(sigma);
+
+fp sum = 0.0;
+
+	for (int irho = min_irho() ; irho <= max_irho() ; ++irho)
+	{
+	const fp rho = rho_of_irho(irho);
+	const fp r = ghosted_gridfn(ghosted_radius_gfn, irho,isigma);
+	const fp partial_surface_r_wrt_rho
+		= partial_rho  (ghosted_radius_gfn, irho,isigma);
+	const fp partial_surface_r_wrt_sigma
+		= partial_sigma(ghosted_radius_gfn, irho,isigma);
+
+	const fp g_xx = gridfn(g_xx_gfn, irho,isigma);
+	const fp g_xy = gridfn(g_xy_gfn, irho,isigma);
+	const fp g_xz = gridfn(g_xz_gfn, irho,isigma);
+	const fp g_yy = gridfn(g_yy_gfn, irho,isigma);
+	const fp g_yz = gridfn(g_yz_gfn, irho,isigma);
+	const fp g_zz = gridfn(g_zz_gfn, irho,isigma);
+
+	fp g_rho_rho, g_rho_sigma, g_sigma_sigma;
+	rho_sigma_metric(r, rho, sigma,
+			 partial_surface_r_wrt_rho,
+			 partial_surface_r_wrt_sigma,
+			 g_xx, g_xy, g_xz,
+			       g_yy, g_yz,
+				     g_zz,
+			 g_rho_rho, g_rho_sigma,
+				    g_sigma_sigma);
+
+	const fp coeff = integration_coeff(method,
+					   max_irho()-min_irho(),
+					   irho      -min_irho());
+
+	sum += coeff * sqrt(g_rho_rho);
+	}
+
+return delta_rho() * sum;
+}
+
+//******************************************************************************
+
+//
+// This function computes an approximation to the arc length of a surface
+// over the patch's nominal bounds along the sigma direction (i.e. in a
+// drho=constant plane where drho is a multiple of 90 degrees)
+//
+// Arguments:
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
+//
+fp patch::sigma_arc_length(int ghosted_radius_gfn,
+			   int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					 int g_yy_gfn, int g_yz_gfn,
+						       int g_zz_gfn,
+			   enum integration_method method)
+	const
+{
+fp drho;
+if	(is_valid_drho(  0.0)) then drho =   0.0;
+else if (is_valid_drho( 90.0)) then drho =  90.0;
+else if (is_valid_drho(180.0)) then drho = 180.0;
+else if (is_valid_drho(-90.0)) then drho = -90.0;
+else	error_exit(PANIC_EXIT,
+"***** patch::sigma_arc_length(): can't find valid drho\n"
+"                                 which is a multiple of 90 degrees!\n"
+"                                 %s patch: [min,max]_drho()=[%g,%g]\n"
+		   ,
+		   name(), min_drho(), max_drho());
+const fp rho = rho_of_drho(drho);
+const int irho = irho_of_rho(rho);
+
+fp sum = 0.0;
+
+	for (int isigma = min_isigma() ; isigma <= max_isigma() ; ++isigma)
+	{
+	const fp sigma = sigma_of_isigma(isigma);
+	const fp r = ghosted_gridfn(ghosted_radius_gfn, irho,isigma);
+	const fp partial_surface_r_wrt_rho
+		= partial_rho  (ghosted_radius_gfn, irho,isigma);
+	const fp partial_surface_r_wrt_sigma
+		= partial_sigma(ghosted_radius_gfn, irho,isigma);
+
+	const fp g_xx = gridfn(g_xx_gfn, irho,isigma);
+	const fp g_xy = gridfn(g_xy_gfn, irho,isigma);
+	const fp g_xz = gridfn(g_xz_gfn, irho,isigma);
+	const fp g_yy = gridfn(g_yy_gfn, irho,isigma);
+	const fp g_yz = gridfn(g_yz_gfn, irho,isigma);
+	const fp g_zz = gridfn(g_zz_gfn, irho,isigma);
+
+	fp g_rho_rho, g_rho_sigma, g_sigma_sigma;
+	rho_sigma_metric(r, rho, sigma,
+			 partial_surface_r_wrt_rho,
+			 partial_surface_r_wrt_sigma,
+			 g_xx, g_xy, g_xz,
+			       g_yy, g_yz,
+				     g_zz,
+			 g_rho_rho, g_rho_sigma,
+				    g_sigma_sigma);
+
+	const fp coeff = integration_coeff(method,
+					   max_isigma()-min_isigma(),
+					   isigma      -min_isigma());
+
+	sum += coeff * sqrt(g_sigma_sigma);
+	}
+
+return delta_sigma() * sum;
+}
+
+//******************************************************************************
+
+//
+// This function computes the arc length of a surface over a z patch's
+// nominal bounds along the coord=0 plane (coord = 'x' or 'y')
+//
+// Arguments:
+// coord = (in) The coordinate plane in which to compute the arc length.
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
+//
+fp z_patch::plane_arc_length(char coord,
+			     int ghosted_radius_gfn,
+			     int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					   int g_yy_gfn, int g_yz_gfn,
+							 int g_zz_gfn,
+			     enum integration_method method)
+	const
+{
+switch	(coord)
+	{
+case 'x':
+	return rho_arc_length(ghosted_radius_gfn,
+			      g_xx_gfn, g_xy_gfn, g_xz_gfn,
+					g_yy_gfn, g_yz_gfn,
+						  g_zz_gfn,
+			      method);
+case 'y':
+	return sigma_arc_length(ghosted_radius_gfn,
+				g_xx_gfn, g_xy_gfn, g_xz_gfn,
+					  g_yy_gfn, g_yz_gfn,
+						    g_zz_gfn,
+				method);
+case 'z':
+	error_exit(ERROR_EXIT,
+"***** z_patch::plane_arc_length(): %s patch, coord='z', but\n"
+"                                   +/-z patch doesn't contain z=0 plane!\n"
+,
+		   name());					/*NOTREACHED*/
+default:
+	error_exit(PANIC_EXIT,
+"***** z_patch::plane_arc_length(): bad coord=(int)%d!\n",
+		   int(coord));					/*NOTREACHED*/
+	}
+}
+
+//******************************************************************************
+
+//
+// This function computes the arc length of a surface over an x patch's
+// nominal bounds along the coord=0 plane (coord = 'y' or 'z')
+//
+// Arguments:
+// coord = (in) The coordinate plane in which to compute the arc length.
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
+//
+fp x_patch::plane_arc_length(char coord,
+			     int ghosted_radius_gfn,
+			     int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					   int g_yy_gfn, int g_yz_gfn,
+							 int g_zz_gfn,
+			     enum integration_method method)
+	const
+{
+switch	(coord)
+	{
+case 'y':
+	return rho_arc_length(ghosted_radius_gfn,
+			      g_xx_gfn, g_xy_gfn, g_xz_gfn,
+					g_yy_gfn, g_yz_gfn,
+						  g_zz_gfn,
+			      method);
+case 'z':
+	return sigma_arc_length(ghosted_radius_gfn,
+				g_xx_gfn, g_xy_gfn, g_xz_gfn,
+					  g_yy_gfn, g_yz_gfn,
+						    g_zz_gfn,
+				method);
+case 'x':
+	error_exit(ERROR_EXIT,
+"***** x_patch::plane_arc_length(): %s patch, coord='x', but\n"
+"                                   +/-x patch doesn't contain x=0 plane!\n"
+,
+		   name());					/*NOTREACHED*/
+default:
+	error_exit(PANIC_EXIT,
+"***** x_patch::plane_arc_length(): bad coord=(int)%d!\n",
+		   int(coord));					/*NOTREACHED*/
+	}
+}
+
+//******************************************************************************
+
+//
+// This function computes the arc length of a surface over an y patch's
+// nominal bounds along the coord=0 plane (coord = 'x' or 'z')
+//
+// Arguments:
+// coord = (in) The coordinate plane in which to compute the arc length.
+// ghosted_radius_gfn = (in) The surface radius.
+// g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
+// method = (in) Selects the integration scheme.
+//
+fp y_patch::plane_arc_length(char coord,
+			     int ghosted_radius_gfn,
+			     int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					   int g_yy_gfn, int g_yz_gfn,
+							 int g_zz_gfn,
+			     enum integration_method method)
+	const
+{
+switch	(coord)
+	{
+case 'x':
+	return rho_arc_length(ghosted_radius_gfn,
+			      g_xx_gfn, g_xy_gfn, g_xz_gfn,
+					g_yy_gfn, g_yz_gfn,
+						  g_zz_gfn,
+			      method);
+case 'z':
+	return sigma_arc_length(ghosted_radius_gfn,
+				g_xx_gfn, g_xy_gfn, g_xz_gfn,
+					  g_yy_gfn, g_yz_gfn,
+						    g_zz_gfn,
+				method);
+case 'y':
+	error_exit(ERROR_EXIT,
+"***** y_patch::plane_arc_length(): %s patch, coord='y', but\n"
+"                                   +/-y patch doesn't contain y=0 plane!\n"
+,
+		   name());					/*NOTREACHED*/
+default:
+	error_exit(PANIC_EXIT,
+"***** y_patch::plane_arc_length(): bad coord=(int)%d!\n",
+		   int(coord));					/*NOTREACHED*/
+	}
+}
+
+//******************************************************************************
+
+//
 // This function computes an approximation to the (surface) integral of
 // a gridfn over the patch's nominal area,
 //	$\int f(\rho,\sigma) \, dA$
@@ -277,10 +568,6 @@ else	error_exit(ERROR_EXIT,
 // g_{xx,xy,xz,yy,yz,zz}_gfn = (in) The xyz 3-metric components.
 // method = (in) Selects the integration scheme.
 //
-// Bugs:
-// - The integration coefficients are computed in a rather inefficient
-//   manner.
-//
 fp patch::integrate_gridfn(int unknown_src_gfn,
 			   int ghosted_radius_gfn,
 			   int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
@@ -314,12 +601,15 @@ fp sum = 0.0;
 	const fp g_yz = gridfn(g_yz_gfn, irho,isigma);
 	const fp g_zz = gridfn(g_zz_gfn, irho,isigma);
 
-	const fp Jac = Jacobian_xyz_wrt_rho_sigma(r, rho, sigma,
-						  partial_surface_r_wrt_rho,
-						  partial_surface_r_wrt_sigma,
-						  g_xx, g_xy, g_xz,
-							g_yy, g_yz,
-							      g_zz);
+	fp g_rho_rho, g_rho_sigma, g_sigma_sigma;
+	const fp Jac = rho_sigma_metric(r, rho, sigma,
+					partial_surface_r_wrt_rho,
+					partial_surface_r_wrt_sigma,
+					g_xx, g_xy, g_xz,
+					      g_yy, g_yz,
+						    g_zz,
+					g_rho_rho, g_rho_sigma,
+						   g_sigma_sigma);
 
 	const fp coeff_rho   = integration_coeff(method,
 						 max_irho()-min_irho(),
@@ -328,9 +618,10 @@ fp sum = 0.0;
 						 max_isigma()-min_isigma(),
 						 isigma      -min_isigma());
 
-	sum += fn * sqrt(jtutil::abs(Jac)) * coeff_rho*coeff_sigma;
+	sum += coeff_rho*coeff_sigma * fn * sqrt(jtutil::abs(Jac));
 	}
 	}
+
 return delta_rho() * delta_sigma() * sum;
 }
 
diff --git a/src/patch/patch.hh b/src/patch/patch.hh
index 25a71a8..daaf608 100644
--- a/src/patch/patch.hh
+++ b/src/patch/patch.hh
@@ -150,6 +150,7 @@ public:
 
 	// are we a +[xyz] or -[xyz] patch?
 	bool is_plus() const { return is_plus_; }
+
 	// ... values for the  is_plus_in  constructor argument
 	//     FIXME: these should really be bool, but then we couldn't
 	//            use the "enum hack" for in-class constants
@@ -219,17 +220,6 @@ public:
 					 fp& xcos, fp& ycos, fp& zcos)
 		const = 0;
 
-	// Jacobian of (x,y,z) with respect to (rho,sigma),
-	// as per p.33 of my apparent horizon finding notes
-	// det [g_{ij} (partial x^i/partial y^u) (partial y^j/partial y^v)]
-	fp Jacobian_xyz_wrt_rho_sigma(fp r, fp rho, fp sigma,
-				      fp partial_surface_r_wrt_rho,
-				      fp partial_surface_r_wrt_sigma,
-				      fp g_xx, fp g_xy, fp g_xz,
-					       fp g_yy, fp g_yz,
-							fp g_zz)
-		const;
-
 	// partial (x,y,z) / partial (rho,sigma)
 	virtual void partial_xyz_wrt_r_rho_sigma
 	   (fp r, fp rho, fp sigma,
@@ -260,6 +250,19 @@ public:
 	virtual fp partial2_sigma_wrt_yz(fp x, fp y, fp z) const = 0;
 	virtual fp partial2_sigma_wrt_zz(fp x, fp y, fp z) const = 0;
 
+	// compute (rho,sigma) 2-D induced metric from 3-D xyz metric
+	// as per p.33 of my apparent horizon finding notes
+	// ... returns Jacobian of (rho,sigma) 2-D induced metric
+	fp rho_sigma_metric(fp r, fp rho, fp sigma,
+			    fp partial_surface_r_wrt_rho,
+			    fp partial_surface_r_wrt_sigma,
+			    fp g_xx, fp g_xy, fp g_xz,
+				     fp g_yy, fp g_yz,
+					      fp g_zz,
+			    fp& g_rho_rho, fp& g_rho_sigma,
+					   fp& g_sigma_sigma)
+		const;
+
 	// plotting coordinates (dpx,dpy)
 	// ... character string describing how (dpx,dpy) are
 	//     defined in terms of (mu,nu,phi), eg "90 - drho = 90 - dphi"
@@ -272,7 +275,7 @@ public:
 
 
 	//
-	// ***** gridfn operations *****
+	// ***** line/surface integrals *****
 	//
 public:
 
@@ -363,6 +366,34 @@ public:
 	  enum integration_method
 	    decode_integration_method(const char method_string[]);
 
+	// compute the arc length of a surface over the patch's nominal bounds
+	// ... along the coord=0 plane (coord = 'x', 'y', or 'z')
+	//     (error_exit() if the patch doesn't contain that coordinate plane)
+	virtual fp plane_arc_length(char coord,
+				    int ghosted_radius_gfn,
+				    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+						  int g_yy_gfn, int g_yz_gfn,
+								int g_zz_gfn,
+				    enum integration_method method)
+		const = 0;
+
+	// ... along the rho direction (i.e. in a dsigma=constant plane
+	//     where dsigma is a multiple of 90 degrees)
+	fp rho_arc_length(int ghosted_radius_gfn,
+			  int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					int g_yy_gfn, int g_yz_gfn,
+						      int g_zz_gfn,
+			  enum integration_method method)
+		const;
+	// ... along the sigma direction (i.e. in a drho=constant plane
+	//     where drho is a multiple of 90 degrees)
+	fp sigma_arc_length(int ghosted_radius_gfn,
+			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					  int g_yy_gfn, int g_yz_gfn,
+							int g_zz_gfn,
+			    enum integration_method method)
+		const;
+
 	// compute the surface integral of a gridfn over the patch's
 	// nominal area,
 	//	$\int f(\rho,\sigma) \, dA$
@@ -704,6 +735,16 @@ public:
 		return is_plus() ? drho : 180.0 - drho;
 		}
 
+	// compute the arc length of a surface over the patch's nominal bounds
+	// along the coord=0 plane (coord = 'x' or 'y')
+	fp plane_arc_length(char coord,
+			    int ghosted_radius_gfn,
+			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					  int g_yy_gfn, int g_yz_gfn,
+							int g_zz_gfn,
+			    enum integration_method method)
+		const;
+
 	// constructor, destructor
 	z_patch(patch_system &my_patch_system_in, int patch_number_in,
 		const char* name_in, bool is_plus_in,
@@ -840,6 +881,16 @@ public:
 		return is_plus() ? dsigma : 180.0 - dsigma;
 		}
 
+	// compute the arc length of a surface over the patch's nominal bounds
+	// along the coord=0 plane (coord = 'y' or 'z')
+	fp plane_arc_length(char coord,
+			    int ghosted_radius_gfn,
+			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					  int g_yy_gfn, int g_yz_gfn,
+							int g_zz_gfn,
+			    enum integration_method method)
+		const;
+
 	// constructor, destructor
 	x_patch(patch_system &my_patch_system_in, int patch_number_in,
 		const char* name_in, bool is_plus_in,
@@ -976,6 +1027,15 @@ public:
 	fp dpy_of_rho_sigma(fp rho, fp sigma) const
 		{ return jtutil::degrees_of_radians(rho); }
 
+	// compute the arc length of a surface over the patch's nominal bounds
+	// along the coord=0 plane (coord = 'x' or 'z')
+	fp plane_arc_length(char coord,
+			    int ghosted_radius_gfn,
+			    int g_xx_gfn, int g_xy_gfn, int g_xz_gfn,
+					  int g_yy_gfn, int g_yz_gfn,
+							int g_zz_gfn,
+			    enum integration_method method)
+		const;
 
 	// constructor, destructor
 	y_patch(patch_system &my_patch_system_in, int patch_number_in,