      SUBROUTINE Set_version_hfodd_wigner
      !
      USE hfodd_sizes
      !
      CHARACTER NAMMOD*16
      INTEGER   MODUVE,MODSET
      !
      COMMON /VERMOD/ NAMMOD(NDMODU),MODUVE(NDMODU),MODSET(NDMODU)
      !
      NAMMOD(16) = 'hfodd_wigner   '
      MODUVE(16) = 6
      MODSET(16) = 16
      !
      END SUBROUTINE Set_version_hfodd_wigner
      !
module hfodd_wigner

      SAVE D_WIGN,IWIIMK,BETANG,ISWIGN

      ALLOCATABLE D_WIGN(:,:)
      ALLOCATABLE IWIIMK(:,:,:)
      ALLOCATABLE BETANG(:)
      ALLOCATABLE ISWIGN(:,:)

end module hfodd_wigner
!     CONTAINS
!
!=======================================================================
!         BELOW ARE IMPLEMENTATIONS COPIED FROM THE SUPPLEMENTAL MATERIAL OF
!         https://journals.aps.org/pre/abstract/10.1103/PhysRevE.92.043307
!         A FEW MODIFCATIONS ARE MARKED BY "!JD 190316"
!=======================================================================
!************************************************************************************************************************************************************
!  FUNCTIONS: Wigner's d-matrix
!  Methods: (a) Complex Fourier-series expansion of the d-matrix;
!           (b) the Fourier coefficients are obtained by numerical diagonalizing the angular-momentum operator Jy, using the ZHBEV subroutine of LAPACK.
!
!  update:  (y/m/d) 2015/09/30 by G.R.J
!*************************************************************************************************************************************************************
!
!   [1] The defination of the term X_{m}=<j,m-1|J_{-}|j,m>, where J_{-} is the ladder operator,
!       obeying J_{-}|j,m>=X_{m} |j,m-1>, and {|j,m>} eigenvectors of the angular-momentum operator Jz.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    function X(j,n)
!JD 190316 vvvv
!       real*8 X,n,j
        real   X,n,j
!JD 190316 ^^^^
          X=dsqrt((j+n)*(j-n+1.d0))
        end

!****************************************************************************************************
!   [2]  The Hermitian matrix of J_{y}
!******************************************************************************************************
        subroutine coeffi(A,Ndim)
!        external X
!JD 190316 vvvv
!       implicit real*8 (a-h,o-z)
!       real*8 jmq
!       COMPLEX*16 A(Ndim,Ndim), im
        implicit real   (a-h,o-z)
        real   jmq
        COMPLEX    A(Ndim,Ndim), im
!JD 190316 ^^^^
        integer Ndim,m
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    im=(0.d0,1.d0)
    jmq=(dble(Ndim)-1.d0)/2.d0
        do i=1,Ndim
                do j=1,Ndim
            A(i,j)=0.d0

           if(i.eq.1.and.Ndim.gt.1) then
                A(i,i+1)=X(jmq,dble(i)-jmq)*im/(2.d0)
            endif
            if(i.eq.Ndim.and.Ndim.gt.1) then
                A(i,i-1)=-X(jmq,-dble(i)+jmq+2)*im/(2.d0)
            endif

               if((i.gt.1).and.(i.lt.Ndim))then
                A(i,i+1)=X(jmq,dble(i)-jmq)*im/(2.d0)
                    A(i,i-1)=-X(jmq,-dble(i)+jmq+2)*im/(2.d0)
               endif
            enddo
    enddo
        return
    end

!**************************************************************************************************************************
!   [3] Calculation of the eigenvalues and the right eigenvectors for the Hermitian matrix J_{y},
!       using ZHBEV subroutine of LAPACK.
!****************************************************************************************************************************
    subroutine eigen(A,Eigenvalue,Eigenvec,Ndim)
    character *1 JOBZ,UPLO
    integer Ndim,KD,LDAB,LDZ,INFO
!JD 190316 vvvv
!   double precision Eigenvalue(Ndim)
!   double precision,allocatable:: RWORK(:)
!   double complex,allocatable:: WORK(:),AB(:,:)
!   double complex Eigenvec(Ndim,Ndim),A(Ndim,Ndim)
    real Eigenvalue(Ndim)
    real,allocatable:: RWORK(:)
    complex,allocatable:: WORK(:),AB(:,:)
    complex Eigenvec(Ndim,Ndim),A(Ndim,Ndim)
!JD 190316 ^^^^
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   JOBZ='V'
   UPLO='U'
   KD=1
   LDAB=KD+1
   LDZ=max(1,Ndim)
   INFO=0
   allocate(AB(LDAB,Ndim),WORK(Ndim),RWORK(max(1,3*Ndim-2)))

   call coeffi(A,Ndim)
        do i=1,Ndim
                  do j=i,Ndim
                   if(max(1,j-kd).le.i) then
            AB(KD+1+i-j,j)=A(i,j)
                   end if
          enddo
        enddo

   call ZHBEV (JOBZ, UPLO, Ndim, KD, AB, LDAB, Eigenvalue, Eigenvec, LDZ, WORK, RWORK,INFO)

   deallocate(RWORK,WORK,AB)

   end


!**************************************************************************************************************************
!   [4] The subroutine to calculate the Wigner-d matrix and its first-order derivative with various mv and nv for a given {j,theta}
!****************************************************************************************************************************
!JD 190316 vvvv
!       subroutine Wigner_dmatrix(jmq,theta,wd_matrix,diffwd_matrix)
        subroutine Wigner_dmatrix(jmq,theta,wd_matrix)
!=======================================================================
      USE hfodd_sizes
!=======================================================================
!JD 190316 ^^^^
        integer Ndim,ixx,iyy,mu
!JD 190316 vvvv
!       real *8 jmq, mvar, nvar,theta, wd, wdderivative, wd_matrix(int(2*jmq+1),int(2*jmq+1)), diffwd_matrix(int(2*jmq+1),int(2*jmq+1))
!       double complex  inum
!       double precision,allocatable :: Evalue(:),Eigenvalue(:)
!       double complex,allocatable :: coeffimatrix(:,:), Eigenvector(:,:)
        real    jmq, mvar, nvar,theta, wd
        dimension wd_matrix(int(2*jmq+1),int(2*jmq+1))
        complex  inum
        real,allocatable :: Evalue(:),Eigenvalue(:)
        complex,allocatable :: coeffimatrix(:,:), Eigenvector(:,:)
!JD 190316 ^^^^
! +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
        inum=(0.d0,1.d0)
    Ndim=int(2.d0*jmq+1.d0)
    allocate(Evalue(Ndim),coeffimatrix(Ndim,Ndim),Eigenvector(Ndim,Ndim),Eigenvalue(Ndim))
    Evalue=0.0d0
    coeffimatrix=0.0d0
    Eigenvector=0.0d0
    Eigenvalue=0.0d0
!    call coeffi(coeffimatrix,Ndim)
     call eigen(coeffimatrix,Eigenvalue,Eigenvector,Ndim)

     do ixx=1,Ndim
           do iyy=1,Ndim

              wd=0.d0
!JD 190316 vvvv
!                 wdderivative=0.d0
!JD 190316 ^^^^

              do mu=1,Ndim
                     if(dble(int(jmq)).eq.jmq) then
                Evalue(mu)=dble(floor(Real(Eigenvalue(mu))+0.5d0))
                         else
                Evalue(mu)=dble(floor(Real(Eigenvalue(mu)))+0.5d0)
                         endif
                wd=wd+exp(-inum*Evalue(mu)*theta)*Eigenvector(ixx,mu)*dconjg(Eigenvector(iyy,mu))
!JD 190316 vvvv
!               wdderivative=wdderivative+(-inum*Evalue(mu))*exp(-inum*Evalue(mu)*theta)*Eigenvector(ixx,mu)*dconjg(Eigenvector(iyy,mu))
!JD 190316 ^^^^
                  enddo

                  wd_matrix(ixx,iyy)=wd
!JD 190316 vvvv
!                 diffwd_matrix(ixx,iyy)=wdderivative
!JD 190316 ^^^^
                end do
         end do

     deallocate(Evalue, coeffimatrix, Eigenvector,Eigenvalue)
        return
    end

