add 2nd derivatives partial2_{mu,nu,phi}_wrt_{xx,xy,xz,yy,yz,zz}

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

2 files changed, 62 insertions, 11 deletions

diff --git a/src/patch/coords.cc b/src/patch/coords.cc
index b33a6af..0def972 100644
--- a/src/patch/coords.cc
+++ b/src/patch/coords.cc
@@ -8,6 +8,9 @@
 // (r,(mu,nu,phi)) <--> (x,y,z)
 // ((mu,nu,phi)) --> the 3rd
 //
+// partial_{mu,nu,phi}_wrt_{x,y,z}
+// partial2_{mu,nu,phi}_wrt_{xx,xy,xz,yy,yz,zz}
+//
 // usual polar spherical (r,theta,phi) <--> (x,y,z)
 // ((mu,nu,phi)) <--> usual polar spherical (theta,phi)
 //
@@ -249,23 +252,21 @@ fp phi_of_xy(fp x, fp y) { return arctan_xy(x,y); }
 	  }
 
 //******************************************************************************
+//******************************************************************************
+//******************************************************************************
 
 //
 // these functions compute the partial derivatives
 //	partial {mu,nu,phi} / partial {x,y,z}
-// as per the formulas on p.5 of my apparent horizon finding notes
+// as computed by the maple file "coord_derivs.{maple,out}" in this directory
 //
-fp partial_mu_wrt_x(fp x, fp y, fp z) { return 0.0; }
-fp partial_mu_wrt_y(fp x, fp y, fp z) { return  z / (y*y + z*z); }
-fp partial_mu_wrt_z(fp x, fp y, fp z) { return -y / (y*y + z*z); }
+fp  partial_mu_wrt_x(fp x, fp y, fp z) { return 0.0; }
+fp  partial_mu_wrt_y(fp x, fp y, fp z) { return  z / (y*y + z*z); }
+fp  partial_mu_wrt_z(fp x, fp y, fp z) { return -y / (y*y + z*z); }
 
-//**************************************
-
-fp partial_nu_wrt_x(fp x, fp y, fp z) { return  z / (x*x + z*z); }
-fp partial_nu_wrt_y(fp x, fp y, fp z) { return 0.0; }
-fp partial_nu_wrt_z(fp x, fp y, fp z) { return -x / (x*x + z*z); }
-
-//**************************************
+fp  partial_nu_wrt_x(fp x, fp y, fp z) { return  z / (x*x + z*z); }
+fp  partial_nu_wrt_y(fp x, fp y, fp z) { return 0.0; }
+fp  partial_nu_wrt_z(fp x, fp y, fp z) { return -x / (x*x + z*z); }
 
 fp partial_phi_wrt_x(fp x, fp y, fp z) { return -y / (x*x + y*y); }
 fp partial_phi_wrt_y(fp x, fp y, fp z) { return  x / (x*x + y*y); }
@@ -274,6 +275,36 @@ fp partial_phi_wrt_z(fp x, fp y, fp z) { return 0.0; }
 //******************************************************************************
 
 //
+// these functions compute the 2nd partial derivatives
+//	partial {mu,nu,phi} / partial {xx,xy,xz,yy,yz,zz}
+// as computed by the maple file "coord_derivs.{maple,out}" in this directory
+//
+fp  partial2_mu_wrt_xx(fp x, fp y, fp z) {return    0.0;}
+fp  partial2_mu_wrt_xy(fp x, fp y, fp z) {return    0.0;}
+fp  partial2_mu_wrt_xz(fp x, fp y, fp z) {return    0.0;}
+fp  partial2_mu_wrt_yy(fp x, fp y, fp z) {return   -2.0*y*z  / pow2(y*y + z*z);}
+fp  partial2_mu_wrt_yz(fp x, fp y, fp z) {return (y*y - z*z) / pow2(y*y + z*z);}
+fp  partial2_mu_wrt_zz(fp x, fp y, fp z) {return    2.0*y*z  / pow2(y*y + z*z);}
+
+fp  partial2_nu_wrt_xx(fp x, fp y, fp z) {return   -2.0*x*z  / pow2(x*x + z*z);}
+fp  partial2_nu_wrt_xy(fp x, fp y, fp z) {return    0.0;}
+fp  partial2_nu_wrt_xz(fp x, fp y, fp z) {return (x*x - z*z) / pow2(x*x + z*z);}
+fp  partial2_nu_wrt_yy(fp x, fp y, fp z) {return    0.0;}
+fp  partial2_nu_wrt_yz(fp x, fp y, fp z) {return    0.0;}
+fp  partial2_nu_wrt_zz(fp x, fp y, fp z) {return    2.0*x*z  / pow2(x*x + z*z);}
+
+fp partial2_phi_wrt_xx(fp x, fp y, fp z) {return    2.0*x*y  / pow2(x*x + y*y);}
+fp partial2_phi_wrt_xy(fp x, fp y, fp z) {return (y*y - x*x) / pow2(x*x + y*y);}
+fp partial2_phi_wrt_xz(fp x, fp y, fp z) {return    0.0;}
+fp partial2_phi_wrt_yy(fp x, fp y, fp z) {return   -2.0*x*y  / pow2(x*x + y*y);}
+fp partial2_phi_wrt_yz(fp x, fp y, fp z) {return    0.0;}
+fp partial2_phi_wrt_zz(fp x, fp y, fp z) {return    0.0;}
+
+//******************************************************************************
+//******************************************************************************
+//******************************************************************************
+
+//
 // these functions convert the usual polar spherical (theta,phi) <--> (x,y,z)
 //
 
diff --git a/src/patch/coords.hh b/src/patch/coords.hh
index ecf1124..12a6baa 100644
--- a/src/patch/coords.hh
+++ b/src/patch/coords.hh
@@ -148,6 +148,26 @@ fp partial_phi_wrt_x(fp x, fp y, fp z);
 fp partial_phi_wrt_y(fp x, fp y, fp z);
 fp partial_phi_wrt_z(fp x, fp y, fp z);
 
+// partial^2 {mu,nu,phi} / partial {x,y,z}{x,y,z}
+fp partial2_mu_wrt_xx(fp x, fp y, fp z);
+fp partial2_mu_wrt_xy(fp x, fp y, fp z);
+fp partial2_mu_wrt_xz(fp x, fp y, fp z);
+fp partial2_mu_wrt_yy(fp x, fp y, fp z);
+fp partial2_mu_wrt_yz(fp x, fp y, fp z);
+fp partial2_mu_wrt_zz(fp x, fp y, fp z);
+fp partial2_nu_wrt_xx(fp x, fp y, fp z);
+fp partial2_nu_wrt_xy(fp x, fp y, fp z);
+fp partial2_nu_wrt_xz(fp x, fp y, fp z);
+fp partial2_nu_wrt_yy(fp x, fp y, fp z);
+fp partial2_nu_wrt_yz(fp x, fp y, fp z);
+fp partial2_nu_wrt_zz(fp x, fp y, fp z);
+fp partial2_phi_wrt_xx(fp x, fp y, fp z);
+fp partial2_phi_wrt_xy(fp x, fp y, fp z);
+fp partial2_phi_wrt_xz(fp x, fp y, fp z);
+fp partial2_phi_wrt_yy(fp x, fp y, fp z);
+fp partial2_phi_wrt_yz(fp x, fp y, fp z);
+fp partial2_phi_wrt_zz(fp x, fp y, fp z);
+
 // usual polar spherical (r,theta,phi) <--> (x,y,z)
 void xyz_of_r_theta_phi(fp r, fp theta, fp phi,   fp& x, fp& y, fp& z);
 // ... already have r_of_xyz()