************************************************************************
* 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 77 version 1.0: July 1994
*            version 1.1: March 1995
*            version 1.2: August 1995
*            version 2.0: April 2000
*
* E-mail: info@nn-online.org
* WWW: http://nn-online.org
*-----------------------------------------------------------------------
* IN:  real*8        r         = distance in fermi
*      character*2   type      = reaction, 'PP', 'NP', 'PN', or 'NN'
*      character*3   name      = name of the partial wave (see below),
*                                maximum total angular momentum j = 9
* OUT: real*8        vpot(2,2) = potential matrix in MeV on LSJ-basis
*-----------------------------------------------------------------------
* 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 - epsilon_1 - 3D1 channel
*       3C2 denotes   3P2 - epsilon_2 - 3F2 channel
*       ...
************************************************************************

      subroutine rreid93(r,name,type,vpot)

      implicit none

      real*8  r,vpot(2,2)
      character*2  type
      character*3  name
      character*12  large,small
      integer  l,spin,j,nchan,iso
      logical  first
      real*8  ri,rj,vc,vt,vls,vtpi,x,x0,xc,mpi,vspi,vspis
      real*8  r1,r2,r3,r4,r6
      real*8  f0pi,fcpi,hbc,mpi0,mpic,mpis
      real*8  a(5,5),b(5,5),parspp(5,5),parsnp(5,5)
      real*8  r93_y,r93_yp,r93_w,r93_z

      save a,b

      data f0pi /0.075d0/
      data fcpi /0.075d0/
      data hbc /197.327053d0/
      data mpi0 /134.9739d0/
      data mpic /139.5675d0/
      data mpis /139.5675d0/
      data r1 /1.0d0/
      data r2 /2.0d0/
      data r3 /3.0d0/
      data r4 /4.0d0/
      data r6 /6.0d0/
      data first /.true./
      data large /'CSPDFGHIKLMN'/
      data small /'cspdfghiklmn'/
      data parspp
     ./ 0.1756084d0,-0.1414234d2, 0.1518489d3,-0.6868230d3, 0.1104157d4,
     . -0.4224976d2, 0.2072246d3,-0.3354364d3,-0.1989250d1,-0.6178469d2,
     .  0.2912845d2, 0.1511690d3, 0.8151964d1, 0.5832103d2,-0.2074743d2,
     . -0.5840566d0,-0.1029310d2, 0.2263391d2, 0.2316915d2,-0.1959172d1,
     . -0.2608488d1, 0.1090858d2,-0.4374212d0,-0.2148862d2,-0.6584788d0/
      data parsnp
     ./-0.2234989d2, 0.2551761d3,-0.1063549d4, 0.1609196d4,-0.3505968d1,
     . -0.4248612d1,-0.5352001d1, 0.1827642d3,-0.3927086d3, 0.5812273d2,
     . -0.2904577d1, 0.3802497d2, 0.3395927d0, 0.8318097d0, 0.1923895d1,
     .  0.0913746d0,-0.1274773d2, 0.1458600d3,-0.6432461d3, 0.1022217d4,
     . -0.0461640d0, 0.7950192d1,-0.1925573d1, 0.5066234d2, 0.8359896d1/

*     The parameters
      if (first) then
          do j=1,5
          do l=1,5
              a(j,l) = parspp(l,j)
              b(j,l) = parsnp(l,j)
          enddo
          enddo
          first = .false.
      endif

*     Determine the quantumnumbers
      do l=1,12
          if (name(2:2).eq.small(l:l)) name(2:2) = large(l:l)
      enddo
      nchan = 1
      if (name(2:2).eq.'C') nchan = 2
      if (name(1:1).eq.'1') spin = 0
      if (name(1:1).eq.'3') spin = 1
      read(name,'(2x,i1)') j
      l = j
      if (name.eq.'3P0') l = 1
      if (nchan.eq.2) l = j - 1
      iso = mod(spin+l+1,2)

      ri = real(iso)
      rj = real(j)

*     OPE potential
      x0    = mpi0/hbc * r
      vspis = f0pi*(mpi0/mpis)**2*mpi0/r3*r93_yp(1,8,x0)
      vspi  = f0pi*(mpi0/mpis)**2*mpi0/r3*r93_y(1,8,x0)
      vtpi  = f0pi*(mpi0/mpis)**2*mpi0/r3*r93_z(1,8,x0)
      if (type.eq.'NP' .or. type.eq.'PN') then
          xc    = mpic/hbc * r
          vspis = (r4*ri-r2)*fcpi*(mpic/mpis)**2*mpic/r3*r93_yp(1,8,xc)
     .            - vspis
          vspi  = (r4*ri-r2)*fcpi*(mpic/mpis)**2*mpic/r3*r93_y(1,8,xc)
     .            - vspi
          vtpi  = (r4*ri-r2)*fcpi*(mpic/mpis)**2*mpic/r3*r93_z(1,8,xc)
     .            - vtpi
      endif

*     Start clean
      vpot(1,1) = 0.0d0
      vpot(1,2) = 0.0d0
      vpot(2,1) = 0.0d0
      vpot(2,2) = 0.0d0

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

      if (name.eq.'1S0') then
          if (type.eq.'PP' .or. type.eq.'NN') then
              vpot(1,1) =   a(1,1)*r93_y(2,8,x)
     .                    + a(1,2)*r93_y(3,8,x)
     .                    + a(1,3)*r93_y(4,8,x)
     .                    + a(1,4)*r93_y(5,8,x)
     .                    + a(1,5)*r93_y(6,8,x)
          elseif (type.eq.'NP' .or. type.eq.'PN') then
              vpot(1,1) =   b(1,1)*r93_y(3,8,x)
     .                    + b(1,2)*r93_y(4,8,x)
     .                    + b(1,3)*r93_y(5,8,x)
     .                    + b(1,4)*r93_y(6,8,x)
          endif
      elseif (name.eq.'1D2') then
          vpot(1,1) =   a(2,1)*r93_y(4,8,x)
     .                + a(2,2)*r93_y(5,8,x)
     .                + a(2,3)*r93_y(6,8,x)
      elseif (name.eq.'1G4') then
          vpot(1,1) =   a(2,4)*r93_y(3,8,x)
      elseif (name.eq.'3P0') then
          vpot(1,1) =   a(3,1)*r93_y(3,8,x)
     .                + a(3,2)*r93_y(5,8,x)
     .                + a(2,5)*r93_z(3,8,x)/r3
      elseif (name.eq.'3P1') then
          vpot(1,1) =   a(3,3)*r93_y(3,8,x)
     .                + a(3,4)*r93_y(5,8,x)
     .                + a(3,5)*r93_z(3,8,x)/r3
      elseif (name.eq.'3F3') then
          vpot(1,1) =   a(4,5)*r93_y(3,8,x)
      elseif (name.eq.'3C2' .or. name.eq.'3C4') then
          vc =  a(4,1)*r93_y(3,8,x)
     .        + a(4,2)*r93_y(4,8,x)
     .        + a(4,3)*r93_y(5,8,x)
     .        + a(4,4)*r93_y(6,8,x)
          vt = ( a(5,1)*r93_z(4,8,x) + a(5,2)*r93_z(6,8,x) )/r3
          if (name.eq.'3C2') vls = a(5,3)*r93_w(3,8,x) +
     .                             a(5,4)*r93_w(5,8,x)
          if (name.eq.'3C4') vls = a(5,5)*r93_w(3,8,x)
          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)
      elseif (name.eq.'1P1') then
          vpot(1,1) =   b(2,1)*r93_y(3,8,x)
     .                + b(2,2)*r93_y(4,8,x)
     .                + b(2,3)*r93_y(5,8,x)
     .                + b(2,4)*r93_y(6,8,x)
      elseif (name.eq.'1F3') then
          vpot(1,1) =   b(1,5)*r93_y(3,8,x)
     .                + b(2,5)*r93_y(5,8,x)
      elseif (name.eq.'3D2') then
          vpot(1,1) =   b(3,1)*r93_y(3,8,x)
     .                + b(3,2)*r93_y(5,8,x)
     .                + b(3,3)*r93_z(3,8,x)/r3
      elseif (name.eq.'3G4') then
          vpot(1,1) =   b(3,4)*r93_y(3,8,x)
      elseif (name.eq.'3C1' .or. name.eq.'3C3') then
          vc =  b(4,1)*r93_y(2,8,x)
     .        + b(4,2)*r93_y(3,8,x)
     .        + b(4,3)*r93_y(4,8,x)
     .        + b(4,4)*r93_y(5,8,x)
     .        + b(4,5)*r93_y(6,8,x)
          vt = ( b(3,5)*r93_z(4,8,x) + b(5,5)*r93_z(6,8,x) )/r3
          if (name.eq.'3C1') vls = b(5,1)*r93_w(3,8,x) +
     .                             b(5,2)*r93_w(5,8,x)
          if (name.eq.'3C3') vls = b(5,3)*r93_w(3,8,x) +
     .                             b(5,4)*r93_w(5,8,x)
          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)
      elseif (spin.eq.0) then
*         j >= 5, spin = 0
          if (iso.eq.1) then
              vpot(1,1) =  a(1,1)*r93_y(2,8,x)
     .                   + a(1,2)*r93_y(3,8,x)
     .                   + a(1,3)*r93_y(4,8,x)
     .                   + a(1,4)*r93_y(5,8,x)
     .                   + a(1,5)*r93_y(6,8,x)
          elseif (iso.eq.0) then
              vpot(1,1) =  b(2,1)*r93_y(3,8,x)
     .                   + b(2,2)*r93_y(4,8,x)
     .                   + b(2,3)*r93_y(5,8,x)
     .                   + b(2,4)*r93_y(6,8,x)
          endif
      elseif (spin.eq.1) then
*         j >= 5, spin = 1
          if (iso.eq.1) then
              vc =  a(4,1)*r93_y(3,8,x)
     .            + a(4,2)*r93_y(4,8,x)
     .            + a(4,3)*r93_y(5,8,x)
     .            + a(4,4)*r93_y(6,8,x)
              vt = ( a(5,1)*r93_z(4,8,x) + a(5,2)*r93_z(6,8,x) )/r3
          elseif (iso.eq.0) then
              vc =  b(4,1)*r93_y(2,8,x)
     .            + b(4,2)*r93_y(3,8,x)
     .            + b(4,3)*r93_y(4,8,x)
     .            + b(4,4)*r93_y(5,8,x)
     .            + b(4,5)*r93_y(6,8,x)
              vt = ( b(3,5)*r93_z(4,8,x) + b(5,5)*r93_z(6,8,x) )/r3
          endif
          if (nchan.eq.1) then
              if (l.eq.(j-1)) vpot(1,1) = vc - r2*(rj-r1)/(r2*rj+r1)*vt
              if (l.eq. j)    vpot(1,1) = vc + r2*vt
              if (l.eq.(j+1)) vpot(1,1) = vc - r2*(rj+r2)/(r2*rj+r1)*vt
          elseif (nchan.eq.2) then
              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)
          endif
      endif
      vpot(1,1) = mpi * vpot(1,1)
      vpot(1,2) = mpi * vpot(1,2)
      vpot(2,1) = mpi * vpot(2,1)
      vpot(2,2) = mpi * vpot(2,2)

*     All together
      if (nchan.eq.1) then
          if (spin.eq.0) then
              if (l.eq.0) then
                  vpot(1,1) = vpot(1,1) - r3*vspis
              else
                  vpot(1,1) = vpot(1,1) - r3*vspi
              endif
          elseif (l.eq.j) then
              vpot(1,1) = vpot(1,1) + vspi + r2*vtpi
          elseif (name.eq.'3P0') then
              vpot(1,1) = vpot(1,1) + vspi - r2*(rj+r2)/(r2*rj+r1)*vtpi
          endif
      elseif (nchan.eq.2) then
          if (l.eq.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)
      endif

      end

************************************************************************

      function r93_y(in,im,x)

      implicit none

      integer  in,im
      real*8  x,r93_y
      real*8  n,m,r1,r2
      data r1 /1.0d0/
      data r2 /2.0d0/

      n = real(in)
      m = real(im)

      if (x < 1.0d-4) then
          r93_y = - n + m/r2 + n*n/(r2*m)
      else
          r93_y = exp(-n*x)/x - exp(-m*x)/x*(r1+(m*m-n*n)*x/(r2*m))
      endif

      end

************************************************************************

      function r93_yp(in,im,x)

      implicit none

      integer  in,im
      real*8  x,r93_yp
      real*8  n,m,n2,m2,d,r1,r2
      data r1 /1.0d0/
      data r2 /2.0d0/

      n  = real(in)
      m  = real(im)
      n2 = n*n
      m2 = m*m

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

      end

************************************************************************

      function r93_w(in,im,x)

      implicit none

      integer  in,im
      real*8  x,r93_w
      real*8  n,m,n2,m2,x2,r1,r2,r3,r6,r8
      data r1 /1.0d0/
      data r2 /2.0d0/
      data r3 /3.0d0/
      data r6 /6.0d0/
      data r8 /8.0d0/

      n  = real(in)
      m  = real(im)
      n2 = n*n
      m2 = m*m

      if (x < 1.0d-4) then
          r93_w =    (r2*n - r3*m + m2*m/n2)/r6 +
     .             x*(r2*m2 - n2 - m2*m2/n2)/r8
      else
          x2 = x*x
          r93_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 r93_z(in,im,x)

      implicit none

      integer  in,im
      real*8  r93_z,x
      real*8  n,m,n2,m2,x2,r1,r2,r3,r4,r8
      data r1 /1.0d0/
      data r2 /2.0d0/
      data r3 /3.0d0/
      data r4 /4.0d0/
      data r8 /8.0d0/

      n  = real(in)
      m  = real(im)
      n2 = n*n
      m2 = m*m

      if (x < 1.0d-4) then
          r93_z = x*(n2 + r3*m2*m2/n2 - r4*m2)/r8
      else
          x2 = x*x
          r93_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

************************************************************************
