add functions to compute partial {x,y,z} / partial {mu,nu,phi}

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

2 files changed, 126 insertions, 0 deletions

diff --git a/src/patch/coords.cc b/src/patch/coords.cc
index 5b02c5b..67c41f9 100644
--- a/src/patch/coords.cc
+++ b/src/patch/coords.cc
@@ -8,6 +8,7 @@
 // (r,(mu,nu,phi)) <--> (x,y,z)
 // ((mu,nu,phi)) --> the 3rd
 //
+// partial_{x,y,z}_wrt_{mu,nu,phi}
 // partial_{mu,nu,phi}_wrt_{x,y,z}
 // partial2_{mu,nu,phi}_wrt_{xx,xy,xz,yy,yz,zz}
 //
@@ -263,6 +264,117 @@ fp phi_of_xy(fp x, fp y) { return arctan_xy(x,y); }
 
 //
 // these functions compute the partial derivatives
+//	partial {x,y,z} / partial {mu,nu,phi}
+//
+// Bugs: they're slow :(
+//
+
+namespace local_coords
+	  {
+void partial_xyz_wrt_mu_nu(fp r, fp mu, fp nu,
+			   fp& partial_x_wrt_mu, fp& partial_x_wrt_nu,
+			   fp& partial_y_wrt_mu, fp& partial_y_wrt_nu,
+			   fp& partial_z_wrt_mu, fp& partial_z_wrt_nu)
+{
+const fp tan_mu = tan(mu);
+const fp tan_nu = tan(nu);
+const fp r2 = pow2(r);
+const fp tan2_mu = pow2(tan_mu);
+const fp tan2_nu = pow2(tan_nu);
+
+fp x, y, z;
+xyz_of_r_mu_nu(r,mu,nu, x,y,z);
+
+const fp t = 1 + tan2_mu + tan2_nu;	// = $r^2/z^2$
+const fp t2 = pow2(t);
+const fp partial_t_wrt_mu = 2.0 * tan_mu * (1.0+tan2_mu);
+const fp partial_t_wrt_nu = 2.0 * tan_nu * (1.0+tan2_nu);
+
+partial_z_wrt_mu = -0.5 * (r2/(z*t2)) * partial_t_wrt_mu;
+partial_z_wrt_nu = -0.5 * (r2/(z*t2)) * partial_t_wrt_nu;
+
+partial_x_wrt_mu = tan_nu*partial_z_wrt_mu;
+partial_x_wrt_nu = tan_nu*partial_z_wrt_nu + z*(1.0+tan2_mu);
+partial_y_wrt_mu = tan_mu*partial_z_wrt_mu + z*(1.0+tan2_mu);
+partial_y_wrt_nu = tan_mu*partial_z_wrt_nu;
+}
+	  }
+
+//**************************************
+
+namespace local_coords
+	  {
+void partial_xyz_wrt_mu_phi(fp r, fp mu, fp phi,
+			    fp& partial_x_wrt_mu, fp& partial_x_wrt_phi,
+			    fp& partial_y_wrt_mu, fp& partial_y_wrt_phi,
+			    fp& partial_z_wrt_mu, fp& partial_z_wrt_phi)
+{
+const fp  mu_bar = 0.5*PI - mu ;
+const fp phi_bar = 0.5*PI - phi;
+
+const fp tan_mu_bar  = tan(mu_bar);
+const fp tan_phi_bar = tan(phi_bar);
+const fp r2 = pow2(r);
+const fp tan2_mu_bar  = pow2(tan_mu_bar);
+const fp tan2_phi_bar = pow2(tan_phi_bar);
+
+fp x, y, z;
+xyz_of_r_mu_phi(r,mu,phi, x,y,z);
+
+const fp t = 1 + tan2_mu_bar + tan2_phi_bar;	// = $r^2/y^2$
+const fp t2 = pow2(t);
+const fp partial_t_wrt_mu_bar  = 2.0 * tan_mu_bar  * (1.0+tan2_mu_bar );
+const fp partial_t_wrt_phi_bar = 2.0 * tan_phi_bar * (1.0+tan2_phi_bar);
+
+partial_y_wrt_mu  = 0.5 * (r2/(y*t2)) * partial_t_wrt_mu_bar;
+partial_y_wrt_phi = 0.5 * (r2/(y*t2)) * partial_t_wrt_phi_bar;
+
+partial_x_wrt_mu  = tan_phi_bar*partial_y_wrt_mu;
+partial_x_wrt_phi = tan_phi_bar*partial_y_wrt_phi - y*(1.0+tan2_phi_bar);
+partial_z_wrt_mu  = tan_mu_bar *partial_y_wrt_mu  - y*(1.0+tan2_mu_bar );
+partial_z_wrt_phi = tan_mu_bar *partial_y_wrt_phi;
+}
+	  }
+
+//**************************************
+
+namespace local_coords
+	  {
+void partial_xyz_wrt_nu_phi(fp r, fp nu, fp phi,
+			    fp& partial_x_wrt_nu, fp& partial_x_wrt_phi,
+			    fp& partial_y_wrt_nu, fp& partial_y_wrt_phi,
+			    fp& partial_z_wrt_nu, fp& partial_z_wrt_phi)
+{
+const fp  nu_bar = 0.5*PI - nu ;
+
+const fp tan_nu_bar = tan(nu_bar);
+const fp tan_phi    = tan(phi);
+const fp r2 = pow2(r);
+const fp tan2_nu_bar = pow2(tan_nu_bar);
+const fp tan2_phi    = pow2(tan_phi);
+
+fp x, y, z;
+xyz_of_r_nu_phi(r,nu,phi, x,y,z);
+
+const fp t = 1 + tan2_nu_bar + tan2_phi;	// = $r^2/x^2$
+const fp t2 = pow2(t);
+const fp partial_t_wrt_nu_bar  = 2.0 * tan_nu_bar * (1.0+tan2_nu_bar);
+const fp partial_t_wrt_phi     = 2.0 * tan_phi    * (1.0+tan2_phi   );
+
+partial_x_wrt_nu  =  0.5 * (r2/(x*t2)) * partial_t_wrt_nu_bar;
+partial_x_wrt_phi = -0.5 * (r2/(x*t2)) * partial_t_wrt_phi;
+
+partial_y_wrt_nu  = tan_phi   *partial_x_wrt_nu;
+partial_y_wrt_phi = tan_phi   *partial_x_wrt_phi + x*(1.0+tan2_phi   );
+partial_z_wrt_nu  = tan_nu_bar*partial_x_wrt_nu  - x*(1.0+tan2_nu_bar);
+partial_z_wrt_phi = tan_nu_bar*partial_x_wrt_phi;
+}
+	  }
+
+//******************************************************************************
+
+//
+// these functions compute the partial derivatives
 //	partial {mu,nu,phi} / partial {x,y,z}
 // as computed by the maple file "coord_derivs.{maple,out}" in this directory
 //
diff --git a/src/patch/coords.hh b/src/patch/coords.hh
index 54889da..bad88ea 100644
--- a/src/patch/coords.hh
+++ b/src/patch/coords.hh
@@ -146,6 +146,20 @@ fp phi_of_mu_nu(fp mu, fp nu );
 fp nu_of_mu_phi(fp mu, fp phi);
 fp mu_of_nu_phi(fp nu, fp phi);
 
+// partial {x,y,z} / partial {mu,nu,phi}
+void partial_xyz_wrt_mu_nu(fp r, fp mu, fp nu,
+			   fp& partial_x_wrt_mu, fp& partial_x_wrt_nu,
+			   fp& partial_y_wrt_mu, fp& partial_y_wrt_nu,
+			   fp& partial_z_wrt_mu, fp& partial_z_wrt_nu);
+void partial_xyz_wrt_mu_phi(fp r, fp mu, fp phi,
+			    fp& partial_x_wrt_mu, fp& partial_x_wrt_phi,
+			    fp& partial_y_wrt_mu, fp& partial_y_wrt_phi,
+			    fp& partial_z_wrt_mu, fp& partial_z_wrt_phi);
+void partial_xyz_wrt_nu_phi(fp r, fp nu, fp phi,
+			    fp& partial_x_wrt_nu, fp& partial_x_wrt_phi,
+			    fp& partial_y_wrt_nu, fp& partial_y_wrt_phi,
+			    fp& partial_z_wrt_nu, fp& partial_z_wrt_phi);
+
 // partial {mu,nu,phi} / partial {x,y,z}
 fp partial_mu_wrt_y(fp y, fp z);
 fp partial_mu_wrt_z(fp y, fp z);