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