Continued development on the BSSN constraints

Rename dtalpha and dtbeta to A and B.

Add alternative formulae for the conformal Ricci tensor as comments.

13 files changed, 291 insertions, 510 deletions

diff --git a/ML_ADM/schedule.ccl b/ML_ADM/schedule.ccl
index 0aa9f9d..314d2b5 100644
--- a/ML_ADM/schedule.ccl
+++ b/ML_ADM/schedule.ccl
@@ -78,7 +78,7 @@ schedule ML_ADM_RHS AT analysis
   SYNC: shiftrhs
 } "ML_ADM_RHS"
 
-schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER ADM_ApplyBoundConds
+schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER ML_ADM_ApplyBCs
 {
   LANG: C
   
diff --git a/ML_BSSN/interface.ccl b/ML_BSSN/interface.ccl
index dd11e9f..2e9cc98 100644
--- a/ML_BSSN/interface.ccl
+++ b/ML_BSSN/interface.ccl
@@ -78,15 +78,15 @@ CCTK_REAL curvrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweig
 public:
 CCTK_REAL dtlapserhs type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
-  dtalpharhs
+  Arhs
 } "dtlapserhs"
 
 public:
 CCTK_REAL dtshiftrhs type=GF timelevels=1 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
-  dtbeta1rhs,
-  dtbeta2rhs,
-  dtbeta3rhs
+  B1rhs,
+  B2rhs,
+  B3rhs
 } "dtshiftrhs"
 
 public:
@@ -148,15 +148,15 @@ CCTK_REAL curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=
 public:
 CCTK_REAL dtlapse type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
-  dtalpha
+  A
 } "dtlapse"
 
 public:
 CCTK_REAL dtshift type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
-  dtbeta1,
-  dtbeta2,
-  dtbeta3
+  B1,
+  B2,
+  B3
 } "dtshift"
 
 public:
diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl
index a7f2534..5b5a7d0 100644
--- a/ML_BSSN/param.ccl
+++ b/ML_BSSN/param.ccl
@@ -241,7 +241,7 @@ KEYWORD At33_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD dtalpha_bound "Boundary condition to implement"
+KEYWORD A_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -253,7 +253,7 @@ KEYWORD dtalpha_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD dtbeta1_bound "Boundary condition to implement"
+KEYWORD B1_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -265,7 +265,7 @@ KEYWORD dtbeta1_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD dtbeta2_bound "Boundary condition to implement"
+KEYWORD B2_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -277,7 +277,7 @@ KEYWORD dtbeta2_bound "Boundary condition to implement"
 } "skip"
 
 private:
-KEYWORD dtbeta3_bound "Boundary condition to implement"
+KEYWORD B3_bound "Boundary condition to implement"
 {
   "flat" :: "Flat boundary condition"
   "none" :: "No boundary condition"
@@ -613,25 +613,25 @@ CCTK_REAL At33_bound_speed "characteristic speed at boundary"
 } 1.
 
 private:
-CCTK_REAL dtalpha_bound_speed "characteristic speed at boundary"
+CCTK_REAL A_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
 
 private:
-CCTK_REAL dtbeta1_bound_speed "characteristic speed at boundary"
+CCTK_REAL B1_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
 
 private:
-CCTK_REAL dtbeta2_bound_speed "characteristic speed at boundary"
+CCTK_REAL B2_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
 
 private:
-CCTK_REAL dtbeta3_bound_speed "characteristic speed at boundary"
+CCTK_REAL B3_bound_speed "characteristic speed at boundary"
 {
   "0:*" :: "outgoing characteristic speed > 0"
 } 1.
@@ -817,25 +817,25 @@ CCTK_REAL At33_bound_limit "limit value for r -> infinity"
 } 0.
 
 private:
-CCTK_REAL dtalpha_bound_limit "limit value for r -> infinity"
+CCTK_REAL A_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
 
 private:
-CCTK_REAL dtbeta1_bound_limit "limit value for r -> infinity"
+CCTK_REAL B1_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
 
 private:
-CCTK_REAL dtbeta2_bound_limit "limit value for r -> infinity"
+CCTK_REAL B2_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
 
 private:
-CCTK_REAL dtbeta3_bound_limit "limit value for r -> infinity"
+CCTK_REAL B3_bound_limit "limit value for r -> infinity"
 {
   "*:*" :: "value of limit value is unrestricted"
 } 0.
@@ -1021,25 +1021,25 @@ CCTK_REAL At33_bound_scalar "Dirichlet boundary value"
 } 0.
 
 private:
-CCTK_REAL dtalpha_bound_scalar "Dirichlet boundary value"
+CCTK_REAL A_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
 
 private:
-CCTK_REAL dtbeta1_bound_scalar "Dirichlet boundary value"
+CCTK_REAL B1_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
 
 private:
-CCTK_REAL dtbeta2_bound_scalar "Dirichlet boundary value"
+CCTK_REAL B2_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
 
 private:
-CCTK_REAL dtbeta3_bound_scalar "Dirichlet boundary value"
+CCTK_REAL B3_bound_scalar "Dirichlet boundary value"
 {
   "*:*" :: "unrestricted"
 } 0.
diff --git a/ML_BSSN/schedule.ccl b/ML_BSSN/schedule.ccl
index a185b4d..4be8b1f 100644
--- a/ML_BSSN/schedule.ccl
+++ b/ML_BSSN/schedule.ccl
@@ -120,13 +120,13 @@ schedule ML_BSSN_RHS AT analysis
   SYNC: trace_curvrhs
 } "ML_BSSN_RHS"
 
-schedule ML_BSSN_enforce IN MoL_PostStep AFTER BSSN_ApplyBoundConds
+schedule ML_BSSN_enforce IN MoL_PostStep BEFORE ML_BSSN_BoundConds
 {
   LANG: C
   
 } "ML_BSSN_enforce"
 
-schedule ML_BSSN_convertToADMBase IN MoL_PostStep AFTER ML_BSSN_ApplyBoundConds AFTER ML_BSSN_enforce
+schedule ML_BSSN_convertToADMBase IN MoL_PostStep AFTER ML_BSSN_ApplyBCs AFTER ML_BSSN_enforce
 {
   LANG: C
   
diff --git a/ML_BSSN/src/Boundaries.c b/ML_BSSN/src/Boundaries.c
index 03838fd..eac8cbd 100644
--- a/ML_BSSN/src/Boundaries.c
+++ b/ML_BSSN/src/Boundaries.c
@@ -200,48 +200,48 @@ void ML_BSSN_ApplyBoundConds(CCTK_ARGUMENTS)
        CCTK_WARN(-1, "Failed to register At33_bound BC for ML_BSSN::At33!");
   }
   
-  if (CCTK_EQUALS(dtalpha_bound, "none"  ) ||
-      CCTK_EQUALS(dtalpha_bound, "static") ||
-      CCTK_EQUALS(dtalpha_bound, "flat"  ) ||
-      CCTK_EQUALS(dtalpha_bound, "zero"  ) ) 
+  if (CCTK_EQUALS(A_bound, "none"  ) ||
+      CCTK_EQUALS(A_bound, "static") ||
+      CCTK_EQUALS(A_bound, "flat"  ) ||
+      CCTK_EQUALS(A_bound, "zero"  ) ) 
   {
     ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
-                      "ML_BSSN::dtalpha", dtalpha_bound);
+                      "ML_BSSN::A", A_bound);
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register dtalpha_bound BC for ML_BSSN::dtalpha!");
+       CCTK_WARN(-1, "Failed to register A_bound BC for ML_BSSN::A!");
   }
   
-  if (CCTK_EQUALS(dtbeta1_bound, "none"  ) ||
-      CCTK_EQUALS(dtbeta1_bound, "static") ||
-      CCTK_EQUALS(dtbeta1_bound, "flat"  ) ||
-      CCTK_EQUALS(dtbeta1_bound, "zero"  ) ) 
+  if (CCTK_EQUALS(B1_bound, "none"  ) ||
+      CCTK_EQUALS(B1_bound, "static") ||
+      CCTK_EQUALS(B1_bound, "flat"  ) ||
+      CCTK_EQUALS(B1_bound, "zero"  ) ) 
   {
     ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
-                      "ML_BSSN::dtbeta1", dtbeta1_bound);
+                      "ML_BSSN::B1", B1_bound);
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register dtbeta1_bound BC for ML_BSSN::dtbeta1!");
+       CCTK_WARN(-1, "Failed to register B1_bound BC for ML_BSSN::B1!");
   }
   
-  if (CCTK_EQUALS(dtbeta2_bound, "none"  ) ||
-      CCTK_EQUALS(dtbeta2_bound, "static") ||
-      CCTK_EQUALS(dtbeta2_bound, "flat"  ) ||
-      CCTK_EQUALS(dtbeta2_bound, "zero"  ) ) 
+  if (CCTK_EQUALS(B2_bound, "none"  ) ||
+      CCTK_EQUALS(B2_bound, "static") ||
+      CCTK_EQUALS(B2_bound, "flat"  ) ||
+      CCTK_EQUALS(B2_bound, "zero"  ) ) 
   {
     ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
-                      "ML_BSSN::dtbeta2", dtbeta2_bound);
+                      "ML_BSSN::B2", B2_bound);
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register dtbeta2_bound BC for ML_BSSN::dtbeta2!");
+       CCTK_WARN(-1, "Failed to register B2_bound BC for ML_BSSN::B2!");
   }
   
-  if (CCTK_EQUALS(dtbeta3_bound, "none"  ) ||
-      CCTK_EQUALS(dtbeta3_bound, "static") ||
-      CCTK_EQUALS(dtbeta3_bound, "flat"  ) ||
-      CCTK_EQUALS(dtbeta3_bound, "zero"  ) ) 
+  if (CCTK_EQUALS(B3_bound, "none"  ) ||
+      CCTK_EQUALS(B3_bound, "static") ||
+      CCTK_EQUALS(B3_bound, "flat"  ) ||
+      CCTK_EQUALS(B3_bound, "zero"  ) ) 
   {
     ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, 
-                      "ML_BSSN::dtbeta3", dtbeta3_bound);
+                      "ML_BSSN::B3", B3_bound);
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register dtbeta3_bound BC for ML_BSSN::dtbeta3!");
+       CCTK_WARN(-1, "Failed to register B3_bound BC for ML_BSSN::B3!");
   }
   
   if (CCTK_EQUALS(Xt1_bound, "none"  ) ||
@@ -679,75 +679,75 @@ void ML_BSSN_ApplyBoundConds(CCTK_ARGUMENTS)
   
   }
   
-  if (CCTK_EQUALS(dtalpha_bound, "radiative"))
+  if (CCTK_EQUALS(A_bound, "radiative"))
   {
    /* apply radiation boundary condition */
-    CCTK_INT handle_dtalpha_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtalpha_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtalpha_bound , dtalpha_bound_limit, "LIMIT") < 0)
+    CCTK_INT handle_A_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_A_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_A_bound , A_bound_limit, "LIMIT") < 0)
        CCTK_WARN(-1, "could not set LIMIT value in table!");
-    if (Util_TableSetReal(handle_dtalpha_bound ,dtalpha_bound_speed, "SPEED") < 0)
+    if (Util_TableSetReal(handle_A_bound ,A_bound_speed, "SPEED") < 0)
         CCTK_WARN(-1, "could not set SPEED value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtalpha_bound, 
-                      "ML_BSSN::dtalpha", "Radiation");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_A_bound, 
+                      "ML_BSSN::A", "Radiation");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::dtalpha!");
+       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::A!");
   
   }
   
-  if (CCTK_EQUALS(dtbeta1_bound, "radiative"))
+  if (CCTK_EQUALS(B1_bound, "radiative"))
   {
    /* apply radiation boundary condition */
-    CCTK_INT handle_dtbeta1_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtbeta1_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtbeta1_bound , dtbeta1_bound_limit, "LIMIT") < 0)
+    CCTK_INT handle_B1_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_B1_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_B1_bound , B1_bound_limit, "LIMIT") < 0)
        CCTK_WARN(-1, "could not set LIMIT value in table!");
-    if (Util_TableSetReal(handle_dtbeta1_bound ,dtbeta1_bound_speed, "SPEED") < 0)
+    if (Util_TableSetReal(handle_B1_bound ,B1_bound_speed, "SPEED") < 0)
         CCTK_WARN(-1, "could not set SPEED value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtbeta1_bound, 
-                      "ML_BSSN::dtbeta1", "Radiation");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_B1_bound, 
+                      "ML_BSSN::B1", "Radiation");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::dtbeta1!");
+       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::B1!");
   
   }
   
-  if (CCTK_EQUALS(dtbeta2_bound, "radiative"))
+  if (CCTK_EQUALS(B2_bound, "radiative"))
   {
    /* apply radiation boundary condition */
-    CCTK_INT handle_dtbeta2_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtbeta2_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtbeta2_bound , dtbeta2_bound_limit, "LIMIT") < 0)
+    CCTK_INT handle_B2_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_B2_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_B2_bound , B2_bound_limit, "LIMIT") < 0)
        CCTK_WARN(-1, "could not set LIMIT value in table!");
-    if (Util_TableSetReal(handle_dtbeta2_bound ,dtbeta2_bound_speed, "SPEED") < 0)
+    if (Util_TableSetReal(handle_B2_bound ,B2_bound_speed, "SPEED") < 0)
         CCTK_WARN(-1, "could not set SPEED value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtbeta2_bound, 
-                      "ML_BSSN::dtbeta2", "Radiation");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_B2_bound, 
+                      "ML_BSSN::B2", "Radiation");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::dtbeta2!");
+       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::B2!");
   
   }
   
-  if (CCTK_EQUALS(dtbeta3_bound, "radiative"))
+  if (CCTK_EQUALS(B3_bound, "radiative"))
   {
    /* apply radiation boundary condition */
-    CCTK_INT handle_dtbeta3_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtbeta3_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtbeta3_bound , dtbeta3_bound_limit, "LIMIT") < 0)
+    CCTK_INT handle_B3_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_B3_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_B3_bound , B3_bound_limit, "LIMIT") < 0)
        CCTK_WARN(-1, "could not set LIMIT value in table!");
-    if (Util_TableSetReal(handle_dtbeta3_bound ,dtbeta3_bound_speed, "SPEED") < 0)
+    if (Util_TableSetReal(handle_B3_bound ,B3_bound_speed, "SPEED") < 0)
         CCTK_WARN(-1, "could not set SPEED value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtbeta3_bound, 
-                      "ML_BSSN::dtbeta3", "Radiation");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_B3_bound, 
+                      "ML_BSSN::B3", "Radiation");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::dtbeta3!");
+       CCTK_WARN(-1, "Failed to register Radiation BC for ML_BSSN::B3!");
   
   }
   
@@ -1261,67 +1261,67 @@ void ML_BSSN_ApplyBoundConds(CCTK_ARGUMENTS)
   
   }
   
-  if (CCTK_EQUALS(dtalpha_bound, "scalar"))
+  if (CCTK_EQUALS(A_bound, "scalar"))
   {
    /* apply scalar boundary condition */
-    CCTK_INT handle_dtalpha_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtalpha_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtalpha_bound ,dtalpha_bound_scalar, "SCALAR") < 0)
+    CCTK_INT handle_A_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_A_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_A_bound ,A_bound_scalar, "SCALAR") < 0)
       CCTK_WARN(-1, "could not set SCALAR value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtalpha_bound, 
-                      "ML_BSSN::dtalpha", "scalar");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_A_bound, 
+                      "ML_BSSN::A", "scalar");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::dtalpha!");
+       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::A!");
   
   }
   
-  if (CCTK_EQUALS(dtbeta1_bound, "scalar"))
+  if (CCTK_EQUALS(B1_bound, "scalar"))
   {
    /* apply scalar boundary condition */
-    CCTK_INT handle_dtbeta1_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtbeta1_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtbeta1_bound ,dtbeta1_bound_scalar, "SCALAR") < 0)
+    CCTK_INT handle_B1_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_B1_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_B1_bound ,B1_bound_scalar, "SCALAR") < 0)
       CCTK_WARN(-1, "could not set SCALAR value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtbeta1_bound, 
-                      "ML_BSSN::dtbeta1", "scalar");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_B1_bound, 
+                      "ML_BSSN::B1", "scalar");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::dtbeta1!");
+       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::B1!");
   
   }
   
-  if (CCTK_EQUALS(dtbeta2_bound, "scalar"))
+  if (CCTK_EQUALS(B2_bound, "scalar"))
   {
    /* apply scalar boundary condition */
-    CCTK_INT handle_dtbeta2_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtbeta2_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtbeta2_bound ,dtbeta2_bound_scalar, "SCALAR") < 0)
+    CCTK_INT handle_B2_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_B2_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_B2_bound ,B2_bound_scalar, "SCALAR") < 0)
       CCTK_WARN(-1, "could not set SCALAR value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtbeta2_bound, 
-                      "ML_BSSN::dtbeta2", "scalar");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_B2_bound, 
+                      "ML_BSSN::B2", "scalar");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::dtbeta2!");
+       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::B2!");
   
   }
   
-  if (CCTK_EQUALS(dtbeta3_bound, "scalar"))
+  if (CCTK_EQUALS(B3_bound, "scalar"))
   {
    /* apply scalar boundary condition */
-    CCTK_INT handle_dtbeta3_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
-    if (handle_dtbeta3_bound < 0) CCTK_WARN(-1, "could not create table!");
-    if (Util_TableSetReal(handle_dtbeta3_bound ,dtbeta3_bound_scalar, "SCALAR") < 0)
+    CCTK_INT handle_B3_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE);
+    if (handle_B3_bound < 0) CCTK_WARN(-1, "could not create table!");
+    if (Util_TableSetReal(handle_B3_bound ,B3_bound_scalar, "SCALAR") < 0)
       CCTK_WARN(-1, "could not set SCALAR value in table!");
   
-    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_dtbeta3_bound, 
-                      "ML_BSSN::dtbeta3", "scalar");
+    ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, handle_B3_bound, 
+                      "ML_BSSN::B3", "scalar");
   
     if (ierr < 0)
-       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::dtbeta3!");
+       CCTK_WARN(-1, "Error in registering Scalar BC for ML_BSSN::B3!");
   
   }
   
@@ -1645,25 +1645,25 @@ void ML_BSSN_ApplyBoundConds(CCTK_ARGUMENTS)
 #$bound$#ML_BSSN::At33_bound_limit = 0.0
 #$bound$#ML_BSSN::At33_bound_scalar = 0.0
 
-#$bound$#ML_BSSN::dtalpha_bound       = "skip"
-#$bound$#ML_BSSN::dtalpha_bound_speed = 1.0
-#$bound$#ML_BSSN::dtalpha_bound_limit = 0.0
-#$bound$#ML_BSSN::dtalpha_bound_scalar = 0.0
+#$bound$#ML_BSSN::A_bound       = "skip"
+#$bound$#ML_BSSN::A_bound_speed = 1.0
+#$bound$#ML_BSSN::A_bound_limit = 0.0
+#$bound$#ML_BSSN::A_bound_scalar = 0.0
 
-#$bound$#ML_BSSN::dtbeta1_bound       = "skip"
-#$bound$#ML_BSSN::dtbeta1_bound_speed = 1.0
-#$bound$#ML_BSSN::dtbeta1_bound_limit = 0.0
-#$bound$#ML_BSSN::dtbeta1_bound_scalar = 0.0
+#$bound$#ML_BSSN::B1_bound       = "skip"
+#$bound$#ML_BSSN::B1_bound_speed = 1.0
+#$bound$#ML_BSSN::B1_bound_limit = 0.0
+#$bound$#ML_BSSN::B1_bound_scalar = 0.0
 
-#$bound$#ML_BSSN::dtbeta2_bound       = "skip"
-#$bound$#ML_BSSN::dtbeta2_bound_speed = 1.0
-#$bound$#ML_BSSN::dtbeta2_bound_limit = 0.0
-#$bound$#ML_BSSN::dtbeta2_bound_scalar = 0.0
+#$bound$#ML_BSSN::B2_bound       = "skip"
+#$bound$#ML_BSSN::B2_bound_speed = 1.0
+#$bound$#ML_BSSN::B2_bound_limit = 0.0
+#$bound$#ML_BSSN::B2_bound_scalar = 0.0
 
-#$bound$#ML_BSSN::dtbeta3_bound       = "skip"
-#$bound$#ML_BSSN::dtbeta3_bound_speed = 1.0
-#$bound$#ML_BSSN::dtbeta3_bound_limit = 0.0
-#$bound$#ML_BSSN::dtbeta3_bound_scalar = 0.0
+#$bound$#ML_BSSN::B3_bound       = "skip"
+#$bound$#ML_BSSN::B3_bound_speed = 1.0
+#$bound$#ML_BSSN::B3_bound_limit = 0.0
+#$bound$#ML_BSSN::B3_bound_scalar = 0.0
 
 #$bound$#ML_BSSN::Xt1_bound       = "skip"
 #$bound$#ML_BSSN::Xt1_bound_speed = 1.0
diff --git a/ML_BSSN/src/ML_BSSN_Minkowski.c b/ML_BSSN/src/ML_BSSN_Minkowski.c
index 37dd6ad..e110dda 100644
--- a/ML_BSSN/src/ML_BSSN_Minkowski.c
+++ b/ML_BSSN/src/ML_BSSN_Minkowski.c
@@ -96,11 +96,11 @@ void ML_BSSN_Minkowski_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     /* Declare shorthands */
     
     /* Declare local copies of grid functions */
+    CCTK_REAL AL = INITVALUE;
     CCTK_REAL alphaL = INITVALUE;
     CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE;
+    CCTK_REAL B1L = INITVALUE, B2L = INITVALUE, B3L = INITVALUE;
     CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE;
-    CCTK_REAL dtalphaL = INITVALUE;
-    CCTK_REAL dtbeta1L = INITVALUE, dtbeta2L = INITVALUE, dtbeta3L = INITVALUE;
     CCTK_REAL gt11L = INITVALUE, gt12L = INITVALUE, gt13L = INITVALUE, gt22L = INITVALUE, gt23L = INITVALUE, gt33L = INITVALUE;
     CCTK_REAL phiL = INITVALUE;
     CCTK_REAL trKL = INITVALUE;
@@ -156,7 +156,7 @@ void ML_BSSN_Minkowski_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     
     alphaL  =  1;
     
-    dtalphaL  =  0;
+    AL  =  0;
     
     beta1L  =  0;
     
@@ -164,14 +164,15 @@ void ML_BSSN_Minkowski_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     
     beta3L  =  0;
     
-    dtbeta1L  =  0;
+    B1L  =  0;
     
-    dtbeta2L  =  0;
+    B2L  =  0;
     
-    dtbeta3L  =  0;
+    B3L  =  0;
     
     
     /* Copy local copies back to grid functions */
+    A[index] = AL;
     alpha[index] = alphaL;
     At11[index] = At11L;
     At12[index] = At12L;
@@ -179,13 +180,12 @@ void ML_BSSN_Minkowski_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL
     At22[index] = At22L;
     At23[index] = At23L;
     At33[index] = At33L;
+    B1[index] = B1L;
+    B2[index] = B2L;
+    B3[index] = B3L;
     beta1[index] = beta1L;
     beta2[index] = beta2L;
     beta3[index] = beta3L;
-    dtalpha[index] = dtalphaL;
-    dtbeta1[index] = dtbeta1L;
-    dtbeta2[index] = dtbeta2L;
-    dtbeta3[index] = dtbeta3L;
     gt11[index] = gt11L;
     gt12[index] = gt12L;
     gt13[index] = gt13L;
diff --git a/ML_BSSN/src/ML_BSSN_RHS.c b/ML_BSSN/src/ML_BSSN_RHS.c
index e97f2b5..c34c396 100644
--- a/ML_BSSN/src/ML_BSSN_RHS.c
+++ b/ML_BSSN/src/ML_BSSN_RHS.c
@@ -132,12 +132,13 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
     CCTK_REAL Rt11 = INITVALUE, Rt12 = INITVALUE, Rt13 = INITVALUE, Rt22 = INITVALUE, Rt23 = INITVALUE, Rt33 = INITVALUE;
     
     /* Declare local copies of grid functions */
+    CCTK_REAL AL = INITVALUE;
     CCTK_REAL alphaL = INITVALUE, alpharhsL = INITVALUE;
+    CCTK_REAL ArhsL = INITVALUE;
     CCTK_REAL At11L = INITVALUE, At11rhsL = INITVALUE, At12L = INITVALUE, At12rhsL = INITVALUE, At13L = INITVALUE, At13rhsL = INITVALUE;
     CCTK_REAL At22L = INITVALUE, At22rhsL = INITVALUE, At23L = INITVALUE, At23rhsL = INITVALUE, At33L = INITVALUE, At33rhsL = INITVALUE;
+    CCTK_REAL B1L = INITVALUE, B1rhsL = INITVALUE, B2L = INITVALUE, B2rhsL = INITVALUE, B3L = INITVALUE, B3rhsL = INITVALUE;
     CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE;
-    CCTK_REAL dtalphaL = INITVALUE, dtalpharhsL = INITVALUE;
-    CCTK_REAL dtbeta1L = INITVALUE, dtbeta1rhsL = INITVALUE, dtbeta2L = INITVALUE, dtbeta2rhsL = INITVALUE, dtbeta3L = INITVALUE, dtbeta3rhsL = INITVALUE;
     CCTK_REAL gt11L = INITVALUE, gt11rhsL = INITVALUE, gt12L = INITVALUE, gt12rhsL = INITVALUE, gt13L = INITVALUE, gt13rhsL = INITVALUE;
     CCTK_REAL gt22L = INITVALUE, gt22rhsL = INITVALUE, gt23L = INITVALUE, gt23rhsL = INITVALUE, gt33L = INITVALUE, gt33rhsL = INITVALUE;
     CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE;
@@ -310,6 +311,7 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
     CCTK_REAL PDstandard4th3Xt3 = INITVALUE;
     
     /* Assign local copies of grid functions */
+    AL = A[index];
     alphaL = alpha[index];
     At11L = At11[index];
     At12L = At12[index];
@@ -317,13 +319,12 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
     At22L = At22[index];
     At23L = At23[index];
     At33L = At33[index];
+    B1L = B1[index];
+    B2L = B2[index];
+    B3L = B3[index];
     beta1L = beta1[index];
     beta2L = beta2[index];
     beta3L = beta3[index];
-    dtalphaL = dtalpha[index];
-    dtbeta1L = dtbeta1[index];
-    dtbeta2L = dtbeta2[index];
-    dtbeta3L = dtbeta3[index];
     gt11L = gt11[index];
     gt12L = gt12[index];
     gt13L = gt13[index];
@@ -1431,38 +1432,38 @@ void ML_BSSN_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal
         alphaL*(-2*(At13L*Atm13 + At23L*Atm23 + At33L*Atm33) + At33L*trKL);
     
     alpharhsL  =  beta1L*PDstandard4th1alpha + beta2L*PDstandard4th2alpha + beta3L*PDstandard4th3alpha - 
-        dtalphaL*harmonicF*pow(alphaL,harmonicN);
+        AL*harmonicF*pow(alphaL,harmonicN);
     
-    dtalpharhsL  =  -(AlphaDriver*dtalphaL) + trKrhsL;
+    ArhsL  =  -(AL*AlphaDriver) + trKrhsL;
     
-    beta1rhsL  =  dtbeta1L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
+    beta1rhsL  =  B1L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    beta2rhsL  =  dtbeta2L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
+    beta2rhsL  =  B2L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    beta3rhsL  =  dtbeta3L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
+    beta3rhsL  =  B3L*ShiftGammaCoeff*pow(alphaL,ShiftAlphaPower);
     
-    dtbeta1rhsL  =  -(BetaDriver*dtbeta1L) + Xt1rhsL;
+    B1rhsL  =  -(B1L*BetaDriver) + Xt1rhsL;
     
-    dtbeta2rhsL  =  -(BetaDriver*dtbeta2L) + Xt2rhsL;
+    B2rhsL  =  -(B2L*BetaDriver) + Xt2rhsL;
     
-    dtbeta3rhsL  =  -(BetaDriver*dtbeta3L) + Xt3rhsL;
+    B3rhsL  =  -(B3L*BetaDriver) + Xt3rhsL;
     
     
     /* Copy local copies back to grid functions */
     alpharhs[index] = alpharhsL;
+    Arhs[index] = ArhsL;
     At11rhs[index] = At11rhsL;
     At12rhs[index] = At12rhsL;
     At13rhs[index] = At13rhsL;
     At22rhs[index] = At22rhsL;
     At23rhs[index] = At23rhsL;
     At33rhs[index] = At33rhsL;
+    B1rhs[index] = B1rhsL;
+    B2rhs[index] = B2rhsL;
+    B3rhs[index] = B3rhsL;
     beta1rhs[index] = beta1rhsL;
     beta2rhs[index] = beta2rhsL;
     beta3rhs[index] = beta3rhsL;
-    dtalpharhs[index] = dtalpharhsL;
-    dtbeta1rhs[index] = dtbeta1rhsL;
-    dtbeta2rhs[index] = dtbeta2rhsL;
-    dtbeta3rhs[index] = dtbeta3rhsL;
     gt11rhs[index] = gt11rhsL;
     gt12rhs[index] = gt12rhsL;
     gt13rhs[index] = gt13rhsL;
diff --git a/ML_BSSN/src/ML_BSSN_constraints.c b/ML_BSSN/src/ML_BSSN_constraints.c
index f0f575c..4d6a26d 100644
--- a/ML_BSSN/src/ML_BSSN_constraints.c
+++ b/ML_BSSN/src/ML_BSSN_constraints.c
@@ -159,7 +159,6 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     CCTK_REAL PDstandard4th1gt11 = INITVALUE;
     CCTK_REAL PDstandard4th2gt11 = INITVALUE;
     CCTK_REAL PDstandard4th3gt11 = INITVALUE;
-    CCTK_REAL PDstandard4th11gt11 = INITVALUE;
     CCTK_REAL PDstandard4th22gt11 = INITVALUE;
     CCTK_REAL PDstandard4th33gt11 = INITVALUE;
     CCTK_REAL PDstandard4th12gt11 = INITVALUE;
@@ -171,8 +170,6 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     CCTK_REAL PDstandard4th1gt12 = INITVALUE;
     CCTK_REAL PDstandard4th2gt12 = INITVALUE;
     CCTK_REAL PDstandard4th3gt12 = INITVALUE;
-    CCTK_REAL PDstandard4th11gt12 = INITVALUE;
-    CCTK_REAL PDstandard4th22gt12 = INITVALUE;
     CCTK_REAL PDstandard4th33gt12 = INITVALUE;
     CCTK_REAL PDstandard4th12gt12 = INITVALUE;
     CCTK_REAL PDstandard4th13gt12 = INITVALUE;
@@ -183,9 +180,7 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     CCTK_REAL PDstandard4th1gt13 = INITVALUE;
     CCTK_REAL PDstandard4th2gt13 = INITVALUE;
     CCTK_REAL PDstandard4th3gt13 = INITVALUE;
-    CCTK_REAL PDstandard4th11gt13 = INITVALUE;
     CCTK_REAL PDstandard4th22gt13 = INITVALUE;
-    CCTK_REAL PDstandard4th33gt13 = INITVALUE;
     CCTK_REAL PDstandard4th12gt13 = INITVALUE;
     CCTK_REAL PDstandard4th13gt13 = INITVALUE;
     CCTK_REAL PDstandard4th21gt13 = INITVALUE;
@@ -196,7 +191,6 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     CCTK_REAL PDstandard4th2gt22 = INITVALUE;
     CCTK_REAL PDstandard4th3gt22 = INITVALUE;
     CCTK_REAL PDstandard4th11gt22 = INITVALUE;
-    CCTK_REAL PDstandard4th22gt22 = INITVALUE;
     CCTK_REAL PDstandard4th33gt22 = INITVALUE;
     CCTK_REAL PDstandard4th12gt22 = INITVALUE;
     CCTK_REAL PDstandard4th13gt22 = INITVALUE;
@@ -208,8 +202,6 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     CCTK_REAL PDstandard4th2gt23 = INITVALUE;
     CCTK_REAL PDstandard4th3gt23 = INITVALUE;
     CCTK_REAL PDstandard4th11gt23 = INITVALUE;
-    CCTK_REAL PDstandard4th22gt23 = INITVALUE;
-    CCTK_REAL PDstandard4th33gt23 = INITVALUE;
     CCTK_REAL PDstandard4th12gt23 = INITVALUE;
     CCTK_REAL PDstandard4th13gt23 = INITVALUE;
     CCTK_REAL PDstandard4th21gt23 = INITVALUE;
@@ -221,7 +213,6 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     CCTK_REAL PDstandard4th3gt33 = INITVALUE;
     CCTK_REAL PDstandard4th11gt33 = INITVALUE;
     CCTK_REAL PDstandard4th22gt33 = INITVALUE;
-    CCTK_REAL PDstandard4th33gt33 = INITVALUE;
     CCTK_REAL PDstandard4th12gt33 = INITVALUE;
     CCTK_REAL PDstandard4th13gt33 = INITVALUE;
     CCTK_REAL PDstandard4th21gt33 = INITVALUE;
@@ -243,15 +234,6 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     CCTK_REAL PDstandard4th1trK = INITVALUE;
     CCTK_REAL PDstandard4th2trK = INITVALUE;
     CCTK_REAL PDstandard4th3trK = INITVALUE;
-    CCTK_REAL PDstandard4th1Xt1 = INITVALUE;
-    CCTK_REAL PDstandard4th2Xt1 = INITVALUE;
-    CCTK_REAL PDstandard4th3Xt1 = INITVALUE;
-    CCTK_REAL PDstandard4th1Xt2 = INITVALUE;
-    CCTK_REAL PDstandard4th2Xt2 = INITVALUE;
-    CCTK_REAL PDstandard4th3Xt2 = INITVALUE;
-    CCTK_REAL PDstandard4th1Xt3 = INITVALUE;
-    CCTK_REAL PDstandard4th2Xt3 = INITVALUE;
-    CCTK_REAL PDstandard4th3Xt3 = INITVALUE;
     
     /* Assign local copies of grid functions */
     At11L = At11[index];
@@ -295,17 +277,12 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     PDstandard4th1gt11 = PDstandard4th1(gt11, i, j, k);
     PDstandard4th2gt11 = PDstandard4th2(gt11, i, j, k);
     PDstandard4th3gt11 = PDstandard4th3(gt11, i, j, k);
-    PDstandard4th11gt11 = PDstandard4th11(gt11, i, j, k);
     PDstandard4th22gt11 = PDstandard4th22(gt11, i, j, k);
     PDstandard4th33gt11 = PDstandard4th33(gt11, i, j, k);
-    PDstandard4th12gt11 = PDstandard4th12(gt11, i, j, k);
-    PDstandard4th13gt11 = PDstandard4th13(gt11, i, j, k);
     PDstandard4th23gt11 = PDstandard4th23(gt11, i, j, k);
     PDstandard4th1gt12 = PDstandard4th1(gt12, i, j, k);
     PDstandard4th2gt12 = PDstandard4th2(gt12, i, j, k);
     PDstandard4th3gt12 = PDstandard4th3(gt12, i, j, k);
-    PDstandard4th11gt12 = PDstandard4th11(gt12, i, j, k);
-    PDstandard4th22gt12 = PDstandard4th22(gt12, i, j, k);
     PDstandard4th33gt12 = PDstandard4th33(gt12, i, j, k);
     PDstandard4th12gt12 = PDstandard4th12(gt12, i, j, k);
     PDstandard4th13gt12 = PDstandard4th13(gt12, i, j, k);
@@ -313,9 +290,7 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     PDstandard4th1gt13 = PDstandard4th1(gt13, i, j, k);
     PDstandard4th2gt13 = PDstandard4th2(gt13, i, j, k);
     PDstandard4th3gt13 = PDstandard4th3(gt13, i, j, k);
-    PDstandard4th11gt13 = PDstandard4th11(gt13, i, j, k);
     PDstandard4th22gt13 = PDstandard4th22(gt13, i, j, k);
-    PDstandard4th33gt13 = PDstandard4th33(gt13, i, j, k);
     PDstandard4th12gt13 = PDstandard4th12(gt13, i, j, k);
     PDstandard4th13gt13 = PDstandard4th13(gt13, i, j, k);
     PDstandard4th23gt13 = PDstandard4th23(gt13, i, j, k);
@@ -323,17 +298,12 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     PDstandard4th2gt22 = PDstandard4th2(gt22, i, j, k);
     PDstandard4th3gt22 = PDstandard4th3(gt22, i, j, k);
     PDstandard4th11gt22 = PDstandard4th11(gt22, i, j, k);
-    PDstandard4th22gt22 = PDstandard4th22(gt22, i, j, k);
     PDstandard4th33gt22 = PDstandard4th33(gt22, i, j, k);
-    PDstandard4th12gt22 = PDstandard4th12(gt22, i, j, k);
     PDstandard4th13gt22 = PDstandard4th13(gt22, i, j, k);
-    PDstandard4th23gt22 = PDstandard4th23(gt22, i, j, k);
     PDstandard4th1gt23 = PDstandard4th1(gt23, i, j, k);
     PDstandard4th2gt23 = PDstandard4th2(gt23, i, j, k);
     PDstandard4th3gt23 = PDstandard4th3(gt23, i, j, k);
     PDstandard4th11gt23 = PDstandard4th11(gt23, i, j, k);
-    PDstandard4th22gt23 = PDstandard4th22(gt23, i, j, k);
-    PDstandard4th33gt23 = PDstandard4th33(gt23, i, j, k);
     PDstandard4th12gt23 = PDstandard4th12(gt23, i, j, k);
     PDstandard4th13gt23 = PDstandard4th13(gt23, i, j, k);
     PDstandard4th23gt23 = PDstandard4th23(gt23, i, j, k);
@@ -342,10 +312,7 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     PDstandard4th3gt33 = PDstandard4th3(gt33, i, j, k);
     PDstandard4th11gt33 = PDstandard4th11(gt33, i, j, k);
     PDstandard4th22gt33 = PDstandard4th22(gt33, i, j, k);
-    PDstandard4th33gt33 = PDstandard4th33(gt33, i, j, k);
     PDstandard4th12gt33 = PDstandard4th12(gt33, i, j, k);
-    PDstandard4th13gt33 = PDstandard4th13(gt33, i, j, k);
-    PDstandard4th23gt33 = PDstandard4th23(gt33, i, j, k);
     PDstandard4th1phi = PDstandard4th1(phi, i, j, k);
     PDstandard4th2phi = PDstandard4th2(phi, i, j, k);
     PDstandard4th3phi = PDstandard4th3(phi, i, j, k);
@@ -358,15 +325,6 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     PDstandard4th1trK = PDstandard4th1(trK, i, j, k);
     PDstandard4th2trK = PDstandard4th2(trK, i, j, k);
     PDstandard4th3trK = PDstandard4th3(trK, i, j, k);
-    PDstandard4th1Xt1 = PDstandard4th1(Xt1, i, j, k);
-    PDstandard4th2Xt1 = PDstandard4th2(Xt1, i, j, k);
-    PDstandard4th3Xt1 = PDstandard4th3(Xt1, i, j, k);
-    PDstandard4th1Xt2 = PDstandard4th1(Xt2, i, j, k);
-    PDstandard4th2Xt2 = PDstandard4th2(Xt2, i, j, k);
-    PDstandard4th3Xt2 = PDstandard4th3(Xt2, i, j, k);
-    PDstandard4th1Xt3 = PDstandard4th1(Xt3, i, j, k);
-    PDstandard4th2Xt3 = PDstandard4th2(Xt3, i, j, k);
-    PDstandard4th3Xt3 = PDstandard4th3(Xt3, i, j, k);
     
     /* Precompute derivatives (old style) */
     
@@ -445,250 +403,45 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     Gt333  =  khalf*(-(gtu31*PDstandard4th1gt33) - gtu32*PDstandard4th2gt33 + 2*gtu31*PDstandard4th3gt13 + 
           2*gtu32*PDstandard4th3gt23 + gtu33*PDstandard4th3gt33);
     
-    Rt11  =  -(gtu11*khalf*PDstandard4th11gt11) + gtu21*
-         (2*Gt211*Gt212*gt22L + 6*Gt112*gt13L*Gt311 + 4*Gt113*gt11L*Gt312 + 4*gt13L*Gt312*Gt313 + 4*gt13L*Gt211*Gt322 + 
-           4*gt13L*Gt311*Gt323 + 2*Gt311*Gt312*gt33L - PDstandard4th12gt11) - gtu31*PDstandard4th13gt11 + 
-        2*gt11L*PDstandard4th1Xt1 + gt12L*(6*Gt111*Gt212*gtu21 + 4*Gt211*Gt222*gtu21 + 4*Gt212*Gt222*gtu22 + 
-           6*Gt113*Gt211*gtu31 + 6*Gt113*Gt212*gtu32 + 6*Gt113*Gt213*gtu33 + 2*PDstandard4th1Xt2) + 
-        gt13L*(6*Gt111*Gt312*gtu21 + 4*Gt212*Gt312*gtu21 + 6*Gt112*Gt312*gtu22 + 6*Gt113*Gt311*gtu31 + 
-           6*Gt113*Gt312*gtu32 + 6*Gt113*Gt313*gtu33 + 2*PDstandard4th1Xt3) - gtu22*khalf*PDstandard4th22gt11 + 
-        gtu32*(2*Gt212*Gt213*gt22L + 4*gt12L*Gt233*Gt312 + 2*Gt213*gt23L*Gt312 + 6*Gt112*gt13L*Gt313 + 
-           4*gt13L*Gt213*Gt322 + 4*gt13L*Gt212*Gt323 + 4*gt13L*Gt313*Gt323 + 4*gt13L*Gt312*Gt333 + 2*Gt312*Gt313*gt33L - 
-           PDstandard4th23gt11) - gtu33*khalf*PDstandard4th33gt11 + 
-        Gt111*(10*Gt113*gt11L*gtu31 + 6*gt12L*Gt213*gtu31 + 2*gt11L*Xt1L) + 
-        Gt211*(4*Gt112*gt11L*gtu11 + 6*Gt111*gt12L*gtu11 + 2*gt23L*Gt311*gtu11 + 4*gt13L*Gt312*gtu11 + 
-           6*Gt112*gt12L*gtu21 + 4*gt11L*Gt123*gtu31 + 4*gt12L*Gt223*gtu31 + 2*Gt213*gt22L*gtu31 + 2*gt12L*Xt1L) + 
-        Gt311*(6*Gt111*gt13L*gtu11 + 4*gt12L*Gt213*gtu11 + 4*gt13L*Gt313*gtu11 + 4*gt11L*Gt123*gtu21 + 
-           4*gt12L*Gt223*gtu21 + 2*Gt212*gt23L*gtu21 + 4*gt11L*Gt133*gtu31 + 4*gt12L*Gt233*gtu31 + 2*Gt213*gt23L*gtu31 + 
-           2*Gt313*gt33L*gtu31 + 2*gt13L*Xt1L) + 2*gt12L*Gt212*Xt2L + 
-        Gt112*(10*Gt111*gt11L*gtu21 + 4*gt11L*Gt212*gtu21 + 6*gt12L*Gt212*gtu22 + 4*gt11L*Gt213*gtu31 + 
-           10*Gt113*gt11L*gtu32 + 2*gt11L*Xt2L) + Gt312*
-         (4*gt12L*Gt213*gtu21 + 2*Gt211*gt23L*gtu21 + 4*gt11L*Gt123*gtu22 + 4*gt12L*Gt223*gtu22 + 2*Gt212*gt23L*gtu22 + 
-           4*gt13L*Gt213*gtu31 + 4*gt11L*Gt133*gtu32 + 2*gt13L*Xt2L) + 2*Gt113*gt11L*Xt3L + 
-        Gt213*(4*gt11L*Gt122*gtu32 + 6*Gt112*gt12L*gtu32 + 4*gt11L*Gt123*gtu33 + 2*gt12L*Xt3L) + 
-        Gt313*(4*Gt113*gt11L*gtu31 + 6*Gt111*gt13L*gtu31 + 4*gt12L*Gt213*gtu31 + 2*Gt211*gt23L*gtu31 + 
-           4*gt11L*Gt123*gtu32 + 4*gt12L*Gt223*gtu32 + 2*Gt212*gt23L*gtu32 + 4*gt11L*Gt133*gtu33 + 4*gt12L*Gt233*gtu33 + 
-           2*Gt213*gt23L*gtu33 + 2*gt13L*Xt3L) + 5*gt11L*gtu11*SQR(Gt111) + 5*gt11L*gtu22*SQR(Gt112) + 
-        5*gt11L*gtu33*SQR(Gt113) + gt22L*gtu11*SQR(Gt211) + gt22L*gtu22*SQR(Gt212) + 
-        4*(gt12L*Gt211*Gt212*gtu11 + Gt113*gt11L*Gt311*gtu11 + gt11L*Gt122*Gt211*gtu21 + gt11L*Gt122*Gt212*gtu22 + 
-           gt13L*Gt212*Gt322*gtu22 + gt13L*Gt312*Gt323*gtu22 + gt12L*Gt212*Gt213*gtu31 + gt13L*Gt211*Gt323*gtu31 + 
-           gt13L*Gt311*Gt333*gtu31 + gt11L*Gt123*Gt212*gtu32 + gt12L*Gt213*Gt222*gtu32 + gt12L*Gt212*Gt223*gtu32 + 
-           gt12L*Gt213*Gt223*gtu33 + gt13L*Gt213*Gt323*gtu33 + gt13L*Gt313*Gt333*gtu33 + gt12L*gtu21*SQR(Gt212)) + 
-        gt22L*gtu33*SQR(Gt213) + gt33L*gtu11*SQR(Gt311) + gt33L*gtu22*SQR(Gt312) + 4*gt13L*gtu31*SQR(Gt313) + 
-        gt33L*gtu33*SQR(Gt313);
-    
-    Rt12  =  gt22L*PDstandard4th1Xt2 + gt23L*PDstandard4th1Xt3 + gt11L*PDstandard4th2Xt1 + 
-        gt12L*(PDstandard4th1Xt1 + PDstandard4th2Xt2) + gt13L*PDstandard4th2Xt3 + 
-        khalf*(-(gtu11*PDstandard4th11gt12) - 2*gtu21*PDstandard4th12gt12 - 2*gtu31*PDstandard4th13gt12 - 
-           gtu22*PDstandard4th22gt12 - 2*gtu32*PDstandard4th23gt12 - gtu33*PDstandard4th33gt12) + 
-        (Gt112*gt11L + gt12L*(Gt111 + Gt212) + Gt211*gt22L + gt23L*Gt311 + gt13L*Gt312)*Xt1L + 
-        (gt11L*Gt122 + gt12L*(Gt112 + Gt222) + Gt212*gt22L + gt23L*Gt312 + gt13L*Gt322)*Xt2L + 
-        (gt11L*Gt123 + gt12L*(Gt113 + Gt223) + Gt213*gt22L + gt23L*Gt313 + gt13L*Gt323)*Xt3L + 
-        gtu11*(gt11L*(Gt112*(3*Gt111 + 2*Gt212) + 2*Gt113*Gt312) + 
-           Gt111*(gt12L*Gt212 + 2*(Gt211*gt22L + gt23L*Gt311) + gt13L*Gt312) + 
-           Gt211*(5*Gt112*gt12L + 3*(Gt212*gt22L + gt23L*Gt312)) + 
-           Gt311*(3*Gt112*gt13L + Gt212*gt23L + 2*(Gt113*gt12L + Gt213*gt22L + gt23L*Gt313) + Gt312*gt33L) + 
-           2*(Gt312*(gt12L*Gt213 + gt13L*(Gt212 + Gt313)) + gt12L*(SQR(Gt111) + SQR(Gt212)))) + 
-        gtu22*(Gt112*(3*gt11L*Gt122 + gt12L*Gt222 + 2*(Gt212*gt22L + gt23L*Gt312) + gt13L*Gt322) + 
-           Gt222*(3*Gt212*gt22L + gt23L*Gt312 + 2*(gt11L*Gt122 + gt13L*Gt322)) + 
-           Gt312*(3*Gt122*gt13L + 2*(Gt223*gt22L + gt23L*Gt323)) + 
-           Gt322*(3*Gt212*gt23L + 2*(gt11L*Gt123 + gt13L*Gt323) + Gt312*gt33L) + 
-           gt12L*(5*Gt122*Gt212 + 2*(Gt123*Gt312 + Gt223*Gt322 + SQR(Gt112) + SQR(Gt222)))) + 
-        gtu33*(Gt113*(3*gt11L*Gt123 + gt12L*Gt223 + 2*(Gt213*gt22L + gt23L*Gt313) + gt13L*Gt323) + 
-           Gt223*(3*Gt213*gt22L + gt23L*Gt313 + 2*(gt11L*Gt123 + gt13L*Gt323)) + 
-           Gt313*(3*Gt123*gt13L + 2*(gt22L*Gt233 + gt23L*Gt333)) + 
-           Gt323*(3*Gt213*gt23L + 2*(gt11L*Gt133 + gt13L*Gt333) + Gt313*gt33L) + 
-           gt12L*(5*Gt123*Gt213 + 2*(Gt133*Gt313 + Gt233*Gt323 + SQR(Gt113) + SQR(Gt223)))) + 
-        gtu21*(Gt111*(3*gt11L*Gt122 + gt12L*(2*Gt112 + Gt222) + gt13L*Gt322) + 
-           2*(Gt112*(gt11L*Gt222 + Gt211*gt22L) + Gt223*gt22L*Gt311 + Gt111*(Gt112*gt12L + Gt212*gt22L + gt23L*Gt312) + 
-              (Gt113*gt11L + gt13L*Gt313)*Gt322 + gt13L*(Gt222*Gt312 + Gt212*Gt322)) + 
-           gt12L*(5*Gt122*Gt211 + Gt212*(4*Gt112 + 2*Gt222) + 
-              2*(Gt212*(Gt112 + Gt222) + Gt123*Gt311 + (Gt113 + Gt223)*Gt312 + Gt213*Gt322)) + 
-           Gt312*(2*(Gt213*gt22L + gt23L*Gt313) + gt13L*(4*Gt112 + 2*Gt323)) + 
-           gt23L*(4*Gt212*Gt312 + Gt311*(Gt222 + 2*(Gt112 + Gt323))) + 
-           gt11L*(2*(Gt122*Gt212 + Gt123*Gt312) + 3*SQR(Gt112)) + 
-           3*(Gt122*gt13L*Gt311 + Gt211*(Gt222*gt22L + gt23L*Gt322) + gt22L*SQR(Gt212)) + gt33L*(Gt311*Gt322 + SQR(Gt312)))\
-         + gtu31*(Gt111*(3*gt11L*Gt123 + gt12L*(2*Gt113 + Gt223) + gt13L*Gt323) + 
-           gt11L*(2*(Gt123*Gt212 + Gt133*Gt312) + Gt113*(3*Gt112 + 2*Gt323)) + 
-           3*((Gt212*Gt213 + Gt211*Gt223)*gt22L + Gt123*gt13L*Gt311 + gt23L*(Gt213*Gt312 + Gt211*Gt323)) + 
-           gt12L*(5*Gt123*Gt211 + 3*Gt112*Gt213 + Gt212*(Gt113 + 2*Gt223) + 
-              2*(Gt212*Gt223 + Gt133*Gt311 + Gt233*Gt312 + Gt113*Gt313 + Gt213*(Gt112 + Gt323))) + 
-           gt13L*(3*Gt112*Gt313 + 2*Gt212*Gt323 + Gt312*(Gt113 + 2*Gt333)) + 
-           gt23L*(Gt212*Gt313 + Gt311*(Gt223 + 2*(Gt113 + Gt333))) + (Gt312*Gt313 + Gt311*Gt323)*gt33L + 
-           2*(Gt223*(Gt112*gt11L + gt13L*Gt312) + gt22L*(Gt113*Gt211 + Gt233*Gt311 + Gt213*Gt313) + 
-              Gt111*(Gt113*gt12L + Gt213*gt22L + gt23L*Gt313) + gt13L*Gt313*Gt323 + gt23L*SQR(Gt313))) + 
-        gtu32*(3*((Gt213*Gt222 + Gt212*Gt223)*gt22L + gt13L*(Gt123*Gt312 + Gt122*Gt313) + Gt213*gt23L*Gt322) + 
-           Gt112*(3*gt11L*Gt123 + gt12L*(2*Gt113 + Gt223) + gt13L*Gt323) + 
-           gt12L*(5*(Gt123*Gt212 + Gt122*Gt213) + Gt222*(Gt113 + 4*Gt223) + 
-              2*(Gt133*Gt312 + Gt123*Gt313 + Gt233*Gt322 + Gt223*Gt323)) + 
-           gt11L*(3*Gt113*Gt122 + 2*(Gt133*Gt322 + Gt123*(Gt222 + Gt323))) + 
-           gt13L*(2*Gt222*Gt323 + Gt322*(Gt113 + 2*Gt333)) + 
-           gt23L*(Gt222*Gt313 + 3*Gt212*Gt323 + Gt312*(Gt223 + 2*(Gt113 + Gt333))) + (Gt313*Gt322 + Gt312*Gt323)*gt33L + 
-           2*(gt22L*(Gt113*Gt212 + Gt233*Gt312) + Gt112*(Gt113*gt12L + Gt213*gt22L + gt23L*Gt313) + 
-              Gt223*(gt11L*Gt122 + gt22L*Gt313 + gt13L*Gt322) + gt23L*Gt313*Gt323 + gt13L*SQR(Gt323)));
-    
-    Rt13  =  gt23L*PDstandard4th1Xt2 + gt33L*PDstandard4th1Xt3 + 
-        khalf*(-(gtu11*PDstandard4th11gt13) - 2*gtu21*PDstandard4th12gt13 - 2*gtu31*PDstandard4th13gt13 - 
-           gtu22*PDstandard4th22gt13 - 2*gtu32*PDstandard4th23gt13 - gtu33*PDstandard4th33gt13) + gt11L*PDstandard4th3Xt1 + 
-        gt12L*PDstandard4th3Xt2 + gt13L*(PDstandard4th1Xt1 + PDstandard4th3Xt3) + 
-        (Gt113*gt11L + gt12L*Gt213 + Gt211*gt23L + gt13L*(Gt111 + Gt313) + Gt311*gt33L)*Xt1L + 
-        (gt11L*Gt123 + gt12L*Gt223 + Gt212*gt23L + gt13L*(Gt112 + Gt323) + Gt312*gt33L)*Xt2L + 
-        (gt11L*Gt133 + gt12L*Gt233 + Gt213*gt23L + gt13L*(Gt113 + Gt333) + Gt313*gt33L)*Xt3L + 
-        gtu21*((Gt212*Gt213 + Gt211*Gt223)*gt22L + Gt123*(Gt111*gt11L + gt12L*Gt211 + 3*gt13L*Gt311) + 
-           gt12L*(3*Gt113*Gt212 + Gt112*Gt213 + Gt223*(Gt111 + 2*Gt313)) + 
-           gt11L*(3*Gt112*Gt113 + 2*(Gt122*Gt213 + Gt123*Gt313)) + gt13L*(Gt112*Gt313 + Gt111*Gt323) + 
-           gt23L*(Gt213*Gt312 + Gt212*Gt313 + Gt211*Gt323) + Gt312*Gt313*gt33L + 
-           3*(Gt113*gt13L*Gt312 + Gt311*(Gt223*gt23L + Gt323*gt33L)) + 
-           2*(Gt213*(gt12L*Gt222 + gt13L*Gt322) + gt11L*(Gt111*Gt123 + Gt112*Gt223 + Gt113*Gt323) + 
-              gt12L*(Gt123*Gt211 + Gt212*Gt223 + Gt213*Gt323) + 
-              gt13L*(Gt112*(Gt111 + Gt212) + Gt123*Gt311 + (Gt113 + Gt223)*Gt312 + 2*Gt313*Gt323) + 
-              (Gt112*Gt311 + Gt312*(Gt212 + Gt313))*gt33L + Gt111*(Gt112*gt13L + Gt212*gt23L + Gt312*gt33L) + 
-              Gt211*(Gt122*gt13L + Gt222*gt23L + Gt322*gt33L) + gt23L*(Gt112*Gt211 + Gt213*Gt312 + SQR(Gt212)))) + 
-        gtu32*(gt11L*(Gt123*(3*Gt113 + 2*Gt223) + 2*Gt133*Gt323) + Gt223*(gt23L*(2*Gt212 + 3*Gt313) + 2*gt13L*Gt323) + 
-           Gt233*(Gt212*gt22L + 3*gt23L*Gt312 + 2*(gt11L*Gt122 + gt12L*Gt323)) + 
-           Gt123*(3*gt12L*Gt213 + gt13L*(2*Gt212 + 5*Gt313) + 2*gt11L*Gt333) + 
-           gt13L*(5*Gt133*Gt312 + 2*Gt233*Gt322 + 4*Gt323*Gt333) + 3*(Gt313*Gt323 + Gt312*Gt333)*gt33L + 
-           Gt212*(gt23L*Gt333 + 2*Gt323*gt33L) + Gt113*(gt12L*Gt223 + gt13L*Gt323 + 2*(Gt212*gt23L + Gt312*gt33L)) + 
-           Gt112*(3*gt11L*Gt133 + gt12L*Gt233 + gt13L*(2*Gt113 + Gt333) + 2*(Gt113*gt13L + Gt213*gt23L + Gt313*gt33L)) + 
-           Gt213*(Gt223*gt22L + gt23L*Gt323 + 2*(Gt122*gt13L + Gt222*gt23L + Gt322*gt33L)) + 
-           gt12L*(3*Gt133*Gt212 + 2*(Gt222*Gt233 + Gt223*Gt333 + SQR(Gt223)))) + 
-        gtu31*(Gt133*(Gt111*gt11L + gt12L*Gt211 + 3*gt13L*Gt311) + Gt233*(Gt211*gt22L + 3*gt23L*Gt311) + 
-           gt12L*(4*Gt113*Gt213 + Gt233*(Gt111 + 2*Gt313)) + gt13L*(4*Gt113*Gt313 + 2*Gt213*Gt323 + Gt111*Gt333) + 
-           Gt211*(gt23L*Gt333 + 2*Gt323*gt33L) + gt11L*(2*(Gt123*Gt213 + Gt133*Gt313) + 3*SQR(Gt113)) + gt22L*SQR(Gt213) + 
-           gt33L*(3*Gt311*Gt333 + SQR(Gt313)) + 2*(Gt211*(Gt123*gt13L + Gt223*gt23L) + Gt213*(gt12L*Gt223 + gt23L*Gt313) + 
-              gt23L*(Gt113*Gt211 + Gt213*(Gt212 + Gt313)) + gt11L*(Gt111*Gt133 + Gt112*Gt233 + Gt113*Gt333) + 
-              gt12L*(Gt133*Gt211 + Gt212*Gt233 + Gt213*Gt333) + 
-              gt13L*(Gt112*Gt213 + Gt133*Gt311 + Gt233*Gt312 + Gt113*(Gt111 + Gt313) + 2*Gt313*Gt333) + 
-              Gt111*(Gt113*gt13L + Gt213*gt23L + Gt313*gt33L) + gt33L*(Gt113*Gt311 + Gt213*Gt312 + SQR(Gt313)))) + 
-        gtu11*(Gt211*(3*Gt113*gt12L + 2*Gt112*gt13L + Gt213*gt22L + gt23L*Gt313) + 
-           gt11L*(2*Gt112*Gt213 + Gt113*(3*Gt111 + 2*Gt313)) + 
-           Gt111*(gt12L*Gt213 + gt13L*Gt313 + 2*(Gt211*gt23L + Gt311*gt33L)) + 
-           Gt311*(5*Gt113*gt13L + 3*(Gt213*gt23L + Gt313*gt33L)) + 
-           2*(Gt212*(gt12L*Gt213 + Gt211*gt23L) + Gt213*(gt13L*Gt312 + gt12L*Gt313) + Gt211*Gt312*gt33L + 
-              gt13L*(SQR(Gt111) + SQR(Gt313)))) + gtu22*
-         (Gt212*(3*Gt123*gt12L + 2*Gt122*gt13L + Gt223*gt22L + gt23L*Gt323) + 
-           gt11L*(2*Gt122*Gt223 + Gt123*(3*Gt112 + 2*Gt323)) + 
-           Gt112*(gt12L*Gt223 + gt13L*Gt323 + 2*(Gt212*gt23L + Gt312*gt33L)) + 
-           Gt312*(5*Gt123*gt13L + 3*(Gt223*gt23L + Gt323*gt33L)) + 
-           2*(Gt222*(gt12L*Gt223 + Gt212*gt23L) + Gt223*(gt13L*Gt322 + gt12L*Gt323) + Gt212*Gt322*gt33L + 
-              gt13L*(SQR(Gt112) + SQR(Gt323)))) + gtu33*
-         (gt12L*(3*Gt133*Gt213 + 2*Gt233*(Gt223 + Gt333)) + Gt213*(gt22L*Gt233 + gt23L*(2*Gt223 + Gt333)) + 
-           Gt113*(3*gt11L*Gt133 + gt12L*Gt233 + gt13L*Gt333 + 2*(Gt213*gt23L + Gt313*gt33L)) + 
-           Gt313*(5*Gt133*gt13L + 3*(Gt233*gt23L + Gt333*gt33L)) + 
-           2*(Gt123*(gt13L*Gt213 + gt11L*Gt233) + gt11L*Gt133*Gt333 + Gt323*(gt13L*Gt233 + Gt213*gt33L) + 
-              gt13L*(SQR(Gt113) + SQR(Gt333))));
-    
-    Rt22  =  6*(Gt122*gt12L*Gt212*gtu21 + Gt112*gt12L*Gt222*gtu21 + Gt122*gt12L*Gt222*gtu22 + Gt123*gt12L*Gt212*gtu31 + 
-           Gt123*gt12L*Gt222*gtu32 + Gt123*gt12L*Gt223*gtu33) - gtu11*khalf*PDstandard4th11gt22 + 
-        gtu21*(10*Gt212*Gt222*gt22L + 4*Gt122*gt23L*Gt311 + 2*Gt122*gt13L*Gt312 + 6*Gt222*gt23L*Gt312 + 
-           4*Gt113*gt12L*Gt322 + 4*gt23L*Gt312*Gt323 + 2*Gt312*Gt322*gt33L - PDstandard4th12gt22) + 
-        gtu31*(10*Gt212*Gt223*gt22L + 2*Gt123*gt13L*Gt312 + 4*Gt112*gt23L*Gt313 + 4*Gt113*gt12L*Gt323 + 
-           4*gt23L*Gt312*Gt333 + 2*Gt312*Gt323*gt33L - PDstandard4th13gt22) - gtu22*khalf*PDstandard4th22gt22 + 
-        gtu32*(6*Gt223*gt23L*Gt322 + 2*Gt122*gt13L*Gt323 + 4*gt23L*Gt322*Gt333 + 2*Gt322*Gt323*gt33L - 
-           PDstandard4th23gt22) + gt12L*(4*Gt111*Gt123*gtu31 + 6*Gt112*Gt223*gtu31 + 4*Gt113*Gt122*gtu32 + 
-           4*Gt113*Gt123*gtu33 + 2*PDstandard4th2Xt1) + 
-        gt22L*(4*Gt122*Gt213*gtu32 + 10*Gt222*Gt223*gtu32 + 4*Gt123*Gt213*gtu33 + 2*PDstandard4th2Xt2) + 
-        gt23L*(6*Gt212*Gt322*gtu21 + 4*Gt313*Gt322*gtu21 + 6*Gt222*Gt322*gtu22 + 4*Gt123*Gt311*gtu31 + 
-           6*Gt212*Gt323*gtu31 + 4*Gt313*Gt323*gtu31 + 4*Gt122*Gt313*gtu32 + 6*Gt223*Gt323*gtu33 + 2*PDstandard4th2Xt3) - 
-        gtu33*khalf*PDstandard4th33gt22 + 2*Gt212*gt22L*Xt1L + 
-        Gt112*(4*Gt111*gt12L*gtu11 + 6*gt12L*Gt212*gtu11 + 2*gt11L*Gt122*gtu21 + 4*Gt122*gt12L*gtu22 + 
-           2*gt11L*Gt123*gtu31 + 4*Gt123*gt12L*gtu32 + 2*gt12L*Xt1L) + 
-        Gt312*(4*Gt113*gt12L*gtu11 + 2*Gt112*gt13L*gtu11 + 4*Gt213*gt22L*gtu11 + 6*Gt212*gt23L*gtu11 + 
-           4*gt23L*Gt313*gtu11 + 4*Gt123*gt12L*gtu21 + 4*Gt122*gt23L*gtu22 + 4*gt12L*Gt133*gtu31 + 4*gt22L*Gt233*gtu31 + 
-           6*Gt223*gt23L*gtu31 + 4*Gt123*gt23L*gtu32 + 2*gt23L*Xt1L) + 2*Gt222*gt22L*Xt2L + 
-        Gt122*(2*gt11L*Gt123*gtu32 + 6*gt12L*Gt223*gtu32 + 2*gt12L*Xt2L) + 
-        Gt322*(2*Gt112*gt13L*gtu21 + 4*Gt213*gt22L*gtu21 + 4*Gt123*gt12L*gtu22 + 2*Gt122*gt13L*gtu22 + 
-           4*Gt223*gt22L*gtu22 + 4*gt23L*Gt323*gtu22 + 4*gt12L*Gt133*gtu32 + 2*Gt123*gt13L*gtu32 + 4*gt22L*Gt233*gtu32 + 
-           2*gt23L*Xt2L) + 2*Gt123*gt12L*Xt3L + 2*Gt223*gt22L*Xt3L + 
-        Gt323*(2*Gt112*gt13L*gtu31 + 4*Gt213*gt22L*gtu31 + 4*Gt123*gt12L*gtu32 + 4*Gt223*gt22L*gtu32 + 
-           6*Gt222*gt23L*gtu32 + 4*gt12L*Gt133*gtu33 + 2*Gt123*gt13L*gtu33 + 4*gt22L*Gt233*gtu33 + 4*gt23L*Gt333*gtu33 + 
-           2*gt23L*Xt3L) + gt11L*gtu11*SQR(Gt112) + 4*
-         (Gt112*Gt211*gt22L*gtu11 + Gt112*gt23L*Gt311*gtu11 + Gt111*Gt122*gt12L*gtu21 + Gt122*Gt211*gt22L*gtu21 + 
-           Gt112*Gt212*gt22L*gtu21 + Gt223*gt22L*Gt312*gtu21 + Gt112*gt23L*Gt312*gtu21 + Gt122*Gt212*gt22L*gtu22 + 
-           Gt112*Gt113*gt12L*gtu31 + Gt123*Gt211*gt22L*gtu31 + Gt112*Gt213*gt22L*gtu31 + Gt123*Gt212*gt22L*gtu32 + 
-           Gt123*gt23L*Gt313*gtu33 + gt12L*gtu21*SQR(Gt112)) + gt11L*gtu22*SQR(Gt122) + gt11L*gtu33*SQR(Gt123) + 
-        5*gt22L*gtu11*SQR(Gt212) + 5*gt22L*gtu22*SQR(Gt222) + 5*gt22L*gtu33*SQR(Gt223) + gt33L*gtu11*SQR(Gt312) + 
-        gt33L*gtu22*SQR(Gt322) + 4*gt23L*gtu32*SQR(Gt323) + gt33L*gtu33*SQR(Gt323);
-    
-    Rt23  =  gt13L*PDstandard4th2Xt1 + gt33L*PDstandard4th2Xt3 + 
-        khalf*(-(gtu11*PDstandard4th11gt23) - 2*gtu21*PDstandard4th12gt23 - 2*gtu31*PDstandard4th13gt23 - 
-           gtu22*PDstandard4th22gt23 - 2*gtu32*PDstandard4th23gt23 - gtu33*PDstandard4th33gt23) + gt12L*PDstandard4th3Xt1 + 
-        gt22L*PDstandard4th3Xt2 + gt23L*(PDstandard4th2Xt2 + PDstandard4th3Xt3) + 
-        (Gt113*gt12L + Gt112*gt13L + Gt213*gt22L + gt23L*(Gt212 + Gt313) + Gt312*gt33L)*Xt1L + 
-        (Gt123*gt12L + Gt122*gt13L + Gt223*gt22L + gt23L*(Gt222 + Gt323) + Gt322*gt33L)*Xt2L + 
-        (gt12L*Gt133 + Gt123*gt13L + gt22L*Gt233 + gt23L*(Gt223 + Gt333) + Gt323*gt33L)*Xt3L + 
-        gtu21*(gt11L*(Gt113*Gt122 + Gt112*Gt123) + gt12L*
-            (Gt113*Gt222 + 3*(Gt122*Gt213 + Gt112*Gt223) + Gt123*(Gt212 + 2*Gt313)) + 
-           gt23L*(2*Gt123*Gt311 + Gt222*(4*Gt212 + Gt313) + 5*(Gt223*Gt312 + Gt213*Gt322) + (Gt212 + 4*Gt313)*Gt323) + 
-           gt22L*(Gt212*(2*Gt113 + 3*Gt223) + 2*Gt223*Gt313 + Gt213*(3*Gt222 + 2*Gt323)) + 
-           3*(Gt313*Gt322 + Gt312*Gt323)*gt33L + 2*(Gt111*(Gt123*gt12L + Gt122*gt13L) + gt13L*(Gt122*Gt212 + Gt112*Gt222) + 
-              Gt211*(Gt123*gt22L + Gt122*gt23L) + Gt113*(gt23L*Gt312 + gt12L*Gt323) + 
-              (Gt122*Gt311 + Gt222*Gt312 + Gt212*Gt322)*gt33L + Gt112*(Gt113*gt12L + Gt212*gt23L + Gt312*gt33L)) + 
-           gt13L*(Gt122*Gt313 + 3*(Gt123*Gt312 + Gt113*Gt322) + Gt112*Gt323 + 2*SQR(Gt112))) + 
-        gtu31*(Gt133*(Gt112*gt11L + gt12L*Gt212 + 2*(Gt211*gt22L + gt23L*Gt311)) + 
-           gt23L*(5*Gt233*Gt312 + Gt223*(4*Gt212 + Gt313) + Gt213*(2*Gt112 + 5*Gt323) + Gt212*Gt333) + 
-           Gt113*(gt11L*Gt123 + 2*(Gt112*gt13L + Gt213*gt22L + gt23L*Gt313) + gt13L*Gt323 + gt12L*(Gt223 + 2*Gt333)) + 
-           Gt123*(3*gt12L*Gt213 + gt13L*(2*Gt212 + Gt313) + 2*Gt311*gt33L) + 
-           Gt333*(Gt112*gt13L + 4*gt23L*Gt313 + 3*Gt312*gt33L) + 
-           3*(Gt112*gt12L*Gt233 + gt22L*(Gt213*Gt223 + Gt212*Gt233) + Gt133*gt13L*Gt312 + Gt313*Gt323*gt33L) + 
-           2*(Gt111*(gt12L*Gt133 + Gt123*gt13L) + Gt123*Gt211*gt23L + gt13L*(Gt112*Gt223 + Gt113*Gt323) + 
-              Gt213*gt22L*Gt333 + (Gt223*Gt312 + Gt212*Gt323)*gt33L + Gt313*(gt12L*Gt133 + gt22L*Gt233 + Gt112*gt33L) + 
-              gt12L*SQR(Gt113))) + gtu11*(Gt113*(Gt112*gt11L + gt12L*Gt212 + 2*(Gt211*gt22L + gt23L*Gt311)) + 
-           (Gt112*gt13L + Gt212*gt23L)*Gt313 + Gt213*(5*gt23L*Gt312 + 2*gt22L*Gt313) + 
-           3*(Gt213*(Gt112*gt12L + Gt212*gt22L) + Gt312*(Gt113*gt13L + Gt313*gt33L)) + 
-           2*(Gt111*(Gt113*gt12L + Gt112*gt13L) + Gt113*gt12L*Gt313 + Gt212*Gt312*gt33L + 
-              Gt112*(gt13L*Gt212 + Gt211*gt23L + Gt311*gt33L) + gt23L*(SQR(Gt212) + SQR(Gt313)))) + 
-        gtu32*(Gt133*(gt11L*Gt122 + gt12L*Gt222 + 3*gt13L*Gt322) + gt12L*(4*Gt123*Gt223 + Gt122*Gt233 + 2*Gt133*Gt323) + 
-           Gt233*(3*gt23L*Gt322 + 2*gt22L*Gt323) + gt23L*Gt323*(4*Gt223 + 2*Gt333) + Gt122*(gt13L*Gt333 + 2*Gt313*gt33L) + 
-           Gt222*(gt22L*Gt233 + gt23L*(2*Gt223 + Gt333) + 2*Gt323*gt33L) + gt11L*SQR(Gt123) + 
-           3*(Gt322*Gt333*gt33L + gt22L*SQR(Gt223)) + gt33L*SQR(Gt323) + 
-           2*(Gt113*(Gt123*gt12L + Gt122*gt13L) + Gt133*(Gt112*gt12L + Gt212*gt22L + gt23L*Gt312) + 
-              Gt233*(Gt122*gt12L + Gt222*gt22L + gt23L*Gt322) + gt13L*(Gt122*Gt223 + Gt123*(Gt112 + Gt323)) + 
-              Gt223*gt22L*Gt333 + Gt123*(Gt213*gt22L + gt23L*Gt313 + gt13L*(Gt222 + Gt323) + gt12L*Gt333) + 
-              gt23L*(Gt123*Gt212 + Gt122*Gt213 + Gt223*(Gt222 + Gt323) + Gt323*Gt333) + 
-              gt33L*(Gt123*Gt312 + Gt223*Gt322 + SQR(Gt323)))) + 
-        gtu22*(Gt123*(gt11L*Gt122 + gt12L*Gt222 + 2*(Gt212*gt22L + gt23L*Gt312)) + (Gt122*gt13L + Gt222*gt23L)*Gt323 + 
-           Gt223*(5*gt23L*Gt322 + 2*gt22L*Gt323) + 3*
-            (Gt223*(Gt122*gt12L + Gt222*gt22L) + Gt322*(Gt123*gt13L + Gt323*gt33L)) + 
-           2*(Gt112*(Gt123*gt12L + Gt122*gt13L) + Gt123*gt12L*Gt323 + Gt222*Gt322*gt33L + 
-              Gt122*(gt13L*Gt222 + Gt212*gt23L + Gt312*gt33L) + gt23L*(SQR(Gt222) + SQR(Gt323)))) + 
-        gtu33*(Gt133*(gt11L*Gt123 + gt12L*Gt223 + 2*(Gt213*gt22L + gt23L*Gt313)) + (Gt123*gt13L + Gt223*gt23L)*Gt333 + 
-           Gt233*(5*gt23L*Gt323 + 2*gt22L*Gt333) + 3*
-            ((Gt123*gt12L + Gt223*gt22L)*Gt233 + Gt323*(Gt133*gt13L + Gt333*gt33L)) + 
-           2*(Gt113*(gt12L*Gt133 + Gt123*gt13L) + gt12L*Gt133*Gt333 + Gt223*Gt323*gt33L + 
-              Gt123*(gt13L*Gt223 + Gt213*gt23L + Gt313*gt33L) + gt23L*(SQR(Gt223) + SQR(Gt333))));
-    
-    Rt33  =  6*(Gt133*gt13L*Gt313*gtu31 + Gt233*gt23L*Gt313*gtu31 + Gt113*gt13L*Gt333*gtu31 + Gt213*gt23L*Gt333*gtu31 + 
-           Gt133*gt13L*Gt323*gtu32 + Gt133*gt13L*Gt333*gtu33) + 
-        gtu21*(4*Gt122*gt13L*Gt213 + 2*Gt113*gt12L*Gt223 + 4*Gt112*gt13L*Gt223 + 2*Gt213*Gt223*gt22L + 
-           4*Gt123*Gt211*gt23L + 4*Gt212*Gt223*gt23L + 6*Gt123*gt13L*Gt313 + 6*Gt223*gt23L*Gt313 + 6*Gt113*gt13L*Gt323 + 
-           6*Gt213*gt23L*Gt323 + 4*Gt123*Gt311*gt33L - PDstandard4th12gt33) - gtu31*PDstandard4th13gt33 + 
-        gtu22*(4*Gt222*Gt223*gt23L + 6*Gt123*gt13L*Gt323 + 6*Gt223*gt23L*Gt323 + 4*Gt223*Gt322*gt33L - 
-           khalf*PDstandard4th22gt33) + gtu32*(2*Gt223*gt22L*Gt233 + 4*Gt133*Gt212*gt23L + 6*Gt233*gt23L*Gt323 + 
-           6*Gt123*gt13L*Gt333 + 6*Gt223*gt23L*Gt333 + 4*Gt123*Gt313*gt33L + 10*Gt323*Gt333*gt33L - PDstandard4th23gt33) - 
-        gtu33*khalf*PDstandard4th33gt33 + 2*(gt12L*Gt133*Gt213*gtu31 + Gt113*gt12L*Gt233*gtu31 + gt12L*Gt133*Gt223*gtu32 + 
-           gt13L*PDstandard4th3Xt1) + 2*gt23L*PDstandard4th3Xt2 + 
-        gt33L*(4*Gt213*Gt322*gtu21 + 10*Gt313*Gt323*gtu21 + 4*Gt123*Gt312*gtu22 + 4*Gt133*Gt311*gtu31 + 
-           4*Gt213*Gt323*gtu31 + 10*Gt313*Gt333*gtu31 + 4*Gt133*Gt312*gtu32 + 4*Gt133*Gt313*gtu33 + 2*PDstandard4th3Xt3) + 
-        2*Gt213*gt23L*Xt1L + 2*Gt313*gt33L*Xt1L + Gt113*
-         (4*Gt111*gt13L*gtu11 + 2*gt12L*Gt213*gtu11 + 2*gt11L*Gt123*gtu21 + 2*gt11L*Gt133*gtu31 + 4*Gt123*gt13L*gtu32 + 
-           4*Gt133*gt13L*gtu33 + 2*gt13L*Xt1L) + 2*Gt223*gt23L*Xt2L + 2*Gt323*gt33L*Xt2L + 
-        Gt123*(4*Gt111*gt13L*gtu21 + 2*gt12L*Gt213*gtu21 + 4*Gt112*gt13L*gtu22 + 2*gt12L*Gt223*gtu22 + 
-           4*Gt212*gt23L*gtu22 + 4*gt13L*Gt213*gtu31 + 2*gt11L*Gt133*gtu32 + 4*gt13L*Gt233*gtu33 + 2*gt13L*Xt2L) + 
-        2*Gt133*gt13L*Xt3L + 2*Gt333*gt33L*Xt3L + Gt233*
-         (4*Gt112*gt13L*gtu31 + 2*Gt213*gt22L*gtu31 + 2*Gt123*gt12L*gtu32 + 2*gt12L*Gt133*gtu33 + 4*Gt223*gt23L*gtu33 + 
-           6*gt23L*Gt333*gtu33 + 4*Gt323*gt33L*gtu33 + 2*gt23L*Xt3L) + 
-        gtu11*(4*Gt212*Gt213*gt23L + 6*Gt113*gt13L*Gt313 + 6*Gt213*gt23L*Gt313 + 4*Gt113*Gt311*gt33L + 
-           4*Gt213*Gt312*gt33L - khalf*PDstandard4th11gt33 + gt11L*SQR(Gt113)) + 
-        4*(Gt112*gt13L*Gt213*gtu11 + Gt113*Gt211*gt23L*gtu11 + Gt112*Gt113*gt13L*gtu21 + Gt113*Gt212*gt23L*gtu21 + 
-           Gt213*Gt222*gt23L*gtu21 + Gt113*Gt312*gt33L*gtu21 + Gt223*Gt312*gt33L*gtu21 + Gt122*gt13L*Gt223*gtu22 + 
-           Gt111*Gt133*gt13L*gtu31 + Gt133*Gt211*gt23L*gtu31 + Gt113*Gt213*gt23L*gtu31 + Gt213*Gt223*gt23L*gtu31 + 
-           Gt212*Gt233*gt23L*gtu31 + Gt233*Gt312*gt33L*gtu31 + Gt113*Gt313*gt33L*gtu31 + Gt112*Gt133*gt13L*gtu32 + 
-           Gt123*gt13L*Gt223*gtu32 + Gt122*gt13L*Gt233*gtu32 + Gt123*Gt213*gt23L*gtu32 + Gt222*Gt233*gt23L*gtu32 + 
-           Gt233*Gt322*gt33L*gtu32 + Gt223*Gt323*gt33L*gtu32 + Gt133*Gt213*gt23L*gtu33 + gt13L*gtu31*SQR(Gt113)) + 
-        gt11L*gtu22*SQR(Gt123) + gt11L*gtu33*SQR(Gt133) + gt22L*gtu11*SQR(Gt213) + gt22L*gtu22*SQR(Gt223) + 
-        4*gt23L*gtu32*SQR(Gt223) + gt22L*gtu33*SQR(Gt233) + 5*gt33L*gtu11*SQR(Gt313) + 5*gt33L*gtu22*SQR(Gt323) + 
-        5*gt33L*gtu33*SQR(Gt333);
+    Rt11  =  khalf*(-(gtu22*PDstandard4th11gt22) - 2*(Gt111*Gt122 + Gt111*Gt133 + Gt211*Gt222 + Gt211*Gt233 + Gt311*Gt322 + 
+             Gt311*Gt333 + gtu32*PDstandard4th11gt23) - gtu33*PDstandard4th11gt33 + 2*gtu22*PDstandard4th12gt12 + 
+          2*gtu32*PDstandard4th12gt13 + 2*gtu32*PDstandard4th13gt12 + 2*gtu33*PDstandard4th13gt13 - 
+          gtu22*PDstandard4th22gt11 - 2*gtu32*PDstandard4th23gt11 - gtu33*PDstandard4th33gt11 + 2*SQR(Gt112) + 
+          2*SQR(Gt113) + 2*SQR(Gt212) + 2*SQR(Gt213) + 2*SQR(Gt312) + 2*SQR(Gt313));
+    
+    Rt12  =  khalf*(2*(Gt113*Gt123 + Gt213*Gt223 + Gt313*Gt323) - 
+          2*(Gt112*Gt133 + Gt212*Gt233 + Gt312*Gt333 + gtu21*PDstandard4th12gt12) - gtu32*PDstandard4th12gt23 - 
+          gtu33*PDstandard4th12gt33 + gtu31*(PDstandard4th11gt23 - PDstandard4th12gt13 - PDstandard4th13gt12) + 
+          gtu32*PDstandard4th13gt22 + gtu33*PDstandard4th13gt23 + gtu21*(PDstandard4th11gt22 + PDstandard4th22gt11) + 
+          gtu32*PDstandard4th22gt13 + gtu31*PDstandard4th23gt11 - gtu32*PDstandard4th23gt12 + gtu33*PDstandard4th23gt13 - 
+          gtu33*PDstandard4th33gt12);
+    
+    Rt13  =  khalf*(2*(Gt112*Gt123 + Gt212*Gt223 + Gt312*Gt323) - 
+          2*(Gt113*Gt122 + Gt213*Gt222 + Gt313*Gt322 + gtu31*PDstandard4th13gt13) + 
+          gtu21*(PDstandard4th11gt23 - PDstandard4th12gt13 - PDstandard4th13gt12 + PDstandard4th23gt11) + 
+          gtu22*(PDstandard4th12gt23 - PDstandard4th13gt22 - PDstandard4th22gt13 + PDstandard4th23gt12) + 
+          gtu31*(PDstandard4th11gt33 + PDstandard4th33gt11) + 
+          gtu32*(PDstandard4th12gt33 - PDstandard4th13gt23 - PDstandard4th23gt13 + PDstandard4th33gt12));
+    
+    Rt22  =  khalf*(-2*(Gt122*(Gt111 + Gt133) + Gt222*(Gt211 + Gt233) + Gt322*(Gt311 + Gt333) + gtu31*PDstandard4th13gt22) + 
+          gtu11*(-PDstandard4th11gt22 + 2*PDstandard4th12gt12 - PDstandard4th22gt11) + 
+          gtu31*(-2*PDstandard4th22gt13 + 2*(PDstandard4th12gt23 + PDstandard4th23gt12)) + 
+          gtu33*(-PDstandard4th22gt33 + 2*PDstandard4th23gt23 - PDstandard4th33gt22) + 
+          2*(SQR(Gt112) + SQR(Gt123) + SQR(Gt212) + SQR(Gt223) + SQR(Gt312) + SQR(Gt323)));
+    
+    Rt23  =  khalf*(2*(Gt112*Gt113 + Gt212*Gt213 + Gt312*Gt313) + 
+          gtu11*(-PDstandard4th11gt23 + PDstandard4th12gt13 + PDstandard4th13gt12 - PDstandard4th23gt11) + 
+          gtu21*(-PDstandard4th12gt23 + PDstandard4th13gt22 + PDstandard4th22gt13 - PDstandard4th23gt12) - 
+          2*(Gt111*Gt123 + Gt211*Gt223 + Gt311*Gt323 + gtu32*PDstandard4th23gt23) + 
+          gtu31*(PDstandard4th12gt33 - PDstandard4th13gt23 - PDstandard4th23gt13 + PDstandard4th33gt12) + 
+          gtu32*(PDstandard4th22gt33 + PDstandard4th33gt22));
+    
+    Rt33  =  khalf*(gtu11*(-PDstandard4th11gt33 + 2*PDstandard4th13gt13 - PDstandard4th33gt11) - 
+          2*((Gt111 + Gt122)*Gt133 + (Gt211 + Gt222)*Gt233 + (Gt311 + Gt322)*Gt333 + 
+             gtu21*(PDstandard4th12gt33 + PDstandard4th33gt12)) + 
+          gtu22*(-PDstandard4th22gt33 + 2*PDstandard4th23gt23 - PDstandard4th33gt22) + 
+          2*(gtu21*(PDstandard4th13gt23 + PDstandard4th23gt13) + SQR(Gt113) + SQR(Gt123) + SQR(Gt213) + SQR(Gt223) + 
+             SQR(Gt313) + SQR(Gt323)));
     
     Rphi11  =  2*(-PDstandard4th11phi - gt11L*gtu11*PDstandard4th11phi - 2*gt11L*gtu21*PDstandard4th12phi - 
           2*gt11L*gtu31*PDstandard4th13phi - gt11L*gtu22*PDstandard4th22phi - 2*gt11L*gtu32*PDstandard4th23phi - 
@@ -927,7 +680,7 @@ void ML_BSSN_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REA
     M3L  =  (-gK113 + gK131)*gu11 + (-2*gK123 + gK132 + gK231)*gu21 + (-gK223 + gK232)*gu22 + (-gK133 + gK331)*gu31 + 
         (-gK233 + gK332)*gu32;
     
-    cSL  =  Log(detgt);
+    cSL  =  trR;
     
     cXt1L  =  Gt111*gtu11 + Gt122*gtu22 + 2*(Gt112*gtu21 + Gt113*gtu31 + Gt123*gtu32) + Gt133*gtu33 - Xt1L;
     
diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBase.c b/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
index 973cf21..c1a4db7 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBase.c
@@ -101,16 +101,16 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     CCTK_REAL K11 = INITVALUE, K12 = INITVALUE, K13 = INITVALUE, K22 = INITVALUE, K23 = INITVALUE, K33 = INITVALUE;
     
     /* Declare local copies of grid functions */
+    CCTK_REAL AL = INITVALUE;
     CCTK_REAL alpL = INITVALUE;
     CCTK_REAL alphaL = INITVALUE;
     CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE;
+    CCTK_REAL B1L = INITVALUE, B2L = INITVALUE, B3L = INITVALUE;
     CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE;
     CCTK_REAL betaxL = INITVALUE;
     CCTK_REAL betayL = INITVALUE;
     CCTK_REAL betazL = INITVALUE;
     CCTK_REAL dtalpL = INITVALUE;
-    CCTK_REAL dtalphaL = INITVALUE;
-    CCTK_REAL dtbeta1L = INITVALUE, dtbeta2L = INITVALUE, dtbeta3L = INITVALUE;
     CCTK_REAL dtbetaxL = INITVALUE;
     CCTK_REAL dtbetayL = INITVALUE;
     CCTK_REAL dtbetazL = INITVALUE;
@@ -236,7 +236,7 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     
     alphaL  =  alpL;
     
-    dtalphaL  =  dtalpL;
+    AL  =  dtalpL;
     
     beta1L  =  betaxL;
     
@@ -244,14 +244,15 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     
     beta3L  =  betazL;
     
-    dtbeta1L  =  dtbetaxL;
+    B1L  =  dtbetaxL;
     
-    dtbeta2L  =  dtbetayL;
+    B2L  =  dtbetayL;
     
-    dtbeta3L  =  dtbetazL;
+    B3L  =  dtbetazL;
     
     
     /* Copy local copies back to grid functions */
+    A[index] = AL;
     alpha[index] = alphaL;
     At11[index] = At11L;
     At12[index] = At12L;
@@ -259,13 +260,12 @@ void ML_BSSN_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C
     At22[index] = At22L;
     At23[index] = At23L;
     At33[index] = At33L;
+    B1[index] = B1L;
+    B2[index] = B2L;
+    B3[index] = B3L;
     beta1[index] = beta1L;
     beta2[index] = beta2L;
     beta3[index] = beta3L;
-    dtalpha[index] = dtalphaL;
-    dtbeta1[index] = dtbeta1L;
-    dtbeta2[index] = dtbeta2L;
-    dtbeta3[index] = dtbeta3L;
     gt11[index] = gt11L;
     gt12[index] = gt12L;
     gt13[index] = gt13L;
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBase.c b/ML_BSSN/src/ML_BSSN_convertToADMBase.c
index f5df769..d942c00 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBase.c
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBase.c
@@ -99,16 +99,16 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     CCTK_REAL K11 = INITVALUE, K12 = INITVALUE, K13 = INITVALUE, K22 = INITVALUE, K23 = INITVALUE, K33 = INITVALUE;
     
     /* Declare local copies of grid functions */
+    CCTK_REAL AL = INITVALUE;
     CCTK_REAL alpL = INITVALUE;
     CCTK_REAL alphaL = INITVALUE;
     CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE;
+    CCTK_REAL B1L = INITVALUE, B2L = INITVALUE, B3L = INITVALUE;
     CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE;
     CCTK_REAL betaxL = INITVALUE;
     CCTK_REAL betayL = INITVALUE;
     CCTK_REAL betazL = INITVALUE;
     CCTK_REAL dtalpL = INITVALUE;
-    CCTK_REAL dtalphaL = INITVALUE;
-    CCTK_REAL dtbeta1L = INITVALUE, dtbeta2L = INITVALUE, dtbeta3L = INITVALUE;
     CCTK_REAL dtbetaxL = INITVALUE;
     CCTK_REAL dtbetayL = INITVALUE;
     CCTK_REAL dtbetazL = INITVALUE;
@@ -132,6 +132,7 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     /* Declare derivatives */
     
     /* Assign local copies of grid functions */
+    AL = A[index];
     alphaL = alpha[index];
     At11L = At11[index];
     At12L = At12[index];
@@ -139,13 +140,12 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     At22L = At22[index];
     At23L = At23[index];
     At33L = At33[index];
+    B1L = B1[index];
+    B2L = B2[index];
+    B3L = B3[index];
     beta1L = beta1[index];
     beta2L = beta2[index];
     beta3L = beta3[index];
-    dtalphaL = dtalpha[index];
-    dtbeta1L = dtbeta1[index];
-    dtbeta2L = dtbeta2[index];
-    dtbeta3L = dtbeta3[index];
     gt11L = gt11[index];
     gt12L = gt12[index];
     gt13L = gt13[index];
@@ -216,7 +216,7 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     
     alpL  =  alphaL;
     
-    dtalpL  =  dtalphaL;
+    dtalpL  =  AL;
     
     betaxL  =  beta1L;
     
@@ -224,11 +224,11 @@ void ML_BSSN_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCT
     
     betazL  =  beta3L;
     
-    dtbetaxL  =  dtbeta1L;
+    dtbetaxL  =  B1L;
     
-    dtbetayL  =  dtbeta2L;
+    dtbetayL  =  B2L;
     
-    dtbetazL  =  dtbeta3L;
+    dtbetazL  =  B3L;
     
     
     /* Copy local copies back to grid functions */
diff --git a/ML_BSSN/src/RegisterMoL.c b/ML_BSSN/src/RegisterMoL.c
index aaea9c7..0a5e78e 100644
--- a/ML_BSSN/src/RegisterMoL.c
+++ b/ML_BSSN/src/RegisterMoL.c
@@ -21,10 +21,10 @@ void ML_BSSN_RegisterVars(CCTK_ARGUMENTS)
   ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::At22"),  CCTK_VarIndex("ML_BSSN::At22rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::At23"),  CCTK_VarIndex("ML_BSSN::At23rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::At33"),  CCTK_VarIndex("ML_BSSN::At33rhs"));
-  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::dtalpha"),  CCTK_VarIndex("ML_BSSN::dtalpharhs"));
-  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::dtbeta1"),  CCTK_VarIndex("ML_BSSN::dtbeta1rhs"));
-  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::dtbeta2"),  CCTK_VarIndex("ML_BSSN::dtbeta2rhs"));
-  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::dtbeta3"),  CCTK_VarIndex("ML_BSSN::dtbeta3rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::A"),  CCTK_VarIndex("ML_BSSN::Arhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::B1"),  CCTK_VarIndex("ML_BSSN::B1rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::B2"),  CCTK_VarIndex("ML_BSSN::B2rhs"));
+  ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::B3"),  CCTK_VarIndex("ML_BSSN::B3rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::Xt1"),  CCTK_VarIndex("ML_BSSN::Xt1rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::Xt2"),  CCTK_VarIndex("ML_BSSN::Xt2rhs"));
   ierr += MoLRegisterEvolved(CCTK_VarIndex("ML_BSSN::Xt3"),  CCTK_VarIndex("ML_BSSN::Xt3rhs"));
diff --git a/ML_BSSN/src/RegisterSymmetries.c b/ML_BSSN/src/RegisterSymmetries.c
index 2cf86ef..5bcde48 100644
--- a/ML_BSSN/src/RegisterSymmetries.c
+++ b/ML_BSSN/src/RegisterSymmetries.c
@@ -52,22 +52,22 @@ void ML_BSSN_RegisterSymmetries(CCTK_ARGUMENTS)
   sym[0] = 1;
   sym[1] = 1;
   sym[2] = 1;
-  SetCartSymVN(cctkGH, sym, "ML_BSSN::dtalpha");
+  SetCartSymVN(cctkGH, sym, "ML_BSSN::A");
   
   sym[0] = -1;
   sym[1] = 1;
   sym[2] = 1;
-  SetCartSymVN(cctkGH, sym, "ML_BSSN::dtbeta1");
+  SetCartSymVN(cctkGH, sym, "ML_BSSN::B1");
   
   sym[0] = 1;
   sym[1] = -1;
   sym[2] = 1;
-  SetCartSymVN(cctkGH, sym, "ML_BSSN::dtbeta2");
+  SetCartSymVN(cctkGH, sym, "ML_BSSN::B2");
   
   sym[0] = 1;
   sym[1] = 1;
   sym[2] = -1;
-  SetCartSymVN(cctkGH, sym, "ML_BSSN::dtbeta3");
+  SetCartSymVN(cctkGH, sym, "ML_BSSN::B3");
   
   sym[0] = -1;
   sym[1] = 1;
diff --git a/m/McLachlan.m b/m/McLachlan.m
index 299a221..0bb1b34 100644
--- a/m/McLachlan.m
+++ b/m/McLachlan.m
@@ -43,7 +43,7 @@ KD = KroneckerDelta;
 (* Register the tensor quantities with the TensorTools package *)
 Map [DefineTensor,
      {g, K, alpha, beta, H, M, detg, gu, G, R, trR, Km, trK,
-      phi, gt, At, Xt, dtalpha, dtbeta, Atm, Atu, trA, cXt, cS, cA,
+      phi, gt, At, Xt, A, B, Atm, Atu, trA, cXt, cS, cA,
       e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gt, Rt, Rphi, gK}];
 
 (* NOTE: It seems as if Lie[.,.] did not take these tensor weights
@@ -108,15 +108,15 @@ declaredGroups = Join [evolvedGroups, evaluatedGroups];
 declaredGroupNames = Map [First, declaredGroups];
 
 evolvedGroupsBSSN =
-  {SetGroupName [CreateGroupFromTensor [phi       ], "log_confac"],
-   SetGroupName [CreateGroupFromTensor [gt[la,lb] ], "metric"    ],
-   SetGroupName [CreateGroupFromTensor [Xt[ua]    ], "Gamma"     ],
-   SetGroupName [CreateGroupFromTensor [trK       ], "trace_curv"],
-   SetGroupName [CreateGroupFromTensor [At[la,lb] ], "curv"      ],
-   SetGroupName [CreateGroupFromTensor [alpha     ], "lapse"     ],
-   SetGroupName [CreateGroupFromTensor [dtalpha   ], "dtlapse"   ],
-   SetGroupName [CreateGroupFromTensor [beta[ua]  ], "shift"     ],
-   SetGroupName [CreateGroupFromTensor [dtbeta[ua]], "dtshift"   ]};
+  {SetGroupName [CreateGroupFromTensor [phi      ], "log_confac"],
+   SetGroupName [CreateGroupFromTensor [gt[la,lb]], "metric"    ],
+   SetGroupName [CreateGroupFromTensor [Xt[ua]   ], "Gamma"     ],
+   SetGroupName [CreateGroupFromTensor [trK      ], "trace_curv"],
+   SetGroupName [CreateGroupFromTensor [At[la,lb]], "curv"      ],
+   SetGroupName [CreateGroupFromTensor [alpha    ], "lapse"     ],
+   SetGroupName [CreateGroupFromTensor [A        ], "dtlapse"   ],
+   SetGroupName [CreateGroupFromTensor [beta[ua] ], "shift"     ],
+   SetGroupName [CreateGroupFromTensor [B[ua]    ], "dtshift"   ]};
 evaluatedGroupsBSSN =
   {SetGroupName [CreateGroupFromTensor [H      ], "Ham"],
    SetGroupName [CreateGroupFromTensor [M[la]  ], "mom"],
@@ -171,15 +171,15 @@ initialCalcBSSN =
   (* Where -> Interior, *)
   Equations -> 
   {
-    phi        -> 0,
-    gt[la,lb]  -> KD[la,lb],
-    trK        -> 0,
-    At[la,lb]  -> 0,
-    Xt[ua]     -> 0,
-    alpha      -> 1,
-    dtalpha    -> 0,
-    beta[ua]   -> 0,
-    dtbeta[ua] -> 0
+    phi       -> 0,
+    gt[la,lb] -> KD[la,lb],
+    trK       -> 0,
+    At[la,lb] -> 0,
+    Xt[ua]    -> 0,
+    alpha     -> 1,
+    A         -> 0,
+    beta[ua]  -> 0,
+    B[ua]     -> 0
   }
 }
 
@@ -246,17 +246,17 @@ convertFromADMBaseCalcBSSN =
     trK       -> gu[ua,ub] K[la,lb],
     At[la,lb] -> em4phi (K[la,lb] - (1/3) g[la,lb] trK),
     
-    alpha   -> alp,
+    alpha -> alp,
     (* TODO: this is wrong *)
-    dtalpha -> dtalp,
+    A     -> dtalp,
     
-    beta1   -> betax,
-    beta2   -> betay,
-    beta3   -> betaz,
+    beta1 -> betax,
+    beta2 -> betay,
+    beta3 -> betaz,
     (* TODO: this is wrong *)
-    dtbeta1 -> dtbetax,
-    dtbeta2 -> dtbetay,
-    dtbeta3 -> dtbetaz
+    B1    -> dtbetax,
+    B2    -> dtbetay,
+    B3    -> dtbetaz
   }
 }
 
@@ -284,7 +284,7 @@ convertFromADMBaseCalcBSSNGamma =
 convertToADMBaseCalc =
 {
   Name -> "ML_ADM_convertToADMBase",
-  Schedule -> {"IN MoL_PostStep AFTER ADM_ApplyBoundConds"},
+  Schedule -> {"IN MoL_PostStep AFTER ML_ADM_ApplyBCs"},
   Equations -> 
   {
     gxx     -> g11,
@@ -314,7 +314,7 @@ convertToADMBaseCalc =
 convertToADMBaseCalcBSSN =
 {
   Name -> "ML_BSSN_convertToADMBase",
-  Schedule -> {"IN MoL_PostStep AFTER ML_BSSN_ApplyBoundConds AFTER ML_BSSN_enforce"},
+  Schedule -> {"IN MoL_PostStep AFTER ML_BSSN_ApplyBCs AFTER ML_BSSN_enforce"},
   Shorthands -> {e4phi, g[la,lb], K[la,lb]},
   Equations -> 
   {
@@ -336,14 +336,14 @@ convertToADMBaseCalcBSSN =
     alp      -> alpha,
     (* TODO: this is wrong *)
     (* TODO: rename dtalp->A, dtbeta->B *)
-    dtalp    -> dtalpha,
+    dtalp    -> A,
     betax    -> beta1,
     betay    -> beta2,
     betaz    -> beta3,
     (* TODO: this is wrong *)
-    dtbetax  -> dtbeta1,
-    dtbetay  -> dtbeta2,
-    dtbetaz  -> dtbeta3
+    dtbetax  -> B1,
+    dtbetay  -> B2,
+    dtbetaz  -> B3
   }
 }
 
@@ -466,19 +466,19 @@ evolCalcBSSN =
                       + Lie[At[la,lb], beta] - (2/3) At[la,lb] PD[beta[uc],lc],
     
     (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *)
-    dot[alpha]      -> - harmonicF alpha^harmonicN dtalpha
+    dot[alpha] -> - harmonicF alpha^harmonicN A
                        + Lie[alpha, beta],
-    dot[dtalpha]    -> dot[trK] - AlphaDriver dtalpha,
+    dot[A]    -> dot[trK] - AlphaDriver A,
     (* dot[beta[ua]] -> eta Xt[ua], *)
-    dot[beta[ua]]   -> ShiftGammaCoeff alpha^ShiftAlphaPower dtbeta[ua],
-    dot[dtbeta[ua]] -> dot[Xt[ua]] - BetaDriver dtbeta[ua]
+    dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua],
+    dot[B[ua]]    -> dot[Xt[ua]] - BetaDriver B[ua]
   }
 }
 
 enforceCalcBSSN =
 {
   Name -> "ML_BSSN_enforce",
-  Schedule -> {"IN MoL_PostStep AFTER BSSN_ApplyBoundConds"},
+  Schedule -> {"IN MoL_PostStep BEFORE ML_BSSN_BoundConds"},
   Shorthands -> {detgt, gtu[ua,ub], trA},
   Equations -> 
   {
@@ -537,12 +537,37 @@ constraintsCalcBSSN =
                     (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]),
     
     (* PRD 62, 044034 (2000), eqn. (18) *)
+    (* Note: This differs from the Goddard formulation,
+             which is e.g. described in PRD 70 (2004) 124025, eqn. (6).
+             Goddard has a Gamma^k replaced by its definition in terms
+             of metric derivatives.  *)
+(*
     Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm]
                  + gt[lk,li] PD[Xt[uk],lj] + gt[lk,lj] PD[Xt[uk],li]
                  + Xt[uk] gt[li,ln] Gt[un,lj,lk] + Xt[uk] gt[lj,ln] Gt[un,li,lk]
                  + gtu[ul,um] (+ 2 Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,lm]
                                + 2 Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,lm]
                                + Gt[uk,li,lm] gt[lk,ln] Gt[un,ll,lj]),
+*)
+
+    (* From the long turducken paper.
+       This expression seems to give the same result as the one from 044034.  *)
+(*
+    Rt[li,lj] -> - (1/2) gtu[uk,ul] PD[gt[li,lj],lk,ll]
+                 + gt[lk,li] PD[Xt[uk],lj] + gt[lk,lj] PD[Xt[uk],li]
+                 + gt[li,ln] Gt[un,lj,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm]
+                 + gt[lj,ln] Gt[un,li,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm]
+                 + gtu[ul,us] (+ 2 Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,ls]
+                               + 2 Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,ls]
+                               + Gt[uk,li,ls] gt[lk,ln] Gt[un,ll,lj]),
+*)
+
+    (* Below would be a straightforward calculation,
+       without taking any Gamma^i into account.
+       This expression gives a different answer!  *)
+    Rt[la,lb] -> + Gt[u1,l2,la] Gt[l1,lb,u2] - Gt[u1,la,lb] Gt[l1,l2,u2]
+                 + 1/2 gtu[u1,u2] (- PD[gt[l1,l2],la,lb] + PD[gt[l1,la],l2,lb]
+                                   - PD[gt[la,lb],l1,l2] + PD[gt[l2,lb],l1,la]),
     (* PRD 62, 044034 (2000), eqn. (15) *)
     Rphi[li,lj] -> - 2 CDt[phi,lj,li]
                    - 2 gt[li,lj] gtu[ul,un] CDt[phi,ll,ln]
@@ -558,12 +583,13 @@ constraintsCalcBSSN =
     gu[ua,ub]   -> em4phi gtu[ua,ub],
     (* ddetg[la] -> PD[e4phi detg,la], *)
     ddetg[la]   -> e4phi ddetgt[la] + 4 detgt e4phi PD[phi,la],
+    (* TODO: check this equation, maybe simplify it by omitting ddetg *)
     G[ua,lb,lc] -> Gt[ua,lb,lc]
                    + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb]
                                  - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]),
     
     R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
-    trR -> gu[ua,ub] R[la,lb],
+    trR      -> gu[ua,ub] R[la,lb],
     
     (* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *)
     (* Km[ua,lb] -> gu[ua,uc] K[lc,lb], *)
@@ -578,7 +604,8 @@ constraintsCalcBSSN =
     M[la] -> gu[ub,uc] (gK[lc,la,lb] - gK[lc,lb,la]),
     
     (* det gamma-tilde *)
-    cS -> Log [detgt],
+    (* TODO cS -> Log [detgt], *)
+    cS -> trR,
     
     (* Gamma constraint *)
     cXt[ua] -> gtu[ub,uc] Gt[ua,lb,lc] - Xt[ua],