!*******************************************************************************
! REID93: Updated Reid potential, regularized with a dipole form factor
! of 8 pion masses, in momentum space on LSJ basis.
! Reference: V.G.J. Stoks et al., Phys. Rev. C 49, 2950 (1994).
!-------------------------------------------------------------------------------
! This is f90 version 1.1, 20 November 1999
! WWW: http://nn-online.org
! Email: info@nn-online.org
!-------------------------------------------------------------------------------
! IN:  real(wp)     :: qi        = center of mass momentum initial state in MeV
!      real(wp)     :: qf        = center of mass momentum final   state in MeV
!      character(2) :: type      = reaction; 'PP', 'NN', 'NP', or 'PN'
!      character(3) :: name      = name of the partial wave (see below),
!                                  maximum total angular momentum j = 9
! OUT: real(wp)     :: vpot(2,2) = potential matrix in MeV^{-2}
!-------------------------------------------------------------------------------
! Integer variable 'wp' defines the precision of 'real' variables. It is by
! default put to double precision. This code is made with double-precision
! of 8 bytes in mind.  We do not know about the behavior of the code for
! other precisions, but do not expect anything anomalous.
!-------------------------------------------------------------------------------
! Defining the K-matrix as   2i*mu*q*K = (1-S)(1+S)^-1
! (i.e. for singlet channel   tan(delta) = -2*mu*q*K )
! the partial-wave Lippmann-Schwinger equation reads
!     K(q',q) = V(q',q) + 2/pi int dk k^2 V(q',k) G(q,k) K(k,q)
! with
!     G(q,k) = P / (E(q) - k^2/2/mu)
!     V(q',k) = 1 / (4*pi) * vpot(qi=k,qf=q')
! S.Y.M. convention is used, so additional minus sign in off-diagonal tensor
! potential (shows up in mixing parameter).
!-------------------------------------------------------------------------------
! Potential decomposition in momentum space plane-wave basis:
! V(qf,qi) =    Vc
!             + Vs   (sigma_1.sigma_2)        (only in one-pion-exchange)
!             + Vt   [(sigma_1.k)(sigma_1.k)-k2/3(sigma_1.sigma_1)]
!             + Vls  (i/2)(sigma_1+sigma_1).n
!
! k = qf - qi ,   q = (qf+qi)/2 ,   n = qi x qf = q x k
!-------------------------------------------------------------------------------
! The OPE-potential distinguishes between neutral and charged pion masses,
! and has coupling constants f0pi**2 = fcpi**2 = 0.075. The delta-function
! (smeared out due to the form factor) only contributes to the S waves.
!-------------------------------------------------------------------------------
! The variable 'name' contains the name of the partial wave in spectral
! notation:  singlets:                 1S0  1P1  1D2  1F3  1G4 ...
!            triplets uncoupled:       3P0  3P1  3D2  3F3  3G4 ...
!            triplets coupled:              3C1  3C2  3C3  3C4 ...
! where 3C1 denotes   3S1 - EPS1 - 3D1 channel
!       3C2 denotes   3P2 - EPS2 - 3F2 channel
!       ...
! Only capitals are recognized!
!*******************************************************************************

      subroutine reid93(qi,qf,type,name,vpot)

      implicit none

!     integer, parameter        :: wp = selected_real_kind(15)
      integer, parameter        :: wp = kind(1.d0)       ! double precision

      real(wp), intent(in)      :: qi,qf
      character(2), intent(in)  :: type
      character(3), intent(in)  :: name
      real(wp), intent(out)     :: vpot(2,2)

      integer, save             :: nchan,l,spin,j,iso
      integer, save             :: jmm,jm,jp,jpp,jmax
      real(wp), save            :: tjmm,tjm,tj,tjp,tjpp,tjj
      character(3), save        :: name0 = '***'
      integer                   :: ll,ie,im
      real(wp)                  :: eln(0:12)
      real(wp)                  :: vpis(-1:1),vpit(-2:2)
      real(wp)                  :: vc(-1:1),vl(-2:2),vt(-2:2)
      real(wp)                  :: fac,vten,vpisl0,l2,m2,x,y
      real(wp)                  :: qiqf,qi2,qf2,qi2f2,c2psi,s2psi,cpsi2,spsi2

      real(wp), parameter       :: f0pi = 0.075_wp
      real(wp), parameter       :: fcpi = 0.075_wp
      real(wp), parameter       :: pi = 3.1415926535897932384626433832795028_wp
      real(wp), parameter       :: m_pi0 = 134.9739_wp
      real(wp), parameter       :: m_pic = 139.5675_wp
      real(wp), parameter       :: m_pis = m_pic
      real(wp), parameter       :: m_pis2 = m_pis**2
      real(wp), parameter       :: m_pi = (m_pi0+2.0*m_pic)/3.0
      real(wp), parameter       :: m_pi2 = m_pi**2
      real(wp), parameter       :: zero = 0.0_wp
      real(wp), parameter       :: half = 0.5_wp
      real(wp), parameter       :: one = 1.0_wp
      real(wp), parameter       :: two = 2.0_wp
      real(wp), parameter       :: three = 3.0_wp
      real(wp), parameter       :: four = 4.0_wp
      real(wp), parameter       :: five = 5.0_wp

      real(wp)                  :: ul(6,-2:12) = zero
      real(wp)                  :: ulc(-2:12) = zero

      real(wp), save            :: a(5,5), b(5,5)
      logical, save             :: first = .true.

! pp parameters (isospin = 1)
      real(wp), parameter       :: aa(25) =                                 &
      (/  0.1756084e0_wp,-0.1414234e2_wp, 0.1518489e3_wp,-0.6868230e3_wp,   &
          0.1104157e4_wp,-0.4224976e2_wp, 0.2072246e3_wp,-0.3354364e3_wp,   &
         -0.1989250e1_wp,-0.6178469e2_wp, 0.2912845e2_wp, 0.1511690e3_wp,   &
          0.8151964e1_wp, 0.5832103e2_wp,-0.2074743e2_wp,-0.5840566e0_wp,   &
         -0.1029310e2_wp, 0.2263391e2_wp, 0.2316915e2_wp,-0.1959172e1_wp,   &
         -0.2608488e1_wp, 0.1090858e2_wp,-0.4374212e0_wp,-0.2148862e2_wp,   &
         -0.6584788e0_wp                                                    &
      /)

! np parameters (isospin = 0)
      real(wp), parameter       :: bb(25) =                                 &
      (/ -0.2234989e2_wp, 0.2551761e3_wp,-0.1063549e4_wp, 0.1609196e4_wp,   &
         -0.3505968e1_wp,-0.4248612e1_wp,-0.5352001e1_wp, 0.1827642e3_wp,   &
         -0.3927086e3_wp, 0.5812273e2_wp,-0.2904577e1_wp, 0.3802497e2_wp,   &
          0.3395927e0_wp, 0.8318097e0_wp, 0.1923895e1_wp, 0.0913746e0_wp,   &
         -0.1274773e2_wp, 0.1458600e3_wp,-0.6432461e3_wp, 0.1022217e4_wp,   &
         -0.0461640e0_wp, 0.7950192e1_wp,-0.1925573e1_wp, 0.5066234e2_wp,   &
          0.8359896e1_wp                                                    &
      /)

      if (first) then
          a = reshape(aa, (/5,5/), order=(/2,1/) )
          b = reshape(bb, (/5,5/), order=(/2,1/) )
          first = .false.
      endif

      if (name /= name0) then
          name0 = name
          nchan = 1
          if (name(2:2) == 'C') nchan = 2
          if (name(1:1) == '1') spin = 0
          if (name(1:1) == '3') spin = 1
          read(name,'(2x,i1)') j
          if (j > 9) stop 'REID93: j > 9 not allowed'
          l = j
          if (name == '3P0') l = 1
          if (nchan == 2) l = j-1
          iso = mod(spin+l+1,2)
          jmm  = j-2
          jm   = j-1
          jp   = j+1
          jpp  = j+2
          tjmm = two*j-three
          tjm  = two*j-one
          tj   = two*j+one
          tjp  = two*j+three
          tjpp = two*j+five
          tjj  = sqrt(j*(j+one))
          jmax = j+2
      endif

      qi2 = qi**2
      qf2 = qf**2
      qiqf = qi*qf
      qi2f2 = qi2+qf2
      s2psi = two*qiqf/qi2f2
      spsi2 = qf2/qi2f2
      cpsi2 = qi2/qi2f2

!     Neutral OPE
      m2 = m_pi0**2
      l2 = 64.0*m2
      x = half*(qi2f2+m2)/qiqf
      y = half*(qi2f2+l2)/qiqf
      call xdip(x,y,jmax,eln)
      ul(1,0:jmax) = eln(0:jmax)
      fac = two*pi*f0pi/(m_pis2*qiqf)
      do ll=-1,1
          vpis(ll) = fac*ul(1,j+ll)*m_pi0**2/three
      enddo
      vpisl0 = four*pi*f0pi/(m_pis2*three) * (y-x)**2/(one-y**2)
      do ll=-2,2
          vpit(ll) = -fac*ul(1,j+ll)
      enddo

!     Charged OPE
      select case(type)
      case('NP','PN','np','pn')
          m2 = m_pic**2
          l2 = 64.0*m2
          x = half*(qi2f2+m2)/qiqf
          y = half*(qi2f2+l2)/qiqf
          call xdip(x,y,jmax,eln)
          ulc(0:jmax) = eln(0:jmax)
          fac = (four*iso-two)*two*pi*fcpi/(m_pis2*qiqf)
          do ll=-1,1
              vpis(ll) = fac*ulc(j+ll)*m_pic**2/three - vpis(ll)
          enddo
          vpisl0 = four*pi*fcpi/(m_pis2*three) *                             &
              (four*iso-two)*(y-x)**2/(one-y**2) - vpisl0
          do ll=-2,2
              vpit(ll) = -fac*ulc(j+ll) - vpit(ll)
          enddo
      end select

!     Other Yukawa's with multiples of average pion mass
      l2 = 64.0*m_pi2
      y = half*(qi2f2+l2)/qiqf
      do im=2,6
          x = half*(qi2f2+im**2*m_pi2)/qiqf
          call xdip(x,y,jmax,eln)
          ul(im,0:jmax) = eln(0:jmax)
      enddo

      fac = two*pi/qiqf

!     Potential for each partial wave separately
      select case(name)
      case('1S0')
          select case(type)
          case('PP','NN','pp','nn')
              vpot(1,1) = fac*( a(1,1)*ul(2,j)+a(1,2)*ul(3,j)+a(1,3)*ul(4,j) &
                               +a(1,4)*ul(5,j)+a(1,5)*ul(6,j) )
          case('NP','PN','np','pn')
              vpot(1,1) = fac*( b(1,1)*ul(3,j)+b(1,2)*ul(4,j)+b(1,3)*ul(5,j) &
                               +b(1,4)*ul(6,j) )
          end select
      case('1D2')
          vpot(1,1) = fac*(a(2,1)*ul(4,j)+a(2,2)*ul(5,j)+a(2,3)*ul(6,j) )
      case('1G4')
          vpot(1,1) = fac*a(2,4)*ul(3,j)
      case('3P0')
          vten = ul(3,jp)-half*s2psi*(tjpp*ul(3,j)+tj*ul(3,jpp))/tjp
          vpot(1,1) = fac*( a(3,1)*ul(3,jp)+a(3,2)*ul(5,jp)                  &
                           +a(2,5)*(-vten/9.0/m_pi2)*qi2f2/three )
      case('3F3')
          vpot(1,1) = fac*a(4,5)*ul(3,j)
      case('3P1')
          vten = ul(3,j)-half*s2psi*(tjp*ul(3,jm)+tjm*ul(3,jp))/tj
          vpot(1,1) = fac*( a(3,3)*ul(3,j)+a(3,4)*ul(5,j)                    &
                           +a(3,5)*(-vten/9.0/m_pi2)*qi2f2/three )
      case('3P3')
          vpot(1,1) = fac*a(4,5)*ul(3,j)
      case('1P1','1H5','1K7','1M9')
          vpot(1,1) = fac*( b(2,1)*ul(3,j)+b(2,2)*ul(4,j)+b(2,3)*ul(5,j)     &
                           +b(2,4)*ul(6,j) )
      case('1F3')
          vpot(1,1) = fac*(b(1,5)*ul(3,j)+b(2,5)*ul(5,j) )
      case('3D2')
          vten = ul(3,j)-half*s2psi*(tjp*ul(3,jm)+tjm*ul(3,jp))/tj
          vpot(1,1) = fac*( b(3,1)*ul(3,j)+b(3,2)*ul(5,j)                    &
                           +b(3,3)*(-vten/9.0/m_pi2)*qi2f2/three )
      case('3G4')
          vpot(1,1) = fac*b(3,4)*ul(3,j)
      case('1I6','1L8')     ! 1S0 PP parameters used
          vpot(1,1) = fac*( a(1,1)*ul(2,j)+a(1,2)*ul(3,j)+a(1,3)*ul(4,j)     &
                           +a(1,4)*ul(5,j)+a(1,5)*ul(6,j) )
      case default       ! spin=1, uncoupled j>5 and coupled
          select case(iso)
          case(1)
              do ll=-1,1
                  vc(ll) = fac*( a(4,1)*ul(3,j+ll)+a(4,2)*ul(4,j+ll)     &
                                +a(4,3)*ul(5,j+ll)+a(4,4)*ul(6,j+ll) )
              enddo
              do ll=-2,2
                  vt(ll) = -fac*( a(5,1)*ul(4,j+ll)/16.0/m_pi2           &
                                 +a(5,2)*ul(6,j+ll)/36.0/m_pi2 )
              enddo
          case(0)
              do ll=-1,1
                  vc(ll) = fac*( b(4,1)*ul(2,j+ll)+b(4,2)*ul(3,j+ll)     &
                                +b(4,3)*ul(4,j+ll)+b(4,4)*ul(5,j+ll)     &
                                +b(4,5)*ul(6,j+ll) )
              enddo
              do ll=-2,2
                  vt(ll) = -fac*( b(3,5)*ul(4,j+ll)/16.0/m_pi2           &
                                 +b(5,5)*ul(6,j+ll)/36.0/m_pi2 )
              enddo
          end select

          select case(nchan)
          case(1)
              vl(:) = zero
              vpot(1,1) = vc(0) - qiqf*(vl(-1)-vl(1))/tj +               &
                          two/three*qi2f2*                               &
                          (vt(0)-half*s2psi*(tjp*vt(-1)+tjm*vt(1))/tj)
              vpot(1,2) = zero
              vpot(2,1) = zero
              vpot(2,2) = zero
          case(2)
              select case(name)
              case('3C1')
                  do ll=-2,2
                      vl(ll) = fac*( b(5,1)*ul(3,j+ll)/ 9.0/m_pi2        &
                                    +b(5,2)*ul(5,j+ll)/25.0/m_pi2 )
                  enddo
              case('3C2')
                  do ll=-2,2
                      vl(ll) = fac*( a(5,3)*ul(3,j+ll)/ 9.0/m_pi2        &
                                    +a(5,4)*ul(5,j+ll)/25.0/m_pi2 )
                  enddo
              case('3C3')
                  do ll=-2,2
                      vl(ll) = fac*( b(5,3)*ul(3,j+ll)/ 9.0/m_pi2        &
                                    +b(5,4)*ul(5,j+ll)/25.0/m_pi2 )
                  enddo
              case('3C4')
                  do ll=-2,2
                      vl(ll) = fac*a(5,5)*ul(3,j+ll)/ 9.0/m_pi2
                  enddo
              case default
                  vl(:) = zero
              end select

              vpot(1,1) = vc(-1) + qiqf*jm/tjm*(vl(-2)-vl(0)) +          &
                          two/three*qi2f2*jm/tj*                         &
                          (-vt(-1)+half*s2psi*(tjmm*vt(0)+tj*vt(-2))/tjm)
              vpot(1,2) = -two*qi2f2*tjj/tj*                             &
                          (-s2psi*vt(0)+cpsi2*vt(-1)+spsi2*vt(1))
              vpot(2,1) = -two*qi2f2*tjj/tj*                             &
                          (-s2psi*vt(0)+spsi2*vt(-1)+cpsi2*vt(1))
              vpot(2,2) = vc(1) - qiqf*jpp/tjp*(vl(0)-vl(2)) +           &
                          two/three*qi2f2*jpp/tj*                        &
                          (-vt(1)+half*s2psi*(tjpp*vt(0)+tj*vt(2))/tjp)
          end select
      end select

!     Add OPE potential
      select case(nchan)
      case(1)
          if (spin == 0) then
              vpot(1,1) = vpot(1,1) - three*vpis(0)
              if (l == 0) vpot(1,1) = vpot(1,1) - three*vpisl0
          elseif (l == j) then
              vpot(1,1) = vpot(1,1) + vpis(0) + two/three*qi2f2*             &
                          (vpit(0)-half*s2psi*(tjp*vpit(-1)+tjm*vpit(1))/tj)
          elseif (name == '3P0') then
              vpot(1,1) = vpot(1,1) + vpis(1) + two/three*qi2f2*jpp/tj*      &
                          (-vpit(1)+half*s2psi*(tjpp*vpit(0)+tj*vpit(2))/tjp)
          endif
      case(2)
          vpot(1,1) = vpot(1,1) + vpis(-1) + two/three*qi2f2*jm/tj*          &
                      (-vpit(-1)+half*s2psi*(tjmm*vpit(0)+tj*vpit(-2))/tjm)
          if (l == 0) vpot(1,1) = vpot(1,1) + vpisl0
          vpot(1,2) = vpot(1,2) - two*qi2f2*tjj/tj*                          &
                      (-s2psi*vpit(0)+cpsi2*vpit(-1)+spsi2*vpit(1))
          vpot(2,1) = vpot(2,1) - two*qi2f2*tjj/tj*                          &
                      (-s2psi*vpit(0)+spsi2*vpit(-1)+cpsi2*vpit(1))
          vpot(2,2) = vpot(2,2) + vpis(1) + two/three*qi2f2*jpp/tj*          &
                      (-vpit(1)+half*s2psi*(tjpp*vpit(0)+tj*vpit(2))/tjp)
      end select

      return

      contains

!*******************************************************************************

      subroutine xdip(x,y,jmax,eln)

      implicit none

!     integer, parameter    :: wp = kind(1.d0)
      real(wp), intent(in)  :: x,y
      logical, save         :: first = .true.
      integer, intent(in)   :: jmax
      real(wp)              :: p(12)
      real(wp), intent(out) :: eln(0:12)
      integer, parameter    :: nogp = 64
      real(wp), save        :: zz(nogp),wz(nogp)
      real(wp)              :: qx0,qy0,form
      integer               :: iz,l

      qx0 = 0.5*log((x+1.0)/(x-1.0))
      qy0 = 0.5*log((y+1.0)/(y-1.0))
      eln(0) = qx0-qy0+(y-x)/(1.0-y**2)
      if (jmax == 0) return
      if (first) then
          call gauss(-1.0_wp,1.0_wp,zz,wz,nogp)
          first = .false.
      endif

      eln(1:) = 0.0_wp
      do iz=1,nogp
          form = (y-x)/(y-zz(iz))
          form = form**2/(x-zz(iz))
!         Calculate Legendre polynomials first kind and integrate
          p(1) = zz(iz)
          p(2) = 1.5*zz(iz)*zz(iz)-0.5
          eln(1) = eln(1)+0.5*wz(iz)*form*p(1)
          eln(2) = eln(2)+0.5*wz(iz)*form*p(2)
          do l=3,jmax
              p(l) = ( (2*l-1)*zz(iz)*p(l-1) - (l-1)*p(l-2) ) / l
              eln(l) = eln(l)+0.5*wz(iz)*form*p(l)
          enddo
      enddo

      end subroutine xdip

!*******************************************************************************
! Calculate abscissas and weigths for Gauss-Legendre integration
! (n-points Gauss-Legendre quadrature).
!
! IN:    x1,x2     begin and end of integration range
!        n         number of integration points
! OUT:   x(1:n)    abscissas
!        w(1:n)    weights
!
! Literature: * Numerical Recipes in Fortran, paragraph 4.5 (1st edition)
!             * Abramowitz and Stegun, 25.4.29 and 25.4.30
!             * Mathews and Walker, paragraph 13.2
!             * J.J. de Swart, Lecture notes MTF (THEF-NIJM 84.04) par. F4.8
!
! Only for integrals with finite integration range. For improper integrals:
! perform a mapping first to an integral with a finite integration range.

      subroutine gauss(x1,x2,x,w,n)

      implicit none
!     integer, parameter      :: wp = kind(1.d0)
      real(wp), intent(in)    :: x1,x2
      integer, intent(in)     :: n
      real(wp), intent(out)   :: x(n),w(n)
      integer                 :: i,j,m
      real(wp)                :: xm,xl,z,z1,p1,p2,p3,pp
      real(wp), parameter     :: eps = 3e-14_wp
      real(wp), parameter     :: pi = 3.1415926535897932384626433832795028_wp
      real(wp), parameter     :: r0 = 0.0_wp
      real(wp), parameter     :: r025 = 0.25_wp
      real(wp), parameter     :: r05 = 0.5_wp
      real(wp), parameter     :: r1 = 1.0_wp
      real(wp), parameter     :: r2 = 2.0_wp

      m  = (n+1)/2
      xm = r05*(x2+x1)
      xl = r05*(x2-x1)
      do i=1,m
          z = cos(pi*(i-r025)/(n+r05))
          do
              p1 = r1
              p2 = r0
              do j=1,n
                  p3 = p2
                  p2 = p1
                  p1 = ((2*j-1)*z*p2-(j-1)*p3)/j
              enddo
              pp = n*(z*p1-p2)/(z*z-r1)
              z1 = z
              z  = z1-p1/pp
              if (abs(z-z1) <= eps) exit
          enddo
          x(i) = xm-xl*z
          x(n+1-i) = xm+xl*z
          w(i) = r2*xl/((r1-z*z)*pp*pp)
          w(n+1-i) = w(i)
      enddo

      end subroutine gauss

!*******************************************************************************

      end subroutine reid93
