1 files changed, 127 insertions, 0 deletions
diff --git a/src/tensor.F90 b/src/tensor.F90
new file mode 100644
index 0000000..0fef831
--- /dev/null
+++ b/src/tensor.F90
@@ -0,0 +1,127 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+module tensor
+  implicit none
+  DECLARE_CCTK_PARAMETERS
+  private
+  
+  public calc_trace
+  
+  public calc_det
+  public calc_detderiv
+  public calc_detdot
+  
+  public calc_inv
+  public calc_invderiv
+  public calc_invdot
+  
+contains
+  
+  subroutine calc_trace (kk, gu, trk)
+    CCTK_REAL, intent(in)  :: kk(3,3)
+    CCTK_REAL, intent(in)  :: gu(3,3)
+    CCTK_REAL, intent(out) :: trk
+    integer :: i,j
+    trk = 0
+    do i=1,3
+       do j=1,3
+          trk = trk + gu(i,j) * kk(i,j)
+       end do
+    end do
+  end subroutine calc_trace
+  
+  
+  
+  subroutine calc_det (gg, dtg)
+    CCTK_REAL, intent(in)  :: gg(3,3)
+    CCTK_REAL, intent(out) :: dtg
+    dtg =      gg(1,1) * gg(2,2) * gg(3,3) &
+         & + 2*gg(1,2) * gg(1,3) * gg(2,3) &
+         & -   gg(1,1) * gg(2,3)**2        &
+         & -   gg(2,2) * gg(1,3)**2        &
+         & -   gg(3,3) * gg(1,2)**2
+  end subroutine calc_det
+  
+  subroutine calc_pdet (gg, pgg, pdtg)
+    CCTK_REAL, intent(in)  :: gg(3,3), pgg(3,3)
+    CCTK_REAL, intent(out) :: pdtg
+    pdtg =     pgg(1,1) *   gg(2,2) *  gg(3,3) &
+         & +    gg(1,1) *  pgg(2,2) *  gg(3,3) &
+         & +    gg(1,1) *   gg(2,2) * pgg(3,3) &
+         & + 2*pgg(1,2) *   gg(1,3) *  gg(2,3) &
+         & + 2* gg(1,2) *  pgg(1,3) *  gg(2,3) &
+         & + 2* gg(1,2) *   gg(1,3) * pgg(2,3) &
+         & -   pgg(1,1) *   gg(2,3)**2         &
+         & -    gg(1,1) * 2*gg(2,3) * pgg(2,3) &
+         & -   pgg(2,2) *   gg(1,3)**2         &
+         & -    gg(2,2) * 2*gg(1,3) * pgg(1,3) &
+         & -   pgg(3,3) *   gg(1,2)**2         &
+         & -    gg(3,3) * 2*gg(1,2) * pgg(1,2)
+  end subroutine calc_pdet
+  
+  subroutine calc_detderiv (gg, dgg, ddtg)
+    CCTK_REAL, intent(in)  :: gg(3,3), dgg(3,3,3)
+    CCTK_REAL, intent(out) :: ddtg(3)
+    integer :: i
+    do i=1,3
+       call calc_pdet (gg, dgg(:,:,i), ddtg(i))
+    end do
+  end subroutine calc_detderiv
+  
+  subroutine calc_detdot (gg, gg_dot, dtg_dot)
+    CCTK_REAL, intent(in)  :: gg(3,3), gg_dot(3,3)
+    CCTK_REAL, intent(out) :: dtg_dot
+    call calc_pdet (gg, gg_dot, dtg_dot)
+  end subroutine calc_detdot
+  
+  
+  
+  subroutine calc_inv (gg, dtg, gu)
+    CCTK_REAL, intent(in)  :: gg(3,3), dtg
+    CCTK_REAL, intent(out) :: gu(3,3)
+    gu(1,1) = (gg(2,2) * gg(3,3) - gg(2,3) ** 2     ) / dtg
+    gu(2,2) = (gg(1,1) * gg(3,3) - gg(1,3) ** 2     ) / dtg
+    gu(3,3) = (gg(1,1) * gg(2,2) - gg(1,2) ** 2     ) / dtg
+    gu(1,2) = (gg(1,3) * gg(2,3) - gg(1,2) * gg(3,3)) / dtg
+    gu(1,3) = (gg(1,2) * gg(2,3) - gg(1,3) * gg(2,2)) / dtg
+    gu(2,3) = (gg(1,3) * gg(1,2) - gg(2,3) * gg(1,1)) / dtg
+    gu(2,1) = gu(1,2)
+    gu(3,1) = gu(1,3)
+    gu(3,2) = gu(2,3)
+  end subroutine calc_inv
+  
+  subroutine calc_pinv (gu, pgg, pgu)
+    CCTK_REAL, intent(in)  :: gu(3,3), pgg(3,3)
+    CCTK_REAL, intent(out) :: pgu(3,3)
+    integer :: i,j,k,l
+    do i=1,3
+       do j=1,3
+          pgu(i,j) = 0
+          do k=1,3
+             do l=1,3
+                pgu(i,j) = pgu(i,j) - gu(i,k) * gu(j,l) * pgg(k,l)
+             end do
+          end do
+       end do
+    end do
+  end subroutine calc_pinv
+  
+  subroutine calc_invderiv (gu, dgg, dgu)
+    CCTK_REAL, intent(in)  :: gu(3,3), dgg(3,3,3)
+    CCTK_REAL, intent(out) :: dgu(3,3,3)
+    integer :: i
+    do i=1,3
+       call calc_pinv (gu, dgg(:,:,i), dgu(:,:,i))
+    end do
+  end subroutine calc_invderiv
+  
+  subroutine calc_invdot (gu, gg_dot, gu_dot)
+    CCTK_REAL, intent(in)  :: gu(3,3), gg_dot(3,3)
+    CCTK_REAL, intent(out) :: gu_dot(3,3)
+    call calc_pinv (gu, gg_dot, gu_dot)
+  end subroutine calc_invdot
+  
+end module tensor