1 files changed, 112 insertions, 0 deletions
diff --git a/src/qlm_3determinant.F90 b/src/qlm_3determinant.F90
new file mode 100644
index 0000000..d263009
--- /dev/null
+++ b/src/qlm_3determinant.F90
@@ -0,0 +1,112 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+
+
+SUBROUTINE qlm_calc_3determinant (CCTK_ARGUMENTS, hn)
+  USE adm_metric
+  USE cctk
+  USE qlm_boundary
+  USE qlm_derivs
+  USE qlm_variables
+  USE tensor
+  
+  IMPLICIT NONE
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+  INTEGER :: hn
+  
+  CCTK_REAL    :: gg(3,3)
+  CCTK_REAL    :: alfa   
+  CCTK_REAL    :: beta(3)
+  CCTK_REAL    :: g4(0:3,0:3)
+  CCTK_COMPLEX :: mm(0:3)
+  CCTK_REAL    :: rr(3), rr_p(3), rr_p_p(3)
+  CCTK_REAL    :: Ttilde(0:3)
+  CCTK_REAL    :: a
+  CCTK_COMPLEX :: d
+
+  CCTK_REAL    :: t0, t1, t2
+  logical      :: ce0, ce1, ce2
+  CCTK_REAL    :: delta_space(2)
+
+  INTEGER      :: i, j
+  CCTK_REAL    :: theta, phi
+  
+  IF (veryverbose/=0) THEN
+     CALL CCTK_INFO ("Computing 3-Volume Element of Horizon")
+  END IF
+  
+  t0 = qlm_time(hn)
+  t1 = qlm_time_p(hn)
+  t2 = qlm_time_p_p(hn)
+  
+  ce0 = qlm_have_valid_data(hn) == 0
+  ce1 = qlm_have_valid_data_p(hn) == 0
+  ce2 = qlm_have_valid_data_p_p(hn) == 0
+  
+  delta_space(:) = (/ qlm_delta_theta(hn), qlm_delta_phi(hn) /)
+  
+  ! Calculate the coordinates
+  DO j = 1+qlm_nghostsphi(hn), qlm_nphi(hn)-qlm_nghostsphi(hn)
+     DO i = 1+qlm_nghoststheta(hn), qlm_ntheta(hn)-qlm_nghoststheta(hn)
+        theta = qlm_origin_theta(hn) + (i-1)*qlm_delta_theta(hn)
+        phi   = qlm_origin_phi(hn)   + (j-1)*qlm_delta_phi(hn)
+        
+        ! Get the stuff from the arrays
+        gg(1,1) = qlm_gxx(i,j)
+        gg(1,2) = qlm_gxy(i,j)
+        gg(1,3) = qlm_gxz(i,j)
+        gg(2,2) = qlm_gyy(i,j)
+        gg(2,3) = qlm_gyz(i,j)
+        gg(3,3) = qlm_gzz(i,j)
+        gg(2,1) = gg(1,2)
+        gg(3,1) = gg(1,3)
+        gg(3,2) = gg(2,3)
+        
+        alfa = qlm_alpha(i,j)
+        
+        beta(1) = qlm_betax(i,j)
+        beta(2) = qlm_betay(i,j)
+        beta(3) = qlm_betaz(i,j)
+        
+        mm(0) = qlm_m0(i,j,hn)
+        mm(1) = qlm_m1(i,j,hn)
+        mm(2) = qlm_m2(i,j,hn)
+        mm(3) = qlm_m3(i,j,hn)
+        
+        ! Build the four-metric
+        CALL calc_4metric (gg,alfa,beta, g4)
+        
+        ! Find a 3rd vector of triad
+        rr(1) = qlm_x(i,j,hn)
+        rr(2) = qlm_y(i,j,hn)
+        rr(3) = qlm_z(i,j,hn)
+        
+        rr_p(1) = qlm_x_p(i,j,hn)
+        rr_p(2) = qlm_y_p(i,j,hn)
+        rr_p(3) = qlm_z_p(i,j,hn)
+        
+        rr_p_p(1) = qlm_x_p_p(i,j,hn)
+        rr_p_p(2) = qlm_y_p_p(i,j,hn)
+        rr_p_p(3) = qlm_z_p_p(i,j,hn)
+        
+        Ttilde(0)   = 1
+        Ttilde(1:3) = timederiv (rr, rr_p, rr_p_p, t0,t1,t2, ce0,ce1,ce2)
+        
+        ! Compute some scalar products
+        a=SUM(MATMUL(g4, Ttilde)*Ttilde)
+        d=SUM(MATMUL(g4, mm)*Ttilde)
+        
+        ! This is the determinant
+        qlm_3det(i,j,hn)=a-2*CONJG(d)*d
+        
+     END DO
+  END DO
+
+  CALL set_boundary (CCTK_PASS_FTOF, hn, qlm_3det(:,:,hn), +1)
+
+END SUBROUTINE qlm_calc_3determinant