1 files changed, 130 insertions, 0 deletions
diff --git a/src/qlm_killing_transport.F90 b/src/qlm_killing_transport.F90
new file mode 100644
index 0000000..3ec4f7e
--- /dev/null
+++ b/src/qlm_killing_transport.F90
@@ -0,0 +1,130 @@
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+
+
+subroutine qlm_killing_transport (CCTK_ARGUMENTS, hn)
+  use cctk
+  use constants
+  use qlm_boundary
+  use qlm_killing_transportation
+  use lapack
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+  integer :: hn
+  
+  interface
+     integer function TAT_isnan (x)
+       implicit none
+       CCTK_REAL x
+     end function TAT_isnan
+  end interface
+  
+  CCTK_REAL :: xi(2,3), chi(3)
+  CCTK_REAL :: vec(3,3)
+  CCTK_REAL :: wr(3), wi(3), vl(3,3), vr(3,3)
+  integer   :: i0, j0
+  integer   :: n
+  integer   :: info
+  character :: msg*1000
+  
+  integer, parameter :: lwork = 100
+  CCTK_REAL :: work(lwork)
+  
+  
+  
+  if (veryverbose/=0) then
+     call CCTK_INFO ("Transporting Killing vector field")
+  end if
+  
+  
+  
+  ! latitude of "equator"
+  i0 = (qlm_ntheta(hn)+1)/2
+  ! longitude of zero meridian
+  j0 = 1+qlm_nghostsphi(hn)
+  
+  vec = delta3
+  
+  if (veryverbose/=0) then
+     write (msg, '("Initial vectors (xi^t xi^p chi):")')
+     call CCTK_INFO (msg)
+     do n=1,3
+        write (msg, '(i2,3g18.8)') n, vec(:,n)
+        call CCTK_INFO (msg)
+     end do
+  end if
+  
+  xi(1,:) = vec(1,:)
+  xi(2,:) = vec(2,:)
+  chi(:)  = vec(3,:)
+  
+  do n=1,3
+     call transport_along_equator (CCTK_PASS_FTOF, hn, i0, xi(:,n), chi(n))
+  end do
+  
+  vec(1,:) = xi(1,:)
+  vec(2,:) = xi(2,:)
+  vec(3,:) = chi(:)
+  
+  if (veryverbose/=0) then
+     write (msg, '("Final vectors (xi^t xi^p chi):")')
+     call CCTK_INFO (msg)
+     do n=1,3
+        write (msg, '(i2,3g18.8)') n, vec(:,n)
+        call CCTK_INFO (msg)
+     end do
+  end if
+  
+  if (TAT_isnan(sum(vec)) /= 0) then
+     ! qlm_calc_error(hn) = 1
+     qlm_have_killing_vector(hn) = 0
+     call CCTK_WARN (3, "There are nans in the final vectors")
+     goto 9999
+  end if
+  
+  call geev ('n', 'v', 3, vec, 3, wr, wi, vl, 3, vr, 3, work, lwork, info)
+  if (info/=0) then
+     ! qlm_calc_error(hn) = 1
+     qlm_have_killing_vector(hn) = 0
+     write (msg, '("Error in call to GEEV, info=",i2)') info
+     call CCTK_WARN (3, msg)
+     goto 9999
+  end if
+  
+  if (veryverbose/=0) then
+     write (msg, '("Eigenvectors and eigenvalues:")')
+     call CCTK_INFO (msg)
+     do n=1,3
+        write (msg, '(i2,3g18.8,(" (",g18.8,",",g18.8,")"))') &
+             n, vr(:,n), wr(n), wi(n)
+        call CCTK_INFO (msg)
+     end do
+  end if
+  
+  xi(1,:) = vr(1,:)
+  xi(2,:) = vr(2,:)
+  chi(:)  = vr(3,:)
+  
+  ! TODO:
+  ! This transport scheme is not ideal.
+  ! It leads to large fluctuations in the phi direction.
+  n=3
+  qlm_killing_eigenvalue_re(hn) = wr(n)
+  qlm_killing_eigenvalue_im(hn) = wi(n)
+  if (abs(cmplx(wr(n),wi(n),kind(wr)) - (1,0)) > 1.0d-4) then
+     call CCTK_WARN (3, "Did not manage to find an eigenvector with the eigenvalue 1")
+  end if
+  call transport_along_equator (CCTK_PASS_FTOF, hn, i0, xi(:,n), chi(n))
+  call transport_along_meridians (CCTK_PASS_FTOF, hn, i0)
+  
+  call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_t(:,:,hn), -1)
+  call set_boundary (CCTK_PASS_FTOF, hn, qlm_xi_p(:,:,hn), -1)
+  call set_boundary (CCTK_PASS_FTOF, hn, qlm_chi (:,:,hn), +1)
+  
+9999 continue
+end subroutine qlm_killing_transport