gaussp.f 1.52 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
      MODULE m_gaussp
!**************************************************************
!     generates gaussian points to exactly integrate spherical
!     harmonics up to lmax, i.e., (lm|l'm') for l,l'<=lmax
!     number of points = (2*lmax+1)*(lmax+1 + mod(lmax+1,2))
!**************************************************************
      CONTAINS 
      SUBROUTINE gaussp(
     >                  lmax,
     <                  vgauss,wt)

      USE m_grule
      USE m_constants
      IMPLICIT NONE

      INTEGER, INTENT (IN)  :: lmax
      REAL,    INTENT (OUT) :: vgauss(3,*),wt(*) ! points/weights

      INTEGER :: ngpt,nphi,i,j,k
      REAL    :: delphi,phi,rxy
      REAL    :: xx(lmax/2+1),w(lmax/2+1)

!   determine the number of points cos(theta); ngpt always even
      ngpt= lmax+1 + mod(lmax+1,2)
      CALL grule(            ! outputs ngpt/2 points
     >           ngpt,
     <           xx,w)

!  in phi, use nyquist frequency, i.e.,  2*(lmax+1)
      nphi = 2*lmax+1
      delphi = 2.0*pi_const/nphi

      j = 0
      DO i = 1, ngpt/2
         rxy=sqrt(1.0-xx(i)*xx(i))
         DO k=1,nphi
            phi=k*delphi
            j=j+1
            vgauss(1,j) = rxy*cos(phi)
            vgauss(2,j) = rxy*sin(phi)
            vgauss(3,j) = xx(i)
            wt(j) = w(i)*delphi
            j=j+1
            vgauss(1,j) = vgauss(1,j-1)
            vgauss(2,j) = vgauss(2,j-1)
            vgauss(3,j) = -xx(i)
            wt(j) = w(i)*delphi
         ENDDO
      ENDDO

      RETURN
      END SUBROUTINE gaussp
      END MODULE m_gaussp