SUBROUTINE anmomthe(ttheta,pttheta,ntheta,nthetmax,ltheta,angtot,thetaj,debug)
USE Numeric_Kinds_Module
USE FileUnits_Module
IMPLICIT NONE
!==============================================================================
!     SUBROUTINE angmom calculates the DVR angular momentum matrices
!     j_{theta}^2 and puts the result in angtot
!     for use in calculating the full 2-d hamiltonian
!     thetaj is a scratch array for intermediate results.
!     note that the full angular momentum array is calculated rather than
!     just calculating the elements for the truncated DVR. Some increase
!     in speed may be achieved by only calculating the elements of angtot
!     for those angles that have gaussians on them (as defined in npttheta).
!
!     ttheta=dvr transformation matrix from tdvf2b
!     pttheta=dvr points
!     ntheta=number of angular basis functions
!     nthetmax=maximum number of angular basis functions
!     angtot=transformed j^2 matrix
!     thetaj=scratch array of dimension nthetmax
!==============================================================================
INTEGER  :: i, j, k, ntheta, jtheta, ltheta, nthetmax
REAL(dp) :: temp
REAL(dp) :: ttheta(nthetmax,ntheta), thetaj(ntheta)
REAL(dp) :: pttheta(ntheta), angtot(nthetmax,ntheta)
LOGICAL  :: debug

!==============================================================================
! Begin by calculating FBR j_{theta}^2 (which is diagonal) and placing
!   results in thetaj

DO j=1,ntheta
   jtheta=j-1+ltheta
   thetaj(j)=16.d0*(jtheta*(jtheta+1))
ENDDO

!==============================================================================
! Now transform the j^2 matrix using ttheta

DO i=1,ntheta
   DO j=1,ntheta
      temp=0.0d0
      DO k=1,ntheta
         temp=temp+ttheta(i,k)*thetaj(k)*ttheta(j,k)
      ENDDO
      angtot(i,j)=temp
  ENDDO
ENDDO

RETURN
END