


















































































































 


































































































































































 



































































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                





































module NullEvol_sYlm







































  implicit none
  REAL*8 :: dummy_var
  INTEGER*4, parameter :: wp = kind(dummy_var)

  REAL*8, private, parameter ::&
          pi = 3.1415926535897932384626433832795_wp
  logical, save :: FirstTimeCalled  = .true.
  REAL*8, dimension(0:20), save :: fact_
  COMPLEX*16, dimension(:,:,:), allocatable, save :: zp
  COMPLEX*16, dimension(:,:), allocatable,  save :: tm
  REAL*8, dimension(:,:,:),  allocatable, save :: P, P_d

contains

  subroutine rsYlm(s, l, m, nnx, nny, zeta, out)
    implicit none
    INTEGER*4, intent(in) :: s, l, m, nnx, nny
    COMPLEX*16, dimension(nnx,nny), intent(in) :: zeta
    COMPLEX*16, dimension(nnx,nny,2), intent(out) :: out

    COMPLEX*16, dimension(nnx,nny,2) :: tmp
    COMPLEX*16 :: ii = (0, 1.0)

    select case (m)
    case(:-1)
      call sYlm(s,l,-m,nnx,nny,zeta,out)
      call sYlm(s,l,m,nnx,nny,zeta,tmp)
      out = -ii/ dsqrt(2.0d0)*(out -(-1)**m * tmp)
    case(0)
      call sYlm(s,l,m,nnx,nny,zeta,out)
    case(1:)
      call sYlm(s,l,m,nnx,nny,zeta,out)
      call sYlm(s,l,-m,nnx,nny,zeta,tmp)
      out = 1.0d0/ dsqrt(2.0d0)*(out +(-1)**m * tmp)
    case DEFAULT
      call CCTK_Warn(0,79,"NullEvol_sYlm.F90","NullEvolve","bad m")
    end select
    call sYlm(s,l,m,nnx,nny,zeta,tmp)
  end subroutine rsYlm

  subroutine sYlm(s, l, m, nnx, nny, zeta, F_out)
    implicit none
    INTEGER*4, intent(in) :: s, l, m, nnx, nny
    COMPLEX*16, dimension(nnx,nny), intent(in) :: zeta
    COMPLEX*16, dimension(nnx,nny,2), intent(out) :: F_out

    integer :: ss, mm, n

    if (l < 0 ) then
       call CCTK_Warn(0,93,"NullEvol_sYlm.F90","NullEvolve", "sYlm ERROR: l < 0 ")
    endif
    if (abs(m) > l)  then
      call CCTK_Warn(0,96,"NullEvol_sYlm.F90","NullEvolve", "sYlm ERROR: abs(m) > l")
    endif
    if (abs(s) > l)  then
      call CCTK_Warn(0,99,"NullEvol_sYlm.F90","NullEvolve", "sYlm ERROR: abs(s) > l")
    endif

    if (FirstTimeCalled) then
      FirstTimeCalled = .false.
      fact_(0) = 1.0_wp
      do n = 1, 20
        fact_(n) = n*fact_(n-1)
      end do
      allocate(zp(nnx,nny,0:20), P(nnx,nny,0:20),&
           P_d(nnx,nny,0:20), tm(nnx,nny))

      zp(:,:,0) = 1.0_wp   
      P(:,:,0)  = 1.0_wp
      P_d(:,:,0)  = 1.0_wp
      do n = 1, 20
        zp(:,:,n) = zeta**n
        P(:,:,n) = (1.0_wp + abs(zeta*conjg(zeta)))**n
        P_d(:,:,n) = 1.0_wp / (1.0_wp + abs(zeta*conjg(zeta)))**n
      end do
    endif 

    ss = s
    mm = m
    if (ss< 0) ss= -ss
    if (mm< 0) mm= -mm

    if ( (l + ss) .le. 20 .and. (l+mm) .le. 20 ) then

      call sYlm_north_opt(s, l, m, nnx, nny,F_out(:,:,1))
      call sYlm_south_opt(s, l, m, nnx, nny,F_out(:,:,2))
    else

      call sYlm_north(s, l, m, nnx, nny, zeta, F_out(:,:,1))
      call sYlm_south(s, l, m, nnx, nny, zeta, F_out(:,:,2))
    endif
  end subroutine sYlm

  subroutine sYlm_north_opt(s, l, m, nnx,nny, north)
    implicit none  
    INTEGER*4, intent(in) :: s, l, m, nnx, nny
    COMPLEX*16, dimension(nnx,nny), intent(out) :: north
    INTEGER*4 :: n

    north = 0.0

    tm = 0.0 
    do n  = max(0, s+m), min(l+s, l+m)
      tm = tm + (-1)**n * zp(:,:,n) * conjg(zp(:,:,n-m-s)) / &
              (fact_(n)*fact_(l+m-n)*fact_(l+s-n)*fact_(n-m-s))
    end do
    north = dsqrt(1.0_wp + 2.0_wp * l)*P_d(:,:,l)  * &
            dsqrt(fact_(l-m)*fact_(l+m)*fact_(l-s)*fact_(l+s)) * tm /&
            (2.0_wp * dsqrt(pi))
  end subroutine sYlm_north_opt


  subroutine sYlm_north(s, l, m, nnx, nny, zeta, north)
    implicit none  
    INTEGER*4, intent(in) :: s, l, m, nnx, nny
    COMPLEX*16, dimension(nnx,nny), intent(in) :: zeta
    COMPLEX*16, dimension(nnx,nny), intent(out) :: north

    COMPLEX*16, dimension(nnx,nny) :: zetabar, tmp
    INTEGER*4 :: n

    zetabar = conjg(zeta)

    if (abs(m) > l)  then
      call CCTK_Warn(0,168,"NullEvol_sYlm.F90","NullEvolve", "sYlm_north ERROR: abs(m) > l")
    endif
    if (abs(s) > l)  then
      call CCTK_Warn(0,171,"NullEvol_sYlm.F90","NullEvolve", "sYlm_north ERROR: abs(s) > l")
    endif

    north = 0.0

    tmp = 0.0 
    do n  = max(0, s+m), min(l+s, l+m)
      tmp = tmp + (-1)**n * zeta**n * zetabar**(n-m-s) / &
              (fact(n)*fact(l+m-n)*fact(l+s-n)*fact(n-m-s))
    end do
    north = dsqrt(1.0_wp + 2.0_wp * l) / (  &
                      (1.0_wp + zeta * zetabar )**l ) * &
            dsqrt(fact(l-m)*fact(l+m)*fact(l-s)*fact(l+s)) * tmp /&
            (2.0_wp * dsqrt(pi))
  end subroutine sYlm_north 

  subroutine sYlm_south_opt(s, l, m, nnx, nny, south)
    implicit none  
    INTEGER*4, intent(in) :: s, l, m, nnx, nny
    COMPLEX*16, dimension(nnx,nny), intent(out) :: south
    INTEGER*4 :: n

    south = 0.0

    tm = 0.0 
    do n  = max(0, s+m), min(l+s, l+m)
      tm = tm  + (-1)**n * zp(:,:,l+s-n) * conjg(zp(:,:,l+m-n)) /&
              (fact_(n)*fact_(l+m-n)*fact_(l+s-n)*fact_(n-m-s))
    end do

    south = (-1)**s * dsqrt(1.0_wp + 2.0_wp * l) * P_d(:,:,l) * &
            dsqrt(fact_(l-m)*fact_(l+m)*fact_(l-s)*fact_(l+s)) * tm /&
            (2.0_wp * dsqrt(pi))


  end subroutine sYlm_south_opt


  subroutine sYlm_south(s, l, m, nnx,nny, zeta, south)
    implicit none  
    INTEGER*4, intent(in) :: s, l, m, nnx,nny
    COMPLEX*16, dimension(nnx,nny), intent(in) :: zeta
    COMPLEX*16, dimension(nnx,nny), intent(out) :: south

    COMPLEX*16, dimension(nnx,nny) :: zetabar, tmp
    INTEGER*4 :: n

    zetabar = conjg(zeta)

    if (abs(m) > l)  then
      call CCTK_Warn(0,221,"NullEvol_sYlm.F90","NullEvolve", "sYlm_south ERROR: abs(m) > l")
    endif
    if (abs(s) > l)  then
      call CCTK_Warn(0,224,"NullEvol_sYlm.F90","NullEvolve", "sYlm_south ERROR: abs(s) > l")
    endif

    south = 0.0

    tmp = 0.0 
    do n  = max(0, s+m), min(l+s, l+m)
      tmp = tmp  + (-1)**n * zeta**(l+s-n) * zetabar**(l+m-n) / &
              (fact(n)*fact(l+m-n)*fact(l+s-n)*fact(n-m-s))
    end do

    south = (-1)**s * dsqrt(1.0_wp + 2.0_wp * l) / &
                      (1.0_wp + zeta * zetabar )**l * &
            dsqrt(fact(l-m)*fact(l+m)*fact(l-s)*fact(l+s)) * tmp /&
            (2.0_wp * dsqrt(pi))

  end subroutine sYlm_south 

  REAL*8 function fact(j) 
    implicit none
    INTEGER*4, intent(in) :: j

    INTEGER*4 :: i
    REAL*8 :: res

    if (j < 0 ) then 
      call CCTK_Warn(0,250,"NullEvol_sYlm.F90","NullEvolve", "fact called with a negative argument")
    endif
    if (j > 180 ) then
      call CCTK_Warn(0,253,"NullEvol_sYlm.F90","NullEvolve", "fact cannot caclulate n! for n > 180")
    endif
    if (j == 0 .OR. j == 1) then
      fact = 1.0_wp
    else
      res = 1.0_wp
      do i = 2, j
        res = res * i
      end do
      fact = res
    end if
    return
  end function fact 

end module NullEvol_sYlm
