Line data    Source code

       1             :       MODULE m_gaussp
       2             : !**************************************************************
       3             : !     generates gaussian points to exactly integrate spherical
       4             : !     harmonics up to lmax, i.e., (lm|l'm') for l,l'<=lmax
       5             : !     number of points = (2*lmax+1)*(lmax+1 + mod(lmax+1,2))
       6             : !**************************************************************
       7             :       CONTAINS 
       8        1104 :       SUBROUTINE gaussp(
       9             :      >                  lmax,
      10             :      <                  vgauss,wt)
      11             : 
      12             :       USE m_grule
      13             :       USE m_constants
      14             :       IMPLICIT NONE
      15             : 
      16             :       INTEGER, INTENT (IN)  :: lmax
      17             :       REAL,    INTENT (OUT) :: vgauss(3,*),wt(*) ! points/weights
      18             : 
      19             :       INTEGER :: ngpt,nphi,i,j,k
      20             :       REAL    :: delphi,phi,rxy
      21        1104 :       REAL    :: xx(lmax/2+1),w(lmax/2+1)
      22             : 
      23             : !   determine the number of points cos(theta); ngpt always even
      24        1104 :       ngpt= lmax+1 + mod(lmax+1,2)
      25             :       CALL grule(            ! outputs ngpt/2 points
      26             :      >           ngpt,
      27        1104 :      <           xx,w)
      28             : 
      29             : !  in phi, use nyquist frequency, i.e.,  2*(lmax+1)
      30        1104 :       nphi = 2*lmax+1
      31        1104 :       delphi = 2.0*pi_const/nphi
      32             : 
      33        1104 :       j = 0
      34        6868 :       DO i = 1, ngpt/2
      35        5764 :          rxy=sqrt(1.0-xx(i)*xx(i))
      36      111928 :          DO k=1,nphi
      37      105060 :             phi=k*delphi
      38      105060 :             j=j+1
      39      105060 :             vgauss(1,j) = rxy*cos(phi)
      40      105060 :             vgauss(2,j) = rxy*sin(phi)
      41      105060 :             vgauss(3,j) = xx(i)
      42      105060 :             wt(j) = w(i)*delphi
      43      105060 :             j=j+1
      44      105060 :             vgauss(1,j) = vgauss(1,j-1)
      45      105060 :             vgauss(2,j) = vgauss(2,j-1)
      46      105060 :             vgauss(3,j) = -xx(i)
      47      110824 :             wt(j) = w(i)*delphi
      48             :          ENDDO
      49             :       ENDDO
      50             : 
      51        1104 :       RETURN
      52             :       END SUBROUTINE gaussp
      53             :       END MODULE m_gaussp