!*******************************************************************************
! Reid93 soft core phenomenological potential, updated version,
! including one-pion-exchange with neutral-pion and charged-pion
! masses; coupling f^2=0.075. Tensor potential is regularized
! to equal zero at r=0
! Reference: V.G.J. Stoks et al., Phys. Rev. C 49, 2950 (1994).
!-------------------------------------------------------------------------------
! Fortran 90 version 2.1: August 2001
!
! E-mail: info@nn-online.org
! WWW: http://nn-online.org
!-------------------------------------------------------------------------------
! 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.
!-------------------------------------------------------------------------------
! IN:  real(wp)          :: r         = distance in fermi
!      character(2)      :: type      = reaction, 'PP', 'NP', 'PN', or 'NN'
!      character(3)      :: pname     = name of the partial wave (see below),
!                                       maximum total angular momentum j = 9
! OUT: real(wp)          :: vpot(2,2) = potential matrix in MeV on LSJ-basis
!-------------------------------------------------------------------------------
! The variable 'pname' 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 - epsilon_1 - 3D1 channel
!       3C2 denotes   3P2 - epsilon_2 - 3F2 channel
!       ...
!*******************************************************************************

      subroutine rreid93(r,pname,type,vpot)

      implicit none

      integer, parameter          :: wp = kind(1.0d0)

      real(wp), intent(in)        :: r
      character(2), intent(in)    :: type
      character(3), intent(in)    :: pname
      real(wp), intent(out)       :: vpot(2,2)
      real(wp)                    :: ri,rj,vc,vt,vls,vtpi,f,x,mpi,vspi,vspis
      integer                     :: l,spin,j,nchan,iso
      character(3)                :: name
      character(*), parameter     :: large = 'CSPDFGHIKLMN'
      character(*), parameter     :: small = 'cspdfghiklmn'
      real(wp), parameter         :: r1 = 1.0_wp
      real(wp), parameter         :: r2 = 2.0_wp
      real(wp), parameter         :: r3 = 3.0_wp
      real(wp), parameter         :: r4 = 4.0_wp
      real(wp), parameter         :: r5 = 5.0_wp
      real(wp), parameter         :: r6 = 6.0_wp
      real(wp), parameter         :: r8 = 8.0_wp
      real(wp), parameter         :: f0pi = 0.075_wp
      real(wp), parameter         :: fcpi = 0.075_wp
      real(wp), parameter         :: hbc = 197.327053_wp
      real(wp), parameter         :: mpi0 = 134.9739_wp
      real(wp), parameter         :: mpic = 139.5675_wp
      real(wp), parameter         :: mpis = mpic
      real(wp), parameter         :: a(5,5) = reshape( (/                    &
            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 /) , (/5,5/) , order=(/2,1/) )
      real(wp), parameter         :: b(5,5) = reshape( (/                    &
           -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 /) , (/5,5/) , order=(/2,1/) )

!     Determine the quantumnumbers
      name = pname
      do l=1,12
          if (name(2:2) == small(l:l)) name(2:2) = large(l:l)
      enddo
      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
      l = j
      if (name == '3P0') l = 1
      if (nchan == 2) l = j - 1
      iso = mod(spin+l+1,2)

      ri = real(iso,kind=wp)
      rj = real(j,kind=wp)

!     OPE potential
      x     = mpi0/hbc * r
      f     = f0pi*(mpi0/mpis)**2*mpi0/r3
      vspis = f*yp(r1,r8,x)
      vspi  = f*y(r1,r8,x)
      vtpi  = f*z(r1,r8,x)
      select case(type)
      case('NP','PN','np','pn')
          x     = mpic/hbc * r
          f     = (r4*ri-r2)*fcpi*(mpic/mpis)**2*mpic/r3
          vspis = f*yp(r1,r8,x) - vspis
          vspi  = f*y(r1,r8,x)  - vspi
          vtpi  = f*z(r1,r8,x)  - vtpi
      endselect

!     Start clean
      vpot(:,:) = 0.0_wp

!     Non-OPE contribution to the potential
      mpi = (mpi0 + r2*mpic)/r3
      x = mpi/hbc * r

      select case(name)
      case('1S0')
          select case(type)
          case('PP','NN','pp','nn')
              vpot(1,1) =   a(1,1)*y(r2,r8,x) &
                          + a(1,2)*y(r3,r8,x) &
                          + a(1,3)*y(r4,r8,x) &
                          + a(1,4)*y(r5,r8,x) &
                          + a(1,5)*y(r6,r8,x)
          case('NP','PN','pn','np')
              vpot(1,1) =   b(1,1)*y(r3,r8,x) &
                          + b(1,2)*y(r4,r8,x) &
                          + b(1,3)*y(r5,r8,x) &
                          + b(1,4)*y(r6,r8,x)
          endselect
      case('1D2')
          vpot(1,1) =   a(2,1)*y(r4,r8,x) &
                      + a(2,2)*y(r5,r8,x) &
                      + a(2,3)*y(r6,r8,x)
      case('1G4')
          vpot(1,1) =   a(2,4)*y(r3,r8,x)
      case('3P0')
          vpot(1,1) =   a(3,1)*y(r3,r8,x) &
                      + a(3,2)*y(r5,r8,x) &
                      + a(2,5)*z(r3,r8,x)/r3
      case('3P1')
          vpot(1,1) =   a(3,3)*y(r3,r8,x) &
                      + a(3,4)*y(r5,r8,x) &
                      + a(3,5)*z(r3,r8,x)/r3
      case('3F3')
          vpot(1,1) =   a(4,5)*y(r3,r8,x)
      case('3C2','3C4')
          vc =  a(4,1)*y(r3,r8,x) &
              + a(4,2)*y(r4,r8,x) &
              + a(4,3)*y(r5,r8,x) &
              + a(4,4)*y(r6,r8,x)
          vt = (   a(5,1)*z(r4,r8,x) &
                 + a(5,2)*z(r6,r8,x) )/r3
          select case(name)
          case('3C2')
              vls = a(5,3)*w(r3,r8,x) + a(5,4)*w(r5,r8,x)
          case('3C4')
              vls = a(5,5)*w(r3,r8,x)
          endselect
          vpot(1,1) = vc + (rj-r1)*vls - r2*(rj-r1)/(r2*rj+r1)*vt
          vpot(2,2) = vc - (rj+r2)*vls - r2*(rj+r2)/(r2*rj+r1)*vt
          vpot(1,2) = r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vt
          vpot(2,1) = vpot(1,2)
      case('1P1')
          vpot(1,1) =   b(2,1)*y(r3,r8,x) &
                      + b(2,2)*y(r4,r8,x) &
                      + b(2,3)*y(r5,r8,x) &
                      + b(2,4)*y(r6,r8,x)
      case('1F3')
          vpot(1,1) =   b(1,5)*y(r3,r8,x) &
                      + b(2,5)*y(r5,r8,x)
      case('3D2')
          vpot(1,1) =   b(3,1)*y(r3,r8,x) &
                      + b(3,2)*y(r5,r8,x) &
                      + b(3,3)*z(r3,r8,x)/r3
      case('3G4')
          vpot(1,1) =   b(3,4)*y(r3,r8,x)
      case('3C1','3C3')
          vc =  b(4,1)*y(r2,r8,x) &
              + b(4,2)*y(r3,r8,x) &
              + b(4,3)*y(r4,r8,x) &
              + b(4,4)*y(r5,r8,x) &
              + b(4,5)*y(r6,r8,x)
          vt = ( b(3,5)*z(r4,r8,x) + b(5,5)*z(r6,r8,x) )/r3
          select case(name)
          case('3C1')
              vls = b(5,1)*w(r3,r8,x) + b(5,2)*w(r5,r8,x)
          case('3C3')
              vls = b(5,3)*w(r3,r8,x) + b(5,4)*w(r5,r8,x)
          endselect
          vpot(1,1) = vc + (rj-r1)*vls - r2*(rj-r1)/(r2*rj+r1)*vt
          vpot(2,2) = vc - (rj+r2)*vls - r2*(rj+r2)/(r2*rj+r1)*vt
          vpot(1,2) = r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vt
          vpot(2,1) = vpot(1,2)
      case default
!         j > 5
          select case(spin)
          case(0)
              select case(iso)
              case(1)
                  vpot(1,1) =  a(1,1)*y(r2,r8,x) &
                             + a(1,2)*y(r3,r8,x) &
                             + a(1,3)*y(r4,r8,x) &
                             + a(1,4)*y(r5,r8,x) &
                             + a(1,5)*y(r6,r8,x)
              case(0)
                  vpot(1,1) =  b(2,1)*y(r3,r8,x) &
                             + b(2,2)*y(r4,r8,x) &
                             + b(2,3)*y(r5,r8,x) &
                             + b(2,4)*y(r6,r8,x)
              endselect
          case(1)
              select case(iso)
              case(1)
                  vc =  a(4,1)*y(r3,r8,x) &
                      + a(4,2)*y(r4,r8,x) &
                      + a(4,3)*y(r5,r8,x) &
                      + a(4,4)*y(r6,r8,x)
                  vt = ( a(5,1)*z(r4,r8,x) + a(5,2)*z(r6,r8,x) )/r3
              case(0)
                  vc =  b(4,1)*y(r2,r8,x) &
                      + b(4,2)*y(r3,r8,x) &
                      + b(4,3)*y(r4,r8,x) &
                      + b(4,4)*y(r5,r8,x) &
                      + b(4,5)*y(r6,r8,x)
                  vt =  ( b(3,5)*z(r4,r8,x) + b(5,5)*z(r6,r8,x) )/r3
              endselect
              select case(nchan)
              case(1)
                  if (l == (j-1)) vpot(1,1) = vc - r2*(rj-r1)/(r2*rj+r1)*vt
                  if (l == j)     vpot(1,1) = vc + r2*vt
                  if (l == (j+1)) vpot(1,1) = vc - r2*(rj+r2)/(r2*rj+r1)*vt
              case(2)
                  vpot(1,1) = vc - r2*(rj-r1)/(r2*rj+r1)*vt
                  vpot(2,2) = vc - r2*(rj+r2)/(r2*rj+r1)*vt
                  vpot(1,2) = r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vt
                  vpot(2,1) = vpot(1,2)
              endselect
          endselect
      endselect
      vpot(:,:) = mpi * vpot(:,:)

!     All together
      select case(nchan)
      case(1)
          if (spin == 0) then
              if (l == 0) then
                  vpot(1,1) = vpot(1,1) - r3*vspis
              else
                  vpot(1,1) = vpot(1,1) - r3*vspi
              endif
          elseif (l == j) then
              vpot(1,1) = vpot(1,1) + vspi + r2*vtpi
          elseif (name == '3P0') then
              vpot(1,1) = vpot(1,1) + vspi - r2*(rj+r2)/(r2*rj+r1)*vtpi
          endif
      case(2)
          if (l == 0) then
              vpot(1,1) = vpot(1,1) + vspis - r2*(rj-r1)/(r2*rj+r1)*vtpi
          else
              vpot(1,1) = vpot(1,1) + vspi - r2*(rj-r1)/(r2*rj+r1)*vtpi
          endif
          vpot(2,2) = vpot(2,2) + vspi - r2*(rj+r2)/(r2*rj+r1)*vtpi
          vpot(1,2) = vpot(1,2) + r6*sqrt(rj*(rj+r1))/(r2*rj+r1)*vtpi
          vpot(2,1) = vpot(1,2)
      endselect

      contains

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

      function y(n,m,x)

      implicit none
      real(wp), intent(in) :: n,m
      real(wp), intent(in) :: x
      real(wp)             :: y
      real(wp), parameter  :: r1 = 1.0_wp
      real(wp), parameter  :: r2 = 2.0_wp


      if (x < 1.0e-4_wp) then
          y = - n + m/r2 + n**2/(r2*m)
      else
          y = exp(-n*x)/x - exp(-m*x)/x*(r1+(m**2-n**2)*x/(r2*m))
      endif

      end function y

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

      function yp(n,m,x)

      implicit none
      real(wp), intent(in) :: n,m
      real(wp), intent(in) :: x
      real(wp)             :: yp
      real(wp), parameter  :: r1 = 1.0_wp
      real(wp), parameter  :: r2 = 2.0_wp
      real(wp)             :: n2,m2,d


      n2 = n**2
      m2 = m**2

      if (x < 1.0e-4_wp) then
          d  = m*(m2-n2)/(r2*n2)
          yp = - n + m - d + x*(n2/r2 + m2/r2 + m*d)
      else
          yp = exp(-n*x)/x - exp(-m*x)/x*(r1+(m2-n2)*m*x/(r2*n2))
      endif

      end function yp

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

      function w(n,m,x)

      implicit none
      real(wp), intent(in) :: n,m
      real(wp), intent(in) :: x
      real(wp)             :: w
      real(wp)             :: n2,m2,x2
      real(wp), parameter  :: r1 = 1.0_wp
      real(wp), parameter  :: r2 = 2.0_wp
      real(wp), parameter  :: r3 = 3.0_wp
      real(wp), parameter  :: r6 = 6.0_wp
      real(wp), parameter  :: r8 = 8.0_wp


      n2 = n**2
      m2 = m**2

      if (x < 1.0e-4_wp) then
          w = (r2*n - r3*m + m2*m/n2)/r6 + x*(r2*m2 - n2 - m2*m2/n2)/r8
      else
          x2 = x*x
          w =  exp(-n*x)/x*(r1/(n*x)+r1/(n2*x2))          &
             - exp(-m*x)/x*(r1/(m*x)+r1/(m2*x2))*m2/n2    &
             - exp(-m*x)/x*(m2-n2)/(r2*n2)
      endif

      end function w

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

      function z(n,m,x)

      implicit none
      real(wp), intent(in) :: n,m
      real(wp), intent(in) :: x
      real(wp)             :: z
      real(wp)             :: n2,m2,x2
      real(wp), parameter  :: r1 = 1.0_wp
      real(wp), parameter  :: r2 = 2.0_wp
      real(wp), parameter  :: r3 = 3.0_wp
      real(wp), parameter  :: r4 = 4.0_wp
      real(wp), parameter  :: r8 = 8.0_wp


      n2 = n**2
      m2 = m**2

      if (x < 1.0e-4_wp) then
          z  = x*(n2 + r3*m2*m2/n2 - r4*m2)/r8
      else
          x2 = x*x
          z  =  exp(-n*x)/x*(r1+r3/(n*x)+r3/(n2*x2))           &
              - exp(-m*x)/x*(r1+r3/(m*x)+r3/(m2*x2))*m2/n2     &
              - exp(-m*x)*(r1+r1/(m*x))*m*(m2-n2)/(r2*n2)
      endif

      end function z

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

      end subroutine rreid93
