      SUBROUTINE REID93P(QI,QF,TYPE,VPOT)
************************************************************************
**    Version 1.0: June 1994                                          **
**    Version 1.1: March 1995                                         **
**    Version 1.2: August 1995                                        **
**    E-mail: info@nn-online.org
**    Reference: Stoks et al. Phys.Rev. C49 (1994) June               **
**                                                                    **
**    Updated Reid potential, regularized with a dipole form factor   **
**    of 8 pion masses, in momentum space on LSJ basis.               **
**                                                                    **
**    INPUT :  QI    center of mass momentum initial state in MeV     **
**    -----    QF    center of mass momentum final   state in MeV     **
**             TYPE  'PP', 'NN', 'NP', or 'PN' (character*2)          **
**             Name partial wave via COMMON/EMANHP/PHNAME (see below) **
**             Maximum total angular momentum J=10 !!                 **
**                                                                    **
**    OUTPUT:  This subroutine returns a 2x2 potential matrix VPOT    **
**    ------   in MeV**-2 which is the partial-wave momentum-space    **
**             potential for the partial wave PHNAME (see below)      **
**--------------------------------------------------------------------**
**    Defining the K-matrix as   2i*mu*q*K = (1-S)(1+S)^-1            **
**    (so 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(kq)        **
**    with                                                            **
**       G(q,k) = P / (E(q) - k^2/2/mu)                               **
**       V(q'k) = 1 / (4*pi) * VPOT(QI=k,QF=q')                       **
**--------------------------------------------------------------------**
**    Potential decomposition in momentum space plane-wave basis:     **
**    V(QF,QI) = VC                                                   **
**             + VS   (SIG1.SIG2)        (only in one-pion-exchnage)  **
**             + VT   [(SIG1.K)(SIG2.K)-K2/3(SIG1.SIG2)]              **
**             + VLS  (i/2)(SIG1+SIG2).N                              **
**                                                                    **
**          K = QF - QI ,   Q = (QF+QI)/2 ,   N = QI X QF = Q X K     **
**                                                                    **
**    NOTE: In the partial wave decomposition we used the             **
**          SYM-convention.                                           **
**          If you use another convention in your Lippmann-Schwinger  **
**          programm, you may need an extra minus sign for the        **
**          the off-diagonal tensor potential VPOT(1,2) and VPOT(2,1) **
**                                                                    **
**    One-pion-exchange part distinguishes between neutral and charged**
**    pion masses, and has coupling constants F0PI=FCPI=0.075         **
**    The delta-function (smeared out due to the form factor) ONLY    **
**    contributes to the S waves.                                     **
**                                                                    **
**    COMMON-block which has to be filled beforehand:                 **
**    +  COMMON/EMANHP/PHNAME                                         **
**              PHNAME is character*3 and contains the name of the    **
**              partial wave in the 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 ...   **
**                                                                    **
************************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      INTEGER SPIN
      CHARACTER TYPE*2, PHNAME*3, PHNAM0*3, CAPS*12, SMAL*12
      REAL*8 VPOT(2,2), ELN(0:12), UL(6,-2:12),ULC(-2:12)
      REAL*8 VPIS(-1:1),VPIT(-2:2), VC(-1:1),VL(-2:2),VT(-2:2)
      REAL*8 PARSPP(5,5),PARSNP(5,5), A(5,5),B(5,5)
      COMMON/EMANHP/PHNAME
      DATA CAPS/'CSPDFGHIKLMN'/, SMAL/'cspdfghiklmn'/
      DATA F0PI/0.075D0/, FCPI/0.075D0/, PI/3.14159265358979D0/
      DATA PIOM,PIOMC,PIOMS/134.9739D0,139.5675D0,139.5675D0/
      DATA ICAL/0/, PHNAM0/'***'/, UL/90*0D0/, ULC/15*0D0/
      DATA PARSPP/
     1  .1756084D0,-.1414234D2, .1518489D3,-.6868230D3, .1104157D4
     2,-.4224976D2, .2072246D3,-.3354364D3,-.1989250D1,-.6178469D2
     3, .2912845D2, .1511690D3, .8151964D1, .5832103D2,-.2074743D2
     4,-.5840566D0,-.1029310D2, .2263391D2, .2316915D2,-.1959172D1
     5,-.2608488D1, .1090858D2,-.4374212D0,-.2148862D2,-.6584788D0/
      DATA PARSNP/
     1 -.2234989D2, .2551761D3,-.1063549D4, .1609196D4,-.3505968D1
     2,-.4248612D1,-.5352001D1, .1827642D3,-.3927086D3, .5812273D2
     3,-.2904577D1, .3802497D2, .3395927D0, .8318097D0, .1923895D1
     4, .0913746D0,-.1274773D2, .1458600D3,-.6432461D3, .1022217D4
     5,-.0461640D0, .7950192D1,-.1925573D1, .5066234D2, .83598955D1/
      SAVE A,B, PIOMS2,PIOMM2, NCHAN,L,SPIN,J,ISO
      SAVE JMM,JM,JP,JPP,TJMM,TJM,TJ,TJP,TJPP,TJJ,JMAX

      VPOT(1,1) = 0D0
      VPOT(1,2) = 0D0
      VPOT(2,1) = 0D0
      VPOT(2,2) = 0D0

      IF(ICAL.EQ.0) THEN
        DO 1 I1=1,5
          DO 1 I2=1,5
            A(I1,I2)=PARSPP(I2,I1)
  1         B(I1,I2)=PARSNP(I2,I1)
        PIOMS2=PIOMS*PIOMS
        PIOMM=(PIOM+2D0*PIOMC)/3D0
        PIOMM2=PIOMM*PIOMM
        ICAL=1
      ENDIF
 
      IF(PHNAME.NE.PHNAM0) THEN
        DO 5 LL=1,12
          IF(PHNAME(2:2).EQ.SMAL(LL:LL)) THEN
            WRITE(PHNAME(2:2),'(A1)') CAPS(LL:LL)
          ENDIF
    5   CONTINUE
        PHNAM0=PHNAME
        NCHAN=1
        IF(PHNAME(2:2).EQ.'C') NCHAN=2
        IF(PHNAME(1:1).EQ.'1') SPIN=0
        IF(PHNAME(1:1).EQ.'3') SPIN=1
        READ(PHNAME,'(2X,I1)') J
        IF(J.GT.10) WRITE(*,*)
     .              '**** Partial wave exceeds allowable maximum J=10'
        IF(J.GT.10) STOP
        L=J
        IF(PHNAME.EQ.'3P0') L=1
        IF(NCHAN.EQ.2) L=J-1
        ISO=MOD(SPIN+L+1,2)
        JMM=J-2
        JM =J-1
        JP =J+1
        JPP=J+2
        TJMM=2D0*J-3D0
        TJM =2D0*J-1D0
        TJ  =2D0*J+1D0
        TJP =2D0*J+3D0
        TJPP=2D0*J+5D0
        TJJ =DSQRT(J*(J+1D0))
        JMAX=J+2
      ENDIF

      IF(TYPE.EQ.'pp') THEN
        TYPE='PP'
      ELSEIF(TYPE.EQ.'np') THEN
        TYPE='NP'
      ELSEIF(TYPE.EQ.'pn') THEN
        TYPE='PN'
      ELSEIF(TYPE.EQ.'nn') THEN
        TYPE='NN'
      ENDIF
 
      QI2=QI*QI
      QF2=QF*QF
      QIQF=QI*QF
      QI2F2=QI2+QF2
      S2PSI=2D0*QIQF/QI2F2
      SPSI2=QF2/QI2F2
      CPSI2=QI2/QI2F2
 
***   Neutral one-pion-exchange part
      AMES2=PIOM*PIOM
      ALAM2=64D0*AMES2
      X=0.5D0*(QI2F2+AMES2)/QIQF
      Y=0.5D0*(QI2F2+ALAM2)/QIQF
      CALL SDIP(X,Y,JMAX,ELN)
      DO 10 IE=0,JMAX
   10   UL(1,IE)=ELN(IE)
      FAC=2D0*PI*F0PI/PIOMS2/QIQF
      DO 11 LL=-1,1
   11   VPIS(LL)= FAC*UL(1,J+LL)*PIOM*PIOM/3D0
      VPISL0=4D0*PI*F0PI/PIOMS2/3D0 * (Y-X)*(Y-X)/(1D0-Y*Y)
      DO 12 LL=-2,2
   12   VPIT(LL)=-FAC*UL(1,J+LL)
 
***   Charged one-pion-exchange part
      IF(TYPE.EQ.'NP' .OR. TYPE.EQ.'PN') THEN
        AMES2=PIOMC*PIOMC
        ALAM2=64D0*AMES2
        X=0.5D0*(QI2F2+AMES2)/QIQF
        Y=0.5D0*(QI2F2+ALAM2)/QIQF
        CALL SDIP(X,Y,JMAX,ELN)
        DO 20 IE=0,JMAX
   20     ULC(IE)=ELN(IE)
        FAC=(4D0*ISO-2D0)*2D0*PI*FCPI/PIOMS2/QIQF
        DO 21 LL=-1,1
   21     VPIS(LL)= FAC*ULC(J+LL)*PIOMC*PIOMC/3D0 - VPIS(LL)
        VPISL0=4D0*PI*FCPI/PIOMS2/3D0 * (4D0*ISO-2D0)
     .        *(Y-X)*(Y-X)/(1D0-Y*Y) - VPISL0
        DO 22 LL=-2,2
   22     VPIT(LL)=-FAC*ULC(J+LL) - VPIT(LL)
      ENDIF
 
***   Other Yukawa's with multiples of mean pion mass
      ALAM2=64D0*PIOMM2
      Y=0.5D0*(QI2F2+ALAM2)/QIQF
      DO 30 IM=2,6
        X=0.5D0*(QI2F2+IM*IM*PIOMM2)/QIQF
        CALL SDIP(X,Y,JMAX,ELN)
        DO 35 IE=0,JMAX
   35     UL(IM,IE)=ELN(IE)
   30 CONTINUE
 
      FAC=2D0*PI/QIQF
 
***   Potential for each partial wave separately
      IF(PHNAME.EQ.'1S0') THEN
        IF(TYPE.EQ.'PP' .OR. TYPE.EQ.'NN') THEN
          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) )
        ELSEIF(TYPE.EQ.'NP' .OR. TYPE.EQ.'PN') THEN
          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) )
        ENDIF
      ELSEIF(PHNAME.EQ.'1D2') THEN
        VPOT(1,1)=FAC*(A(2,1)*UL(4,J)+A(2,2)*UL(5,J)+A(2,3)*UL(6,J) )
      ELSEIF(PHNAME.EQ.'1G4') THEN
        VPOT(1,1)=FAC* A(2,4)*UL(3,J)
      ELSEIF(PHNAME.EQ.'3P0') THEN
        VTEN=UL(3,JP)-0.5D0*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/9D0/PIOMM2)*QI2F2/3D0 )
      ELSEIF(PHNAME.EQ.'3P1') THEN
        VTEN=UL(3,J)-0.5D0*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/9D0/PIOMM2)*QI2F2/3D0 )
      ELSEIF(PHNAME.EQ.'3F3') THEN
        VPOT(1,1)=FAC* A(4,5)*UL(3,J)
      ELSEIF(ISO.EQ.1 .AND. SPIN.EQ.1) THEN
        DO 101 LL=-1,1
  101     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) )
        DO 102 LL=-2,2
          VT(LL)=-FAC*(A(5,1)*UL(4,J+LL)/16D0/PIOMM2
     .                +A(5,2)*UL(6,J+LL)/36D0/PIOMM2 )
          IF(PHNAME.EQ.'3C2') THEN
            VL(LL)=FAC*(A(5,3)*UL(3,J+LL)/ 9D0/PIOMM2
     .                 +A(5,4)*UL(5,J+LL)/25D0/PIOMM2 )
          ELSEIF(PHNAME.EQ.'3C4') THEN
            VL(LL)=FAC* A(5,5)*UL(3,J+LL)/ 9D0/PIOMM2
          ELSE
            VL(LL)=0D0
          ENDIF
  102   CONTINUE
        IF(NCHAN.EQ.1)
     .    VPOT(1,1)=VC(0) - QIQF*(VL(-1)-VL(1))/TJ + 2D0/3D0*QI2F2*
     .              (VT(0)-0.5D0*S2PSI*(TJP*VT(-1)+TJM*VT(1))/TJ)
      ELSEIF(PHNAME.EQ.'1P1') THEN
        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) )
      ELSEIF(PHNAME.EQ.'1F3') THEN
        VPOT(1,1)=FAC*(B(1,5)*UL(3,J)+B(2,5)*UL(5,J) )
      ELSEIF(PHNAME.EQ.'3D2') THEN
        VTEN=UL(3,J)-0.5D0*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/9D0/PIOMM2)*QI2F2/3D0 )
      ELSEIF(PHNAME.EQ.'3G4') THEN
        VPOT(1,1)=FAC* B(3,4)*UL(3,J)
      ELSEIF(ISO.EQ.0 .AND. SPIN.EQ.1) THEN
        DO 201 LL=-1,1
  201     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) )
        DO 202 LL=-2,2
          VT(LL)=-FAC*(B(3,5)*UL(4,J+LL)/16D0/PIOMM2
     .                +B(5,5)*UL(6,J+LL)/36D0/PIOMM2 )
          IF(PHNAME.EQ.'3C1') THEN
            VL(LL)=FAC*(B(5,1)*UL(3,J+LL)/ 9D0/PIOMM2
     .                 +B(5,2)*UL(5,J+LL)/25D0/PIOMM2 )
          ELSEIF(PHNAME.EQ.'3C3') THEN
            VL(LL)=FAC*(B(5,3)*UL(3,J+LL)/ 9D0/PIOMM2
     .                 +B(5,4)*UL(5,J+LL)/25D0/PIOMM2 )
          ELSE
            VL(LL)=0D0
          ENDIF
  202   CONTINUE
        IF(NCHAN.EQ.1)
     .    VPOT(1,1)=VC(0) - QIQF*(VL(-1)-VL(1))/TJ + 2D0/3D0*QI2F2*
     .              (VT(0)-0.5D0*S2PSI*(TJP*VT(-1)+TJM*VT(1))/TJ)
      ELSEIF(J.GE.5 .AND. SPIN.EQ.0) THEN
        IF(ISO.EQ.1) THEN
          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) )
        ELSEIF(ISO.EQ.0) THEN
          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) )
        ENDIF
      ENDIF
 
      IF(NCHAN.EQ.2) THEN
        VPOT(1,1)=VC(-1) + QIQF*JM/TJM*(VL(-2)-VL(0)) +
     .            2D0/3D0*QI2F2*JM/TJ*
     .            (-VT(-1)+0.5D0*S2PSI*(TJMM*VT(0)+TJ*VT(-2))/TJM)
        VPOT(1,2)=-2D0*QI2F2*TJJ/TJ*
     .            (-S2PSI*VT(0)+CPSI2*VT(-1)+SPSI2*VT(1))
        VPOT(2,1)=-2D0*QI2F2*TJJ/TJ*
     .            (-S2PSI*VT(0)+SPSI2*VT(-1)+CPSI2*VT(1))
        VPOT(2,2)=VC(1) - QIQF*JPP/TJP*(VL(0)-VL(2)) +
     .            2D0/3D0*QI2F2*JPP/TJ*
     .            (-VT(1)+0.5D0*S2PSI*(TJPP*VT(0)+TJ*VT(2))/TJP)
      ENDIF
 
***   Add one-pion-exchange part
      IF(NCHAN.EQ.1) THEN
        IF(SPIN.EQ.0) THEN
          VPOT(1,1) = VPOT(1,1) - 3D0*VPIS(0)
          IF(L.EQ.0) VPOT(1,1) = VPOT(1,1) - 3D0*VPISL0
        ELSEIF(L.EQ.J) THEN
          VPOT(1,1) = VPOT(1,1) + VPIS(0) + 2D0/3D0*QI2F2*
     .        (VPIT(0)-0.5D0*S2PSI*(TJP*VPIT(-1)+TJM*VPIT(1))/TJ)
        ELSEIF(PHNAME.EQ.'3P0') THEN
          VPOT(1,1) = VPOT(1,1) + VPIS(1) + 2D0/3D0*QI2F2*JPP/TJ*
     .        (-VPIT(1)+0.5D0*S2PSI*(TJPP*VPIT(0)+TJ*VPIT(2))/TJP)
        ENDIF
      ELSE
        VPOT(1,1) = VPOT(1,1) + VPIS(-1) + 2D0/3D0*QI2F2*JM/TJ*
     .      (-VPIT(-1)+0.5D0*S2PSI*(TJMM*VPIT(0)+TJ*VPIT(-2))/TJM)
        IF(L.EQ.0) VPOT(1,1) = VPOT(1,1) + VPISL0
        VPOT(1,2) = VPOT(1,2) - 2D0*QI2F2*TJJ/TJ*
     .      (-S2PSI*VPIT(0)+CPSI2*VPIT(-1)+SPSI2*VPIT(1))
        VPOT(2,1) = VPOT(2,1) - 2D0*QI2F2*TJJ/TJ*
     .      (-S2PSI*VPIT(0)+SPSI2*VPIT(-1)+CPSI2*VPIT(1))
        VPOT(2,2) = VPOT(2,2) + VPIS(1) + 2D0/3D0*QI2F2*JPP/TJ*
     .      (-VPIT(1)+0.5D0*S2PSI*(TJPP*VPIT(0)+TJ*VPIT(2))/TJP)
      ENDIF
 
      RETURN
      END
*******************************************************************
      SUBROUTINE SDIP(X,Y,JMAX,ELN)
      IMPLICIT REAL*8 (A-H,O-Z)
      PARAMETER (NOGP=64)
      DIMENSION P(12), ELN(0:12), ZZ(NOGP), WZ(NOGP)
      DATA ICAL/0/
      SAVE ZZ,WZ
 
      QX0 = 0.5D0*DLOG((X+1D0)/(X-1D0))
      QY0 = 0.5D0*DLOG((Y+1D0)/(Y-1D0))
      ELN(0)=QX0-QY0+(Y-X)/(1D0-Y*Y)
      IF(JMAX.EQ.0) RETURN
      IF(ICAL.EQ.0) THEN
        CALL GAUS(NOGP/2,ZZ,WZ,2)
        ICAL=1
      ENDIF
      DO 1 L=1,JMAX
    1   ELN(L)=0D0
      DO 2 IZ=1,NOGP
        FORM=(Y-X)/(Y-ZZ(IZ))
        FORM=FORM*FORM/(X-ZZ(IZ))
C*      Calculate Legendre polynomials first kind and integrate
        P(1)=ZZ(IZ)
        P(2)=1.5D0*ZZ(IZ)*ZZ(IZ)-0.5D0
        ELN(1)=ELN(1)+0.5D0*WZ(IZ)*FORM*P(1)
        ELN(2)=ELN(2)+0.5D0*WZ(IZ)*FORM*P(2)
        DO 3 L=3,JMAX
          P(L)=( (2*L-1)*ZZ(IZ)*P(L-1) - (L-1)*P(L-2) ) / L
          ELN(L)=ELN(L)+0.5D0*WZ(IZ)*FORM*P(L)
    3   CONTINUE
    2 CONTINUE
      RETURN
      END
*******************************************************************
      SUBROUTINE GAUS(N,X,W,NT)
C*      The abcissas and weights for Gaussian integration,
C*      N=4,8,12,16,20,32,48 only:
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION WGAUSS(61), WEULRA(32),WEULRB(48)
      DIMENSION X(1),W(1)
      DATA WGAUSS
     1/.33998104358486,.65214515486255,.86113631159405,.34785484513745
     *,.18343464249565,.36268378337836,.52553240991633,.31370664587789
     *,.79666647741362,.22238103445337,.96028985649753,.10122853629038
     *,.12523340851147,.24914704581340,.36783149899818,.23349253653835
     *,.58731795428661,.20316742672307,.76990267419430,.16007832854335
     *,.90411725637048,.10693932599532,.98156063424672,.04717533638651
     *,.09501250983764,.18945061045507,.28160355077926,.18260341504492
     *,.45801677765722,.16915651939500,.61787624440264,.14959598881658
     *,.75540440835500,.12462897125553,.86563120238783,.09515851168249
     *,.94457502307323,.06225352393865,.98940093499165,.02715245941175
     *,.07652652113350,.15275338713072,.22778585114164,.14917298647260
     *,.37370608871542,.14209610931838,.51086700195082,.13168863844918
     *,.63605368072652,.11819453196152,.74633190646015,.10193011981724
     *,.83911697182222,.08327674157670,.91223442825132,.06267204833411
     *,.96397192727791,.04060142980039,.99312859918509,.01761400713915
     *,.00/

C  32 POINTS:
      DATA WEULRA
     1/.04830766568774,.09654008851473,.14447196158280,.09563872007927
     *,.23928736225214,.09384439908080,.33186860228213,.09117387869576
     *,.42135127613064,.08765209300440,.50689990893223,.08331192422695
     *,.58771575724076,.07819389578707,.66304426693022,.07234579410885
     *,.73218211874029,.06582222277636,.79448379596794,.05868409347854
     *,.84936761373257,.05099805926238,.89632115576605,.04283589802223
     *,.93490607593774,.03427386291302,.96476225558751,.02539206530926
     *,.98561151154527,.01627439473091,.99726386184948,.00701861000947
     */

C  48 POINTS:
      DATA WEULRB
     1/.03238017096287,.06473769681268,.09700469920946,.06446616443595
     *,.16122235606889,.06392423858465,.22476379039469,.06311419228625
     *,.28736248735546,.06203942315989,.34875588629216,.06070443916589
     *,.40868648199071,.05911483969840,.46690290475096,.05727729210040
     *,.52316097472223,.05519950369998,.57722472608397,.05289018948519
     *,.62886739677651,.05035903555385,.67787237963266,.04761665849249
     *,.72403413092381,.04467456085669,.76715903251574,.04154508294346
     *,.80706620402944,.03824135106583,.84358826162439,.03477722256477
     *,.87657202027425,.03116722783280,.90587913671557,.02742650970836
     *,.93138669070655,.02357076083932,.95298770316043,.01961616045736
     *,.97059159254625,.01557931572294,.98412458372283,.01147723457923
     *,.99353017226635,.00732755390128,.99877100725243,.00315334605231
     */

      IF(MOD(N,4).NE.0) THEN
C*      Simpson's rule: only odd N
        H=2D0/(N+1D0)
        X(1)=-1D0+H
        W23=2D0*H/3D0
        W43=2D0*W23
        W(1)=W43
        NN=(N-1)/2
        DO 10 I=1,NN
          J=2*I
          X(J)=X(J-1)+H
          X(J+1)=X(J)+H
          W(J)=W23
          W(J+1)=W43
   10   CONTINUE
        RETURN
      ENDIF

      IF(N.LE.20) THEN
        K=N/4-1
        INIT=2*K*(K+1)+1
        M=N/2
        DO 1 I=1,M
          J=INIT-2+I*2
          X(M+I)=WGAUSS(J)
          X(M+1-I)=-WGAUSS(J)
          W(M+I)=WGAUSS(J+1)
          W(M+1-I)=WGAUSS(J+1)
    1   CONTINUE
      ENDIF
      IF(N.EQ.32) THEN
        M=N/2
        DO 2 I=1,M
          J=2*I-1
          X(M+I)=WEULRA(J)
          X(M+1-I)=-WEULRA(J)
          W(M+I)=WEULRA(J+1)
          W(M+1-I)=WEULRA(J+1)
    2   CONTINUE
      ENDIF
      IF(N.EQ.48) THEN
        M=N/2
        DO 3 I=1,M
          J=2*I-1
          X(M+I)=WEULRB(J)
          X(M+1-I)=-WEULRB(J)
          W(M+I)=WEULRB(J+1)
          W(M+1-I)=WEULRB(J+1)
    3   CONTINUE
      ENDIF
      IF(NT.EQ.1) RETURN

      DO 5 I=1,N
        X(N+I)=(X(I)+1D0)/2D0
        X(I)=(X(I)-1D0)/2D0
        W(N+I)=W(I)/2D0
        W(I)=W(N+I)
    5 CONTINUE

      RETURN
      END
