module sec_ham
! generates sector Hamiltonians
! Time-stamp: "2003-10-03 17:36:23 cjn"
  use precisn, only: wp
  use io_units, only: fo
  use error_out, only: error_check

  implicit none

  integer, save               :: nx       ! quadradure order
  real(wp), allocatable, save :: xi(:)    ! quadrature nodes
  real(wp), allocatable, save :: wi(:)    ! quadrature weights
  real(wp), allocatable, save :: pl(:,:)  ! normalized Legendre polynomials
  real(wp), allocatable, save :: vbl(:)   ! sector left bdry amplitudes
  real(wp), allocatable, save :: vbr(:)   ! sector right bdry amplitudesx

  real(wp), save              :: r12, rsq, r1b, r2b

  integer, save               :: ncq

  private
  public def_rb_vals, h_el, gauleg, h_reset, vbl, vbr, legndr, xi
  public def_ncq

contains

  subroutine def_ncq (ncp)
    integer, intent(in)   :: ncp   ! spin partition offset
    ncq = ncp - 1
  end subroutine def_ncq

  subroutine def_rb_vals (ral, rar)
! define sector geometry parameters
    use hfile_data, only: bbloch
    real(wp), intent(in)    :: ral   ! left endpoint of sector
    real(wp), intent(in)    :: rar   ! right endpoint of sector

    r12 = 1.0_wp / (rar - ral)
    if (bbloch == 0.0_wp) then
       r1b = 0.0_wp
       r2b = 0.0_wp
    else
       r1b = r12 * bbloch / ral
       r2b = r12 * bbloch / rar
    end if
    rsq = 2.0_wp * r12**2
  end subroutine def_rb_vals

  function h_el (iglob, jglob)
! Hamiltonian matrix element for global indices iglob, jglob
! should be called for iglob <= jglob
    use sector_potl, only: nc, v
    use hfile_data, only: ethr
    use def_sectors, only: nl
    real(wp)              :: h_el
    integer, intent(in)   :: iglob, jglob  ! coordinates of element
    integer               :: i, j, in, jn, ic, jc
    real(wp)              :: rmn, ss
    integer               :: ij

    j = jglob
    jc = 1 + (j-1)/nl     ! channel index
    jn = j - (jc-1)*nl    ! basis index (inner loop)
    i = iglob
    ic = 1 + (i-1)/nl     ! channel index
    in = i - (ic-1)*nl    ! basis index

    ss = 0.0_wp
    if (ic == jc) then
       rmn = SQRT(REAL((2*jn-1)*(2*in-1),wp))
       ss = - rmn * (r2b - r1b * (-1)**MOD(in+jn,2))
       if (in == jn) ss = ss + ethr(ncq+jc)
       if (MOD(in+jn,2) == 0) ss = ss + rsq * rmn * REAL(in*(in-1),wp)
    end if
    h_el = ss + SUM(pl(:,in) * wi(:) * pl(:,jn) * v(:,ic,jc))
  end function h_el

  subroutine legndr
! Legndr: normalized legendre polynomials at quadrature points
    use def_sectors, only: nl, nx, w_sect
    real(wp)               :: sgn
    integer                :: k, status

    allocate (pl(nx,nl), vbl(nl), vbr(nl), stat=status)
    call error_check (status, 'legndr: allocation')

    pl(:,1) = SQRT(0.5_wp)
    vbr(1) = SQRT(1.0_wp/w_sect)
    vbl(1) = vbr(1)
    if (nl == 1) return
    pl(:,2) = SQRT(1.5_wp) * xi
    vbr(2) = SQRT(3.0_wp/w_sect)
    vbl(2) = - vbr(2)
    sgn = -1.0_wp
    do k = 3, nl
       pl(:,k) = (SQRT(REAL(2*k-3,wp)) * xi * pl(:,k-1) - REAL(k-2,wp) * &
            pl(:,k-2) / SQRT(REAL(2*k-5,wp))) * SQRT(REAL(2*k-1,wp)) / &
            REAL(k-1,wp)
       vbr(k) = SQRT(REAL(2*k-1,wp)/w_sect)
       sgn = - sgn
       vbl(k) = sgn * vbr(k)
    end do
  end subroutine legndr

  subroutine gauleg
!  Gauss-Legendre quadrature abscissas and weights
    use def_sectors, only: nx
    integer                 :: i, ierr
    real(wp)                :: b(nx-1), d(nx), e(nx)     ! automatic
    real(wp)                :: z(nx,nx)                 ! automatic
    real(wp)                :: work(MAX(1,2*nx-2))     ! automatic
    integer                 :: status

    allocate (xi(nx), wi(nx), stat=status)
    call error_check (status, 'gauleg: allocation')

    do i = 1, nx-1
       b(i) = REAL(i*i, wp) / REAL((2*i+1) * (2*i-1), wp)
    end do
    d = 0.0_wp                  ! diagonal elements
    e(:nx-1) = - SQRT(b(1:nx-1))   ! sub-diagonal elements
    call dsteqr ('i', nx, d, e(:nx-1), z, nx, work, ierr) ! diagonalization
    if (ierr /= 0) then
       write (fo,'(a,i6)') 'GAULEG: error code from dsteqr:', ierr
    end if
    e = 2.0_wp * z(1,:)**2
    xi = d
    wi = e
  end subroutine gauleg

  subroutine h_reset
! deallocate memory
    integer    :: status
    deallocate (xi, wi, vbr, vbl, pl, stat=status)
    call error_check (status, 'h_reset: deallocation')
  end subroutine h_reset
end module sec_ham
