SUBROUTINE rotational(j, l, rotdim, rot, io, frame)

!=========================================================================================
! Written by Jeff Crawford
!
! This routine calculates the body fixed or space fixed rotational eigenfunctions in 
!    Jacobi coordinates and maps them to APH coordinates.
!
! Need to update for jtotal .ne. zero
!=========================================================================================
USE Numeric_Kinds_Module
USE CommonInfo_Module
USE Jacobi_Module
USE APH_Module
USE QuantumNumber_Module
USE SurfaceJacobi_Module
USE Factrl_Module
IMPLICIT NONE

!=========================================================================================
!     I N P U T 
CHARACTER(LEN=1)    :: io
CHARACTER(LEN=2)    :: frame
INTEGER,INTENT(IN)  :: j, rotdim, l

!=========================================================================================
!     O U T P U T 
REAL(dp),INTENT(OUT) :: rot(rotdim)
!COMPLEX(dp),INTENT(OUT) :: rot(rotdim)

!=========================================================================================
!     I N T E R N A L S 
INTEGER  :: irho, itheta, ichi, ispt, ipt, iomega, omega, rhostart, rhoend
REAL(dp) :: cgcoef, gcoord, realpart, wigner_d, costheta, beta, coeff
REAL(dp) :: shfac
!=========================================================================================
!     A L L O C A T A B L E
REAL(dp),ALLOCATABLE :: plmn(:)

!=========================================================================================
!     F U N C T I O N S
INTEGER  :: aphindex
REAL(dp) :: cleb
!===============================================================================
! Calculate the assoaiated legendre polynomial a function of capital theta
! Need to add option to be space fixed or body fixed.

  IF (io.eq.'I') THEN
    rhostart=1
    rhoend=nrho
  ELSEIF (io.eq.'O') THEN
    rhostart=ind_rho_infty
    rhoend=rhostart
  ELSE
    PRINT*,'Incorrect I/O designation in rotational'
    STOP'rotational'
  ENDIF

! Changed  jmax to j for all further commands
ALLOCATE(plmn(0:jmax))

IF (frame.eq.'SF') THEN

  IF (jtotal.eq.0) THEN
    
    ipt=0

    wigner_d=1.d0

    DO irho=rhostart,rhoend
      DO itheta=1,ntheta 
        DO ichi=1,nchi

          ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)
          ipt=ipt+1

          gcoord=thetacap(ispt)
          costheta=Cos(gcoord)
          cgcoef=cleb(j,j,0,0,0,0)
          Call plm(jmax, 0, costheta, plmn)
          !rot(ipt)=plmn(j)*wigner_d
          !shfac=Sqrt((2.d0*j+1.d0)*ftl(j))/Sqrt(2.d0*ftl(j))
          rot(ipt)=Sqrt(2.d0*j+1.d0)*cgcoef*plmn(j)*wigner_d
          ! add in phase factor sqrt((2*l+1)*(l-m)!)/(4*pi*(l+m))) 
          !rot(ipt)=Cmplx(realpart,0.d0,dp)

        ENDDO
      ENDDO 
    ENDDO 

  ELSE
    PRINT*,'Finish rotational routine'
    STOP'rotational'

    rot=0.d0
  
 
    ! SUM OVER OMEGA
     omega=MIN(jtotal,j)
     coeff=Sqrt((2.d0*l+1.d0)/(2.d0*jtotal+1.d0))
     DO iomega=-omega,omega
        cgcoef=cleb(j,l,jtotal,omega,0,omega)
        ipt=0
        DO irho=rhostart,rhoend
          DO itheta=1,ntheta
            DO ichi=1,nchi

               ispt=aphindex(ntheta,nchi,irho,itheta,ichi,.false.)
               ipt=ipt+1
               gcoord=thetacap(ispt)
               costheta=Cos(gcoord)
               beta=betaqf(ispt)

               Call plm(j, l, costheta, plmn)
               Call wigner(jtotal, l, iomega, beta, wigner_d)
               rot(ipt)=cgcoef*plmn(j)*wigner_d
         !      rot(ipt)=rot(ipt)+Cmplx(realpart,0.d0,dp)
         !      CALL aphdef() and djmk() (modify both)

            ENDDO ! chi
          ENDDO ! theta
        ENDDO ! rho
     ENDDO ! iomega

     rot=coeff*rot

  ENDIF

ELSE

  PRINT*,'ERROR!'
  PRINT*,'INVALID FRAME OPTION'
  STOP'rotational'

ENDIF

!===============================================================================
! Check Normalization by interpolating to a uniform grid.
! Would have to sort it to ascending order too. Time consuming!

DEALLOCATE(plmn)

Return
End Subroutine rotational