WeylScal4:

	Remove variable modifications in Kranc script.
	Fix tensor parity tags.
	Add old-style BC compatibility.
	Add intelligence to runmath.sh in finding Kranc.
	Change bash script permissions.
	Add testsuites.
	Expand README description.



git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/WeylScal4/trunk@52 4f5cb9a8-4dd8-4c2d-9bbd-173fa4467843

4 files changed, 313 insertions, 87 deletions

diff --git a/src/Boundaries.c b/src/Boundaries.c
index ef1a009..b6e6cd0 100644
--- a/src/Boundaries.c
+++ b/src/Boundaries.c
@@ -31,11 +31,219 @@ void WeylScal4_SelectBoundConds(CCTK_ARGUMENTS)
   DECLARE_CCTK_PARAMETERS;
   
   CCTK_INT ierr = 0;
+  
+  if (CCTK_EQUALS(Psi4i_group_bound, "none"  ) ||
+      CCTK_EQUALS(Psi4i_group_bound, "static") ||
+      CCTK_EQUALS(Psi4i_group_bound, "flat"  ) ||
+      CCTK_EQUALS(Psi4i_group_bound, "zero"  ) ) 
+  {
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
+                      "WeylScal4::Psi4i_group", Psi4i_group_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Psi4i_group_bound BC for WeylScal4::Psi4i_group!");
+  }
+  
+  if (CCTK_EQUALS(Psi4r_group_bound, "none"  ) ||
+      CCTK_EQUALS(Psi4r_group_bound, "static") ||
+      CCTK_EQUALS(Psi4r_group_bound, "flat"  ) ||
+      CCTK_EQUALS(Psi4r_group_bound, "zero"  ) ) 
+  {
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
+                      "WeylScal4::Psi4r_group", Psi4r_group_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Psi4r_group_bound BC for WeylScal4::Psi4r_group!");
+  }
+  
+  if (CCTK_EQUALS(Psi4i_bound, "none"  ) ||
+      CCTK_EQUALS(Psi4i_bound, "static") ||
+      CCTK_EQUALS(Psi4i_bound, "flat"  ) ||
+      CCTK_EQUALS(Psi4i_bound, "zero"  ) ) 
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
+                      "WeylScal4::Psi4i", Psi4i_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Psi4i_bound BC for WeylScal4::Psi4i!");
+  }
+  
+  if (CCTK_EQUALS(Psi4r_bound, "none"  ) ||
+      CCTK_EQUALS(Psi4r_bound, "static") ||
+      CCTK_EQUALS(Psi4r_bound, "flat"  ) ||
+      CCTK_EQUALS(Psi4r_bound, "zero"  ) ) 
+  {
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
+                      "WeylScal4::Psi4r", Psi4r_bound);
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Psi4r_bound BC for WeylScal4::Psi4r!");
+  }
+  
+  if (CCTK_EQUALS(Psi4i_group_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_Psi4i_group_bound = -1;
+    if (handle_Psi4i_group_bound < 0) handle_Psi4i_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4i_group_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4i_group_bound , Psi4i_group_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_Psi4i_group_bound ,Psi4i_group_bound_speed, "SPEED") < 0)
+       CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_group_bound, 
+                      "WeylScal4::Psi4i_group", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4i_group!");
+  
+  }
+  
+  if (CCTK_EQUALS(Psi4r_group_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_Psi4r_group_bound = -1;
+    if (handle_Psi4r_group_bound < 0) handle_Psi4r_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4r_group_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4r_group_bound , Psi4r_group_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_Psi4r_group_bound ,Psi4r_group_bound_speed, "SPEED") < 0)
+       CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_group_bound, 
+                      "WeylScal4::Psi4r_group", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4r_group!");
+  
+  }
+  
+  if (CCTK_EQUALS(Psi4i_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_Psi4i_bound = -1;
+    if (handle_Psi4i_bound < 0) handle_Psi4i_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4i_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4i_bound , Psi4i_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_Psi4i_bound ,Psi4i_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_bound, 
+                      "WeylScal4::Psi4i", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4i!");
+  
+  }
+  
+  if (CCTK_EQUALS(Psi4r_bound, "radiative"))
+  {
+   /* select radiation boundary condition */
+    static CCTK_INT handle_Psi4r_bound = -1;
+    if (handle_Psi4r_bound < 0) handle_Psi4r_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4r_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4r_bound , Psi4r_bound_limit, "LIMIT") < 0)
+       CCTK_WARN(0, "could not set LIMIT value in table!");
+    if (Util_TableSetReal(handle_Psi4r_bound ,Psi4r_bound_speed, "SPEED") < 0)
+        CCTK_WARN(0, "could not set SPEED value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_bound, 
+                      "WeylScal4::Psi4r", "Radiation");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Radiation BC for WeylScal4::Psi4r!");
+  
+  }
+  
+  if (CCTK_EQUALS(Psi4i_group_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_Psi4i_group_bound = -1;
+    if (handle_Psi4i_group_bound < 0) handle_Psi4i_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4i_group_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4i_group_bound ,Psi4i_group_bound_scalar, "SCALAR") < 0)
+        CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_group_bound, 
+                      "WeylScal4::Psi4i_group", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Scalar BC for WeylScal4::Psi4i_group!");
+  
+  }
+  
+  if (CCTK_EQUALS(Psi4r_group_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_Psi4r_group_bound = -1;
+    if (handle_Psi4r_group_bound < 0) handle_Psi4r_group_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4r_group_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4r_group_bound ,Psi4r_group_bound_scalar, "SCALAR") < 0)
+        CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_group_bound, 
+                      "WeylScal4::Psi4r_group", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Failed to register Scalar BC for WeylScal4::Psi4r_group!");
+  
+  }
+  
+  if (CCTK_EQUALS(Psi4i_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_Psi4i_bound = -1;
+    if (handle_Psi4i_bound < 0) handle_Psi4i_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4i_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4i_bound ,Psi4i_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4i_bound, 
+                      "WeylScal4::Psi4i", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for WeylScal4::Psi4i!");
+  
+  }
+  
+  if (CCTK_EQUALS(Psi4r_bound, "scalar"))
+  {
+   /* select scalar boundary condition */
+    static CCTK_INT handle_Psi4r_bound = -1;
+    if (handle_Psi4r_bound < 0) handle_Psi4r_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_Psi4r_bound < 0) CCTK_WARN(0, "could not create table!");
+    if (Util_TableSetReal(handle_Psi4r_bound ,Psi4r_bound_scalar, "SCALAR") < 0)
+      CCTK_WARN(0, "could not set SCALAR value in table!");
+  
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_Psi4r_bound, 
+                      "WeylScal4::Psi4r", "scalar");
+  
+    if (ierr < 0)
+       CCTK_WARN(0, "Error in registering Scalar BC for WeylScal4::Psi4r!");
+  
+  }
   return;
 }
 
 
 
 /* template for entries in parameter file:
+#$bound$#WeylScal4::Psi4i_group_bound       = "skip"
+#$bound$#WeylScal4::Psi4i_group_bound_speed = 1.0
+#$bound$#WeylScal4::Psi4i_group_bound_limit = 0.0
+#$bound$#WeylScal4::Psi4i_group_bound_scalar = 0.0
+
+#$bound$#WeylScal4::Psi4r_group_bound       = "skip"
+#$bound$#WeylScal4::Psi4r_group_bound_speed = 1.0
+#$bound$#WeylScal4::Psi4r_group_bound_limit = 0.0
+#$bound$#WeylScal4::Psi4r_group_bound_scalar = 0.0
+
+#$bound$#WeylScal4::Psi4i_bound       = "skip"
+#$bound$#WeylScal4::Psi4i_bound_speed = 1.0
+#$bound$#WeylScal4::Psi4i_bound_limit = 0.0
+#$bound$#WeylScal4::Psi4i_bound_scalar = 0.0
+
+#$bound$#WeylScal4::Psi4r_bound       = "skip"
+#$bound$#WeylScal4::Psi4r_bound_speed = 1.0
+#$bound$#WeylScal4::Psi4r_bound_limit = 0.0
+#$bound$#WeylScal4::Psi4r_bound_scalar = 0.0
+
 */
 
diff --git a/src/RegisterSymmetries.c b/src/RegisterSymmetries.c
index 20511fd..df4787a 100644
--- a/src/RegisterSymmetries.c
+++ b/src/RegisterSymmetries.c
@@ -16,9 +16,9 @@ void WeylScal4_RegisterSymmetries(CCTK_ARGUMENTS)
   
   
   /* Register symmetries of grid functions */
-  sym[0] = 1;
-  sym[1] = 1;
-  sym[2] = 1;
+  sym[0] = -1;
+  sym[1] = -1;
+  sym[2] = -1;
   SetCartSymVN(cctkGH, sym, "WeylScal4::Psi4i");
   
   sym[0] = 1;
diff --git a/src/psis_calc_2nd.c b/src/psis_calc_2nd.c
index f97fa78..193c552 100644
--- a/src/psis_calc_2nd.c
+++ b/src/psis_calc_2nd.c
@@ -353,85 +353,94 @@ void psis_calc_2nd_Body(cGH const * restrict const cctkGH, int const dir, int co
         CCTK_REAL vc3 = ((-(gInv33*va2) + gInv32*va3)*vb1 + (gInv33*va1 - 
           gInv31*va3)*vb2 + (-(gInv32*va1) + gInv31*va2)*vb3)*pow(detg,0.5);
         
-        CCTK_REAL omega11 = 2*(gyzL*va2*va3 + va1*(gxyL*va2 + gxzL*va3)) + 
-          gxxL*SQR(va1) + gyyL*SQR(va2) + gzzL*SQR(va3);
+        CCTK_REAL wa1 = va1;
         
-        va1 = va1*pow(omega11,-khalf);
+        CCTK_REAL wa2 = va2;
         
-        va2 = va2*pow(omega11,-khalf);
+        CCTK_REAL wa3 = va3;
         
-        va3 = va3*pow(omega11,-khalf);
+        CCTK_REAL omega11 = 2*(gyzL*wa2*wa3 + wa1*(gxyL*wa2 + gxzL*wa3)) + 
+          gxxL*SQR(wa1) + gyyL*SQR(wa2) + gzzL*SQR(wa3);
         
-        CCTK_REAL omega12 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vb1 + (gxyL*va1 + 
-          gyyL*va2 + gyzL*va3)*vb2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vb3;
+        CCTK_REAL ea1 = wa1*pow(omega11,-khalf);
         
-        vb1 = -(omega12*va1) + vb1;
+        CCTK_REAL ea2 = wa2*pow(omega11,-khalf);
         
-        vb2 = -(omega12*va2) + vb2;
+        CCTK_REAL ea3 = wa3*pow(omega11,-khalf);
         
-        vb3 = -(omega12*va3) + vb3;
+        CCTK_REAL omega12 = ea1*(gxxL*vb1 + gxyL*vb2 + gxzL*vb3) + 
+          ea2*(gxyL*vb1 + gyyL*vb2 + gyzL*vb3) + ea3*(gxzL*vb1 + gyzL*vb2 + 
+          gzzL*vb3);
         
-        CCTK_REAL omega22 = 2*(gyzL*vb2*vb3 + vb1*(gxyL*vb2 + gxzL*vb3)) + 
-          gxxL*SQR(vb1) + gyyL*SQR(vb2) + gzzL*SQR(vb3);
+        CCTK_REAL wb1 = -(ea1*omega12) + vb1;
         
-        vb1 = vb1*pow(omega22,-khalf);
+        CCTK_REAL wb2 = -(ea2*omega12) + vb2;
         
-        vb2 = vb2*pow(omega22,-khalf);
+        CCTK_REAL wb3 = -(ea3*omega12) + vb3;
         
-        vb3 = vb3*pow(omega22,-khalf);
+        CCTK_REAL omega22 = 2*(gyzL*wb2*wb3 + wb1*(gxyL*wb2 + gxzL*wb3)) + 
+          gxxL*SQR(wb1) + gyyL*SQR(wb2) + gzzL*SQR(wb3);
         
-        CCTK_REAL omega13 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vc1 + (gxyL*va1 + 
-          gyyL*va2 + gyzL*va3)*vc2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vc3;
+        CCTK_REAL eb1 = wb1*pow(omega22,-khalf);
         
-        CCTK_REAL omega23 = (gxxL*vb1 + gxyL*vb2 + gxzL*vb3)*vc1 + (gxyL*vb1 + 
-          gyyL*vb2 + gyzL*vb3)*vc2 + (gxzL*vb1 + gyzL*vb2 + gzzL*vb3)*vc3;
+        CCTK_REAL eb2 = wb2*pow(omega22,-khalf);
         
-        vc1 = -(omega13*va1) - omega23*vb1 + vc1;
+        CCTK_REAL eb3 = wb3*pow(omega22,-khalf);
         
-        vc2 = -(omega13*va2) - omega23*vb2 + vc2;
+        CCTK_REAL omega13 = ea1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + 
+          ea2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + ea3*(gxzL*vc1 + gyzL*vc2 + 
+          gzzL*vc3);
         
-        vc3 = -(omega13*va3) - omega23*vb3 + vc3;
+        CCTK_REAL omega23 = eb1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + 
+          eb2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + eb3*(gxzL*vc1 + gyzL*vc2 + 
+          gzzL*vc3);
         
-        CCTK_REAL omega33 = 2*(gyzL*vc2*vc3 + vc1*(gxyL*vc2 + gxzL*vc3)) + 
-          gxxL*SQR(vc1) + gyyL*SQR(vc2) + gzzL*SQR(vc3);
+        CCTK_REAL wc1 = -(ea1*omega13) - eb1*omega23 + vc1;
         
-        vc1 = vc1*pow(omega33,-khalf);
+        CCTK_REAL wc2 = -(ea2*omega13) - eb2*omega23 + vc2;
         
-        vc2 = vc2*pow(omega33,-khalf);
+        CCTK_REAL wc3 = -(ea3*omega13) - eb3*omega23 + vc3;
         
-        vc3 = vc3*pow(omega33,-khalf);
+        CCTK_REAL omega33 = 2*(gyzL*wc2*wc3 + wc1*(gxyL*wc2 + gxzL*wc3)) + 
+          gxxL*SQR(wc1) + gyyL*SQR(wc2) + gzzL*SQR(wc3);
+        
+        CCTK_REAL ec1 = wc1*pow(omega33,-khalf);
+        
+        CCTK_REAL ec2 = wc2*pow(omega33,-khalf);
+        
+        CCTK_REAL ec3 = wc3*pow(omega33,-khalf);
         
         CCTK_REAL isqrt2 = 0.707106781186547524;
         
-        CCTK_REAL n1 = -(isqrt2*vb1);
+        CCTK_REAL n1 = -(eb1*isqrt2);
         
-        CCTK_REAL n2 = -(isqrt2*vb2);
+        CCTK_REAL n2 = -(eb2*isqrt2);
         
-        CCTK_REAL n3 = -(isqrt2*vb3);
+        CCTK_REAL n3 = -(eb3*isqrt2);
         
-        CCTK_REAL rm1 = isqrt2*vc1;
+        CCTK_REAL rm1 = ec1*isqrt2;
         
-        CCTK_REAL rm2 = isqrt2*vc2;
+        CCTK_REAL rm2 = ec2*isqrt2;
         
-        CCTK_REAL rm3 = isqrt2*vc3;
+        CCTK_REAL rm3 = ec3*isqrt2;
         
-        CCTK_REAL im1 = isqrt2*va1;
+        CCTK_REAL im1 = ea1*isqrt2;
         
-        CCTK_REAL im2 = isqrt2*va2;
+        CCTK_REAL im2 = ea2*isqrt2;
         
-        CCTK_REAL im3 = isqrt2*va3;
+        CCTK_REAL im3 = ea3*isqrt2;
         
-        CCTK_REAL rmbar1 = isqrt2*vc1;
+        CCTK_REAL rmbar1 = ec1*isqrt2;
         
-        CCTK_REAL rmbar2 = isqrt2*vc2;
+        CCTK_REAL rmbar2 = ec2*isqrt2;
         
-        CCTK_REAL rmbar3 = isqrt2*vc3;
+        CCTK_REAL rmbar3 = ec3*isqrt2;
         
-        CCTK_REAL imbar1 = -(isqrt2*va1);
+        CCTK_REAL imbar1 = -(ea1*isqrt2);
         
-        CCTK_REAL imbar2 = -(isqrt2*va2);
+        CCTK_REAL imbar2 = -(ea2*isqrt2);
         
-        CCTK_REAL imbar3 = -(isqrt2*va3);
+        CCTK_REAL imbar3 = -(ea3*isqrt2);
         
         CCTK_REAL nn = isqrt2;
         
diff --git a/src/psis_calc_4th.c b/src/psis_calc_4th.c
index 3be80e2..32fd7c2 100644
--- a/src/psis_calc_4th.c
+++ b/src/psis_calc_4th.c
@@ -353,85 +353,94 @@ void psis_calc_4th_Body(cGH const * restrict const cctkGH, int const dir, int co
         CCTK_REAL vc3 = ((-(gInv33*va2) + gInv32*va3)*vb1 + (gInv33*va1 - 
           gInv31*va3)*vb2 + (-(gInv32*va1) + gInv31*va2)*vb3)*pow(detg,0.5);
         
-        CCTK_REAL omega11 = 2*(gyzL*va2*va3 + va1*(gxyL*va2 + gxzL*va3)) + 
-          gxxL*SQR(va1) + gyyL*SQR(va2) + gzzL*SQR(va3);
+        CCTK_REAL wa1 = va1;
         
-        va1 = va1*pow(omega11,-khalf);
+        CCTK_REAL wa2 = va2;
         
-        va2 = va2*pow(omega11,-khalf);
+        CCTK_REAL wa3 = va3;
         
-        va3 = va3*pow(omega11,-khalf);
+        CCTK_REAL omega11 = 2*(gyzL*wa2*wa3 + wa1*(gxyL*wa2 + gxzL*wa3)) + 
+          gxxL*SQR(wa1) + gyyL*SQR(wa2) + gzzL*SQR(wa3);
         
-        CCTK_REAL omega12 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vb1 + (gxyL*va1 + 
-          gyyL*va2 + gyzL*va3)*vb2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vb3;
+        CCTK_REAL ea1 = wa1*pow(omega11,-khalf);
         
-        vb1 = -(omega12*va1) + vb1;
+        CCTK_REAL ea2 = wa2*pow(omega11,-khalf);
         
-        vb2 = -(omega12*va2) + vb2;
+        CCTK_REAL ea3 = wa3*pow(omega11,-khalf);
         
-        vb3 = -(omega12*va3) + vb3;
+        CCTK_REAL omega12 = ea1*(gxxL*vb1 + gxyL*vb2 + gxzL*vb3) + 
+          ea2*(gxyL*vb1 + gyyL*vb2 + gyzL*vb3) + ea3*(gxzL*vb1 + gyzL*vb2 + 
+          gzzL*vb3);
         
-        CCTK_REAL omega22 = 2*(gyzL*vb2*vb3 + vb1*(gxyL*vb2 + gxzL*vb3)) + 
-          gxxL*SQR(vb1) + gyyL*SQR(vb2) + gzzL*SQR(vb3);
+        CCTK_REAL wb1 = -(ea1*omega12) + vb1;
         
-        vb1 = vb1*pow(omega22,-khalf);
+        CCTK_REAL wb2 = -(ea2*omega12) + vb2;
         
-        vb2 = vb2*pow(omega22,-khalf);
+        CCTK_REAL wb3 = -(ea3*omega12) + vb3;
         
-        vb3 = vb3*pow(omega22,-khalf);
+        CCTK_REAL omega22 = 2*(gyzL*wb2*wb3 + wb1*(gxyL*wb2 + gxzL*wb3)) + 
+          gxxL*SQR(wb1) + gyyL*SQR(wb2) + gzzL*SQR(wb3);
         
-        CCTK_REAL omega13 = (gxxL*va1 + gxyL*va2 + gxzL*va3)*vc1 + (gxyL*va1 + 
-          gyyL*va2 + gyzL*va3)*vc2 + (gxzL*va1 + gyzL*va2 + gzzL*va3)*vc3;
+        CCTK_REAL eb1 = wb1*pow(omega22,-khalf);
         
-        CCTK_REAL omega23 = (gxxL*vb1 + gxyL*vb2 + gxzL*vb3)*vc1 + (gxyL*vb1 + 
-          gyyL*vb2 + gyzL*vb3)*vc2 + (gxzL*vb1 + gyzL*vb2 + gzzL*vb3)*vc3;
+        CCTK_REAL eb2 = wb2*pow(omega22,-khalf);
         
-        vc1 = -(omega13*va1) - omega23*vb1 + vc1;
+        CCTK_REAL eb3 = wb3*pow(omega22,-khalf);
         
-        vc2 = -(omega13*va2) - omega23*vb2 + vc2;
+        CCTK_REAL omega13 = ea1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + 
+          ea2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + ea3*(gxzL*vc1 + gyzL*vc2 + 
+          gzzL*vc3);
         
-        vc3 = -(omega13*va3) - omega23*vb3 + vc3;
+        CCTK_REAL omega23 = eb1*(gxxL*vc1 + gxyL*vc2 + gxzL*vc3) + 
+          eb2*(gxyL*vc1 + gyyL*vc2 + gyzL*vc3) + eb3*(gxzL*vc1 + gyzL*vc2 + 
+          gzzL*vc3);
         
-        CCTK_REAL omega33 = 2*(gyzL*vc2*vc3 + vc1*(gxyL*vc2 + gxzL*vc3)) + 
-          gxxL*SQR(vc1) + gyyL*SQR(vc2) + gzzL*SQR(vc3);
+        CCTK_REAL wc1 = -(ea1*omega13) - eb1*omega23 + vc1;
         
-        vc1 = vc1*pow(omega33,-khalf);
+        CCTK_REAL wc2 = -(ea2*omega13) - eb2*omega23 + vc2;
         
-        vc2 = vc2*pow(omega33,-khalf);
+        CCTK_REAL wc3 = -(ea3*omega13) - eb3*omega23 + vc3;
         
-        vc3 = vc3*pow(omega33,-khalf);
+        CCTK_REAL omega33 = 2*(gyzL*wc2*wc3 + wc1*(gxyL*wc2 + gxzL*wc3)) + 
+          gxxL*SQR(wc1) + gyyL*SQR(wc2) + gzzL*SQR(wc3);
+        
+        CCTK_REAL ec1 = wc1*pow(omega33,-khalf);
+        
+        CCTK_REAL ec2 = wc2*pow(omega33,-khalf);
+        
+        CCTK_REAL ec3 = wc3*pow(omega33,-khalf);
         
         CCTK_REAL isqrt2 = 0.707106781186547524;
         
-        CCTK_REAL n1 = -(isqrt2*vb1);
+        CCTK_REAL n1 = -(eb1*isqrt2);
         
-        CCTK_REAL n2 = -(isqrt2*vb2);
+        CCTK_REAL n2 = -(eb2*isqrt2);
         
-        CCTK_REAL n3 = -(isqrt2*vb3);
+        CCTK_REAL n3 = -(eb3*isqrt2);
         
-        CCTK_REAL rm1 = isqrt2*vc1;
+        CCTK_REAL rm1 = ec1*isqrt2;
         
-        CCTK_REAL rm2 = isqrt2*vc2;
+        CCTK_REAL rm2 = ec2*isqrt2;
         
-        CCTK_REAL rm3 = isqrt2*vc3;
+        CCTK_REAL rm3 = ec3*isqrt2;
         
-        CCTK_REAL im1 = isqrt2*va1;
+        CCTK_REAL im1 = ea1*isqrt2;
         
-        CCTK_REAL im2 = isqrt2*va2;
+        CCTK_REAL im2 = ea2*isqrt2;
         
-        CCTK_REAL im3 = isqrt2*va3;
+        CCTK_REAL im3 = ea3*isqrt2;
         
-        CCTK_REAL rmbar1 = isqrt2*vc1;
+        CCTK_REAL rmbar1 = ec1*isqrt2;
         
-        CCTK_REAL rmbar2 = isqrt2*vc2;
+        CCTK_REAL rmbar2 = ec2*isqrt2;
         
-        CCTK_REAL rmbar3 = isqrt2*vc3;
+        CCTK_REAL rmbar3 = ec3*isqrt2;
         
-        CCTK_REAL imbar1 = -(isqrt2*va1);
+        CCTK_REAL imbar1 = -(ea1*isqrt2);
         
-        CCTK_REAL imbar2 = -(isqrt2*va2);
+        CCTK_REAL imbar2 = -(ea2*isqrt2);
         
-        CCTK_REAL imbar3 = -(isqrt2*va3);
+        CCTK_REAL imbar3 = -(ea3*isqrt2);
         
         CCTK_REAL nn = isqrt2;