[PATCH 2/4] Separate calculation into a part for psi4

Separate calculation into a part for psi4, a part for the other psis
and a part for the invariants. Add parameter for enabling each part.


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

1 files changed, 201 insertions, 181 deletions

diff --git a/m/WeylScal4.m b/m/WeylScal4.m
index 171a38e..c6f8c91 100644
--- a/m/WeylScal4.m
+++ b/m/WeylScal4.m
@@ -105,206 +105,213 @@ g11=gxx; g21=gxy; g22=gyy; g31=gxz; g32=gyz; g33=gzz;
 
 realParameters = {{Name -> offset, Default -> 10^(-15)},xorig,yorig,zorig};
 
-PsisCalc[fdOrder_, PD_] := 
-{
-  Name -> "psis_calc_" <> fdOrder,
-  Where -> Interior,
-  After -> "ADMBase_SetADMVars",
-  ConditionalOnKeyword -> {"fd_order", fdOrder},
-  Shorthands -> shorthands,
-  Equations ->
-  {
-        detg -> gDet,
-        invdetg -> 1 / detg,
-        gInv[ua,ub] -> invdetg gDet MatrixInverse[g[ua,ub]],
-        gamma[ua, lb, lc] -> 1/2 gInv[ua,ud] (PD[g[lb,ld], lc] + 
-          PD[g[lc,ld], lb] - PD[g[lb,lc],ld]),
-
-(****************************************************************************
-  Offset the origin
- ****************************************************************************)
-
-	xmoved -> x - xorig,
-	ymoved -> y - yorig,
-	zmoved -> z - zorig,
-
-(****************************************************************************
-  Compute the local tetrad
- ****************************************************************************)
-
-        (* azmuthal *)
-	va1 -> -ymoved, va2 -> xmoved+offset, va3 -> 0,
-
-        (* radial *)
-	vb1 -> xmoved+offset, vb2 -> ymoved, vb3 -> zmoved,
-
-        (* polar *)
-	vc[ua] -> Sqrt[detg] gInv[ua,ud] Eps[ld,lb,lc] va[ub] vb[uc],
-
-	(* Orthonormalize using Gram-Schmidt*)
-
-        (* Orthonormalize in the order phi, r, theta *)
-        wa[ua] -> va[ua],
-        omega11 -> wa[ua] wa[ub] g[la,lb],
-        ea[ua] -> wa[ua] / Sqrt[omega11],
-
-        omega12 -> ea[ua] vb[ub] g[la,lb],
-        wb[ua] -> vb[ua] - omega12 ea[ua],
-        omega22 -> wb[ua] wb[ub] g[la,lb],
-        eb[ua] -> wb[ua]/Sqrt[omega22],
-
-        omega13 -> ea[ua] vc[ub] g[la,lb],
-        omega23 -> eb[ua] vc[ub] g[la,lb],
-        wc[ua] -> vc[ua] - omega13 ea[ua] - omega23 eb[ua],
-        omega33 -> wc[ua] wc[ub] g[la,lb],
-        ec[ua] -> wc[ua]/Sqrt[omega33],
-
-	(* Create Spatial Portion of Null Tetrad *)
-	isqrt2  ->  0.7071067811865475244,
-	ltet[ua] -> isqrt2 eb[ua],
-	n[ua] -> - isqrt2 eb[ua],
-	rm[ua] -> isqrt2 ec[ua],
-	im[ua] -> isqrt2 ea[ua],
-	rmbar[ua] -> isqrt2 ec[ua],
-	imbar[ua] -> -isqrt2 ea[ua],
-
-	(* nn here is the projection of both l^a and n^a with u^a (the time-like unit
-	   vector normal to the hypersurface).  We do NOT save the t component of the 
-	   tetrads in this code to avoid unnecessary factors of lapse and shift. *)
-	nn -> isqrt2,
-
-
-(****************************************************************************
-  Compute the NP pseudoscalars
- ****************************************************************************)
-
-	(* Calculate the relevant Riemann Quantities *)
+psi4Eqs[PD_] := Flatten@{
+  detg -> gDet,
+  invdetg -> 1 / detg,
+  gInv[ua,ub] -> invdetg gDet MatrixInverse[g[ua,ub]],
+  gamma[ua, lb, lc] -> 1/2 gInv[ua,ud] (PD[g[lb,ld], lc] + PD[g[lc,ld], lb] - PD[g[lb,lc],ld]),
+
+  (****************************************************************************
+    Offset the origin
+   ****************************************************************************)
+
+  xmoved -> x - xorig,
+  ymoved -> y - yorig,
+  zmoved -> z - zorig,
+
+  (****************************************************************************
+    Compute the local tetrad
+   ****************************************************************************)
+  (* azmuthal *)
+  va1 -> -ymoved, va2 -> xmoved+offset, va3 -> 0,
+
+  (* radial *)
+  vb1 -> xmoved+offset, vb2 -> ymoved, vb3 -> zmoved,
+
+  (* polar *)
+  vc[ua] -> Sqrt[detg] gInv[ua,ud] Eps[ld,lb,lc] va[ub] vb[uc],
+
+  (* Orthonormalize using Gram-Schmidt*)
+  (* Orthonormalize in the order phi, r, theta *)
+  wa[ua] -> va[ua],
+  omega11 -> wa[ua] wa[ub] g[la,lb],
+  ea[ua] -> wa[ua] / Sqrt[omega11],
+
+  omega12 -> ea[ua] vb[ub] g[la,lb],
+  wb[ua] -> vb[ua] - omega12 ea[ua],
+  omega22 -> wb[ua] wb[ub] g[la,lb],
+  eb[ua] -> wb[ua]/Sqrt[omega22],
+
+  omega13 -> ea[ua] vc[ub] g[la,lb],
+  omega23 -> eb[ua] vc[ub] g[la,lb],
+  wc[ua] -> vc[ua] - omega13 ea[ua] - omega23 eb[ua],
+  omega33 -> wc[ua] wc[ub] g[la,lb],
+  ec[ua] -> wc[ua]/Sqrt[omega33],
+
+  (* Create Spatial Portion of Null Tetrad *)
+  isqrt2  ->  0.7071067811865475244,
+  ltet[ua] -> isqrt2 eb[ua],
+  n[ua] -> - isqrt2 eb[ua],
+  rm[ua] -> isqrt2 ec[ua],
+  im[ua] -> isqrt2 ea[ua],
+  rmbar[ua] -> isqrt2 ec[ua],
+  imbar[ua] -> -isqrt2 ea[ua],
+
+  (* nn here is the projection of both l^a and n^a with u^a (the time-like unit
+     vector normal to the hypersurface).  We do NOT save the t component of the
+     tetrads in this code to avoid unnecessary factors of lapse and shift. *)
+  nn -> isqrt2,
+
+
+  (****************************************************************************
+    Compute the NP pseudoscalars
+   ****************************************************************************)
+
+  (* Calculate the relevant Riemann Quantities *)
 		
-        (* The 3-Riemann *)                
-	R[la,lb,lc,ld] -> 1/2 ( PD[g[la,ld],lc,lb] + PD[g[lb,lc],ld,la] )
-			- 1/2 ( PD[g[la,lc],lb,ld] + PD[g[lb,ld],la,lc] )
-			+ g[lj,le] gamma[uj,lb,lc] gamma[ue,la,ld] 	
-			- g[lj,le] gamma[uj,lb,ld] gamma[ue,la,lc],
-
-        (* The 4-Riemann projected into the slice on all its indices. 
-           The Gauss equation. *)
-        R4p[li,lj,lk,ll] -> R[li,lj,lk,ll] + 
-          2 AntiSymmetrize[k[li,lk] k[ll,lj], lk, ll],
-
-        (* The 4-Riemann projected in the unit normal direction on one 
-           index, then into the slice on the remaining indices. The Codazzi
-           equation. *)
-	Ro[lj,lk,ll] ->  - 2 AntiSymmetrize[ PD[k[lj,lk],ll], lk,ll]
-			 - 2 AntiSymmetrize[ gamma[up,lj,lk] k[ll,lp], lk,ll],
-
-        (* The 4-Riemann projected in the unit normal direction on two
-           indices, and into the slice on the remaining two. *)
-	Rojo[lj,ll] ->  gInv[uc,ud] (R[lj,lc,ll,ld] )   
-                        - k[lj,lp] gInv[up,ud] k[ld,ll]
-			+ k[lc,ld] gInv[uc,ud] k[lj,ll],
-
-	(* Calculate End Quantities
-		NOTE: In writing this, I assume m[0]=0!! to save lots of work *)
-
-	Psi4r -> R4p[li,lj,lk,ll] n[ui] n[uk] 
-                   ( rmbar[uj] rmbar[ul] - imbar[uj] imbar[ul] )
-		 + 2 Ro[lj,lk,ll] n[uk] nn 
-                   ( rmbar[uj] rmbar[ul] - imbar[uj] imbar[ul] )
+  (* The 3-Riemann *)
+  R[la,lb,lc,ld] -> 1/2 ( PD[g[la,ld],lc,lb] + PD[g[lb,lc],ld,la] )
+			      - 1/2 ( PD[g[la,lc],lb,ld] + PD[g[lb,ld],la,lc] )
+			      + g[lj,le] gamma[uj,lb,lc] gamma[ue,la,ld]
+			      - g[lj,le] gamma[uj,lb,ld] gamma[ue,la,lc],
+
+  (* The 4-Riemann projected into the slice on all its indices. The Gauss equation. *)
+  R4p[li,lj,lk,ll] -> R[li,lj,lk,ll] + 2 AntiSymmetrize[k[li,lk] k[ll,lj], lk, ll],
+
+  (* The 4-Riemann projected in the unit normal direction on one
+     index, then into the slice on the remaining indices. The Codazzi equation. *)
+  Ro[lj,lk,ll] -> - 2 AntiSymmetrize[ PD[k[lj,lk],ll], lk,ll]
+			      - 2 AntiSymmetrize[ gamma[up,lj,lk] k[ll,lp], lk,ll],
+
+  (* The 4-Riemann projected in the unit normal direction on two
+     indices, and into the slice on the remaining two. *)
+  Rojo[lj,ll] ->  gInv[uc,ud] (R[lj,lc,ll,ld] ) - k[lj,lp] gInv[up,ud] k[ld,ll]
+			      + k[lc,ld] gInv[uc,ud] k[lj,ll],
+
+  (* Calculate End Quantities
+     NOTE: In writing this, I assume m[0]=0!! to save lots of work *)
+
+  Psi4r -> R4p[li,lj,lk,ll] n[ui] n[uk] ( rmbar[uj] rmbar[ul] - imbar[uj] imbar[ul] )
+		 + 2 Ro[lj,lk,ll] n[uk] nn ( rmbar[uj] rmbar[ul] - imbar[uj] imbar[ul] )
 		 + Rojo[lj,ll] nn nn ( rmbar[uj] rmbar[ul] - imbar[uj] imbar[ul] 
-              (* + terms in mbar^0 == 0*) ),
+         (* + terms in mbar^0 == 0 *) ),
 
-	Psi4i -> R4p[la,lb,lc,ld] n[ua] n[uc] ( - rm[ub] im[ud] - im[ub] rm[ud] )
+  Psi4i -> R4p[la,lb,lc,ld] n[ua] n[uc] ( - rm[ub] im[ud] - im[ub] rm[ud] )
 		 + 2 Ro[la,lb,lc] n[ub] nn ( - rm[ua] im[uc] - im[ua] rm[uc] )
-		 + Rojo[la,lb] nn nn ( - rm[ua] im[ub] - im[ua] rm[ub] )  ,
-
+		 + Rojo[la,lb] nn nn ( - rm[ua] im[ub] - im[ua] rm[ub] )
+};
 
-	Psi3r -> R4p[la,lb,lc,ld] ltet[ua] n[ub] rm[uc] n[ud]
-		 + Ro[la,lb,lc] ( nn (n[ua]-ltet[ua]) rm[ub] n[uc]
-                                      - nn rm[ua] ltet[ub] n[uc] )
+otherPsiEqs = {
+  Psi3r -> R4p[la,lb,lc,ld] ltet[ua] n[ub] rm[uc] n[ud]
+		 + Ro[la,lb,lc] ( nn (n[ua]-ltet[ua]) rm[ub] n[uc] - nn rm[ua] ltet[ub] n[uc] )
 		 - Rojo[la,lb] nn (n[ua]-ltet[ua]) nn rm[ub],
 
-	Psi3i -> - R4p[la,lb,lc,ld] ltet[ua] n[ub] im[uc] n[ud]
-		 - Ro[la,lb,lc] ( nn (n[ua]-ltet[ua]) im[ub] n[uc] 
-                                  - nn im[ua] ltet[ub] n[uc] )
+  Psi3i -> - R4p[la,lb,lc,ld] ltet[ua] n[ub] im[uc] n[ud]
+		 - Ro[la,lb,lc] ( nn (n[ua]-ltet[ua]) im[ub] n[uc] - nn im[ua] ltet[ub] n[uc] )
 		 + Rojo[la,lb] nn (n[ua]-ltet[ua]) nn im[ub],
 
-	Psi2r -> R4p[la,lb,lc,ld] ltet[ua] n[ud] (rm[ub] rm[uc] + im[ub] im[uc])
-		 + Ro[la,lb,lc] nn ( n[uc] (rm[ua] rm[ub] + im[ua] im[ub]) 
-                                     - ltet[ub] (rm[ua] rm[uc] + im[ua] im[uc]) )
+  Psi2r -> R4p[la,lb,lc,ld] ltet[ua] n[ud] (rm[ub] rm[uc] + im[ub] im[uc])
+		 + Ro[la,lb,lc] nn ( n[uc] (rm[ua] rm[ub] + im[ua] im[ub]) - ltet[ub] (rm[ua] rm[uc] + im[ua] im[uc]) )
 		 - Rojo[la,lb] nn nn (rm[ua] rm[ub] + im[ua] im[ub]),
 
-	Psi2i -> R4p[la,lb,lc,ld] ltet[ua] n[ud] (im[ub] rm[uc] - rm[ub] im[uc])
-		 + Ro[la,lb,lc] nn ( n[uc] (im[ua] rm[ub] - rm[ua] im[ub]) 
-                                     - ltet[ub] (rm[ua] im[uc] - im[ua] rm[uc]) )
+  Psi2i -> R4p[la,lb,lc,ld] ltet[ua] n[ud] (im[ub] rm[uc] - rm[ub] im[uc])
+		 + Ro[la,lb,lc] nn ( n[uc] (im[ua] rm[ub] - rm[ua] im[ub]) - ltet[ub] (rm[ua] im[uc] - im[ua] rm[uc]) )
 		 - Rojo[la,lb] nn nn (im[ua] rm[ub] - rm[ua] im[ub]),
 
-	Psi1r -> R4p[la,lb,lc,ld] n[ua] ltet[ub] rm[uc] ltet[ud]
-		 + Ro[la,lb,lc] ( nn ltet[ua] rm[ub] ltet[uc] 
-                                  - nn rm[ua] n[ub] ltet[uc] 
-                                  - nn n[ua] rm[ub] ltet[uc] )
+  Psi1r -> R4p[la,lb,lc,ld] n[ua] ltet[ub] rm[uc] ltet[ud]
+		 + Ro[la,lb,lc] ( nn ltet[ua] rm[ub] ltet[uc] - nn rm[ua] n[ub] ltet[uc] - nn n[ua] rm[ub] ltet[uc] )
 		 + Rojo[la,lb] nn nn ( n[ua] rm[ub] - ltet[ua] rm[ub] ),
 
-	Psi1i -> R4p[la,lb,lc,ld] n[ua] ltet[ub] im[uc] ltet[ud]
-		 + Ro[la,lb,lc] ( nn ltet[ua] im[ub] ltet[uc] 
-                                  - nn im[ua] n[ub] ltet[uc] 
-                                  - nn n[ua] im[ub] ltet[uc] )
+  Psi1i -> R4p[la,lb,lc,ld] n[ua] ltet[ub] im[uc] ltet[ud]
+		 + Ro[la,lb,lc] ( nn ltet[ua] im[ub] ltet[uc] - nn im[ua] n[ub] ltet[uc] - nn n[ua] im[ub] ltet[uc] )
 		 + Rojo[la,lb] nn nn ( n[ua] im[ub] - ltet[ua] im[ub] ),
 
-	Psi0r -> R4p[la,lb,lc,ld] ltet[ua] ltet[uc] (rm[ub] rm[ud] - im[ub] im[ud])
+  Psi0r -> R4p[la,lb,lc,ld] ltet[ua] ltet[uc] (rm[ub] rm[ud] - im[ub] im[ud])
 		 + 2 Ro[la,lb,lc] nn ltet[ub] (rm[ua] rm[uc] - im[ua] im[uc])
 		 + Rojo[la,lb] nn nn (rm[ua] rm[ub] - im[ua] im[ub]),
 
-	Psi0i -> R4p[la,lb,lc,ld] ltet[ua] ltet[uc] (rm[ub] im[ud] + im[ub] rm[ud])
+  Psi0i -> R4p[la,lb,lc,ld] ltet[ua] ltet[uc] (rm[ub] im[ud] + im[ub] rm[ud])
 		 + 2 Ro[la,lb,lc] nn ltet[ub] (rm[ua] im[uc] + im[ua] rm[uc])
-		 + Rojo[la,lb] nn nn (rm[ua] im[ub] + im[ua] rm[ub]),
-
-        (* Scalar invariants I and J as defined in (2.2a) and (2.2b) of arXiv:gr-qc/0407013 *)
-        curvIr -> ComplexExpand[Re[3 (Psi2r+I Psi2i)^2 - 4 (Psi1r+I Psi1i) (Psi3r + I Psi3i) + (Psi4r + I Psi4i) (Psi0r + I Psi0i)]],
-        curvIi -> ComplexExpand[Im[3 (Psi2r+I Psi2i)^2 - 4 (Psi1r+I Psi1i) (Psi3r + I Psi3i) + (Psi4r + I Psi4i) (Psi0r + I Psi0i)]],
-        curvJr -> ComplexExpand[Re[Det[{{Psi4r+I Psi4i,Psi3r+I Psi3i,Psi2r+I Psi2i},
-                                        {Psi3r+I Psi3i,Psi2r+I Psi2i,Psi1r+I Psi1i},
-                                        {Psi2r+I Psi2i,Psi1r+I Psi1i,Psi0r+I Psi0i}}]]],
-        curvJi -> ComplexExpand[Im[Det[{{Psi4r+I Psi4i,Psi3r+I Psi3i,Psi2r+I Psi2i},
-                                        {Psi3r+I Psi3i,Psi2r+I Psi2i,Psi1r+I Psi1i},
-                                        {Psi2r+I Psi2i,Psi1r+I Psi1i,Psi0r+I Psi0i}}]]],
-  
-        (* Scalar invariants J1, J2, J3 and J4 of the Narlikar and Karmarkar basis as defined in B5-B8 of arXiv:0704.1756.
-           Computed from Weyl tensor expressions using xTensor. *)
-        curvJ1 -> -16(3 Psi2i^2-3 Psi2r^2-4 Psi1i Psi3i+4 Psi1r Psi3r+Psi0i Psi4i-Psi0r Psi4r),
-        curvJ2 -> 96(-3 Psi2i^2 Psi2r+Psi2r^3+2 Psi1r Psi2i Psi3i+2 Psi1i Psi2r Psi3i-Psi0r Psi3i^2+2 Psi1i Psi2i Psi3r-2 Psi1r Psi2r Psi3r
-          -2 Psi0i Psi3i Psi3r+Psi0r Psi3r^2-2 Psi1i Psi1r Psi4i+Psi0r Psi2i Psi4i+Psi0i Psi2r Psi4i-Psi1i^2 Psi4r+Psi1r^2 Psi4r
-          +Psi0i Psi2i Psi4r-Psi0r Psi2r Psi4r),
-        curvJ3 -> 64(9 Psi2i^4-54 Psi2i^2 Psi2r^2+9 Psi2r^4-24 Psi1i Psi2i^2 Psi3i+48 Psi1r Psi2i Psi2r Psi3i+24 Psi1i Psi2r^2 Psi3i
-          +16 Psi1i^2 Psi3i^2-16 Psi1r^2 Psi3i^2+24 Psi1r Psi2i^2 Psi3r+48 Psi1i Psi2i Psi2r Psi3r-24 Psi1r Psi2r^2 Psi3r
-          -64 Psi1i Psi1r Psi3i Psi3r-16 Psi1i^2 Psi3r^2+16 Psi1r^2 Psi3r^2+6 Psi0i Psi2i^2 Psi4i-12 Psi0r Psi2i Psi2r Psi4i
-          -6 Psi0i Psi2r^2 Psi4i-8 Psi0i Psi1i Psi3i Psi4i+8 Psi0r Psi1r Psi3i Psi4i+8 Psi0r Psi1i Psi3r Psi4i
-          +8 Psi0i Psi1r Psi3r Psi4i+Psi0i^2 Psi4i^2-Psi0r^2 Psi4i^2-6 Psi0r Psi2i^2 Psi4r-12 Psi0i Psi2i Psi2r Psi4r+6 Psi0r Psi2r^2 Psi4r
-          +8 Psi0r Psi1i Psi3i Psi4r+8 Psi0i Psi1r Psi3i Psi4r+8 Psi0i Psi1i Psi3r Psi4r-8 Psi0r Psi1r Psi3r Psi4r-4 Psi0i Psi0r Psi4i Psi4r-Psi0i^2 Psi4r^2+Psi0r^2 Psi4r^2),
-        curvJ4 -> -640(-15 Psi2i^4 Psi2r+30 Psi2i^2 Psi2r^3-3 Psi2r^5+10 Psi1r Psi2i^3 Psi3i+30 Psi1i Psi2i^2 Psi2r Psi3i-30 Psi1r Psi2i Psi2r^2 Psi3i
-          -10 Psi1i Psi2r^3 Psi3i-16 Psi1i Psi1r Psi2i Psi3i^2-3 Psi0r Psi2i^2 Psi3i^2-8 Psi1i^2 Psi2r Psi3i^2+8 Psi1r^2 Psi2r Psi3i^2
-          -6 Psi0i Psi2i Psi2r Psi3i^2+3 Psi0r Psi2r^2 Psi3i^2+4 Psi0r Psi1i Psi3i^3+4 Psi0i Psi1r Psi3i^3+10 Psi1i Psi2i^3 Psi3r
-          -30 Psi1r Psi2i^2 Psi2r Psi3r-30 Psi1i Psi2i Psi2r^2 Psi3r+10 Psi1r Psi2r^3 Psi3r-16 Psi1i^2 Psi2i Psi3i Psi3r
-          +16 Psi1r^2 Psi2i Psi3i Psi3r-6 Psi0i Psi2i^2 Psi3i Psi3r+32 Psi1i Psi1r Psi2r Psi3i Psi3r+12 Psi0r Psi2i Psi2r Psi3i Psi3r
-          +6 Psi0i Psi2r^2 Psi3i Psi3r+12 Psi0i Psi1i Psi3i^2 Psi3r-12 Psi0r Psi1r Psi3i^2 Psi3r+16 Psi1i Psi1r Psi2i Psi3r^2
-          +3 Psi0r Psi2i^2 Psi3r^2+8 Psi1i^2 Psi2r Psi3r^2-8 Psi1r^2 Psi2r Psi3r^2+6 Psi0i Psi2i Psi2r Psi3r^2-3 Psi0r Psi2r^2 Psi3r^2
-          -12 Psi0r Psi1i Psi3i Psi3r^2-12 Psi0i Psi1r Psi3i Psi3r^2-4 Psi0i Psi1i Psi3r^3+4 Psi0r Psi1r Psi3r^3-6 Psi1i Psi1r Psi2i^2 Psi4i
-          +2 Psi0r Psi2i^3 Psi4i-6 Psi1i^2 Psi2i Psi2r Psi4i+6 Psi1r^2 Psi2i Psi2r Psi4i+6 Psi0i Psi2i^2 Psi2r Psi4i
-          +6 Psi1i Psi1r Psi2r^2 Psi4i-6 Psi0r Psi2i Psi2r^2 Psi4i-2 Psi0i Psi2r^3 Psi4i+12 Psi1i^2 Psi1r Psi3i Psi4i-4 Psi1r^3 Psi3i Psi4i
-          -2 Psi0r Psi1i Psi2i Psi3i Psi4i-2 Psi0i Psi1r Psi2i Psi3i Psi4i-2 Psi0i Psi1i Psi2r Psi3i Psi4i
-          +2 Psi0r Psi1r Psi2r Psi3i Psi4i-2 Psi0i Psi0r Psi3i^2 Psi4i+4 Psi1i^3 Psi3r Psi4i-12 Psi1i Psi1r^2 Psi3r Psi4i
-          -2 Psi0i Psi1i Psi2i Psi3r Psi4i+2 Psi0r Psi1r Psi2i Psi3r Psi4i+2 Psi0r Psi1i Psi2r Psi3r Psi4i
-          +2 Psi0i Psi1r Psi2r Psi3r Psi4i-2 Psi0i^2 Psi3i Psi3r Psi4i+2 Psi0r^2 Psi3i Psi3r Psi4i+2 Psi0i Psi0r Psi3r^2 Psi4i
-          -Psi0r Psi1i^2 Psi4i^2-2 Psi0i Psi1i Psi1r Psi4i^2+Psi0r Psi1r^2 Psi4i^2+2 Psi0i Psi0r Psi2i Psi4i^2+Psi0i^2 Psi2r Psi4i^2
-          -Psi0r^2 Psi2r Psi4i^2-3 Psi1i^2 Psi2i^2 Psi4r+3 Psi1r^2 Psi2i^2 Psi4r+2 Psi0i Psi2i^3 Psi4r+12 Psi1i Psi1r Psi2i Psi2r Psi4r
-          -6 Psi0r Psi2i^2 Psi2r Psi4r+3 Psi1i^2 Psi2r^2 Psi4r-3 Psi1r^2 Psi2r^2 Psi4r-6 Psi0i Psi2i Psi2r^2 Psi4r+2 Psi0r Psi2r^3 Psi4r
-          +4 Psi1i^3 Psi3i Psi4r-12 Psi1i Psi1r^2 Psi3i Psi4r-2 Psi0i Psi1i Psi2i Psi3i Psi4r+2 Psi0r Psi1r Psi2i Psi3i Psi4r
-          +2 Psi0r Psi1i Psi2r Psi3i Psi4r+2 Psi0i Psi1r Psi2r Psi3i Psi4r-Psi0i^2 Psi3i^2 Psi4r+Psi0r^2 Psi3i^2 Psi4r
-          -12 Psi1i^2 Psi1r Psi3r Psi4r+4 Psi1r^3 Psi3r Psi4r+2 Psi0r Psi1i Psi2i Psi3r Psi4r+2 Psi0i Psi1r Psi2i Psi3r Psi4r
-          +2 Psi0i Psi1i Psi2r Psi3r Psi4r-2 Psi0r Psi1r Psi2r Psi3r Psi4r+4 Psi0i Psi0r Psi3i Psi3r Psi4r+Psi0i^2 Psi3r^2 Psi4r-Psi0r^2 Psi3r^2 Psi4r
-          -2 Psi0i Psi1i^2 Psi4i Psi4r+4 Psi0r Psi1i Psi1r Psi4i Psi4r+2 Psi0i Psi1r^2 Psi4i Psi4r+2 Psi0i^2 Psi2i Psi4i Psi4r
-          -2 Psi0r^2 Psi2i Psi4i Psi4r-4 Psi0i Psi0r Psi2r Psi4i Psi4r+Psi0r Psi1i^2 Psi4r^2+2 Psi0i Psi1i Psi1r Psi4r^2-Psi0r Psi1r^2 Psi4r^2
-          -2 Psi0i Psi0r Psi2i Psi4r^2-Psi0i^2 Psi2r Psi4r^2+Psi0r^2 Psi2r Psi4r^2)
-  }
+		 + Rojo[la,lb] nn nn (rm[ua] im[ub] + im[ua] rm[ub])
+};
+
+invariantEqs := {
+  (* Scalar invariants I and J as defined in (2.2a) and (2.2b) of arXiv:gr-qc/0407013 *)
+  curvIr -> ComplexExpand[Re[3 (Psi2r+I Psi2i)^2 - 4 (Psi1r+I Psi1i) (Psi3r + I Psi3i) + (Psi4r + I Psi4i) (Psi0r + I Psi0i)]],
+  curvIi -> ComplexExpand[Im[3 (Psi2r+I Psi2i)^2 - 4 (Psi1r+I Psi1i) (Psi3r + I Psi3i) + (Psi4r + I Psi4i) (Psi0r + I Psi0i)]],
+  curvJr -> ComplexExpand[Re[Det[{{Psi4r+I Psi4i,Psi3r+I Psi3i,Psi2r+I Psi2i},
+                                  {Psi3r+I Psi3i,Psi2r+I Psi2i,Psi1r+I Psi1i},
+                                  {Psi2r+I Psi2i,Psi1r+I Psi1i,Psi0r+I Psi0i}}]]],
+  curvJi -> ComplexExpand[Im[Det[{{Psi4r+I Psi4i,Psi3r+I Psi3i,Psi2r+I Psi2i},
+                                  {Psi3r+I Psi3i,Psi2r+I Psi2i,Psi1r+I Psi1i},
+                                  {Psi2r+I Psi2i,Psi1r+I Psi1i,Psi0r+I Psi0i}}]]],
+
+  (* Scalar invariants J1, J2, J3 and J4 of the Narlikar and Karmarkar basis as defined
+     in B5-B8 of arXiv:0704.1756. Computed from Weyl tensor expressions using xTensor. *)
+  curvJ1 -> -16(3 Psi2i^2-3 Psi2r^2-4 Psi1i Psi3i+4 Psi1r Psi3r+Psi0i Psi4i-Psi0r Psi4r),
+  curvJ2 -> 96(-3 Psi2i^2 Psi2r+Psi2r^3+2 Psi1r Psi2i Psi3i+2 Psi1i Psi2r Psi3i-Psi0r Psi3i^2+2 Psi1i Psi2i Psi3r-2 Psi1r Psi2r Psi3r
+    -2 Psi0i Psi3i Psi3r+Psi0r Psi3r^2-2 Psi1i Psi1r Psi4i+Psi0r Psi2i Psi4i+Psi0i Psi2r Psi4i-Psi1i^2 Psi4r+Psi1r^2 Psi4r
+    +Psi0i Psi2i Psi4r-Psi0r Psi2r Psi4r),
+  curvJ3 -> 64(9 Psi2i^4-54 Psi2i^2 Psi2r^2+9 Psi2r^4-24 Psi1i Psi2i^2 Psi3i+48 Psi1r Psi2i Psi2r Psi3i+24 Psi1i Psi2r^2 Psi3i
+    +16 Psi1i^2 Psi3i^2-16 Psi1r^2 Psi3i^2+24 Psi1r Psi2i^2 Psi3r+48 Psi1i Psi2i Psi2r Psi3r-24 Psi1r Psi2r^2 Psi3r
+    -64 Psi1i Psi1r Psi3i Psi3r-16 Psi1i^2 Psi3r^2+16 Psi1r^2 Psi3r^2+6 Psi0i Psi2i^2 Psi4i-12 Psi0r Psi2i Psi2r Psi4i
+    -6 Psi0i Psi2r^2 Psi4i-8 Psi0i Psi1i Psi3i Psi4i+8 Psi0r Psi1r Psi3i Psi4i+8 Psi0r Psi1i Psi3r Psi4i
+    +8 Psi0i Psi1r Psi3r Psi4i+Psi0i^2 Psi4i^2-Psi0r^2 Psi4i^2-6 Psi0r Psi2i^2 Psi4r-12 Psi0i Psi2i Psi2r Psi4r+6 Psi0r Psi2r^2 Psi4r
+    +8 Psi0r Psi1i Psi3i Psi4r+8 Psi0i Psi1r Psi3i Psi4r+8 Psi0i Psi1i Psi3r Psi4r-8 Psi0r Psi1r Psi3r Psi4r-4 Psi0i Psi0r Psi4i Psi4r-Psi0i^2 Psi4r^2+Psi0r^2 Psi4r^2),
+  curvJ4 -> -640(-15 Psi2i^4 Psi2r+30 Psi2i^2 Psi2r^3-3 Psi2r^5+10 Psi1r Psi2i^3 Psi3i+30 Psi1i Psi2i^2 Psi2r Psi3i-30 Psi1r Psi2i Psi2r^2 Psi3i
+    -10 Psi1i Psi2r^3 Psi3i-16 Psi1i Psi1r Psi2i Psi3i^2-3 Psi0r Psi2i^2 Psi3i^2-8 Psi1i^2 Psi2r Psi3i^2+8 Psi1r^2 Psi2r Psi3i^2
+    -6 Psi0i Psi2i Psi2r Psi3i^2+3 Psi0r Psi2r^2 Psi3i^2+4 Psi0r Psi1i Psi3i^3+4 Psi0i Psi1r Psi3i^3+10 Psi1i Psi2i^3 Psi3r
+    -30 Psi1r Psi2i^2 Psi2r Psi3r-30 Psi1i Psi2i Psi2r^2 Psi3r+10 Psi1r Psi2r^3 Psi3r-16 Psi1i^2 Psi2i Psi3i Psi3r
+    +16 Psi1r^2 Psi2i Psi3i Psi3r-6 Psi0i Psi2i^2 Psi3i Psi3r+32 Psi1i Psi1r Psi2r Psi3i Psi3r+12 Psi0r Psi2i Psi2r Psi3i Psi3r
+    +6 Psi0i Psi2r^2 Psi3i Psi3r+12 Psi0i Psi1i Psi3i^2 Psi3r-12 Psi0r Psi1r Psi3i^2 Psi3r+16 Psi1i Psi1r Psi2i Psi3r^2
+    +3 Psi0r Psi2i^2 Psi3r^2+8 Psi1i^2 Psi2r Psi3r^2-8 Psi1r^2 Psi2r Psi3r^2+6 Psi0i Psi2i Psi2r Psi3r^2-3 Psi0r Psi2r^2 Psi3r^2
+    -12 Psi0r Psi1i Psi3i Psi3r^2-12 Psi0i Psi1r Psi3i Psi3r^2-4 Psi0i Psi1i Psi3r^3+4 Psi0r Psi1r Psi3r^3-6 Psi1i Psi1r Psi2i^2 Psi4i
+    +2 Psi0r Psi2i^3 Psi4i-6 Psi1i^2 Psi2i Psi2r Psi4i+6 Psi1r^2 Psi2i Psi2r Psi4i+6 Psi0i Psi2i^2 Psi2r Psi4i
+    +6 Psi1i Psi1r Psi2r^2 Psi4i-6 Psi0r Psi2i Psi2r^2 Psi4i-2 Psi0i Psi2r^3 Psi4i+12 Psi1i^2 Psi1r Psi3i Psi4i-4 Psi1r^3 Psi3i Psi4i
+    -2 Psi0r Psi1i Psi2i Psi3i Psi4i-2 Psi0i Psi1r Psi2i Psi3i Psi4i-2 Psi0i Psi1i Psi2r Psi3i Psi4i
+    +2 Psi0r Psi1r Psi2r Psi3i Psi4i-2 Psi0i Psi0r Psi3i^2 Psi4i+4 Psi1i^3 Psi3r Psi4i-12 Psi1i Psi1r^2 Psi3r Psi4i
+    -2 Psi0i Psi1i Psi2i Psi3r Psi4i+2 Psi0r Psi1r Psi2i Psi3r Psi4i+2 Psi0r Psi1i Psi2r Psi3r Psi4i
+    +2 Psi0i Psi1r Psi2r Psi3r Psi4i-2 Psi0i^2 Psi3i Psi3r Psi4i+2 Psi0r^2 Psi3i Psi3r Psi4i+2 Psi0i Psi0r Psi3r^2 Psi4i
+    -Psi0r Psi1i^2 Psi4i^2-2 Psi0i Psi1i Psi1r Psi4i^2+Psi0r Psi1r^2 Psi4i^2+2 Psi0i Psi0r Psi2i Psi4i^2+Psi0i^2 Psi2r Psi4i^2
+    -Psi0r^2 Psi2r Psi4i^2-3 Psi1i^2 Psi2i^2 Psi4r+3 Psi1r^2 Psi2i^2 Psi4r+2 Psi0i Psi2i^3 Psi4r+12 Psi1i Psi1r Psi2i Psi2r Psi4r
+    -6 Psi0r Psi2i^2 Psi2r Psi4r+3 Psi1i^2 Psi2r^2 Psi4r-3 Psi1r^2 Psi2r^2 Psi4r-6 Psi0i Psi2i Psi2r^2 Psi4r+2 Psi0r Psi2r^3 Psi4r
+    +4 Psi1i^3 Psi3i Psi4r-12 Psi1i Psi1r^2 Psi3i Psi4r-2 Psi0i Psi1i Psi2i Psi3i Psi4r+2 Psi0r Psi1r Psi2i Psi3i Psi4r
+    +2 Psi0r Psi1i Psi2r Psi3i Psi4r+2 Psi0i Psi1r Psi2r Psi3i Psi4r-Psi0i^2 Psi3i^2 Psi4r+Psi0r^2 Psi3i^2 Psi4r
+    -12 Psi1i^2 Psi1r Psi3r Psi4r+4 Psi1r^3 Psi3r Psi4r+2 Psi0r Psi1i Psi2i Psi3r Psi4r+2 Psi0i Psi1r Psi2i Psi3r Psi4r
+    +2 Psi0i Psi1i Psi2r Psi3r Psi4r-2 Psi0r Psi1r Psi2r Psi3r Psi4r+4 Psi0i Psi0r Psi3i Psi3r Psi4r+Psi0i^2 Psi3r^2 Psi4r-Psi0r^2 Psi3r^2 Psi4r
+    -2 Psi0i Psi1i^2 Psi4i Psi4r+4 Psi0r Psi1i Psi1r Psi4i Psi4r+2 Psi0i Psi1r^2 Psi4i Psi4r+2 Psi0i^2 Psi2i Psi4i Psi4r
+    -2 Psi0r^2 Psi2i Psi4i Psi4r-4 Psi0i Psi0r Psi2r Psi4i Psi4r+Psi0r Psi1i^2 Psi4r^2+2 Psi0i Psi1i Psi1r Psi4r^2-Psi0r Psi1r^2 Psi4r^2
+    -2 Psi0i Psi0r Psi2i Psi4r^2-Psi0i^2 Psi2r Psi4r^2+Psi0r^2 Psi2r Psi4r^2)
+};
+
+Psi4Calc[fdOrder_, PD_] :=
+{
+  Name -> "psi4_calc_" <> fdOrder,
+  Where -> Interior,
+  After -> "ADMBase_SetADMVars",
+  ConditionalOnKeywords -> {{"fd_order", fdOrder}, {"calc_scalars", "psi4"}},
+  Shorthands -> shorthands,
+  Equations -> psi4Eqs[PD]
+};
+
+PsisCalc[fdOrder_, PD_] :=
+{
+  Name -> "psis_calc_" <> fdOrder,
+  Where -> Interior,
+  After -> "ADMBase_SetADMVars",
+  ConditionalOnKeywords -> {{"fd_order", fdOrder}, {"calc_scalars", "psis"}},
+  Shorthands -> shorthands,
+  Equations -> Join[psi4Eqs[PD], otherPsiEqs]
+};
+
+InvariantsCalc[fdOrder_, PD_] :=
+{
+  Name -> "invars_calc_" <> fdOrder,
+  Where -> Interior,
+  After -> "ADMBase_SetADMVars",
+  ConditionalOnKeywords -> {{"fd_order", fdOrder}, {"calc_scalars", "psis_and_invariants"}},
+  Shorthands -> shorthands,
+  Equations -> Join[psi4Eqs[PD], otherPsiEqs, invariantEqs]
 };
 
 
@@ -319,9 +326,16 @@ fdOrderParam =
   AllowedValues -> {"Nth", "2nd", "4th"}
 };
 
+calcScalarsParam = {
+  Name -> "calc_scalars",
+  Description -> "Which scalars to calculate",
+  AllowedValues -> {"psi4", "psis", "psis_and_invariants"},
+  Default -> "psi4"
+};
+
 keywordParameters = 
 {
-  fdOrderParam
+  fdOrderParam, calcScalarsParam
 };
 
 intParameters =
@@ -336,9 +350,15 @@ intParameters =
 
 calculations = 
 {
+  Psi4Calc["Nth", PDstandard],
+  Psi4Calc["2nd", PDstandard2nd],
+  Psi4Calc["4th", PDstandard4th],
   PsisCalc["Nth", PDstandard],
   PsisCalc["2nd", PDstandard2nd],
-  PsisCalc["4th", PDstandard4th]
+  PsisCalc["4th", PDstandard4th],
+  InvariantsCalc["Nth", PDstandard],
+  InvariantsCalc["2nd", PDstandard2nd],
+  InvariantsCalc["4th", PDstandard4th]
 };
 
 CreateKrancThornTT[groups, ".", "WeylScal4",