Some superficial changes:

* Rename Tet -> Theta and Zet -> Z
* Introduce an enum to select between CMBSSNphi, CMBSSNW and CMCCZ4.
* Rename gamashift parameter to GammaShift

1 files changed, 155 insertions, 126 deletions

diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index 6943a3e..e6dd50b 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -103,14 +103,15 @@ Map [DefineTensor,
       xx, rr, th, ph,
       admg, admK, admalpha, admdtalpha, admbeta, admdtbeta, H, M,
       g, detg, gu, G, R, trR, Km, trK, cdphi, cdphi2,
-      phi, gt, At, Xt, Xtn, alpha, A, beta, B, Atm, Atu, trA, Ats, trAts,
+      phi, gt, At, Xt, Xtn, Theta, Z,
+      alpha, A, beta, B, Atm, Atu, trA, Ats, trAts,
       dottrK, dotXt,
       cXt, cS, cA,
       e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt,
       Rt, Rphi, gK,
       T00, T0, T, rho, S,
       x, y, z, r,
-      epsdiss, Tet, Zet, dampk1k2}];
+      epsdiss}];
 
 (* NOTE: It seems as if Lie[.,.] did not take these tensor weights
    into account.  Presumably, CD[.,.] and CDt[.,.] don't do this either.  *)
@@ -159,6 +160,12 @@ T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz;
 (* Expressions *)
 (******************************************************************************)
 
+(* enum constants for conformalMethod; these must be consistent
+   with the definition of the Cactus parameter conformalMethod *)
+CMBSSNphi = 0;
+CMBSSNW   = 1;
+CMCCZ4    = 2;
+
 detgExpr  = Det [MatrixOfComponents [g [la,lb]]];
 ddetgExpr[la_] =
   Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la],
@@ -182,13 +189,13 @@ evolvedGroups =
   {SetGroupName [CreateGroupFromTensor [phi      ], prefix <> "log_confac"],
    SetGroupName [CreateGroupFromTensor [gt[la,lb]], prefix <> "metric"    ],
    SetGroupName [CreateGroupFromTensor [Xt[ua]   ], prefix <> "Gamma"     ],
+   SetGroupName [CreateGroupFromTensor [Theta    ], prefix <> "Theta"     ],
    SetGroupName [CreateGroupFromTensor [trK      ], prefix <> "trace_curv"],
    SetGroupName [CreateGroupFromTensor [At[la,lb]], prefix <> "curv"      ],
    SetGroupName [CreateGroupFromTensor [alpha    ], prefix <> "lapse"     ],
    SetGroupName [CreateGroupFromTensor [A        ], prefix <> "dtlapse"   ],
    SetGroupName [CreateGroupFromTensor [beta[ua] ], prefix <> "shift"     ],
-   SetGroupName [CreateGroupFromTensor [B[ua]    ], prefix <> "dtshift" ],
-   SetGroupName [CreateGroupFromTensor [Tet      ], prefix <> "Zenergy"   ]};
+   SetGroupName [CreateGroupFromTensor [B[ua]    ], prefix <> "dtshift"   ]};
 evaluatedGroups =
   {SetGroupName [CreateGroupFromTensor [H      ], prefix <> "Ham"],
    SetGroupName [CreateGroupFromTensor [M[la]  ], prefix <> "mom"],
@@ -229,16 +236,16 @@ initialCalc =
   ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
   Equations -> 
   {
-    phi       -> IfThen [conformalMethod, 1, 0],
+    phi       -> IfThen[conformalMethod==CMBSSNW, 1, 0],
     gt[la,lb] -> KD[la,lb],
     trK       -> 0,
     At[la,lb] -> 0,
     Xt[ua]    -> 0,
+    Theta     -> 0,
     alpha     -> 1,
     A         -> 0,
     beta[ua]  -> 0,
-    B[ua]     -> 0,
-    Tet       -> 0
+    B[ua]     -> 0
   }
 };
 
@@ -285,13 +292,15 @@ convertFromADMBaseCalc =
     detg      -> detgExpr,
     gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
     
-    phi       -> IfThen [conformalMethod, detg^(-1/6), Log[detg]/12],
-    em4phi    -> IfThen [conformalMethod, phi^2, Exp[-4 phi]],
+    phi       -> IfThen[conformalMethod==CMBSSNW, detg^(-1/6), Log[detg]/12],
+    em4phi    -> IfThen[conformalMethod==CMBSSNW, phi^2, Exp[-4 phi]],
     gt[la,lb] -> em4phi g[la,lb],
     
     trK       -> gu[ua,ub] admK[la,lb],
     At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK),
     
+    Theta     -> 0,
+    
     alpha     -> admalpha,
     
     beta[ua]  -> admbeta[ua]
@@ -365,8 +374,7 @@ initGammaCalc =
   {
     Xt[ua] -> 0,
     A      -> 0,
-    B[ua]  -> 0,
-    Tet    -> 0
+    B[ua]  -> 0
   }
 };
 
@@ -384,7 +392,7 @@ convertToADMBaseCalc =
   Shorthands -> {e4phi},
   Equations -> 
   {
-    e4phi       -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]],
+    e4phi       -> IfThen[conformalMethod==CMBSSNW, 1/phi^2, Exp[4 phi]],
     admg[la,lb] -> e4phi gt[la,lb],
     admK[la,lb] -> e4phi At[la,lb] + (1/3) admg[la,lb] trK,
     admalpha    -> alpha,
@@ -413,13 +421,13 @@ convertToADMBaseDtLapseShiftCalc =
     (* see RHS *)
 (*
     admdtalpha -> - harmonicF alpha^harmonicN
-                    ((1 - LapseAdvectionCoeff) A
-                  + LapseAdvectionCoeff (trK - 2 Tet))
+                    ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
                   + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
 *)
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
-                     + (1 - LapseACoeff) (trK - 2 Tet))
+                     + ((1 - LapseACoeff)
+                        (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0])))
                   + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
     admdtbeta[ua] -> IfThen[harmonicShift,
                             - 1/2 gtu[ua,uj] phi alpha
@@ -455,12 +463,12 @@ convertToADMBaseDtLapseShiftBoundaryCalc =
     (* see RHS, but omit derivatives near the boundary *)
 (*
     admdtalpha -> - harmonicF alpha^harmonicN
-                    ((1 - LapseAdvectionCoeff) A
-                    + LapseAdvectionCoeff (trK - 2 Tet)),
+                    ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
 *)
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
-                     + (1 - LapseACoeff) (trK - 2 Tet)),
+                     + ((1 - LapseACoeff)
+                        (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]))),
     admdtbeta[ua] -> IfThen[harmonicShift,
                             0,
                             (* else *)
@@ -492,12 +500,12 @@ convertToADMBaseFakeDtLapseShiftCalc =
        conditions) *)
 (*
     admdtalpha -> - harmonicF alpha^harmonicN
-                    ((1 - LapseAdvectionCoeff) A
-                    + LapseAdvectionCoeff (trK - 2 Tet)),
+                    ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
 *)
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
-                     + (1 - LapseACoeff) (trK - 2 Tet)),
+                     + ((1 - LapseACoeff)
+                        (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]))),
     admdtbeta[ua] -> IfThen[harmonicShift,
                             0,
                             (* else *)
@@ -524,14 +532,14 @@ evolCalc =
      radiative boundary conditions.  *)
   Where -> InteriorNoSync,
   Shorthands -> {dir[ua],
-                 detgt, gtu[ua,ub],
+                 detgt, gtu[ua,ub], Z[ua],
                  Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua],
                  Rt[la,lb], Rphi[la,lb], R[la,lb],
                  Atm[ua,lb], Atu[ua,ub],
                  e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg,
                  gu[ua,ub], Ats[la,lb], trAts, eta, theta,
-                 rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua],
-                 epsdiss[ua], dotTet, Zet[ua], dampk1k2},
+                 rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua], dotTheta,
+                 epsdiss[ua]},
   Equations -> 
   {
     dir[ua] -> Sign[beta[ua]],
@@ -549,7 +557,7 @@ evolCalc =
        used instead of Xt where no derivatives of Xt are taken *)
     Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk],
 
-    e4phi       -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]],
+    e4phi       -> IfThen[conformalMethod==CMBSSNW, 1/phi^2, Exp[4 phi]],
     em4phi      -> 1 / e4phi,
     g[la,lb]    -> e4phi gt[la,lb],
     detg        -> detgExpr,
@@ -557,13 +565,13 @@ evolCalc =
 
     (* The Z quantities *)
     (* gr-qc:1106.2254 (2011), eqn. (23) *)
-    Zet[ud] -> (1/2) gu[ua,ud] (-PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc]),
-
-    (* shorthand for the product of damping terms k1 and k2 in gr-qc:1106.2254 (2011) *)
-    dampk1k2 -> dampk1 dampk2,
+    Z[ud] -> IfThen[conformalMethod==CMCCZ4,
+                    (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc]
+                                     + gt[la,lc] Xt[uc]),
+                    0],
 
     (* PRD 62, 044034 (2000), eqn. (18) *)
-    (* Added Zet term by changing Xtn to Xt *)
+    (* Adding Z term by changing Xtn to Xt *)
     Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm]
                  + (1/2) gt[lk,li] PD[Xt[uk],lj]
                  + (1/2) gt[lk,lj] PD[Xt[uk],li]
@@ -573,9 +581,9 @@ evolCalc =
                     + Gt[uk,lj,ll] Gtlu[li,lk,ul]
                     + Gt[uk,li,ll] Gtlu[lk,lj,ul]),
 
-    fac1 -> IfThen [conformalMethod, -1/(2 phi), 1],
+    fac1 -> IfThen[conformalMethod==CMBSSNW, -1/(2 phi), 1],
     cdphi[la] -> fac1 CDt[phi,la],
-    fac2 -> IfThen [conformalMethod, 1/(2 phi^2), 0],
+    fac2 -> IfThen[conformalMethod==CMBSSNW, 1/(2 phi^2), 0],
     cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb],
 
     (* PRD 62, 044034 (2000), eqn. (15) *)
@@ -587,13 +595,14 @@ evolCalc =
     Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
     Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
     
-    (* Added Zet terms *)
-    R[la,lb] -> Rt[la,lb] + Rphi[la,lb] +
-                  (2/phi) (g[la,lc] Zet[uc] PD[phi,lb] +
-                  g[lb,lc] Zet[uc] PD[phi,la] -
-                  g[la,lb] Zet[uc] PD[phi,lc]) +
-                  e4phi Zet[uc] PD[gt[la,lb],lc],
-
+    R[la,lb] -> + Rt[la,lb] + Rphi[la,lb]
+                + IfThen[conformalMethod==CMCCZ4,
+                         + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb]
+                                    + g[lb,lc] Z[uc] PD[phi,la]
+                                    - g[la,lb] Z[uc] PD[phi,lc])
+                         + e4phi Z[uc] PD[gt[la,lb],lc],
+                         0],
+    
     (* Matter terms *)
     
     (* rho = n^a n^b T_ab *)
@@ -610,20 +619,23 @@ evolCalc =
     
     (* PRD 62, 044034 (2000), eqn. (10) *)
     (* PRD 67 084023 (2003), eqn. (16) and (23) *)
-    dot[phi]       -> IfThen [conformalMethod, 1/3 phi, -1/6] alpha trK
-                      + IfThen [conformalMethod, -1/3 phi, 1/6] PD[beta[ua],la],
+    dot[phi]       -> IfThen[conformalMethod==CMBSSNW, 1/3 phi, -1/6]
+                      (alpha trK - PD[beta[ua],la]),
     
     (* PRD 62, 044034 (2000), eqn. (9) *)
     (* gr-qc:1106.2254 (2011), eqn. (14) *)
-    (* removing trA from Aij insures that detg = 1 *)
-    dot[gt[la,lb]] -> - 2 alpha (At[la,lb]
-                      - (1/3) At[lc,ld] gtu[uc,ud] gt[la,lb])
+    (* removing trA from Aij ensures that detg = 1 *)
+    dot[gt[la,lb]] -> - 2 alpha (+ At[la,lb]
+                                 - IfThen[conformalMethod==CMCCZ4,
+                                          (1/3) At[lc,ld] gtu[uc,ud] gt[la,lb],
+                                          0])
                       + gt[la,lc] PD[beta[uc],lb] + gt[lb,lc] PD[beta[uc],la]
                       - (2/3) gt[la,lb] PD[beta[uc],lc],
     (* PRD 62, 044034 (2000), eqn. (20) *)
     (* PRD 67 084023 (2003), eqn (26) *)
     (* gr-qc:1106.2254 (2011), eqn. (19) *)
-    (* Added Zet terms by changing Xtn to Xt, and extra Zet and Tet terms *)
+    (* Adding Z terms by changing Xtn to Xt,
+       also adding extra Z and Theta terms *)
     dotXt[ui]      -> - 2 Atu[ui,uj] PD[alpha,lj]
                       + 2 alpha (+ Gt[ui,lj,lk] Atu[uk,uj]
                                  - (2/3) gtu[ui,uj] PD[trK,lj]
@@ -632,34 +644,43 @@ evolCalc =
                       + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll]
                       - Xtn[uj] PD[beta[ui],lj] 
                       + (2/3) Xtn[ui] PD[beta[uj],lj]
-                      + gamashift (2 e4phi) (- Zet[uj] PD[beta[ui],lj]
-                      + (2/3) Zet[ui] PD[beta[uj],lj])
-                      - (4/3) alpha e4phi Zet[ui] trK
-                      + 2 gtu[ui,uj] (alpha PD[Tet,lj] - Tet PD[alpha,lj])
-                      - 2 alpha e4phi dampk1 Zet[ui]
+                      + IfThen[conformalMethod==CMCCZ4,
+                               + IfThen[GammaShift,
+                                        2 e4phi (- Z[uj] PD[beta[ui],lj]
+                                                 + (2/3) Z[ui] PD[beta[uj],lj]),
+                                        0]
+                               - (4/3) alpha e4phi Z[ui] trK
+                               + 2 gtu[ui,uj] (+ alpha PD[Theta,lj]
+                                               - Theta PD[alpha,lj])
+                               - 2 alpha e4phi dampk1 Z[ui],
+                               0]
     (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
                       + addMatter (- 16 Pi alpha gtu[ui,uj] S[lj]),
     dot[Xt[ui]]    -> dotXt[ui],
 
-    (* Equation for Theta *)
     (* gr-qc:1106.2254 (2011), eqn. (18) *)
-    dotTet -> - PD[alpha,la] Zet[ua]
-              + (1/2) alpha (gu[ua,ub] R[la,lb]
-              - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - 2 trK Tet)
-              - dampk1 2 alpha Tet - dampk1k2 alpha Tet,
+    dotTheta       -> IfThen[conformalMethod==CMCCZ4,
+                             - PD[alpha,la] Z[ua]
+                             + (1/2) alpha (gu[ua,ub] R[la,lb]
+                                            - Atm[ua,lb] Atm[ub,la]
+                                            + (2/3) trK^2 - 2 trK Theta)
+                             - dampk1 (2 + dampk2) alpha Theta,
+                             0],
 
-    dot[Tet] -> dotTet,
+    dot[Theta] -> dotTheta,
 
     (* PRD 62, 044034 (2000), eqn. (11) *)
     (* gr-qc:1106.2254 (2011), eqn. (17) *)
-    (* Added the rhs of Tet to K, because K_Z4 = K_BSSN + 2 Tet. *)
-    (* And the Z term, as it has to cancel with the one in Tet *)
+    (* Adding the RHS of Theta to K, because K_Z4 = K_BSSN + 2 Theta *)
+    (* Also adding the Z term, as it has to cancel with the one in Theta *)
     dottrK         -> - em4phi ( gtu[ua,ub] ( PD[alpha,la,lb]
                                 + 2 cdphi[la] PD[alpha,lb] )
                                 - Xtn[ua] PD[alpha,la] )
                       + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
-                      + 2 dotTet + 2 PD[alpha,la] Zet[ua]
-                      + dampk1 alpha Tet - dampk1k2 alpha Tet
+                      + IfThen[conformalMethod==CMCCZ4,
+                               + 2 dotTheta + 2 PD[alpha,la] Z[ua]
+                               + dampk1 (1 - dampk2) alpha Theta,
+                               0]
     (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
                       + addMatter (4 Pi alpha (rho + trS)),
     dot[trK]       -> dottrK,
@@ -667,17 +688,19 @@ evolCalc =
     (* PRD 62, 044034 (2000), eqn. (12) *)
     (* TODO: Should we use the Hamiltonian constraint to make Rij tracefree? *)
     (* gr-qc:1106.2254 (2011), eqn. (15) *)
-    (* Added Zet terms in the Ricci and Tet terms *)
+    (* Adding Z terms in the Ricci and Theta terms *)
     Ats[la,lb]     -> - CDt[alpha,la,lb] +
                       + 2 (PD[alpha,la] cdphi[lb] + PD[alpha,lb] cdphi[la] )
                       + alpha R[la,lb],
     trAts          -> gu[ua,ub] Ats[la,lb],
     dot[At[la,lb]] -> + em4phi (+ Ats[la,lb] - (1/3) g[la,lb] trAts )
-                      + alpha ((trK - 2 Tet) At[la,lb] - 2 At[la,lc] Atm[uc,lb])
+                      + alpha (+ ((trK - IfThen[conformalMethod==CMCCZ4,
+                                                2 Theta,
+                                                0])
+                                  At[la,lb])
+                               - 2 At[la,lc] Atm[uc,lb])
                       + At[la,lc] PD[beta[uc],lb] + At[lb,lc] PD[beta[uc],la]
                       - (2/3) At[la,lb] PD[beta[uc],lc]
-                      (* damping term in trA, alternative to trA cleaning *)
-                      (* - (dampA/3) alpha At[lc,ld] gtu[uc,ud] gt[la,lb] *)
     (* Equation (4.23) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
                       + addMatter (- em4phi alpha 8 Pi
                                      (T[la,lb] - (1/3) g[la,lb] trS)),
@@ -686,19 +709,22 @@ evolCalc =
     (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *)
 (*
     dot[alpha] -> - harmonicF alpha^harmonicN (
-                    (1 - LapseAdvectionCoeff) A
-                  + LapseAdvectionCoeff (trK - 2 Tet))
+                    (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
                   + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
 
 
-    dot[A]     -> (1 - LapseAdvectionCoeff) (dottrK - 2 dotTet - AlphaDriver A),
+    dot[A]     -> (1 - LapseAdvectionCoeff) (dottrK - AlphaDriver A),
 
 *)
     dot[alpha] -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
-                     + (1 - LapseACoeff) (trK + AlphaDriver (alpha - 1) - 2 Tet)),
+                     + ((1 - LapseACoeff)
+                        (+ trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]
+                         + AlphaDriver (alpha - 1)))),
 
-    dot[A]     -> + LapseACoeff (dottrK - 2 dotTet - AlphaDriver A),
+    dot[A]     -> + (LapseACoeff
+                     (+ dottrK - IfThen[conformalMethod==CMCCZ4, 2 dotTheta, 0]
+                      - AlphaDriver A)),
     
     eta -> etaExpr,
     theta -> thetaExpr,
@@ -745,7 +771,10 @@ advectCalc =
 
        dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj],
 
-          dot[Tet] -> dot[Tet] + Upwind[beta[ua], Tet, la],
+        dot[Theta] -> dot[Theta]
+                      + IfThen[conformalMethod==CMCCZ4,
+                               Upwind[beta[ua], Theta, la],
+                               0],
 
           dot[trK] -> dot[trK] + Upwind[beta[ua], trK, la],
 
@@ -775,11 +804,11 @@ evolCalc1 = PartialCalculation[evolCalc, "1",
     ConditionalOnKeyword -> {"RHS_calculation", "split"}
   },
   {
-    dot[Tet],
-    dot[trK],
     dot[phi],
     dot[gt[la,lb]],
     dot[Xt[ui]],
+    dot[Theta],
+    dot[trK],
     dot[alpha],
     dot[A],
     dot[beta[ua]],
@@ -807,7 +836,8 @@ dissCalc =
     epsdiss[ua] -> EpsDiss,
     Sequence@@Table[
       dot[var]       -> dot[var] + epsdiss[ux] PDdiss[var,lx],
-      {var, {phi, gt[la,lb], Xt[ui], trK, Tet, At[la,lb], alpha, A, beta[ua], B[ua]}}]
+      {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb],
+             alpha, A, beta[ua], B[ua]}}]
   }
 };
 
@@ -826,7 +856,8 @@ Table[
     dot[var]    -> dot[var] + epsdiss[ux] PDdiss[var,lx]
   }
 },
-  {var, {phi, gt[la,lb], Xt[ui], trK, At[la,lb], alpha, A, beta[ua], B[ua]}}
+  {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb],
+         alpha, A, beta[ua], B[ua]}}
 ];
 
 RHSStaticBoundaryCalc =
@@ -839,10 +870,10 @@ RHSStaticBoundaryCalc =
   {
     dot[phi]       -> 0,
     dot[gt[la,lb]] -> 0,
-    dot[Tet]       -> 0,
     dot[trK]       -> 0,
     dot[At[la,lb]] -> 0,
     dot[Xt[ua]]    -> 0,
+    dot[Theta]     -> 0,
     dot[alpha]     -> 0,
     dot[A]         -> 0,
     dot[beta[ua]]  -> 0,
@@ -862,10 +893,10 @@ initRHSCalc =
   {
     dot[phi]       -> 0,
     dot[gt[la,lb]] -> 0,
-    dot[Tet]       -> 0,
     dot[trK]       -> 0,
     dot[At[la,lb]] -> 0,
     dot[Xt[ua]]    -> 0,
+    dot[Theta]     -> 0,
     dot[alpha]     -> 0,
     dot[A]         -> 0,
     dot[beta[ua]]  -> 0,
@@ -889,7 +920,7 @@ RHSRadiativeBoundaryCalc =
     
     detgt      -> 1 (* detgtExpr *),
     gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
-    em4phi     -> IfThen [conformalMethod, phi^2, Exp[-4 phi]],
+    em4phi     -> IfThen[conformalMethod==CMBSSNW, phi^2, Exp[-4 phi]],
     gu[ua,ub]  -> em4phi gtu[ua,ub],
     
     nn[la] -> Euc[la,lb] normal[ub],
@@ -902,10 +933,10 @@ RHSRadiativeBoundaryCalc =
     
     dot[phi]       -> - vg su[uc] PDo[phi      ,lc],
     dot[gt[la,lb]] -> -    su[uc] PDo[gt[la,lb],lc],
-    dot[Tet]       -> - vg su[uc] PDo[Tet      ,lc],
     dot[trK]       -> - vg su[uc] PDo[trK      ,lc],
     dot[At[la,lb]] -> -    su[uc] PDo[At[la,lb],lc],
     dot[Xt[ua]]    -> -    su[uc] PDo[Xt[ua]   ,lc],
+    dot[Theta]     -> - vg su[uc] PDo[Theta    ,lc],
     dot[alpha]     -> - vg su[uc] PDo[alpha    ,lc],
     dot[A]         -> - vg su[uc] PDo[A        ,lc],
     dot[beta[ua]]  -> -    su[uc] PDo[beta[ua] ,lc],
@@ -952,12 +983,12 @@ boundaryCalc =
   Where -> BoundaryWithGhosts,
   Equations -> 
   {
-    phi       -> IfThen [conformalMethod, 1, 0],
+    phi       -> IfThen[conformalMethod==CMBSSNW, 1, 0],
     gt[la,lb] -> KD[la,lb],
-    Tet       -> 0,
     trK       -> 0,
     At[la,lb] -> 0,
     Xt[ua]    -> 0,
+    Theta     -> 0,
     alpha     -> 1,
     A         -> 0,
     beta[ua]  -> 0,
@@ -975,13 +1006,13 @@ constraintsCalc =
   Schedule -> Automatic,
   After -> "MoL_PostStep",
   Where -> Interior,
-  Shorthands -> {detgt, ddetgt[la], gtu[ua,ub],
+  Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], Z[ua],
                  Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua],
                  e4phi, em4phi,
                  g[la,lb], detg, gu[ua,ub], ddetg[la], G[ua,lb,lc],
                  Rt[la,lb], Rphi[la,lb], R[la,lb], trR, Atm[ua,lb],
                  gK[la,lb,lc], cdphi[la], cdphi2[la,lb],
-                 rho, S[la], fac1, fac2, Zet[ua]},
+                 rho, S[la], fac1, fac2},
   Equations -> 
   {
     detgt        -> 1 (* detgtExpr *),
@@ -998,17 +1029,19 @@ constraintsCalc =
        used instead of Xt where no derivatives of Xt are taken *)
     Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk],
 
-    e4phi       -> IfThen [conformalMethod, 1/phi^2, Exp[4 phi]],
+    e4phi       -> IfThen[conformalMethod==CMBSSNW, 1/phi^2, Exp[4 phi]],
     em4phi      -> 1 / e4phi,
     g[la,lb]    -> e4phi gt[la,lb],
     detg        -> e4phi^3,
     gu[ua,ub]   -> em4phi gtu[ua,ub],
 
     (* The Z quantities *)
-    Zet[ud] -> (1/2) gu[ua,ud] (-PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc]),
-
+    Z[ud] -> IfThen[conformalMethod==CMCCZ4,
+                    (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc]
+                                     + gt[la,lc] Xt[uc]),
+                    0],
+    
     (* PRD 62, 044034 (2000), eqn. (18) *)
-    (* Added Zet term by changing Xtn to Xt *)
     Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm]
                  + (1/2) gt[lk,li] PD[Xt[uk],lj]
                  + (1/2) gt[lk,lj] PD[Xt[uk],li]
@@ -1040,9 +1073,9 @@ constraintsCalc =
                                    - PD[gt[la,lb],l1,l2] + PD[gt[l2,lb],l1,la]),
 *)
 
-    fac1 -> IfThen [conformalMethod, -1/(2 phi), 1],
+    fac1 -> IfThen[conformalMethod==CMBSSNW, -1/(2 phi), 1],
     cdphi[la] -> fac1 CDt[phi,la],
-    fac2 -> IfThen [conformalMethod, 1/(2 phi^2), 0],
+    fac2 -> IfThen[conformalMethod==CMBSSNW, 1/(2 phi^2), 0],
     cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb],
 
     (* PRD 62, 044034 (2000), eqn. (15) *)
@@ -1058,13 +1091,13 @@ constraintsCalc =
                    + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb]
                                  - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]),
     
-    (* Added Zet terms *)
-    R[la,lb] -> Rt[la,lb] + Rphi[la,lb] +
-                  (2/phi) (g[la,lc] Zet[uc] PD[phi,lb] +
-                  g[lb,lc] Zet[uc] PD[phi,la] -
-                  g[la,lb] Zet[uc] PD[phi,lc]) +
-                  e4phi Zet[uc] PD[gt[la,lb],lc],
-
+    R[la,lb] -> + Rt[la,lb] + Rphi[la,lb]
+                + IfThen[conformalMethod==CMCCZ4,
+                         + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb]
+                                    + g[lb,lc] Z[uc] PD[phi,la]
+                                    - g[la,lb] Z[uc] PD[phi,lc])
+                         + e4phi Z[uc] PD[gt[la,lb],lc],
+                         0],
     trR      -> gu[ua,ub] R[la,lb],
     
     (* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *)
@@ -1227,6 +1260,19 @@ keywordParameters =
 intParameters =
 {
   {
+    Name -> conformalMethod,
+    Description -> "Treatment of conformal factor",
+    AllowedValues -> {{Value -> "0", Description -> "BSSN phi method"},
+                      {Value -> "1", Description -> "BSSN W method"},
+                      {Value -> "2", Description -> "CCZ4 method"}},
+    Default -> 0
+  },
+  {
+    Name -> GammaShift,
+    Description -> "Covariant shift term in Gamma",
+    Default -> 0
+  },
+  {
     Name -> harmonicN,
     Description -> "d/dt alpha = - f alpha^n K  (harmonic=2, 1+log=1)",
     Default -> 2
@@ -1236,29 +1282,32 @@ intParameters =
     Default -> 0
   },
   {
-    Name -> conformalMethod,
-    Description -> "Treatment of conformal factor",
-    AllowedValues -> {{Value -> "0", Description -> "phi method"},
-                      {Value -> "1", Description -> "W method"}},
+    Name -> harmonicShift,
+    Description -> "Whether to use the harmonic shift",
+    AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"},
+                      {Value -> "1", Description -> "Harmonic shift"}},
     Default -> 0
   },
   {
     Name -> fdOrder,
     Default -> derivOrder,
     AllowedValues -> {2,4,6,8}
-  },
-  {
-    Name -> harmonicShift,
-    Description -> "Whether to use the harmonic shift",
-    AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"},
-                      {Value -> "1", Description -> "Harmonic shift"}},
-    Default -> 0
   }
 };
 
 realParameters =
 {
   {
+    Name -> dampk1,
+    Description -> "CCZ4 damping term 1 for Theta and Z",
+    Default -> 0
+  },
+  {
+    Name -> dampk2,
+    Description -> "CCZ4 damping term 2 for Theta and Z",
+    Default -> 0
+  },
+  {
     Name -> LapseACoeff,
     Description -> "Whether to evolve A in time",
     Default -> 0
@@ -1317,26 +1366,6 @@ realParameters =
     Description -> "Dissipation strength",
     AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}},
     Default -> 0
-  },
-  {
-    Name -> dampk1,
-    Description -> "Damping term 1 in Tet and Zet",
-    Default -> 0
-  },
-  {
-    Name -> dampk2,
-    Description -> "Damping term 2 in Tet and Zet",
-    Default -> 0
-  },
-  {
-    Name -> gamashift,
-    Description -> "Covariant shift term in Gamma",
-    Default -> 0
-  },
-  {
-    Name -> dampA,
-    Description -> "Damping term trA",
-    Default -> 0
   }
 };