Remove this experimental routine, that turned out to be a dead end.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@131 2a26948c-0e4f-0410-aee8-f1d3e353619c

1 files changed, 0 insertions, 340 deletions

diff --git a/src/EHFinder_Integrate.F90 b/src/EHFinder_Integrate.F90
deleted file mode 100644
index 04a2a46..0000000
--- a/src/EHFinder_Integrate.F90
+++ /dev/null
@@ -1,340 +0,0 @@
-! Integration over the surface
-! $Header$
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-subroutine EHFinder_Integrate(CCTK_ARGUMENTS)
-
-  use EHFinder_mod
-
-  implicit none
-
-  DECLARE_CCTK_PARAMETERS
-  DECLARE_CCTK_ARGUMENTS
-  DECLARE_CCTK_FUNCTIONS
-
-  CCTK_INT :: i, j, k
-  CCTK_REAL :: idx, idy, idz
-  CCTK_REAL :: al, ar, bl, br, cl, cr
-  CCTK_REAL :: alminus, alplus, blminus, blplus, clminus, clplus
-  CCTK_REAL :: arminus, arplus, brminus, brplus, crminus, crplus
-  CCTK_REAL :: dv, local_area, area, eps, eps_inv
-  CCTK_REAL, dimension(3,3) :: ginv
-  CCTK_REAL :: g2xx, g2xy, g2xz, g2yy, g2yz, g2zz
-  CCTK_REAL :: g22xx, g22xy, g22yy
-  CCTK_REAL, dimension(3) :: norm
-  CCTK_REAL :: dfsquare, dfsquare_i
-  CCTK_REAL, dimension(3) :: vec1, vec2
-  CCTK_REAL :: tmp1, tmp2, tmp3
-  CCTK_REAL :: detg, idetg, sqrtdg2
-  CCTK_REAL :: a, b, c, r1, r2, r3, q
-
-  interface
-    CCTK_REAL function dirac ( x, eps, eps_inv )
-      CCTK_REAL, intent(in) :: x, eps, eps_inv
-    end function dirac
-  end interface
-
-# include "include/physical_part.h"
-
-  idx = half / cctk_delta_space(1)
-  idy = half / cctk_delta_space(2)
-  idz = half / cctk_delta_space(3)
-
-  dv = cctk_delta_space(1)*cctk_delta_space(2)*cctk_delta_space(3)
-  eps = minval(cctk_delta_space)
-  eps_inv = one / eps
-  
-  local_area = zero
-
-  if ( CCTK_EQUALS ( metric_type, 'physical' ) ) then
-    do k = kzl, kzr
-      do j = jyl, jyr
-        do i = ixl, ixr
-          if ( abs(f(i,j,k)) .lt. eps ) then
-
-#include "include/upwind_second2.h"
-
-! Calculate the inverse of the determinant of the spatial 3-metric
-
-            tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2
-            tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k)
-            tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k)
-  
-            detg = gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3
-            idetg = one / detg
-
-! Invert the spatial 3-metric
-
-            ginv(1,1) = tmp1 * idetg
-            ginv(1,2) = tmp2 * idetg
-            ginv(1,3) = tmp3 * idetg
-            ginv(2,1) = ginv(1,2)
-            ginv(2,2) = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
-            ginv(2,3) = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) )*idetg
-            ginv(3,1) = ginv(1,3)
-            ginv(3,2) = ginv(2,3)
-            ginv(3,3) = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
-
-            dfsquare = sqrt ( ginv(1,1) * dfx(i,j,k)**2 + &
-                              two * ginv(1,2) * dfx(i,j,k) * dfy(i,j,k) + &
-                              two * ginv(1,3) * dfx(i,j,k) * dfz(i,j,k) + &
-                              ginv(2,2) * dfy(i,j,k)**2 + &
-                              two * ginv(2,3) * dfy(i,j,k)*dfz(i,j,k) + &
-                              ginv(3,3) *  dfz(i,j,k)**2 )
-!            print*,dfsquare
-
-            dfsquare_i = one / dfsquare
-!            print*,dfsquare_i
-!            dfsquare_i = one / sqrt ( dfx(i,j,k)**2 + &
-!                                      dfy(i,j,k)**2 + &
-!                                      dfz(i,j,k)**2 )
-
-            norm(1) = dfx(i,j,k) * dfsquare_i
-            norm(2) = dfy(i,j,k) * dfsquare_i
-            norm(3) = dfz(i,j,k) * dfsquare_i
-!            print*,norm
-
-!            norm(1) = dfx(i,j,k)
-!            norm(2) = dfy(i,j,k)
-!            norm(3) = dfz(i,j,k)
-
-            call find_orthogonal ( norm, vec1, vec2, ginv, idetg )
-!            print*,norm(1),norm(2),norm(3)
-!            print*,vec1
-!            print*,vec2
-!            print*,norm(1)*vec1(1)+norm(2)*vec1(2)+norm(3)*vec1(3)
-!            print*,norm(1)*vec2(1)+norm(2)*vec2(2)+norm(3)*vec2(3)
-!            print*,gxx(i,j,k)*vec1(1)*vec2(1)+gxy(i,j,k)*vec1(1)*vec2(2)+&
-!                   gxz(i,j,k)*vec1(1)*vec2(3)+gxy(i,j,k)*vec1(2)*vec2(1)+&
-!                   gyy(i,j,k)*vec1(2)*vec2(2)+gyz(i,j,k)*vec1(2)*vec2(3)+&
-!                   gxz(i,j,k)*vec1(3)*vec2(1)+gyz(i,j,k)*vec1(3)*vec2(2)+&
-!                   gzz(i,j,k)*vec1(3)*vec2(3)
-!            print*,gxx(i,j,k)*vec1(1)*vec1(1)+gxy(i,j,k)*vec1(1)*vec1(2)+&
-!                   gxz(i,j,k)*vec1(1)*vec1(3)+gxy(i,j,k)*vec1(2)*vec1(1)+&
-!                   gyy(i,j,k)*vec1(2)*vec1(2)+gyz(i,j,k)*vec1(2)*vec1(3)+&
-!                   gxz(i,j,k)*vec1(3)*vec1(1)+gyz(i,j,k)*vec1(3)*vec1(2)+&
-!                   gzz(i,j,k)*vec1(3)*vec1(3)
-!            print*,gxx(i,j,k)*vec2(1)*vec2(1)+gxy(i,j,k)*vec2(1)*vec2(2)+&
-!                   gxz(i,j,k)*vec2(1)*vec2(3)+gxy(i,j,k)*vec2(2)*vec2(1)+&
-!                   gyy(i,j,k)*vec2(2)*vec2(2)+gyz(i,j,k)*vec2(2)*vec2(3)+&
-!                   gxz(i,j,k)*vec2(3)*vec2(1)+gyz(i,j,k)*vec2(3)*vec2(2)+&
-!                   gzz(i,j,k)*vec2(3)*vec2(3)
- 
-!           print*
-!           print*,gxx(i,j,k),gxy(i,j,k),gxz(i,j,k),&
-!                  gyy(i,j,k),gyz(i,j,k),gzz(i,j,k)
-!           print*
-!            pause
-            g2xx = gxx(i,j,k) - norm(1)**2
-            g2xy = gxy(i,j,k) - norm(1) * norm(2)
-            g2xz = gxz(i,j,k) - norm(1) * norm(3)
-            g2yy = gyy(i,j,k) - norm(2)**2
-            g2yz = gyz(i,j,k) - norm(2) * norm(3)
-            g2zz = gzz(i,j,k) - norm(3)**2
-!            print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz
-
-!            a = - ( g2xx + g2yy + g2zz )
-            b = g2xx * g2yy + g2xx * g2zz + g2yy * g2zz - &
-                                  g2xy**2 - g2xz**2 - g2yz**2
-!            c = g2xy**2 * g2zz + g2xz**2 * g2yy + g2yz**2 * g2xx - &
-!                  g2xx * g2yy * g2zz - two * g2xy * g2xz * g2yz
-
-!            q = - half * ( a + sign ( one, a ) * sqrt ( a**2 - 4 * b ) )
-!            print*,'factors : ',a,b,c,q,b/q,sqrt(b)
-!            print*
-!            sqrtdg2 = sqrt ( sqrt (b) )
-            sqrtdg2 = sqrt (b)
-!            call cubic ( a, b, c, r1, r2, r3 )
-!            print*,'roots : ',r1,r2,r3
-!            print*
-!            pause
-!           print*,g2xx,g2xy,g2xz,g2yy,g2yz,g2zz
-!           print*
-
-!            g22xx = vec1(1)**2*ginv(1,1) + two*vec1(1)*vec1(2)*ginv(1,2) +&
-!                    two*vec1(1)*vec1(3)*ginv(1,3) + vec1(2)**2*ginv(2,2) +&
-!                    two*vec1(2)*vec1(3)*ginv(2,3) + vec1(3)**2 *ginv(3,3)
-!            g22xy = vec1(1)*vec2(1)*ginv(1,1) + vec1(1)*vec2(2)*ginv(1,2) +&
-!                    vec1(1)*vec2(3)*ginv(1,3) + vec1(2)*vec2(1)*ginv(2,1) +&
-!                    vec1(2)*vec2(2)*ginv(2,2) + vec1(2)*vec2(3)*ginv(2,3) +&
-!                    vec1(3)*vec2(1)*ginv(3,1) + vec1(3)*vec2(2)*ginv(3,2) +&
-!                    vec1(3)*vec2(3)*ginv(3,3)
-!            g22yy = vec2(1)**2*ginv(1,1) + two*vec2(1)*vec2(2)*ginv(1,2) +&
-!                    two*vec2(1)*vec2(3)*ginv(1,3) + vec2(2)**2*ginv(2,2) +&
-!                    two*vec2(2)*vec2(3)*ginv(2,3) + vec2(3)**2*ginv(3,3)
-
-!            g22xx = vec1(1)**2*gxx(i,j,k) + two*vec1(1)*vec1(2)*gxy(i,j,k) +&
-!                    two*vec1(1)*vec1(3)*gxz(i,j,k) + vec1(2)**2*gyy(i,j,k) +&
-!                    two*vec1(2)*vec1(3)*gyz(i,j,k) + vec1(3)**2 *gzz(i,j,k)
-!            g22xy = vec1(1)*vec2(1)*gxx(i,j,k) + vec1(1)*vec2(2)*gxy(i,j,k) +&
-!                    vec1(1)*vec2(3)*gxz(i,j,k) + vec1(2)*vec2(1)*gxy(i,j,k) +&
-!                    vec1(2)*vec2(2)*gyy(i,j,k) + vec1(2)*vec2(3)*gyz(i,j,k) +&
-!                    vec1(3)*vec2(1)*gxz(i,j,k) + vec1(3)*vec2(2)*gyz(i,j,k) +&
-!                    vec1(3)*vec2(3)*gzz(i,j,k)
-!            g22yy = vec2(1)**2*gxx(i,j,k) + two*vec2(1)*vec2(2)*gxy(i,j,k) +&
-!                    two*vec2(1)*vec2(3)*gxz(i,j,k) + vec2(2)**2*gyy(i,j,k) +&
-!                    two*vec2(2)*vec2(3)*gyz(i,j,k) + vec2(3)**2*gzz(i,j,k)
-
-            g22xx = vec1(1)**2*g2xx + two*vec1(1)*vec1(2)*g2xy +&
-                    two*vec1(1)*vec1(3)*g2xz + vec1(2)**2*g2yy +&
-                    two*vec1(2)*vec1(3)*g2yz + vec1(3)**2*g2zz
-            g22xy = vec1(1)*vec2(1)*g2xx + vec1(1)*vec2(2)*g2xy +&
-                    vec1(1)*vec2(3)*g2xz + vec1(2)*vec2(1)*g2xy +&
-                    vec1(2)*vec2(2)*g2yy + vec1(2)*vec2(3)*g2yz +&
-                    vec1(3)*vec2(1)*g2xz + vec1(3)*vec2(2)*g2yz +&
-                    vec1(3)*vec2(3)*g2zz
-            g22yy = vec2(1)**2*g2xx + two*vec2(1)*vec2(2)*g2xy +&
-                    two*vec2(1)*vec2(3)*g2xz + vec2(2)**2*g2yy +&
-                    two*vec2(2)*vec2(3)*g2yz + vec2(3)**2*g2zz
-!            print*
-!            print*,g22xx,g22xy,g22yy
-!            print*,g22xx*g22yy-g22xy**2
-
-          tmp1 = g2yy*g2zz - g2yz**2
-          tmp2 = g2xz*g2yz - g2xy*g2zz
-          tmp3 = g2xy*g2yz - g2xz*g2yy        
-!          print*,tmp1,tmp2,tmp3
-
-!          print*,tmp1,tmp2,tmp3
-!
-!         print*, g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3
-!          sqrtdg2 = sqrt ( g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3 )
-!            sqrtdg2 = one / sqrt ( g22xx*g22yy-g22xy**2 )
-!            sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 )
-!            sqrtdg2 = sqrt ( g22xx*g22yy-g22xy**2 )
-
-!           print*,sqrtdg2,g2xx*tmp1 + g2xy*tmp2 + g2xz*tmp3
-!           pause
-
-            local_area = local_area + sqrtdg2 * dirac(f(i,j,k),eps,eps_inv) * &
-                                      dfsquare * sqrt(detg)
-!                         sqrt ( dfx(i,j,k)**2 + dfy(i,j,k)**2 + dfz(i,j,k)**2 )
-          end if
-        end do
-      end do
-    end do
-  end if
-  local_area = dv * local_area
-!  print*,'local_area = ',local_area,sqrtdg2
-
-  call CCTK_ReduceLocScalar ( ierr, cctkGH, -1, sum_handle, &
-                                    local_area, area, CCTK_VARIABLE_REAL )
-
-  eh_area = area
-  print*,'Area = ',eh_area
-
-  return
-end subroutine EHFinder_Integrate
-
-
-subroutine find_orthogonal ( n, v1, v2, ginv, idetg )
-
-  use EHFinder_Constants
-
-  CCTK_REAL, dimension(3) :: n, v1, v2
-  CCTK_REAL, dimension(3) :: v1l
-  CCTK_REAL, dimension(3,3) :: ginv
-  CCTK_REAL :: idetg
-  CCTK_REAL, dimension(1) :: maxn
-  CCTK_REAL :: nv1, nv2
-  CCTK_INT, dimension(1) :: maxp
-  CCTK_INT :: first, second, third
-
-  maxn = maxval (n)
-  maxp = maxloc (n)
-
-!  print*,'Orth'
-!  print*,n
-!  print*,maxn
-!  print*,maxp
-  third = maxp(1)
-  if ( third .eq. 1 ) then
-    first = 2
-    second =3
-  else if ( third .eq. 2 ) then
-    first = 3
-    second = 1
-  else
-    first = 1
-    second = 2
-  end if
-
-!  print*,first,second,third
-  v1l(first) = -n(third)
-  v1l(second) = n(first)
-  v1l(third) = - ( ginv(first,first) * n(first) * v1l(first) + &
-                  ginv(first,second) * n(first) * v1l(second) + &
-                  ginv(second,first) * n(second) * v1l(first) + &
-                  ginv(second,second) * n(second) * v1l(second) + &
-                  ginv(third,first) * n(third) * v1l(first) + &
-                  ginv(third,second) * n(third) * v1l(second) ) / &
-                ( ginv(first,third) * n(first) + &
-                  ginv(second,third) * n(second) + &
-                  ginv(third,third) * n(third) )
-
-  nv1 = sqrt ( ginv(1,1) * v1l(1)**2 + &
-               two * ginv(1,2) * v1l(1) * v1l(2) + &
-               two * ginv(1,3) * v1l(1) * v1l(3) + &
-               ginv(2,2) * v1l(2)**2 + &
-               two * ginv(2,3) * v1l(2) * v1l(3) + &
-               ginv(3,3) * v1l(3)**2 )
-
-  v1l = v1l / nv1
-  v1(1) = ginv(1,1) * v1l(1) + ginv(1,2) * v1l(2) + ginv(1,3) * v1l(3)
-  v1(2) = ginv(2,1) * v1l(1) + ginv(2,2) * v1l(2) + ginv(2,3) * v1l(3)
-  v1(3) = ginv(3,1) * v1l(1) + ginv(3,2) * v1l(2) + ginv(3,3) * v1l(3)
-
-!  print*, v1(1)*v1l(1)+v1(2)*v1l(2)+v1(3)*v1l(3)
-!  print*
-
-  v2(1) = n(2) * v1l(3) - n(3) * v1l(2)
-  v2(2) = n(3) * v1l(1) - n(1) * v1l(3)
-  v2(3) = n(1) * v1l(2) - n(2) * v1l(1)
-  v2 = v2 * sqrt(idetg)
-
-!  print*, v2
-!  print*
-!  print*,v1(1)*n(1)+v1(2)*n(2)+v1(3)*n(3)
-!  print*,v2(1)*n(1)+v2(2)*n(2)+v2(3)*n(3)
-!  print*,v2(1)*v1l(1)+v2(2)*v1l(2)+v2(3)*v1l(3)
-!  print*,'Orth end'
-!  pause 
-  return
-end subroutine find_orthogonal
-  
-CCTK_REAL function dirac ( x, eps, eps_inv )
-
-  use EHFinder_Constants
-
-  CCTK_REAL, intent(in) :: x, eps, eps_inv
-
-  if ( abs(x) .lt. eps ) then
-    dirac = half * ( one + cos(pi*x*eps_inv) ) * eps_inv
-  else
-    dirac = zero
-  end if
-end function dirac
-
-subroutine cubic ( a, b, c, x1, x2, x3 )
-
-  CCTK_REAL, intent(in) :: a, b, c
-  CCTK_REAL, intent(out) :: x1, x2, x3
-  CCTK_REAL :: q, r, theta
-  CCTK_REAL, parameter :: third = 1.0d0 / 3.0d0, &
-                          pi2 = 6.2831853071795864769d0
-
-  q = ( a**2 - 3.0d0 * b ) / 9.0d0
-  r = ( 2.0d0 * a**3 - 9.0d0 * a * b + 27.0d0 * c ) / 54.0d0
-  theta = acos ( r / sqrt( q**3 ) )
-
-  x1 = - 2.0d0 * sqrt( q ) * cos( third * theta ) - third * a
-  x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta + pi2 ) ) - third * a 
-  x2 = - 2.0d0 * sqrt( q ) * cos( third * ( theta - pi2 ) ) - third * a 
-
-  print*,a,b,c
-  print*,q,r,q**3-r**2
-!  pause
-end subroutine cubic