c/*@@
c  @file      AHFinder_fun.F
c  @date      April 1998
c  @author    Miguel Alcubierre
c  @desc 
c             Find horizon function.
c  @enddesc 
c@@*/

c Note that including cactus.h will also include AHFinder.h
!#include "cactus.h"
#include "cctk.h"
#include "cctk_parameters.h"
#include "cctk_arguments.h"

      subroutine AHFinder_fun(CCTK_FARGUMENTS)

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

      integer i,j,k
      integer l,m

      CCTK_REAL LEGEN
      CCTK_REAL xp,yp,zp,rp
      CCTK_REAL phi,cost,cosa,sina
      CCTK_REAL zero,half,one,two
      CCTK_REAL pi,halfpi,twopi
      CCTK_REAL aux1,aux2
      INTEGER istat



!     **************************
!     ***   DEFINE NUMBERS   ***
!     **************************

      zero  = 0.0D0
      half  = 0.5D0
      one   = 1.0D0
      two   = 2.0D0

      nx = cctk_lsh(1)
      ny = cctk_lsh(2)
      nz = cctk_lsh(3)

!     *************************************
!     ***   DEFINE  {pi,halfpi,twopi}   ***
!     *************************************

      pi = 3.141592654D0

      halfpi = half*pi
      twopi  = two*pi


!     *********************************
!     ***   FIND HORIZON FUNCTION   ***
!     *********************************

      do k=1,nz
         do j=1,ny
            do i=1,nx

               xp = x(i,j,k) - xc
               yp = y(i,j,k) - yc
               zp = z(i,j,k) - zc

               rp = dsqrt(xp**2 + yp**2 + zp**2)

!              Find phi.

               if (dabs(xp).gt.dabs(yp)) then
                  phi = atan(dabs(yp/xp))
               else if (dabs(xp).lt.dabs(yp)) then
                  phi = halfpi - atan(dabs(xp/yp))
               else
                  phi = 0.25D0*pi
               end if

               if ((xp.eq.zero).and.(yp.eq.zero)) then
                  phi = zero
               else if ((xp.le.zero).and.(yp.ge.zero)) then
                  phi = pi - phi
               else if ((xp.le.zero).and.(yp.le.zero)) then
                  phi = pi + phi
               else if ((xp.ge.zero).and.(yp.le.zero)) then
                  phi = twopi - phi
               end if

!              Monopole term.

               aux1 = c0(0)

!              Axisymmetric terms.

               if (rp.ne.zero) then
                  cost = zp/rp
               else
                  cost = one
               end if

               do l=1+stepz,lmax,1+stepz
                   aux1 = aux1 + c0(l)*LEGEN(l,0,cost)
               end do


!              Non-axisymmetric terms. Notice how the sum over M is first.
!              This will allow me to use the recursion relations to avoid
!              having to start from scratch every time.  Also, I sum over
!              all l's even if I don't want some terms.  This is because
!              in order to use the recursion relations I need all polynomials.

               if (nonaxi) then
                  do m=1,lmax
                     aux2 = dble(m)*phi
                
                     sina = sin(aux2)
                     cosa = cos(aux2)
                     do l=m,lmax
                        aux2 = LEGEN(l,m,cost)
                        aux1 = aux1 + aux2*cc(l,m)*cosa
                        if (.not.refy) then
                           aux1 = aux1 + aux2*cs(l,m)*sina
                        end if
                     end do
                  end do
               end if

               ahfgrid(i,j,k) = rp - aux1

!              Check for NaNs.

               if (ahfgrid(i,j,k).ge.zero) then
               else if (ahfgrid(i,j,k).lt.zero) then
               else

                  if (myproc.eq.0) then
                     write(*,*) 'NaN in ahfgrid at point',i,j,k
                  end if
                  STOP
               end if

            end do
         end do
      end do

!     Synchronize.

      call CCTK_SyncGroup(cctkGH,"ahfinder::ahfindergrid")

!     ***************
!     ***   END   ***
!     ***************

      end subroutine AHFinder_fun