!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Tensor spherical harmonic subroutine
!!$ Version 1.6
!!$
!!$ Computes independent components of tensor spherical harmonics at
!!$    a stereographic coordinate point z
!!$
!!$ l,m - spherical harmonic mode
!!$ south - logical variable to tell which coordinate patch we are in
!!$ ispole - logical variable which tells whether we are at a pole
!!$ v - vector containing 5 independent components of tensor spherical
!!$     harmonics
!!$
!!$ This subroutine uses recursion relations, rather than finite differences
!!$    to calculate derivatives of spherical harmonics


SUBROUTINE harmonics(l,m,z,south,ispole,v)

  INTEGER, INTENT(in) :: l,m
  COMPLEX*16, INTENT(in) ::  z
  COMPLEX*16, DIMENSION(5), INTENT(out) :: v(5)
  LOGICAL, INTENT(in) :: south,ispole

  COMPLEX*16 Y(10),zc,eip,eipb
  REAL*8 C(10)
  INTEGER i

  zc = CONJG(z)

  IF(ispole) THEN
     eip = (1.,0)
  ELSE IF(south) THEN
     eip = zc/SQRT(DBLE(z*zc))
  ELSE
     eip = z/SQRT(DBLE(z*zc))
  END IF

  eipb = CONJG(eip)

  DO i=1,5
     Y(i) = Yster(l,m+i-3,z,south,ispole)
     Y(i+5) = Yster(l-1,m+i-3,z,south,ispole)
  END DO

!!$ Constant coefficients

  C(1) = SQRT(l*(l+1.0)-m*(m+1))
  C(2) = SQRT(l*(l+1.0)-m*(m-1))

  IF (m.GE.l.OR.m.LE.-l) THEN
     C(3) = 0.
     C(4) = 0.
  ELSE
     C(3) = SQRT(l*(l+1.0)-(m+1)*(m+2))
     C(4) = SQRT(l*(l+1.0)-(m-1)*(m-2))
  END IF

  IF(l.EQ.0) THEN
     C(5) = 0.
     C(6) = 0.
     C(7) = 0.
     C(8) = 0.
     C(9) = 0.
     C(10) = 0
  ELSE
     C(5) = SQRT((2*l+1)*(l-m)*(l-m-1)/(2*l-1.0))
     C(6) = SQRT((2*l+1)*(l+m)*(l+m-1)/(2*l-1.0))
     C(7) = SQRT((2*l+1)*(l-m-1)*(l-m-2)/(2*l-1.0))
     C(8) = SQRT((2*l+1)*(l+m+1)*(l+m)/(2*l-1.0))
     C(9) = SQRT((2*l+1)*(l-m+1)*(l-m)/(2*l-1.0))
     C(10) = SQRT((2*l+1)*(l+m-1)*(l+m-2)/(2*l-1.0))
  END IF

!!$ Calculate vector of values

  v(1) = Y(3)

  v(2) = (C(1)*eipb*Y(4) - C(2)*eip*Y(2))/2

  v(3) = (C(1)*C(3)*(eipb**2)*Y(5) + C(2)*C(4)*(eip**2)*Y(1) &
       - (C(1)**2+C(2)**2)*Y(3))/4

  v(4) = (0,-.50)*(C(5)*eipb*Y(9) + C(6)*eip*Y(7))

  v(5) = (0,-.250)*(C(1)*C(7)*(eipb**2)*Y(10) -    &
       C(2)*C(10)*(eip**2)*Y(6) + (C(1)*C(8)-C(2)*C(9))*Y(8))

  RETURN


CONTAINS



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Spherical harmonics in terms of stereographic coordinate values
!!$ Version 1.4
!!$
!!$ Computes Ylm algorithmically at stereographic coordinate point z 
!!$    south is a logical value which is .true. if z is in the 
!!$    southern hemisphere and .false. if z is in the northern hemisphere
!!$ 
!!$ This routine uses a series formula for associated Legendre functions,
!!$    this could be replaced by another subroutine such as plgndr in 
!!$    Numerical Recipes


  COMPLEX*16 FUNCTION Yster(l,m,z,south,ispole)

    IMPLICIT NONE
    INTEGER :: l,m
    COMPLEX*16 :: z
    LOGICAL :: south,ispole

    INTEGER :: lo,hi,n
    REAL*8 :: alpha,Pi,r,rpow
    COMPLEX*16 :: zc,eip

    Pi = 3.141592653589793238460

    zc = CONJG(z)
    r = SQRT(DBLE(z*zc))

    IF(ispole) THEN
       eip = (1,0)
    ELSE IF(south) THEN
       eip = zc/r
    ELSE
       eip = z/r
    END IF

    lo = max0(0,-m)
    hi = min0(l,l-m)

    Yster = 0

    IF(m.GT.l.OR.m.LT.-l.OR.l.LT.0) RETURN

    alpha = (-1)**l * SQRT((2*l+1)/(4*Pi)*fact(l+m)*fact(l-m))*fact(l)/ &
         (1 + r**2)**l

    DO n = lo,hi

       IF(2*(l-n)-m.EQ.0) THEN
          rpow = 1.
       ELSE
          rpow = r**(2*(l-n)-m)
       END IF

       Yster = Yster + (-1)**n * rpow * eip**m / &
            (fact(n)*fact(l-n)*fact(m+n)*fact(l-m-n))

    ENDDO


    Yster = Yster * alpha
    IF(south) Yster = (-1)**(l+m) * Yster

    RETURN

  END FUNCTION Yster



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



!!$ Factorial function
!!$ Version 1.4
!!$
!!$ Calculates factorials for use in the Yster function
!!$
!!$ This is a modified version of a F77 factorial routine used in 
!!$    the Pitt code for 3-D evolution of waves in stereographic coords.


  REAL*8 FUNCTION fact(n)

    IMPLICIT NONE
    INTEGER n,j,ntop
    REAL*8 a(33)
    SAVE ntop,a
    DATA ntop,a(1)/0,1./

    IF(n.LT.0) THEN
       STOP 'negative factorial'
    ELSE IF(n.LE.ntop) THEN
       fact=a(n+1)
    ELSE IF(n.LE.32) THEN
       DO j=ntop+1,n
          a(j+1)=j*a(j)
       END DO
       ntop = n
       fact=a(n+1)
    END IF
    RETURN
  END FUNCTION fact



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


END SUBROUTINE harmonics



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
