1 files changed, 85 insertions, 0 deletions
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