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

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"

      subroutine AHFinder_fun(CCTK_ARGUMENTS)

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      integer i,j,k
      integer l,m

      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
      CCTK_REAL local_legen,local_legen_old,local_legen_old_old 
      CCTK_REAL cc_cosa(lmax,lmax),cs_sina(lmax,lmax)
      CCTK_REAL dsqrt_factor(lmax,lmax),dsqrt_factor2(lmax)
      CCTK_REAL lm_factor1(0:lmax,0:lmax),lm_factor2(0:lmax,0:lmax)
      CCTK_REAL lm_factor3(lmax)
      CCTK_REAL sinphi,cosphi,sinphi1,cosphi1,sinphi2,cosphi2
      CCTK_REAL legen_base


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

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

      pi     = acos(-one)
      halfpi = half*pi
      twopi  = two*pi


!     *******************************************
!     ***   FACTORS FOR SPHERICAL HARMONICS   ***
!     *******************************************

!     Compute the normalization factors only once for all (l,m)
!
!     For M=0 I use the normalization factor:
!
!     sqrt(2*L + 1)
!
!     For M non-zero I use the normalization factor:
!
!     sqrt( 2 (2*L + 1) (L-M)! / (L+M)! )
!
!     The extra factor of sqrt(2) is there because I
!     use a basis of sines and cosines and not the
!     standard complex exponential.
!
!     With this normalization, I ensure that my basis
!     functions f _lm are such that:
!
!     /
!     | f_lm f_l;mp sin(theta) dtheta dphi  =  4 pi delta_mm' delta_ll'
!     /
!
!     Notice the extra factor of 4 pi.  This is there because
!     for a sphere, I want the (0,0) coefficient to correspond
!     to the radius.
!
!     Also compute the multiplying factors in the recursion relations

!     Initialize.

      lm_factor1 = zero
      lm_factor2 = zero
      lm_factor3 = zero

      dsqrt_factor  = zero
      dsqrt_factor2 = zero

!     M=0.

      do l=1,lmax

         dsqrt_factor2(l) = c0(l)*sqrt(dble(2*l+1))

         lm_factor1(l,0) = dble(2*l-1)/dble(l)
         lm_factor2(l,0) = dble(l-1)/dble(l)
         lm_factor3(l)   = dble(2*l-1)

      end do

!     M > 0.

      if (nonaxi) then

         do l=1,lmax

            dsqrt_factor(l,1) = sqrt(2.0D0*dble(2*l+1)/dble(l*(l+1)))

            do m=2,l
               dsqrt_factor(l,m) = dsqrt_factor(l,m-1)
     .                           / sqrt(dble((l-m+1)*(l+m)))
            end do

         end do

         do l=2,lmax

            lm_factor1(l,1) = dble(2*l-1)/dble(l-1)
            lm_factor2(l,1) = dble(l)/dble(l-1)

            do m=2,l-1
               lm_factor1(l,m) = dble(2*l-1)/dble(l-m)
               lm_factor2(l,m) = dble(l+m-1)/dble(l-m)
            end do

         end do

      end if


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

!     Loop over j.

      do j=1,ny

!        Find yp.

         yp = y(1,j,1) - yc

!        Loop over i.

         do i=1,nx

!           Find xp.

            xp = x(i,1,1) - xc

!           Find sines and cosines of phi.

            if (nonaxi) then

!              Find phi.

               if (yp.ne.zero) then
                  phi = atan2(yp,xp)
               else
                  phi = zero
               end if

!              Compute the expensive sines and cosines only once
!              for all (l,m) using the recursion relations:
!
!              cos(m*phi) = 2*cos((m-1)*phi)*cos(phi)-cos((m-2)*phi)
!              sin(m*phi) = 2*sin((m-1)*phi)*cos(phi)-sin((m-2)*phi)

               cosphi  = cos(phi)
               cosphi1 = one
               cosphi2 = cosphi

               if (.not.refy) then
                  sinphi  = sin(phi)
                  sinphi1 = zero
                  sinphi2 = sinphi
               end if

!              M=1.

               do l=1,lmax

                  cc_cosa(l,1) = cc(l,1)*cosphi*dsqrt_factor(l,1)

                  if (.not.refy) then
                     cs_sina(l,1) = cs(l,1)*sinphi*dsqrt_factor(l,1)
                  end if

               end do

!              M>1.

               do m=2,lmax

                  cosa = 2.0D0*cosphi2*cosphi - cosphi1

                  cosphi1 = cosphi2
                  cosphi2 = cosa

                  do l=m,lmax
                     cc_cosa(l,m) = cc(l,m)*cosa*dsqrt_factor(l,m)
                  end do 

                  if (.not.refy) then

                     sina = 2.0D0*sinphi2*cosphi - sinphi1

                     sinphi1 = sinphi2
                     sinphi2 = sina

                     do l=m,lmax
                        cs_sina(l,m) = cs(l,m)*sina*dsqrt_factor(l,m)
                     end do

                  end if

              end do

            end if

!           Loop over k.

            do k=1,nz

!              Find zp.

               zp = z(i,j,k) - zc

!              Find rp.

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

!              Monopole term.

               aux1 = c0(0)

!              Axisymmetric terms.

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

!              Compute the contribution from M=0, L=1,Lmax using
!              the recursion relations from Numerical Recipes.

               if (lmax.gt.0) then

!                 L=1

                  legen_base = one
                  local_legen_old_old = one
                  local_legen_old = cost

                  aux1 = aux1 + dsqrt_factor2(1)*local_legen_old

!                 L=2,lmax

                  do l=2,lmax

                     local_legen = cost*lm_factor1(l,0)*local_legen_old
     .                                - lm_factor2(l,0)*local_legen_old_old
                     local_legen_old_old = local_legen_old
                     local_legen_old = local_legen

                     aux1 = aux1 + dsqrt_factor2(l)*local_legen

                  end do
!                  
!                 Non-axisymmetric terms.
!
                  if (nonaxi) then

                     aux2 = sqrt((one-cost)*(one+cost))

!                    This general loop can only be used until l=lmax-2 since
!                    the recursion relations requires 2 starting values.

                     do m=1,lmax-2
! L=M
                        legen_base = -legen_base*aux2*lm_factor3(m)     

                        aux1 = aux1 + legen_base*cc_cosa(m,m)

                        if (.not.refy) then
                           aux1 = aux1 + legen_base*cs_sina(m,m)
                        end if

                        local_legen_old_old = legen_base
! L=M+1
                        local_legen_old = cost*lm_factor3(m+1)         
     .                                  *legen_base

                        aux1 = aux1 + local_legen_old*cc_cosa(m+1,m)

                        if (.not.refy) then
                           aux1 = aux1 + local_legen_old*cs_sina(m+1,m)
                        end if

! L=M+2,lmax
                        do l=m+2,lmax                             

                           local_legen = cost*lm_factor1(l,m)*local_legen_old
     .                                      - lm_factor2(l,m)*local_legen_old_old
                           local_legen_old_old = local_legen_old
                           local_legen_old = local_legen

                           aux1 = aux1 + local_legen*cc_cosa(l,m)

                           if (.not.refy) then
                              aux1 = aux1 + local_legen*cs_sina(l,m)
                           end if

                        end do
                     end do

!                    Then do M=lmax-1

                     if (lmax.gt.1) then

! L=lmax-1
                        m=lmax-1                                   
                        legen_base = -legen_base*aux2*lm_factor3(m)   
                        aux1 = aux1 + legen_base*cc_cosa(m,m)

                        if (.not.refy) then
                           aux1 = aux1 + legen_base*cs_sina(m,m)
                        end if

! L=lmax
                        local_legen = cost*lm_factor3(lmax)*   
     .                                            legen_base
                        aux1 = aux1 + local_legen*cc_cosa(lmax,m)

                        if (.not.refy) then
                           aux1 = aux1 + local_legen*cs_sina(lmax,m)
                        end if

                     end if

!                    And finally do M=lmax.

! L=lmax
                     m=lmax                                         

                     legen_base = -legen_base*aux2*lm_factor3(m)
                     aux1 = aux1 + legen_base*cc_cosa(m,m)

                     if (.not.refy) then
                        aux1 = aux1 + legen_base*cs_sina(m,m)
                     end if

                  end if

               end if

!              Find horizon function.

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

            end do
         end do
      end do


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

      end subroutine AHFinder_fun

