Fixed an bug in an error criteria that was making the finder think that

it has found a horizon when in fact it was giving up becuase it was too
close to the origin.

While I was at it, I cleaned up a bit the routine as well.


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

1 files changed, 41 insertions, 76 deletions

diff --git a/src/AHFinder_flow.F b/src/AHFinder_flow.F
index 465bfb1..c48582e 100644
--- a/src/AHFinder_flow.F
+++ b/src/AHFinder_flow.F
@@ -27,15 +27,14 @@
 
       logical status,new,update
 
-      integer itmax
       integer i,j,l,m
+      integer itmax
 
       CCTK_REAL rmx,ahftol
       CCTK_REAL A,B
       CCTK_REAL ahfsum,ahfdiff,ahfdiff_old
-      CCTK_REAL maxchange,minchange
       CCTK_REAL reldiff
-      CCTK_REAL spher
+      CCTK_REAL maxchange,minchange
       CCTK_REAL zero,one
       CCTK_REAL dd,aux
 
@@ -43,17 +42,37 @@
 
 !     Description of variables:
 !
-!     i,l,m     Counters.
+!     new         Flag to indicate that stepsize has changed.
+!     update      Flag to indicate if coefficients will be updated.
+!
+!     i,j,l,m     Counters.
+!
+!     itmax       Maximum number of iterations.
+!
+!     ahftol      Tolerance used to decide if we have converged.
 !
-!     ITMAX     Maximum number of iterations.
+!     A,B         Flow parameters.
 !
-!     A,B       Flow parameters.
+!     hw          Weigth of `H' flow (see Carsten's paper).
+!     cw          Weigth of `C' flow (see Carsten's paper).
+!     nw          Weigth of `N' flow (see Carsten's paper).
 !
-!     hw        Weigth of `H' flow (see Carsten's paper).
-!     cw        Weigth of `C' flow (see Carsten's paper).
-!     nw        Weigth of `N' flow (see Carsten's paper).
+!     ahfsum      Square root of the sum of the squares of the
+!                 expansion coefficients.
 !
-!     ahftol   Tolerance used to decide if we have converged.
+!     ahfdiff     Square root of the sum of the squares of the
+!                 change in the coefficients.
+!
+!     rediff      Relative difference between having done one large
+!                 step and two small ones.
+!
+!     maxchange   A value of "reldiff" larger than this forces a
+!                 reduction of the stepsize.
+!
+!     minchange   A value of "reldiff" smaller than this forces an
+!                 increase of the stepsize.
+!
+!     dd          Maximum of (dx,dy,dz).
 
 
 !     ***************************
@@ -63,7 +82,7 @@
       zero = 0.0D0
       one  = 1.0D0
 
-      dd = min(dx,dy,dz)
+      dd = max(dx,dy,dz)
 
 
 !     ***************************
@@ -182,6 +201,16 @@
 !              We can still try to reduce the stepsize and start
 !              the iteration again.  We then return the coefficients
 !              to their original values and, and do not update them.
+!
+!              Since I will not update the coefficients, I need to
+!              assign values to "afhsum" and "ahfdiff". In order
+!              to force a reduction of the step size, I assign
+!              a very large number to "ahfdiff", which will make
+!              "reldiff" close to 1.0, far too large.  I also
+!              assign a value of 1.0 to "ahfsum", so it is much
+!              smaller than "ahfdiff".  This is important in order
+!              to prevent the algorithm from thinking we are done
+!              (see below for termination criteria).
 
             else if (interror) then
 
@@ -192,7 +221,7 @@
                cs = cs_old
 
                ahfsum   = one
-               ahfdiff  = - ahfdiff_old
+               ahfdiff  = 1.0D10
 
                update = .false.
 
@@ -567,67 +596,3 @@
       end subroutine AHFinder_flow
 
 
-
-
-
-
-
-
-
-
-
-
-      CCTK_REAL function spher(theta,phi)
-
-!     Function to reconstruct a function pointwise
-!     from the expansion coefficients {c0,cc,cs}.
-
-      use AHFinder_dat
-
-      implicit none
-
-      integer l,m
-
-      CCTK_REAL LEGEN
-      CCTK_REAL theta,phi
-      CCTK_REAL F,cost,cosa,sina
-      CCTK_REAL aux
-
-
-!     ******************************
-!     ***   Axisymmetric terms   ***
-!     ******************************
-
-      F = 0.0d0
-
-      cost = cos(theta)
-
-      do l=0,lmax
-         F = F + c0(l)*legen(l,0,cost)
-      end do
-
-
-!     **********************************
-!     ***   Non-axisymmetric terms   ***
-!     **********************************
-
-      do m=1,lmax
-
-         aux  = m*phi
-         sina = sin(aux)
-         cosa = cos(aux)
-
-         do l=m,lmax
-            F = F + LEGEN(l,m,cost)*(cc(l,m)*cosa + cs(l,m)*sina)
-         end do
-
-      end do
-
-      spher = F
-
-
-!     ***************
-!     ***   END   ***
-!     ***************
-
-      end function spher