Fixing very subtle bug when finding distance between horizons in the case

when 1 did not exist.  The file was not being opened when it should have,
and the code could crash later when looking for it.  I hope it is fixed
now.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@263 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 39 insertions, 22 deletions

diff --git a/src/AHFinder_dis.F b/src/AHFinder_dis.F
index 3fdcca7..51451d0 100644
--- a/src/AHFinder_dis.F
+++ b/src/AHFinder_dis.F
@@ -83,6 +83,8 @@
          write(*,*)
          write(*,*) 'Neither horizon 1 nor horizon 2 exist,'
          write(*,*) 'giving up.'
+         write(*,*)
+
          d12 = zero
          goto 100
 
@@ -92,6 +94,7 @@
 
          write(*,*)
          write(*,*) 'Only horizon 1 exists.'
+         write(*,*)
 
 !        Check symmetries.
 
@@ -147,8 +150,8 @@
 
             if (CCTK_Equals(bitant_plane,"xy").eq.1) then
 
-               write(*,*) 'Bitant symmetry, xy plane: calculating'
-               write(*,*) 'twice the distance from horizon 1 to xy plane.'
+               write(*,*) 'Bitant symmetry, xy plane: calculating twice'
+               write(*,*) 'the distance from horizon 1 to xy plane.'
 
                factor = two
 
@@ -158,8 +161,8 @@
 
             else if (CCTK_Equals(bitant_plane,"xz").eq.1) then
 
-               write(*,*) 'Bitant symmetry, xz plane: calculating'
-               write(*,*) 'twice the distance from horizon 1 to xz plane.'
+               write(*,*) 'Bitant symmetry, xz plane: calculating twice'
+               write(*,*) 'the distance from horizon 1 to xz plane.'
 
                factor = two
 
@@ -169,8 +172,8 @@
 
             else if (CCTK_Equals(bitant_plane,"yz").eq.1) then
 
-               write(*,*) 'Bitant symmetry, yz plane: calculating'
-               write(*,*) 'twice the distance from horizon 1 to yz plane.'
+               write(*,*) 'Bitant symmetry, yz plane: calculating twice'
+               write(*,*) 'the distance from horizon 1 to yz plane.'
 
                factor = two
 
@@ -196,6 +199,7 @@
 
          write(*,*)
          write(*,*) 'Only horizon 2 exists.'
+         write(*,*)
 
 !        Check symmetries.
 
@@ -251,8 +255,8 @@
 
             if (CCTK_Equals(bitant_plane,"xy").eq.1) then
 
-               write(*,*) 'Bitant symmetry, xy plane: calculating'
-               write(*,*) 'twice the distance from horizon 2 to xy plane.'
+               write(*,*) 'Bitant symmetry, xy plane: calculating twice'
+               write(*,*) 'the distance from horizon 2 to xy plane.'
 
                factor = two
 
@@ -262,8 +266,8 @@
 
             else if (CCTK_Equals(bitant_plane,"xz").eq.1) then
 
-               write(*,*) 'Bitant symmetry, xz plane: calculating'
-               write(*,*) 'twice the distance from horizon 2 to xz plane.'
+               write(*,*) 'Bitant symmetry, xz plane: calculating twice'
+               write(*,*) 'the distance from horizon 2 to xz plane.'
 
                factor = two
 
@@ -273,8 +277,8 @@
 
             else if (CCTK_Equals(bitant_plane,"yz").eq.1) then
 
-               write(*,*) 'Bitant symmetry, yz plane: calculating'
-               write(*,*) 'twice the distance from horizon 2 to yz plane.'
+               write(*,*) 'Bitant symmetry, yz plane: calculating twice'
+               write(*,*) 'the distance from horizon 2 to yz plane.'
 
                factor = two
 
@@ -284,6 +288,14 @@
 
             end if
 
+         else
+
+            write(*,*) 'No grid symmetries: no distance to calculate.'
+            write(*,*) 'giving up.'
+
+            d12 = zero
+            goto 100
+
          end if
 
       end if
@@ -307,7 +319,12 @@
       if (r12.ne.zero) then
          cost = z12/r12
       else
-         return
+         write(*,*)
+         write(*,*) 'Coordinate distance is zero,'
+         write(*,*) 'giving up.'
+         write(*,*)
+         d12 = zero
+         goto 100
       end if
 
       sint = sqrt(one - cost**2)
@@ -470,21 +487,21 @@
       coord_system_handle = -1
 
       if (interpolation_order.eq.1) then
-        call CCTK_InterpHandle(interp_handle, "first-order uniform cartesian")
-      else if (interpolation_order .eq. 2) then
-        call CCTK_InterpHandle(interp_handle, "second-order uniform cartesian")
-      else if (interpolation_order .eq. 3) then
-        call CCTK_InterpHandle(interp_handle, "third-order uniform cartesian")
+        call CCTK_InterpHandle(interp_handle,"first-order uniform cartesian")
+      else if (interpolation_order.eq.2) then
+        call CCTK_InterpHandle(interp_handle,"second-order uniform cartesian")
+      else if (interpolation_order.eq.3) then
+        call CCTK_InterpHandle(interp_handle,"third-order uniform cartesian")
       endif
 
-      if (interp_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
+      if (interp_handle.lt.0) then
+        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system")
       endif
 
       call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
 
-      if (coord_system_handle .lt. 0) then
-        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
+      if (coord_system_handle.lt.0) then
+        call CCTK_WARN(0,"Cannot get handle for interpolation")
       endif
 
 !     Interpolation.