Add some functionality:

Conversion between 4 and 3+1 dimensions
Check for nans (isnan for Fortran)
Gram-Schmidt orthonormalisation
Calculate Riemann, Ricci, and Weyl tensors in 2, 3, and 4 dimensions
Calculate time derivatives


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@21 b716e942-a2de-43ad-8f52-f3dfe468e4e7

9 files changed, 1003 insertions, 47 deletions

diff --git a/src/adm_metric.F90 b/src/adm_metric.F90
new file mode 100644
index 0000000..3df077e
--- /dev/null
+++ b/src/adm_metric.F90
@@ -0,0 +1,277 @@
+! $Header$
+
+#include "cctk.h"
+
+module adm_metric
+  use tensor
+  implicit none
+  private
+  public calc_4metric
+  public calc_4metricderivs
+  public calc_4metricderivs2
+  public calc_3metric
+  public calc_3metricderivs
+  
+  public calc_3metricdot
+  public calc_extcurv
+contains
+  
+  subroutine calc_4metric (gg, alfa, beta, g4)
+    CCTK_REAL, intent(in)  :: gg(3,3), alfa, beta(3)
+    CCTK_REAL, intent(out) :: g4(0:3,0:3)
+    CCTK_REAL :: betal(3)
+    
+    ! ds^2 = -alpha^2 dt^2 + g_ij (dx^i + beta^i dt) (dx^j + beta^j dt)
+    betal = matmul(gg, beta)
+    g4(0  ,0  ) = -alfa**2 + sum(betal*beta)
+    g4(1:3,0  ) = betal
+    g4(0  ,1:3) = betal
+    g4(1:3,1:3) = gg
+  end subroutine calc_4metric
+  
+  subroutine calc_4metricderivs (gg,alfa,beta, dgg,dalfa,dbeta, &
+       gg_dot,alfa_dot,beta_dot, g4,dg4)
+    CCTK_REAL, intent(in)  :: gg(3,3),alfa,beta(3)
+    CCTK_REAL, intent(in)  :: dgg(3,3,3),dalfa(3),dbeta(3,3)
+    CCTK_REAL, intent(in)  :: gg_dot(3,3),alfa_dot,beta_dot(3)
+    CCTK_REAL, intent(out) :: g4(0:3,0:3),dg4(0:3,0:3,0:3)
+    CCTK_REAL :: d4gg(3,3,0:3),d4alfa(0:3),d4beta(3,0:3)
+    CCTK_REAL :: betal(3),d4betal(3,0:3)
+    integer   :: i,j
+    
+    ! 4-metric
+    forall (i=1:3)
+       betal(i) = sum(gg(i,:) * beta(:))
+    end forall
+    g4(0  ,0  ) = -alfa**2 + sum(betal(:) * beta(:))
+    g4(1:3,0  ) = betal
+    g4(0  ,1:3) = betal
+    g4(1:3,1:3) = gg
+    
+    ! derivatives
+    d4gg  (:,:,0  ) = gg_dot(:,:)
+    d4gg  (:,:,1:3) = dgg(:,:,:)
+    d4alfa(    0  ) = alfa_dot
+    d4alfa(    1:3) = dalfa(:)
+    d4beta(:,  0  ) = beta_dot(:)
+    d4beta(:,  1:3) = dbeta(:,:)
+    
+    forall (i=1:3, j=0:3)
+       d4betal(i,j) = sum(d4gg(i,:,j) * beta(:) + gg(i,:) * d4beta(:,j))
+    end forall
+    
+    forall (i=0:3)
+       dg4(0  ,0  ,i) = - 2 * alfa * d4alfa(i) &
+            &           + sum(d4betal(:,i) * beta(:) + betal(:) * d4beta(:,i))
+       dg4(1:3,0  ,i) = d4betal(:,i)
+       dg4(0  ,1:3,i) = d4betal(:,i)
+       dg4(1:3,1:3,i) = d4gg(:,:,i)
+    end forall
+  end subroutine calc_4metricderivs
+  
+  subroutine calc_4metricderivs2 (gg,alfa,beta, dgg,dalfa,dbeta, &
+       ddgg,ddalfa,ddbeta, gg_dot,alfa_dot,beta_dot, &
+       gg_dot2,alfa_dot2,beta_dot2, dgg_dot,dalfa_dot,dbeta_dot, g4,dg4,ddg4)
+    CCTK_REAL, intent(in)  :: gg(3,3),alfa,beta(3)
+    CCTK_REAL, intent(in)  :: dgg(3,3,3),dalfa(3),dbeta(3,3)
+    CCTK_REAL, intent(in)  :: ddgg(3,3,3,3),ddalfa(3,3),ddbeta(3,3,3)
+    CCTK_REAL, intent(in)  :: gg_dot(3,3),alfa_dot,beta_dot(3)
+    CCTK_REAL, intent(in)  :: gg_dot2(3,3),alfa_dot2,beta_dot2(3)
+    CCTK_REAL, intent(in)  :: dgg_dot(3,3,3),dalfa_dot(3),dbeta_dot(3,3)
+    CCTK_REAL, intent(out) :: g4(0:3,0:3),dg4(0:3,0:3,0:3)
+    CCTK_REAL, intent(out) :: ddg4(0:3,0:3,0:3,0:3)
+    CCTK_REAL :: d4gg(3,3,0:3),d4alfa(0:3),d4beta(3,0:3)
+    CCTK_REAL :: dd4gg(3,3,0:3,0:3),dd4alfa(0:3,0:3),dd4beta(3,0:3,0:3)
+    CCTK_REAL :: betal(3),d4betal(3,0:3),dd4betal(3,0:3,0:3)
+    integer   :: i,j,k
+    
+    ! 4-metric
+    forall (i=1:3)
+       betal(i) = sum(gg(i,:) * beta(:))
+    end forall
+    g4(0  ,0  ) = -alfa**2 + sum(betal(:) * beta(:))
+    g4(1:3,0  ) = betal
+    g4(0  ,1:3) = betal
+    g4(1:3,1:3) = gg
+    
+    ! first derivatives
+    d4gg  (:,:,0  ) = gg_dot(:,:)
+    d4gg  (:,:,1:3) = dgg(:,:,:)
+    d4alfa(    0  ) = alfa_dot
+    d4alfa(    1:3) = dalfa(:)
+    d4beta(:,  0  ) = beta_dot(:)
+    d4beta(:,  1:3) = dbeta(:,:)
+    
+    forall (i=1:3, j=0:3)
+       d4betal(i,j) = sum(d4gg(i,:,j) * beta(:) + gg(i,:) * d4beta(:,j))
+    end forall
+    
+    forall (i=0:3)
+       dg4(0  ,0  ,i) = - 2 * alfa * d4alfa(i) &
+            &           + sum(d4betal(:,i) * beta(:) + betal(:) * d4beta(:,i))
+       dg4(1:3,0  ,i) = d4betal(:,i)
+       dg4(0  ,1:3,i) = d4betal(:,i)
+       dg4(1:3,1:3,i) = d4gg(:,:,i)
+    end forall
+    
+    ! second derivatives
+    dd4gg  (:,:,0  ,0  ) = gg_dot2(:,:)
+    dd4gg  (:,:,1:3,0  ) = dgg_dot(:,:,:)
+    dd4gg  (:,:,0  ,1:3) = dgg_dot(:,:,:)
+    dd4gg  (:,:,1:3,1:3) = ddgg(:,:,:,:)
+    dd4alfa(    0  ,0  ) = alfa_dot2
+    dd4alfa(    1:3,0  ) = dalfa_dot(:)
+    dd4alfa(    0  ,1:3) = dalfa_dot(:)
+    dd4alfa(    1:3,1:3) = ddalfa(:,:)
+    dd4beta(:,  0  ,0  ) = beta_dot2(:)
+    dd4beta(:,  1:3,0  ) = dbeta_dot(:,:)
+    dd4beta(:,  0  ,1:3) = dbeta_dot(:,:)
+    dd4beta(:,  1:3,1:3) = ddbeta(:,:,:)
+    
+    ! betal(i) = gg(i,m) * beta(m)
+    ! d4betal(i,j) = d4gg(i,m,j) * beta(m) + gg(i,m) * d4beta(m,j)
+    ! dd4betal(i,j,k) =   dd4gg(i,m,j,k) * beta(m) + d4gg(i,m,j) * d4beta(m,k)
+    !                   + d4gg(i,m,k) * d4beta(m,j) + gg(i,m) * dd4beta(m,j,k)
+    forall (i=1:3, j=0:3, k=0:3)
+       dd4betal(i,j,k) = sum(+ dd4gg(i,:,j,k) *    beta(:)      &
+            &                +  d4gg(i,:,j)   *  d4beta(:,k)    &
+            &                +  d4gg(i,:,k)   *  d4beta(:,j)    &
+            &                +    gg(i,:)     * dd4beta(:,j,k))
+    end forall
+    
+    ! g4(0  ,0  ) = -alfa**2 + sum(betal*beta)
+    ! g4(1:3,0  ) = betal
+    ! g4(0  ,1:3) = betal
+    ! g4(1:3,1:3) = gg
+    
+    ! dg4(0  ,0  ,i) = -2*alfa*d4alfa(i) &
+    !      &           + d4betal(m,i)*beta(m) + betal(m)*d4beta(m,i)
+    ! dg4(1:3,0  ,i) = d4betal(:,i)
+    ! dg4(0  ,1:3,i) = d4betal(:,i)
+    ! dg4(1:3,1:3,i) = d4gg(:,:,i)
+    forall (i=0:3, j=0:3)
+       ddg4(0  ,0  ,i,j) = - 2 * d4alfa(j) * d4alfa(i)        &
+            &              - 2 * alfa * dd4alfa(i,j)          &
+            &              + sum(+ dd4betal(:,i,j) * beta(:)  &
+            &                    + d4betal(:,i) * d4beta(:,j) &
+            &                    + d4betal(:,j) * d4beta(:,i) &
+            &                    + betal(:) * dd4beta(:,i,j))
+       ddg4(1:3,0  ,i,j) = dd4betal(:,i,j)
+       ddg4(0  ,1:3,i,j) = dd4betal(:,i,j)
+       ddg4(1:3,1:3,i,j) = dd4gg(:,:,i,j)
+    end forall
+  end subroutine calc_4metricderivs2
+  
+  
+  
+  subroutine calc_3metric (g4, gg, alfa, beta)
+    CCTK_REAL, intent(in)  :: g4(0:3,0:3)
+    CCTK_REAL, intent(out) :: gg(3,3), alfa, beta(3)
+    CCTK_REAL :: betal(3)
+    CCTK_REAL :: dtg, gu(3,3)
+    
+    ! ds^2 = -alpha^2 dt^2 + g_ij (dx^i + beta^i dt) (dx^j + beta^j dt)
+    betal = g4(1:3,0)
+    gg = g4(1:3,1:3)
+    call calc_det (gg, dtg)
+    call calc_inv (gg, dtg, gu)
+    beta = matmul(gu, betal)
+    alfa = sqrt(sum(betal*beta) - g4(0,0))
+  end subroutine calc_3metric
+  
+  subroutine calc_3metricderivs (g4,dg4, gg,alfa,beta, dgg,dalfa,dbeta, &
+       gg_dot,alfa_dot,beta_dot)
+    CCTK_REAL, intent(in)  :: g4(0:3,0:3),dg4(0:3,0:3,0:3)
+    CCTK_REAL, intent(out) :: gg(3,3),alfa,beta(3)
+    CCTK_REAL, intent(out) :: dgg(3,3,3),dalfa(3),dbeta(3,3)
+    CCTK_REAL, intent(out) :: gg_dot(3,3),alfa_dot,beta_dot(3)
+    CCTK_REAL :: betal(3),d4betal(3,0:3)
+    CCTK_REAL :: dtg,gu(3,3),dgu(3,3,3),gu_dot(3,3)
+    CCTK_REAL :: d4gg(3,3,0:3),d4gu(3,3,0:3)
+    CCTK_REAL :: d4alfa(0:3),d4beta(3,0:3)
+    integer   :: i,j
+    
+    ! ds^2 = -alpha^2 dt^2 + g_ij (dx^i + beta^i dt) (dx^j + beta^j dt)
+    gg = g4(1:3,1:3)
+    call calc_det (gg, dtg)
+    call calc_inv (gg, dtg, gu)
+    betal = g4(1:3,0)
+    beta = matmul(gu, betal)
+    alfa = sqrt(sum(betal*beta) - g4(0,0))
+    
+    forall (i=0:3)
+       d4gg(:,:,i) = dg4(1:3,1:3,i)
+    end forall
+    gg_dot = d4gg(:,:,0)
+    dgg = d4gg(:,:,1:3)
+    
+    call calc_invderiv (gu, dgg, dgu)
+    call calc_invdot (gu, gg_dot, gu_dot)
+    d4gu(:,:,0) = gu_dot
+    d4gu(:,:,1:3) = dgu
+    
+    forall (i=0:3)
+       d4betal(:,i) = dg4(1:3,0,i)
+    end forall
+    forall (i=0:3, j=1:3)
+       d4beta(j,i) = sum(d4gu(j,:,i)*betal) + sum(gu(j,:)*d4betal(:,i))
+    end forall
+    forall (i=0:3)
+       d4alfa(i) = 1/(2*alfa) &
+            * (sum(d4betal(:,i)*beta) + sum(betal*d4beta(:,i)) - dg4(0,0,i))
+    end forall
+    
+    alfa_dot = d4alfa(0)
+    dalfa = d4alfa(1:3)
+    beta_dot = d4beta(:,0)
+    dbeta = d4beta(:,1:3)
+  end subroutine calc_3metricderivs
+  
+  
+  
+  subroutine calc_3metricdot (gg, dgg, kk, alfa, beta, dbeta, gg_dot)
+    CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3)
+    CCTK_REAL, intent(in)  :: kk(3,3)
+    CCTK_REAL, intent(in)  :: alfa
+    CCTK_REAL, intent(in)  :: beta(3), dbeta(3,3)
+    CCTK_REAL, intent(out) :: gg_dot(3,3)
+    integer :: i,j,k
+    
+    ! d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k
+    
+    do i=1,3
+       do j=1,3
+          gg_dot(i,j) = -2*alfa * kk(i,j)
+          do k=1,3
+             gg_dot(i,j) = gg_dot(i,j) &
+                  + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) &
+                  + beta(k) * dgg(i,j,k)
+          end do
+       end do
+    end do
+  end subroutine calc_3metricdot
+  
+  subroutine calc_extcurv (gg, dgg, gg_dot, alfa, beta, dbeta, kk)
+    CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3), gg_dot(3,3)
+    CCTK_REAL, intent(in)  :: alfa
+    CCTK_REAL, intent(in)  :: beta(3), dbeta(3,3)
+    CCTK_REAL, intent(out) :: kk(3,3)
+    integer :: i,j,k
+    
+    ! d/dt g_ij = -2 alpha K_ij + g_kj beta^k,i + g_ik beta^k,j + beta^k g_ij,k
+    
+    do i=1,3
+       do j=1,3
+          kk(i,j) = - gg_dot(i,j)
+          do k=1,3
+             kk(i,j) = kk(i,j) &
+                  + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) &
+                  + beta(k) * dgg(i,j,k)
+          end do
+          kk(i,j) = kk(i,j) / (2*alfa)
+       end do
+    end do
+    
+  end subroutine calc_extcurv
+    
+end module adm_metric
diff --git a/src/depend.pl b/src/depend.pl
deleted file mode 100644
index b41b8f6..0000000
--- a/src/depend.pl
+++ /dev/null
@@ -1,46 +0,0 @@
-#!/usr/bin/perl -w
-# $Header$
-
-# Create dependencies for Fortran 90 "use" and "include" statements
-
-$src = $ARGV[0];
-$srcfile = $ARGV[1];
-$dest = $ARGV[2];
-$srcdir = $ARGV[3];
-$moddir = $ARGV[4];
-$srcsuffix = $ARGV[5];
-$modsuffix = $ARGV[6];
-@otherdirs = @ARGV[7..$#ARGV];
-
-print "# $src\n";
-print "$dest:";
-
-open (FILE, $src);
-while (<FILE>) {
-  if (/^\s*include\s*['"]([^'"]+)['"]/i) {
-    print " $srcdir$1";
-  } elsif (/^\s*use\s+(\w+)/i) {
-    $found = 0;
-    if (! $found) {
-      if (-e "$srcdir$1$srcsuffix") {
-	$found = 1;
-	print " $moddir$1$modsuffix";
-      }
-    }
-    if (! $found) {
-      foreach $dir (@otherdirs) {
-	if (-e "$dir$1$modsuffix") {
-	  $found = 1;
-	  print " $dir$1$modsuffix";
-	  last;
-	}
-      }
-    }
-    if (! $found) {
-      die "\nWhile tracing depencencies:\nFortran module $1 (referenced in file $srcfile) not found.\nAborting.\n\n";
-    }
-  }
-}
-close (FILE);
-
-print "\n";
diff --git a/src/gram_schmidt.F90 b/src/gram_schmidt.F90
new file mode 100644
index 0000000..61e865a
--- /dev/null
+++ b/src/gram_schmidt.F90
@@ -0,0 +1,121 @@
+! $Header$
+
+#include "cctk.h"
+
+module gram_schmidt
+  implicit none
+  private
+  public gram_schmidt_project4
+  public gram_schmidt_normalise4
+  
+contains
+  
+  subroutine gram_schmidt_project4 (g, y, dy, ddy, y2, x, dx, ddx)
+    CCTK_REAL, intent(in)    :: g(0:3,0:3)
+    CCTK_REAL, intent(in)    :: y(0:3), dy(0:3,0:3), ddy(0:3,0:3,0:3), y2
+    CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3), ddx(0:3,0:3,0:3)
+    CCTK_REAL                :: z(0:3), dz(0:3,0:3), ddz(0:3,0:3,0:3)
+    
+    integer :: a, b, c, d, e
+    
+    do a=0,3
+       z(a) = x(a)
+       do d=0,3
+          do e=0,3
+             z(a) = z(a) - g(d,e) * x(d) * y(e) * y(a) / y2
+          end do
+       end do
+    end do
+    
+    do a=0,3
+       do b=0,3
+          dz(a,b) = dx(a,b)
+          do d=0,3
+             do e=0,3
+                dz(a,b) = dz(a,b) - g(d,e) * ((dx(d,b) * y(e) + x(d) * dy(e,b)) * y(a) + x(d) * y(e) * dy(a,b)) / y2
+             end do
+          end do
+       end do
+    end do
+    
+    do a=0,3
+       do b=0,3
+          do c=0,3
+             ddz(a,b,c) = ddx(a,b,c)
+             do d=0,3
+                do e=0,3
+                   ddz(a,b,c) = ddz(a,b,c) - g(d,e) * (  (ddx(d,b,c) * y(e) + dx(d,b) * dy(e,c) + dx(d,c) * dy(e,b) + x(d) * ddy(e,b,c)) * y(a) + (dx(d,b) * y(e) + x(d) * dy(e,b)) * dy(a,c) &
+                        &                              + dx(d,c) * y(e) * dy(a,b) + x(d) * dy(e,c) * dy(a,b) + x(d) * y(e) * ddy(a,b,c)) / y2
+                end do
+             end do
+          end do
+       end do
+    end do
+    
+    x = z
+    dx = dz
+    ddx = ddz
+  end subroutine gram_schmidt_project4
+  
+  subroutine gram_schmidt_normalise4 (g, x, dx, ddx, x2)
+    CCTK_REAL, intent(in)    :: g(0:3,0:3)
+    CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3), ddx(0:3,0:3,0:3)
+    CCTK_REAL, intent(in)    :: x2
+    CCTK_REAL, parameter     :: two=2, half=1/two
+    CCTK_REAL                :: z2, dz2(0:3), ddz2(0:3,0:3)
+    CCTK_REAL                :: z(0:3), dz(0:3,0:3), ddz(0:3,0:3,0:3)
+    
+    integer :: a, b, c, d, e
+    
+    z2 = 0
+    do d=0,3
+       do e=0,3
+          z2 = z2 + g(d,e) * x(d) * x(e)
+       end do
+    end do
+    
+    do a=0,3
+       dz2(a) = 0
+       do d=0,3
+          do e=0,3
+             dz2(a) = dz2(a) + 2 * g(d,e) * dx(d,a) * x(e)
+          end do
+       end do
+    end do
+    
+    do a=0,3
+       do b=0,3
+          ddz2(a,b) = 0
+          do d=0,3
+             do e=0,3
+                ddz2(a,b) = ddz2(a,b) + 2 * g(d,e) * ddx(d,a,b) * x(e) + 2 * g(d,e) * dx(d,a) * dx(e,b)
+             end do
+          end do
+       end do
+    end do
+    
+    do a=0,3
+       z(a) = x(a) / sqrt(z2 / x2)
+    end do
+    
+    do a=0,3
+       do b=0,3
+          dz(a,b) = (dx(a,b) * z2 - half * x(a) * dz2(b) / x2) / sqrt(z2 / x2)**3
+       end do
+    end do
+    
+    do a=0,3
+       do b=0,3
+          do c=0,3
+             ddz(a,b,c) = ((ddx(a,b,c) * z2 + dx(a,b) * dz2(c) / x2 - half * (dx(a,c) * dz2(b) / x2 + x(a) * ddz2(b,c) / x2**2)) * z2 &
+                  &        - 3*half * (dx(a,b) * z2 - half * x(a) * dz2(b) / x2) * dz2(c) / x2) / sqrt(z2 / x2)**5
+          end do
+       end do
+    end do
+    
+    x = z
+    dx = dz
+    ddx = ddz
+  end subroutine gram_schmidt_normalise4
+  
+end module gram_schmidt
diff --git a/src/isnan.c b/src/isnan.c
new file mode 100644
index 0000000..a2feeb6
--- /dev/null
+++ b/src/isnan.c
@@ -0,0 +1,32 @@
+/* $Header$ */
+
+#include <math.h>
+
+#include "cctk.h"
+
+int CCTK_FCALL CCTK_FNAME(TAT_isnan) (const CCTK_REAL * restrict const x)
+{
+#ifdef HAVE_ISNAN
+  return isnan(*x);
+#else
+  return 0;
+#endif
+}
+
+int CCTK_FCALL CCTK_FNAME(TAT_isinf) (const CCTK_REAL * restrict const x)
+{
+#ifdef HAVE_ISINF
+  return isinf(*x);
+#else
+  return 0;
+#endif
+}
+
+int CCTK_FCALL CCTK_FNAME(TAT_finite) (const CCTK_REAL * restrict const x)
+{
+#ifdef HAVE_FINITE
+  return finite(*x);
+#else
+  return 1;
+#endif
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index 69e3078..462985e 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,21 +2,28 @@
 # $Header$
 
 # Source files in this directory
-SRCS =	cactus.F90				\
+SRCS =	adm_metric.F90				\
+	cactus.F90				\
 	constants.F90				\
 	conversion.F90				\
 	covderivs.F90				\
 	covderivs2.F90				\
 	derivs.F90				\
 	derivs2.F90				\
+	gram_schmidt.F90			\
+	isnan.c					\
 	lapack.F90				\
 	matdet.F90				\
 	matexp.F90				\
 	matinv.F90				\
+	ricci.F90				\
+	ricci2.F90				\
+	ricci4.F90				\
 	rotation.F90				\
 	tensor.F90				\
 	tensor2.F90				\
 	tensor4.F90				\
+	timederivs.F90				\
 	pointwise.F90				\
 	pointwise2.F90
 
diff --git a/src/ricci.F90 b/src/ricci.F90
new file mode 100644
index 0000000..ebe01fc
--- /dev/null
+++ b/src/ricci.F90
@@ -0,0 +1,85 @@
+! $Header$
+
+#include "cctk.h"
+
+module ricci
+  implicit none
+  private
+  public calc_connections
+  public calc_connectionderivs
+  public calc_ricci
+  
+contains
+  
+  subroutine calc_connections (gu, dg, gamma)
+    CCTK_REAL, intent(in)  :: gu(3,3), dg(3,3,3)
+    CCTK_REAL, intent(out) :: gamma(3,3,3)
+    integer :: i,j,k,l
+    ! Gamma^i_jk = 1/2 g^il (g_lj,k + g_lk,j - g_jk,l)
+    do i=1,3
+       do j=1,3
+          do k=1,3
+             gamma(i,j,k) = 0
+             do l=1,3
+                gamma(i,j,k) = gamma(i,j,k) + 0.5d0 * gu(i,l) &
+                     * (dg(l,j,k) + dg(l,k,j) - dg(j,k,l))
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_connections
+  
+  subroutine calc_connectionderivs (gu, dgg, dgu, ddgg, dgamma)
+    CCTK_REAL, intent(in)  :: gu(3,3), dgg(3,3,3), dgu(3,3,3), ddgg(3,3,3,3)
+    CCTK_REAL, intent(out) :: dgamma(3,3,3,3)
+    integer   :: i,j,k,l,m
+    ! Gamma^i_jk,l = 1/2 g^im,l (g_mj,k + g_mk,j - g_jk,m)
+    !                + 1/2 g^im (g_mj,kl + g_mk,jl - g_jk,ml)
+    do i=1,3
+       do j=1,3
+          do k=1,3
+             do l=1,3
+                dgamma(i,j,k,l) = 0
+                do m=1,3
+                   dgamma(i,j,k,l) = dgamma(i,j,k,l) &
+                        + 0.5d0 * dgu(i,m,l) * (dgg(m,j,k) + dgg(m,k,j) - dgg(j,k,m)) &
+                        + 0.5d0 * gu(i,m) * (ddgg(m,j,k,l) + ddgg(m,k,j,l) - ddgg(j,k,m,l))
+                end do
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_connectionderivs
+  
+  subroutine calc_ricci (gamma, dgamma, ri)
+    CCTK_REAL, intent(in)  :: gamma(3,3,3), dgamma(3,3,3,3)
+    CCTK_REAL, intent(out) :: ri(3,3)
+    CCTK_REAL :: sum, cnt
+    integer :: i,j,k,l
+    ! R_ij =   Gamma^k_ij,k - Gamma^k_ik,j
+    !        + Gamma^k_lk Gamma^l_ij - Gamma^k_lj Gamma^l_ki
+    do i=1,3
+       do j=1,3
+          ri(i,j) = 0
+          do k=1,3
+             ri(i,j) = ri(i,j) + dgamma(k,i,j,k) - dgamma(k,i,k,j)
+             do l=1,3
+                ri(i,j) = ri(i,j) + gamma(k,l,k) * gamma(l,i,j) &
+                     &            - gamma(k,l,j) * gamma(l,k,i)
+             end do
+          end do
+       end do
+    end do
+    ! check symmetries
+    sum = 0
+    cnt = 0
+    do i=1,3
+       do j=1,3
+          sum = sum + (ri(i,j) - ri(j,i))**2
+          cnt = cnt + ri(i,j)**2
+       end do
+    end do
+    if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "Ricci tensor is not symmetric")
+  end subroutine calc_ricci
+  
+end module ricci
diff --git a/src/ricci2.F90 b/src/ricci2.F90
new file mode 100644
index 0000000..dd2c429
--- /dev/null
+++ b/src/ricci2.F90
@@ -0,0 +1,74 @@
+! $Header$
+
+#include "cctk.h"
+
+module ricci2
+  implicit none
+  private
+  public calc_2connections
+  public calc_2connectionderivs
+  public calc_2ricci
+  
+contains
+  
+  subroutine calc_2connections (gu, dgg, gamma)
+    CCTK_REAL, intent(in)  :: gu(2,2), dgg(2,2,2)
+    CCTK_REAL, intent(out) :: gamma(2,2,2)
+    integer :: i,j,k,l
+    ! Gamma^i_jk = 1/2 g^il (g_lj,k + g_lk,j - g_jk,l)
+    do i=1,2
+       do j=1,2
+          do k=1,2
+             gamma(i,j,k) = 0
+             do l=1,2
+                gamma(i,j,k) = gamma(i,j,k) + 0.5d0 * gu(i,l) &
+                     * (dgg(l,j,k) + dgg(l,k,j) - dgg(j,k,l))
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2connections
+  
+  subroutine calc_2connectionderivs (gu, dgg, dgu, ddgg, dgamma)
+    CCTK_REAL, intent(in)  :: gu(2,2), dgg(2,2,2), dgu(2,2,2), ddgg(2,2,2,2)
+    CCTK_REAL, intent(out) :: dgamma(2,2,2,2)
+    integer   :: i,j,k,l,m
+    ! Gamma^i_jk,l = 1/2 g^im,l (g_mj,k + g_mk,j - g_jk,m)
+    !                + 1/2 g^im (g_mj,kl + g_mk,jl - g_jk,ml)
+    do i=1,2
+       do j=1,2
+          do k=1,2
+             do l=1,2
+                dgamma(i,j,k,l) = 0
+                do m=1,2
+                   dgamma(i,j,k,l) = dgamma(i,j,k,l) &
+                        + 0.5d0 * dgu(i,m,l) * (dgg(m,j,k) + dgg(m,k,j) - dgg(j,k,m)) &
+                        + 0.5d0 * gu(i,m) * (ddgg(m,j,k,l) + ddgg(m,k,j,l) - ddgg(j,k,m,l))
+                end do
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2connectionderivs
+  
+  subroutine calc_2ricci (gamma, dgamma, ri)
+    CCTK_REAL, intent(in)  :: gamma(2,2,2), dgamma(2,2,2,2)
+    CCTK_REAL, intent(out) :: ri(2,2)
+    integer :: i,j,k,l
+    ! R_ij =   Gamma^k_ij,k - Gamma^k_ik,j
+    !        + Gamma^k_lk Gamma^l_ij - Gamma^k_lj Gamma^l_ki
+    do i=1,2
+       do j=1,2
+          ri(i,j) = 0
+          do k=1,2
+             ri(i,j) = ri(i,j) + dgamma(k,i,j,k) - dgamma(k,i,k,j)
+             do l=1,2
+                ri(i,j) = ri(i,j) + gamma(k,l,k) * gamma(l,i,j) &
+                     &            - gamma(k,l,j) * gamma(l,k,i)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_2ricci
+  
+end module ricci2
diff --git a/src/ricci4.F90 b/src/ricci4.F90
new file mode 100644
index 0000000..ada7733
--- /dev/null
+++ b/src/ricci4.F90
@@ -0,0 +1,181 @@
+! $Header$
+
+#include "cctk.h"
+
+module ricci4
+  implicit none
+  private
+  public calc_4connections
+  public calc_4connectionderivs
+  public calc_4ricci
+  public calc_4riemann
+  public calc_4weyl
+  
+contains
+  
+  subroutine calc_4connections (gu, dg, gamma)
+    CCTK_REAL, intent(in)  :: gu(4,4), dg(4,4,4)
+    CCTK_REAL, intent(out) :: gamma(4,4,4)
+    integer :: i,j,k,l
+    ! Gamma^i_jk = 1/2 g^il (g_lj,k + g_lk,j - g_jk,l)
+    do i=1,4
+       do j=1,4
+          do k=1,4
+             gamma(i,j,k) = 0
+             do l=1,4
+                gamma(i,j,k) = gamma(i,j,k) + 0.5d0 * gu(i,l) &
+                     * (dg(l,j,k) + dg(l,k,j) - dg(j,k,l))
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_4connections
+  
+  subroutine calc_4connectionderivs (gu, dgg, dgu, ddgg, dgamma)
+    CCTK_REAL, intent(in)  :: gu(4,4), dgg(4,4,4), dgu(4,4,4), ddgg(4,4,4,4)
+    CCTK_REAL, intent(out) :: dgamma(4,4,4,4)
+    integer   :: i,j,k,l,m
+    ! Gamma^i_jk,l = 1/2 g^im,l (g_mj,k + g_mk,j - g_jk,m)
+    !                + 1/2 g^im (g_mj,kl + g_mk,jl - g_jk,ml)
+    do i=1,4
+       do j=1,4
+          do k=1,4
+             do l=1,4
+                dgamma(i,j,k,l) = 0
+                do m=1,4
+                   dgamma(i,j,k,l) = dgamma(i,j,k,l) &
+                        + 0.5d0 * dgu(i,m,l) * (dgg(m,j,k) + dgg(m,k,j) - dgg(j,k,m)) &
+                        + 0.5d0 * gu(i,m) * (ddgg(m,j,k,l) + ddgg(m,k,j,l) - ddgg(j,k,m,l))
+                end do
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_4connectionderivs
+  
+  subroutine calc_4ricci (gamma, dgamma, ri)
+    CCTK_REAL, intent(in)  :: gamma(4,4,4), dgamma(4,4,4,4)
+    CCTK_REAL, intent(out) :: ri(4,4)
+    CCTK_REAL :: sum, cnt
+    integer :: i,j,k,l
+    ! R_ij =   Gamma^k_ij,k - Gamma^k_ik,j
+    !        + Gamma^k_lk Gamma^l_ij - Gamma^k_lj Gamma^l_ki
+    do i=1,4
+       do j=1,4
+          ri(i,j) = 0
+          do k=1,4
+             ri(i,j) = ri(i,j) + dgamma(k,i,j,k) - dgamma(k,i,k,j)
+             do l=1,4
+                ri(i,j) = ri(i,j) + gamma(k,l,k) * gamma(l,i,j) &
+                     &            - gamma(k,l,j) * gamma(l,k,i)
+             end do
+          end do
+       end do
+    end do
+    ! check symmetries
+    sum = 0
+    cnt = 0
+    do i=1,4
+       do j=1,4
+          sum = sum + (ri(i,j) - ri(j,i))**2
+          cnt = cnt + ri(i,j)**2
+       end do
+    end do
+    if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "4-Ricci tensor is not symmetric")
+  end subroutine calc_4ricci
+  
+  subroutine calc_4riemann (gg, gamma, dgamma, rm)
+    CCTK_REAL, intent(in)  :: gg(4,4), gamma(4,4,4), dgamma(4,4,4,4)
+    CCTK_REAL, intent(out) :: rm(4,4,4,4)
+    CCTK_REAL :: rmu(4,4,4,4)
+    CCTK_REAL :: sum, cnt
+    integer :: i,j,k,l,m
+    ! R^i_jkl =   Gamma^i_jl,k - Gamma^i_jk,l
+    !           + Gamma^m_jl Gamma^i_mk - Gamma^m_jk Gamma^i_ml
+    ! (Wald, 3.4, eqn. (3.4.4), p. 48)
+    do i=1,4
+       do j=1,4
+          do k=1,4
+             do l=1,4
+                rmu(i,j,k,l) = dgamma(i,j,l,k) - dgamma(i,j,k,l)
+                do m=1,4
+                   rmu(i,j,k,l) = rmu(i,j,k,l) &
+                        + gamma(m,j,l) * gamma(i,m,k) &
+                        - gamma(m,j,k) * gamma(i,m,l)
+                end do
+             end do
+          end do
+       end do
+    end do
+    do i=1,4
+       do j=1,4
+          do k=1,4
+             do l=1,4
+                rm(i,j,k,l) = 0
+                do m=1,4
+                   rm(i,j,k,l) = rm(i,j,k,l) + gg(i,m) * rmu(m,j,k,l)
+                end do
+             end do
+          end do
+       end do
+    end do
+    ! check symmetries
+    sum = 0
+    cnt = 0
+    do i=1,4
+       do j=1,4
+          do k=1,4
+             do l=1,4
+                sum = sum + (rm(i,j,k,l) + rm(i,j,l,k))**2
+                sum = sum + (rm(i,j,k,l) - rm(k,l,i,j))**2
+                sum = sum + (rm(i,j,k,l) + rm(j,i,k,l))**2
+                sum = sum + (rm(i,j,k,l) + rm(j,k,i,l) + rm(k,i,j,l))**2
+                cnt = cnt + rm(i,j,k,l)**2
+             end do
+          end do
+       end do
+    end do
+    if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "4-Riemann tensor is not symmetric")
+  end subroutine calc_4riemann
+  
+  subroutine calc_4weyl (gg, rm, ri, rsc, we)
+    CCTK_REAL, parameter :: one=1, half=one/2, third=one/3
+    CCTK_REAL, intent(in)  :: gg(4,4), rm(4,4,4,4), ri(4,4), rsc
+    CCTK_REAL, intent(out) :: we(4,4,4,4)
+    CCTK_REAL :: sum, cnt
+    integer :: i,j,k,l
+    ! R_ijkl = C_ijkl + 2/(n-2) (g_i[k R_l]j - g_j[k R_l]i)
+    !                 - 2/(n-1)(n-2) R g_i[k g_l]j
+    ! (Wald, 3.2, eqn. (3.2.28), p. 40)
+    do i=1,4
+       do j=1,4
+          do k=1,4
+             do l=1,4
+                we(i,j,k,l) = rm(i,j,k,l) &
+                     -        (  (gg(i,k) * ri(l,j) - gg(i,l) * ri(k,j)) &
+                     &         - (gg(j,k) * ri(l,i) - gg(j,l) * ri(k,i))) &
+                     + third * rsc * (gg(i,k) * gg(l,j) - gg(i,l) * gg(k,j))
+             end do
+          end do
+       end do
+    end do
+    ! check symmetries
+    sum = 0
+    cnt = 0
+    do i=1,4
+       do j=1,4
+          do k=1,4
+             do l=1,4
+                sum = sum + (we(i,j,k,l) + we(i,j,l,k))**2
+                sum = sum + (we(i,j,k,l) - we(k,l,i,j))**2
+                sum = sum + (we(i,j,k,l) + we(j,i,k,l))**2
+                sum = sum + (we(i,j,k,l) + we(j,k,i,l) + we(k,i,j,l))**2
+                cnt = cnt + we(i,j,k,l)**2
+             end do
+          end do
+       end do
+    end do
+    if (sum > 1.0e-12 * cnt) call CCTK_WARN (0, "4-Weyl tensor is not symmetric")
+  end subroutine calc_4weyl
+  
+end module ricci4
diff --git a/src/timederivs.F90 b/src/timederivs.F90
new file mode 100644
index 0000000..42e107d
--- /dev/null
+++ b/src/timederivs.F90
@@ -0,0 +1,225 @@
+! $Header$
+
+#include "cctk.h"
+
+module timederivs
+  implicit none
+  private
+  public abs2
+  public timederiv
+  public timederiv2
+  public timederiv_uneven
+  public timederiv2_uneven
+  public operator(.outer.)
+  public operator(.dot.)
+  
+  interface timederiv
+     module procedure rtimederiv
+  end interface
+  
+  interface timederiv2
+     module procedure rtimederiv2
+  end interface
+  
+  interface timederiv_uneven
+     module procedure rtimederiv_uneven
+     module procedure ctimederiv_uneven
+  end interface
+  
+  interface timederiv2_uneven
+     module procedure rtimederiv2_uneven
+     module procedure ctimederiv2_uneven
+  end interface
+  
+  interface operator(.outer.)
+     module procedure router
+     module procedure couter
+  end interface
+  
+  interface operator(.dot.)
+     module procedure rdot
+     module procedure cdot
+  end interface
+  
+contains
+  
+  ! abs(c)**2 for complex c without a square root
+  elemental function abs2 (a)
+    CCTK_COMPLEX, intent(in) :: a
+    CCTK_REAL :: abs2
+    abs2 = a * conjg(a)
+  end function abs2
+  
+  
+  
+  ! Calculate a time derivate from several time levels
+  elemental function rtimederiv (f0, f1, f2, dt) result (fdot)
+    CCTK_REAL, intent(in) :: f0, f1, f2
+    CCTK_REAL, intent(in) :: dt
+    CCTK_REAL             :: fdot
+    CCTK_REAL :: dt1, dt2
+    CCTK_REAL :: fdot1, fdot2
+    
+    dt1 = dt
+    dt2 = 2*dt
+    
+    ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0)
+    ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0)
+    ! f'(0) = [f(dt1) - f(0)] / dt1 - dt1/2 f''(0)
+    ! f'(0) = [f(dt2) - f(0)] / dt2 - dt2/2 f''(0)
+    fdot1 = (f0 - f1) / dt1
+    fdot2 = (f0 - f2) / dt2
+    fdot = (fdot1 * dt2 - fdot2 * dt1) / (dt2 - dt1)
+  end function rtimederiv
+  
+  elemental function rtimederiv2 (f0, f1, f2, dt) result (fdotdot)
+    CCTK_REAL, intent(in) :: f0, f1, f2
+    CCTK_REAL, intent(in) :: dt
+    CCTK_REAL             :: fdotdot
+    CCTK_REAL :: dt1, dt2
+    CCTK_REAL :: fdotdot1, fdotdot2
+    
+    dt1 = dt
+    dt2 = 2*dt
+    
+    ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0)
+    ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0)
+    ! f''(0) = [f(dt1) - f(0)] / [dt1^2/2] - f'(0) / [dt1/2]
+    ! f''(0) = [f(dt2) - f(0)] / [dt2^2/2] - f'(0) / [dt2/2]
+    fdotdot1 = (f1 - f0) / (dt1**2/2)
+    fdotdot2 = (f2 - f0) / (dt2**2/2)
+    fdotdot = (fdotdot1 * dt1 - fdotdot2 * dt2) / (dt1 - dt2)
+  end function rtimederiv2
+  
+  
+  
+  ! Calculate a time derivate from several time levels with uneven spacing
+  elemental function rtimederiv_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdot)
+    CCTK_REAL, intent(in) :: f0, f1, f2
+    CCTK_REAL, intent(in) :: t0, t1, t2
+    CCTK_INT,  intent(in) :: ce0, ce1, ce2
+    CCTK_REAL             :: fdot
+    CCTK_REAL :: dt1, dt2
+    CCTK_REAL :: fdot1, fdot2
+    
+!!$    dt1 = ih_time(hn) - ih_time_p(hn)
+!!$    dt2 = ih_time(hn) - ih_time_p_p(hn)
+    dt1 = t0 - t1
+    dt2 = t0 - t2
+    
+    if (ce1/=0) then
+       fdot = 0
+    else if (ce2/=0) then
+       fdot = (f0 - f1) / dt1
+    else
+       ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0)
+       ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0)
+       ! f'(0) = [f(dt1) - f(0)] / dt1 - dt1/2 f''(0)
+       ! f'(0) = [f(dt2) - f(0)] / dt2 - dt2/2 f''(0)
+       fdot1 = (f0 - f1) / dt1
+       fdot2 = (f0 - f2) / dt2
+       fdot = (fdot1 * dt2 - fdot2 * dt1) / (dt2 - dt1)
+    end if
+  end function rtimederiv_uneven
+  
+  elemental function ctimederiv_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdot)
+    CCTK_COMPLEX, intent(in) :: f0, f1, f2
+    CCTK_REAL,    intent(in) :: t0, t1, t2
+    CCTK_INT,     intent(in) :: ce0, ce1, ce2
+    CCTK_COMPLEX             :: fdot
+    
+    fdot = cmplx(timederiv_uneven(real(f0),real(f1),real(f2), t0,t1,t2, ce0,ce1,ce2), &
+         &       timederiv_uneven(aimag(f0),aimag(f1),aimag(f2), t0,t1,t2, ce0,ce1,ce2), &
+         &       kind(fdot))
+  end function ctimederiv_uneven
+  
+  elemental function rtimederiv2_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdotdot)
+    CCTK_REAL, intent(in) :: f0, f1, f2
+    CCTK_REAL, intent(in) :: t0, t1, t2
+    CCTK_INT,  intent(in) :: ce0, ce1, ce2
+    CCTK_REAL             :: fdotdot
+    CCTK_REAL :: dt1, dt2
+    CCTK_REAL :: fdotdot1, fdotdot2
+    
+!!$    dt1 = ih_time(hn) - ih_time_p(hn)
+!!$    dt2 = ih_time(hn) - ih_time_p_p(hn)
+    dt1 = t0 - t1
+    dt2 = t0 - t2
+    
+    if (ce1/=0) then
+       fdotdot = 0
+    else if (ce2/=0) then
+       fdotdot = 0
+    else
+       ! f(dt1) = f(0) + dt1 f'(0) + dt1^2/2 f''(0)
+       ! f(dt2) = f(0) + dt2 f'(0) + dt2^2/2 f''(0)
+       ! f''(0) = [f(dt1) - f(0)] / [dt1^2/2] - f'(0) / [dt1/2]
+       ! f''(0) = [f(dt2) - f(0)] / [dt2^2/2] - f'(0) / [dt2/2]
+       fdotdot1 = (f1 - f0) / (dt1**2/2)
+       fdotdot2 = (f2 - f0) / (dt2**2/2)
+       fdotdot = (fdotdot1 * dt1 - fdotdot2 * dt2) / (dt1 - dt2)
+    end if
+  end function rtimederiv2_uneven
+  
+  elemental function ctimederiv2_uneven (f0, f1, f2, t0, t1, t2, ce0, ce1, ce2) result (fdotdot)
+    CCTK_COMPLEX, intent(in) :: f0, f1, f2
+    CCTK_REAL,    intent(in) :: t0, t1, t2
+    CCTK_INT,     intent(in) :: ce0, ce1, ce2
+    CCTK_COMPLEX             :: fdotdot
+    
+    fdotdot = cmplx(timederiv2_uneven(real(f0),real(f1),real(f2), t0,t1,t2, ce0,ce1,ce2), &
+         &          timederiv2_uneven(aimag(f0),aimag(f1),aimag(f2), t0,t1,t2, ce0,ce1,ce2), &
+         &          kind(fdotdot))
+  end function ctimederiv2_uneven
+  
+  
+  
+  function router (left, right) result (result)
+    CCTK_REAL, intent(in) :: left(:), right(:)
+    CCTK_REAL             :: result(size(left,1),size(right,1))
+    integer :: i, j
+    forall (i=1:size(left,1), j=1:size(right,1))
+       result(i,j) = left(i) * right(j)
+    end forall
+  end function router
+  
+  function couter (left, right) result (result)
+    CCTK_COMPLEX, intent(in) :: left(:), right(:)
+    CCTK_COMPLEX             :: result(size(left,1),size(right,1))
+    integer :: i, j
+    forall (i=1:size(left,1), j=1:size(right,1))
+       result(i,j) = left(i) * right(j)
+    end forall
+  end function couter
+  
+  
+  
+  function rdot (left, right) result (result)
+    CCTK_REAL, intent(in) :: left(:), right(:)
+    CCTK_REAL             :: result
+    integer :: i
+    if (size(left,1) /= size(right,1)) then
+       call CCTK_WARN (0, "Array sizes must have the same sizes")
+    end if
+!!$    result = sum((/( left(i) * right(i), i=1,size(left,1) )/))
+    result = 0
+    do i=1,size(left,1)
+       result = result + left(i) * right(i)
+    end do
+  end function rdot
+  
+  function cdot (left, right) result (result)
+    CCTK_COMPLEX, intent(in) :: left(:), right(:)
+    CCTK_COMPLEX             :: result
+    integer :: i
+    if (size(left,1) /= size(right,1)) then
+       call CCTK_WARN (0, "Array sizes must have the same sizes")
+    end if
+!!$    result = sum((/( left(i) * right(i), i=1,size(left,1) )/))
+    result = 0
+    do i=1,size(left,1)
+       result = result + left(i) * right(i)
+    end do
+  end function cdot
+  
+end module timederivs