SUBROUTINE PNIJMLSJ(QI,QF,TYPE,VPOT) ************************************************************************ ** NIJMEGEN NUCLEON-NUCLEON POTENTIAL PROGRAM ** ** ------------------------------------------ ** ** Reference: Stoks et al., Phys. Rev. C 49 (1994), 2950 ** ** ** ** Version 1.0: June 1994 ** ** Version 1.1: March 1995 ** ** Version 1.2: August 1995 ** ** Version 1.3: November 1999 ** ** Version 1.4: April 2000 ** ** Version 1.5: February 2001 ** ** E-mail: info@nn-online.org ** WWW: NN-OnLine: http://nn-online.org ** ** ** Refined and extended version of the 1978 Nijmegen NN potential ** ** in momentum space on LSJ basis. ** ** Basic references for formulas and set up of this program can be ** ** found in Phys.Rev.D17 (1978) 768 and Phys.Rev.C40 (1989) 2226. ** ** For momentum-space projection see THEF-NIJM 91.05 preprint ** ** (available from NN-OnLine: http://nn-online.sci.kun.nl) ** **--------------------------------------------------------------------** ** This code assumes local variables are initialized and saved. ** ** For some compilers this is default behaviour, other compilers ** ** require you to specify special compiler switches when compiling.** **--------------------------------------------------------------------** ** 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=9 !! ** ** ** ** OUTPUT: central VC, spin-spin VS, tensor VT, spin-orbit VLS, ** ** ------ asymmetric spin-orbit VLA, quadratic spin-orbit VQ, ** ** and extra quadratic spin-orbit part VQP. ** ** All potentials are calculated via a partial-wave ** ** decomposition from Jmin up to Jmax, communicated via ** ** COMMON/POTMOM/VC(-1:1),VS(-1:1),VT(-2:2),VLS(-2:2), ** ** VLA(-2:2),VQ(-3:3),VQP(-2:2) ** ** The 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) ** ** + VT [(SIG1.K)(SIG2.K)-K2/3(SIG1.SIG2)] ** ** + VLS (i/2)(SIG1+SIG2).N ** ** + VLA (i/2)(SIG1-SIG2).N (NOT USED !!!) ** ** + VQ (SIG1.N)(SIG2.N) ** ** + VQP [(SIG1.Q)(SIG2.Q)-Q2(SIG1.SIG2) ** ** -(1/4)(SIG1.K)(SIG2.K)+(1/4)K2(SIG1.SIG2)] ** ** ** ** 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) ** ** ** ** NQ12 integer which opts the full Fourier transform for the ** ** quadratic spin-orbit Q12 operator ** ** 0: Fourier from (SIG1.N)(SIG2.N) --> Q12 ==> approximation ** ** 1: exact Fourier Q12 --> (SIG1.N)(SIG2.N)+extra ** ** so same potential in coordinate space as in momentum space** ** ( Presently NQ12=1 is included as a DATA statement ) ** ** ** ** COMMON-blocks which have to be filled beforehand: ** ** + COMMON/CHOICE/IDPAR ** ** IDPAR is an integer and denotes the various different ** ** models that can be chosen. ** ** IDPAR=0: nijm93: potential for pp and np together. ** ** including a phenomenological ** ** parameter to give the 1S0 pp and np ** ** phase shift/scattering length ** ** difference ** ** 1: nijmI : Reidlike model, each partial wave has** ** its own parameterset ** ** 2: nijmII: like nijmI, but fully local ** ** + 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 ... ** ** + COMMON/RELKIN/NONREL ** ** NONREL is a logical used in the IDPAR=1 and 2 options.** ** NONREL=.TRUE. gives the deuteron binding energy of ** ** B=2.224575 MeV using non-relativistic ** ** kinematics. ** ** NONREL=.FALSE. gives the deuteron binding energy of ** ** B=2.224575 MeV using relativistic ** ** kinematics. ** ** Model IDPAR=0 gives the deuteron only using ** ** relativistic kinematics ==> NONREL=.FALSE. ** ** ** ** NOTE: ALL potential models use a fixed fpi**2=0.075 for the ** ** ----- pion-nucleon coupling constant at the pion pole, which ** ** is represented by the DATA FPPPI0/0.075D0/ statement. ** ** ** ************************************************************************ IMPLICIT REAL*8 (A-H,O-Z) REAL*8 VPOT(2,2) INTEGER SPIN CHARACTER PHNAME*3, PHNAM0*3, TYPE*2, CAPS*12, SMAL*12 COMMON/EMANHP/PHNAME COMMON/CHOICE/IDPAR COMMON/JFAC/ JMMM,JMM,JM,J,JP,JPP,JPPP COMMON/QFAC/ QI2,QF2,QIQF,QI2F2,QI2F22,QIPF2 COMMON/POTMOM/VC(-1:1),VS(-1:1),VT(-2:2),VLS(-2:2),VLA(-2:2), . VQ(-3:3),VQP(-2:2) DATA HBARC/197.327053D0/ DATA CAPS/'CSPDFGHIKLMN'/, SMAL/'cspdfghiklmn'/ DATA NQ12/1/, PHNAM0/'***'/, PI/3.14159265358979323846D0/ SAVE NCHAN,SPIN,L,ISO, TJMM,TJM,TJ,TJP,TJPP,TJJ SAVE E00,F00,G00,E11,F11,G11,EMM,EPP,GMM,GPP,FMM,FPP,FPM VPOT(1,1) = 0D0 VPOT(1,2) = 0D0 VPOT(2,1) = 0D0 VPOT(2,2) = 0D0 C* Necessary to avoid underflows: IF(DABS(QI-QF).GT.30*HBARC) RETURN NLOC=1 C* Only IDPAR=2 is fully local model IF(IDPAR.EQ.2) NLOC=0 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.9) WRITE(*,*) . '**** Partial wave exceeds allowable maximum J=9' IF(J.GT.9) STOP L=J IF(PHNAME.EQ.'3P0') L=1 IF(NCHAN.EQ.2) L=J-1 ISO=MOD(SPIN+L+1,2) JMMM=J-3 JMM =J-2 JM =J-1 JP =J+1 JPP =J+2 JPPP=J+3 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)) C** J-dependent coefficients for quadratic spin-orbit E00 = J*JM/TJM/TJ F00 =-2D0*(J*J+JM)/TJM/TJP G00 = JP*JPP/TJ/TJP E11 = JM*JPP/TJM/TJ F11 =-2D0*JM*JPP/TJM/TJP G11 = JM*JPP/TJ/TJP EMM =-JM*JMM/TJM/TJMM EPP =-J*(2D0*J*J+7D0*JP)/TJ/TJ/TJP GMM =-(2D0*J*J-3D0*J+2D0)*JP/TJM/TJ/TJ GPP =-JPP*JPPP/TJP/TJPP FMM = 2D0*(2D0*J*J*J-3D0*J*J-2D0*JM)/TJ/TJ/TJMM FPM = 2D0*TJJ/TJ/TJ FPP = 2D0*(2D0*J*J*J+9D0*J*J+10D0*J+1D0)/TJ/TJ/TJPP ENDIF QI2=QI*QI QF2=QF*QF QIQF=QI*QF QI2F2=QI2+QF2 QI2F22=QI2F2*QI2F2 QIPF2=(QI+QF)*(QI+QF) S2PSI=2D0*QIQF/QI2F2 SPSI2=QF2/QI2F2 CPSI2=QI2/QI2F2 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 CALL VMOM(TYPE,NLOC,NQ12,ISO) IF(NCHAN.EQ.1) THEN IF(SPIN.EQ.0) THEN VPOT(1,1) = VC(0) - 3D0*VS(0) + . QI2*QF2*(E00*VQ(-2)+F00*VQ(0)+G00*VQ(2)) ELSEIF(L.EQ.J) THEN VPOT(1,1) = VC(0) + VS(0) - QIQF*(VLS(-1)-VLS(1))/TJ + . 2D0/3D0*QI2F2* . (VT(0)-0.5D0*S2PSI*(TJP*VT(-1)+TJM*VT(1))/TJ) + . QI2*QF2*(E11*VQ(-2)+F11*VQ(0)+G11*VQ(2)) ELSEIF(PHNAME.EQ.'3P0') THEN VPOT(1,1) = VC(1) + VS(1) - QIQF*JPP/TJP*(VLS(0)-VLS(2)) + . 2D0/3D0*QI2F2*JPP/TJ* . (-VT(1)+0.5D0*S2PSI*(TJPP*VT(0)+TJ*VT(2))/TJP) + . QI2*QF2*(EPP*VQ(-1) + FPP*VQ(1) + GPP*VQ(3)) ENDIF ELSE VPOT(1,1) = VC(-1) + VS(-1) + QIQF*JM/TJM*(VLS(-2)-VLS(0)) + . 2D0/3D0*QI2F2*JM/TJ* . (-VT(-1)+0.5D0*S2PSI*(TJMM*VT(0)+TJ*VT(-2))/TJM) + . QI2*QF2*(EMM*VQ(-3) + FMM*VQ(-1) + GMM*VQ(1)) VPOT(1,2) =-2D0*QI2F2*TJJ/TJ* . (-S2PSI*VT(0)+CPSI2*VT(-1)+SPSI2*VT(1)) - . QI2*QF2*FPM*(VQ(1)-VQ(-1)) VPOT(2,1) =-2D0*QI2F2*TJJ/TJ* . (-S2PSI*VT(0)+SPSI2*VT(-1)+CPSI2*VT(1)) - . QI2*QF2*FPM*(VQ(1)-VQ(-1)) VPOT(2,2) = VC(1) + VS(1) - QIQF*JPP/TJP*(VLS(0)-VLS(2)) + . 2D0/3D0*QI2F2*JPP/TJ* . (-VT(1)+0.5D0*S2PSI*(TJPP*VT(0)+TJ*VT(2))/TJP) + . QI2*QF2*(EPP*VQ(-1) + FPP*VQ(1) + GPP*VQ(3)) ENDIF IF(NQ12.EQ.1) THEN C** Extra contribution from inverse Fourier transform of Q12 C** so momentum space exactly equivalent to configuration space IF(NCHAN.EQ.1) THEN IF(SPIN.EQ.0) THEN VPOT(1,1) = VPOT(1,1) - 2D0*QIQF/TJ*(J*VQP(-1)+JP*VQP(1)) ELSEIF(L.EQ.J) THEN VPOT(1,1) = VPOT(1,1) + QIQF*(-VQP(-1)+VQP(1))/TJ ELSEIF(PHNAME.EQ.'3P0') THEN VPOT(1,1) = VPOT(1,1) + QIQF* . ((2D0*J*J+5D0*J+4D0)/TJ*VQP(0)+JPP*VQP(2))/TJP ENDIF ELSE VPOT(1,1) = VPOT(1,1) + QIQF* . (JM*VQP(-2)+(2D0*J*J-JM)/TJ*VQP(0))/TJM VPOT(1,2) = VPOT(1,2) + 2D0*QIQF*TJJ/TJ*VQP(0) VPOT(2,1) = VPOT(2,1) + 2D0*QIQF*TJJ/TJ*VQP(0) VPOT(2,2) = VPOT(2,2) + QIQF* . ((2D0*J*J+5D0*J+4D0)/TJ*VQP(0)+JPP*VQP(2))/TJP ENDIF ENDIF RETURN END ************************************************************************ SUBROUTINE VMOM(TYPE,NLOC,NQ12,ISO) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 NEUTM CHARACTER PHNAME*3,PHNAM0*3, TYPE*2,TYP0*2 LOGICAL FIRST COMMON/EMANHP/PHNAME COMMON/CHOICE/IDPAR COMMON/POTMOM/VC(-1:1),VS(-1:1),VT(-2:2),VLS(-2:2),VLA(-2:2), . VQ(-3:3),VQP(-2:2) COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0, . AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC COMMON/PARAMS/PAR(6,5) COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2, . ALVE,THV ,GV1, GRO,GOM,GFI, GRO2,GOM2,GFI2, . ALVM, FV1, FRO,FOM,FFI, FRO2,FOM2,FFI2, . ALGS,THGS,GS1, GA0,GEP,GF0, GA02,GEP2,GF02, . GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI, . FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2 COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2, . AVSC1,AVSC2, AEPS,AME1,BEPS,AME2 COMMON/BROADM/GAMRO,THRRO,GAMRC,THRRC,GAMEP,THRE0,THREC,NRO,NEP COMMON/SCALIN/AMT,AMPV COMMON/CUTOFF/ALAM,ALAMP,ALAMV,ALAMS COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4, . GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4 COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2, . AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2, . AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2 DATA FPPPI0/0.075D0/, FIRST/.FALSE./ DATA TYP0/'XX'/, PHNAM0/'***'/ DATA PI/3.14159265358979323846D0/ SAVE ISIGN,CONV,PROTM,NEUTM,HBARC IF(FIRST) GOTO 10 FIRST = .TRUE. CONV = PI/180D0 ** Nucleon and meson masses (Particle Data Group 1990) HBARC = 197.327053D0 PROTM = 938.27231D0 NEUTM = 939.56563D0 AMT = 938.27231D0 AMPV = 139.5675D0 AMPIC = 139.5675D0 AMPI = 134.9739D0 AMETA = 548.8D0 AMETP = 957.5D0 AMROC = 768.3D0 AMRO = 768.7D0 AMOM = 781.95D0 AMFI =1019.412D0 AMA0C = 983.3D0 AMA0 = 983.3D0 AMF0 = 975.6D0 AVSC = ((NEUTM-PROTM)/AMROC)**2 FAC = 2D0*DSQRT(PROTM*NEUTM)/(PROTM+NEUTM) AWPIC = FAC*DSQRT(AMPIC**2 - (NEUTM-PROTM)**2) AWROC = FAC*DSQRT(AMROC**2 - (NEUTM-PROTM)**2) AWA0C = FAC*DSQRT(AMA0C**2 - (NEUTM-PROTM)**2) ** Broad rho-meson: spectral density to two effective Yukawa's ** Yukawa's fitted to PHI0C 0.001 - 2 fm (steps 0.002, ALAMV=825 MeV) GAMRO = 152.4D0 THRRO = AMPIC+AMPIC GAMRC = 149.1D0 THRRC = AMPIC+AMPI NRO = 1 ARO = 0.2655205D0 BRO = 0.5607493D0 AMR1 = 645.3772D0 AMR2 = 878.3667D0 AROC = 0.3875515D0 BROC = 0.4508341D0 AMRC1= 674.1521D0 AMRC2= 929.9742D0 AWRC1 = FAC*DSQRT(AMRC1**2 - (NEUTM-PROTM)**2) AWRC2 = FAC*DSQRT(AMRC2**2 - (NEUTM-PROTM)**2) AVSC1 = ((NEUTM-PROTM)/AMRC1)**2 AVSC2 = ((NEUTM-PROTM)/AMRC2)**2 10 CONTINUE IF(TYPE.EQ.TYP0) THEN IF((IDPAR.EQ.1 .OR. IDPAR.EQ.2) .AND. PHNAME.NE.PHNAM0) GOTO 15 GOTO 20 ENDIF TYP0=TYPE IF(TYPE.EQ.'PP') THEN AMY = PROTM AMN = PROTM I3Y = 1 I3N = 1 ELSEIF(TYPE.EQ.'NP') THEN AMY = NEUTM AMN = PROTM I3Y =-1 I3N = 1 ELSEIF(TYPE.EQ.'PN') THEN AMY = PROTM AMN = NEUTM I3Y = 1 I3N =-1 ELSEIF(TYPE.EQ.'NN') THEN AMY = NEUTM AMN = NEUTM I3Y =-1 I3N =-1 ENDIF ISIGN = I3Y*I3N CALL AMFACS 15 PHNAM0=PHNAME CALL NYMPAR(ISIGN) C pseudovec vector tensor scalar pomeron C FPI GRO FRO GA0 GPOM C PV1 GV1 FOM GEP GA2D C ALPV ALVE FFI GF0 AMPOM C THPV THV RHOC% THGS AMA2D C Cut-off ALAMP ALAMV ALAM ALAMS AMEP C 2 Yukawa AEPS AME1 BEPS AME2 GAMEP ** (joined) cutoffs for pseudoscalar, vector, scalar ALAMP = PAR(5,1) ALAMV = PAR(5,2) ALAMS = PAR(5,4) ALAM = PAR(5,3) IF(ALAMP.EQ.0D0) ALAMP=ALAM IF(ALAMV.EQ.0D0) ALAMV=ALAM IF(ALAMS.EQ.0D0) ALAMS=ALAM ** pseudovector couplings C FPI = PAR(1,1) FPI = DSQRT(FPPPI0*FDEXP(-(AMPI/ALAMP)**2)) PV1 = PAR(2,1) ALPV = PAR(3,1) THPV = PAR(4,1)*CONV ** vector couplings GRO = PAR(1,2) GV1 = PAR(2,2) ALVE = PAR(3,2) THV = PAR(4,2)*CONV ** tensor couplings FRO = PAR(1,3) FOM = PAR(2,3) FFI = PAR(3,3) ** scalar couplings GA0 = PAR(1,4) GEP = PAR(2,4) GF0 = PAR(3,4) THGS = PAR(4,4)*CONV ** diffractive contribution GPOM = PAR(1,5) GA2D = PAR(2,5) AMPOM = PAR(3,5) AMA2D = PAR(4,5) IF(AMA2D.EQ.0D0) AMA2D=AMPOM AMPOM2= AMPOM*AMPOM AMPOM4= AMPOM2*AMPOM2 AMA2D2= AMA2D*AMA2D AMA2D4= AMA2D2*AMA2D2 CALL NYMCOP(ISIGN) ** Broad epsilon-meson: spectral density to two effective Yukawa's ** Yukawa's fitted to PHI0C from 0.001 to 2 fm (steps 0.005) AMEP = PAR(5,5) GAMEP = PAR(6,5) THRE0 = AMPI+AMPI THREC = AMPIC+AMPIC NEP = 0 AEPS = PAR(6,1) AME1 = PAR(6,2) BEPS = PAR(6,3) AME2 = PAR(6,4) ** isospin dependence 20 FACISO= 4D0*ISO-2D0 GA0C2 = FACISO*GA0*GA0 FPIC2 = FACISO*FPI*FPI/(AMPV*AMPV) GROC2 = FACISO*GRO*GRO FROC2 = FACISO*FRO*FRO/(4D0*AMT*AMT) GFPROC= FACISO*(GRO*FRO+FRO*GRO)/(2D0*AMT) GFMROC= FACISO*(GRO*FRO-FRO*GRO)/(2D0*AMT) IF(PHNAME.EQ.'1S0' .AND. (TYPE.EQ.'NP'.OR.TYPE.EQ.'PN')) THEN ** Phenomenological phase-shift difference in PP and NP 1S0 GROC2=GROC2*(1D0+PAR(4,3)) FROC2=FROC2*(1D0+PAR(4,3)) ENDIF FPOM2 = (GPOM/AMT)**2 FA2D2 = (GA2D/AMT)**2 DO 50 KK=-1,1 VC(KK) = 0D0 50 VS(KK) = 0D0 DO 51 KK=-2,2 VT(KK) = 0D0 VLS(KK) = 0D0 VLA(KK) = 0D0 VQ(KK) = 0D0 51 VQP(KK) = 0D0 VQ(-3) = 0D0 VQ( 3) = 0D0 CALL PSVECT(TYPE,NLOC) CALL VECTOR(TYPE,NLOC,NQ12) CALL SCALAR(TYPE,NLOC,NQ12) CALL DIFRAC(TYPE,NLOC,NQ12,ISIGN,FACISO) RETURN END ************************************************************************ SUBROUTINE AMFACS IMPLICIT REAL*8 (A-H,O-Z) COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2, . AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2, . AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2 AMY2 = AMY*AMY AMYI = 1D0/AMY AMYI2 = AMYI*AMYI AMN2 = AMN*AMN AMNI = 1D0/AMN AMNI2 = AMNI*AMNI AMYPN = AMY+AMN AMYMN = AMY-AMN AMYPN2= AMYPN*AMYPN AMYMN2= AMYMN*AMYMN AY2PN2= AMY2+AMN2 AY2MN2= AMY2-AMN2 AMYN = AMY*AMN AMYNI = AMYI*AMNI AMYNI2= AMYNI*AMYNI AYPNI2= AMYI2+AMNI2 AYMNI2= AMYI2-AMNI2 REDM = AMYN/AMYPN ALF = 4D0*REDM/AMYPN IF(AMY.EQ.AMN) ALF = 1D0 RETURN END ************************************************************************ SUBROUTINE NYMPAR(ISIGN) IMPLICIT REAL*8 (A-H,O-Z) COMMON/PARAMS/PAR(6,5) COMMON/CHOICE/IDPAR DIMENSION PARTR(5,6) C pseudovec vector tensor scalar pomeron C FPI GRO FRO GA0 GPOM C PV1 GV1 FOM GEP GA2D C ALPV ALVE FFI GF0 AMPOM C THPV THV RHOC% THGS AMA2D C Cut-off ALAMP ALAMV ALAM ALAMS AMEP C 2 Yukawa AEPS AME1 BEPS AME2 GAMEP C---------------------------------------------------------------------- ** Parameter-set with separate cutoffs, fitted to all data Nijm93 DATA PARTR/ 1 .2720668D+00,.9209319D+00,.3770582D+01,.1384689D+01,.5228672D+01 2,.1595311D+00,.2594356D+01,.5816365D+00,.5310001D+01,.2204600D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.3484505D+01,.2081618D+03 4,-.230000D+02,.3750000D+02,.4371144D-01,.3790000D+02,.0000000D+00 5,.1177107D+04,.9045040D+03,.0000000D+00,.5544013D+03,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ IF(IDPAR.EQ.0) THEN DO 1 I=1,6 DO 1 J=1,5 1 PAR(I,J)=PARTR(J,I) ELSEIF(IDPAR.EQ.1) THEN CALL PHSRDL(ISIGN) ELSEIF(IDPAR.EQ.2) THEN CALL PHSLOC(ISIGN) ENDIF WRITE(*,*) ' NIJMEGEN POTENTIAL PARAMETERS ARE:' DO 4 I=1,6 4 WRITE(*,5) (PAR(I,J),J=1,5) 5 FORMAT(5(1X,D15.7)) RETURN END ************************************************************************ C This subroutine reads the parameters of the 0-350 MeV fitted C Reidlike potential, Jan 1993, chi**2/data=1.03 Nijm I C----------------------------------------------------------------------- SUBROUTINE PHSRDL(ISIGN) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER PHNAME*3 LOGICAL NONREL COMMON/PARAMS/PAR(6,5) COMMON/EMANHP/PHNAME COMMON/RELKIN/NONREL DIMENSION APP1S0(5,6),ANP1S0(5,6),PAR1P1(5,6),PAR1D2(5,6), . PAR1F3(5,6),PAR1G4(5,6),PAR3P1(5,6),PAR3D2(5,6), . PAR3F3(5,6),PAR3G4(5,6),PAR3P0(5,6),PAR3C1(5,6), . PAR3C2(5,6),PAR3C3(5,6),PAR3C4(5,6),PARRST(5,6), . PNR3C1(5,6) DATA APP1S0/ 1 .2702427D+00,.6723103D+00,.3697581D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4846532D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4744299D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA ANP1S0/ 1 .2702427D+00,.6723103D+00,.3423514D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4942984D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1P1/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.3290519D+01,.5473693D-01,.4784616D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1D2/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.0000000D+00,.5473693D-01,.2348432D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.3639503D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1F3/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.4812761D+01,.5473693D-01,.6316315D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1G4/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.4142482D+01,.5473693D-01,.4859055D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3P1/ 1 .2702427D+00,.6723103D+00,.0000000D+00,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4980991D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4923498D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3D2/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.5880449D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.1693534D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3F3/ 1 .2702427D+00,.6723103D+00,.8961465D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.5458122D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3G4/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.3201615D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3P0/ 1 .2702427D+00,.6723103D+00,.3160965D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.3726932D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C1/ 1 .2706937D+00,.2598577D+01,.2125742D+01,.8334762D+00,.2751792D+01 2,.0000000D+00,.2536422D+01,.5473693D-01,.4906761D+01,.44371888D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.8848614D+03,.6465221D+03,.0000000D+00,.6990612D+03,.7600000D+03 6,.1690008D+00,.5831763D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C2/ 1 .2702427D+00,.6723103D+00,.5188648D+01,.8334762D+00,.3668014D+01 2,.1341669D+00,.0000000D+00,.5473693D-01,.3995761D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.8275346D+03,.0000000D+00,.5831699D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C3/ 1 .2702427D+00,.6723103D+00,.3877361D+01,.8334762D+00,.4723766D+01 2,.2811080D+00,.0000000D+00,.5473693D-01,.4818122D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C4/ 1 .2702427D+00,.6723103D+00,.7551377D+01,.8334762D+00,.2360920D+01 2,.2871299D+00,.0000000D+00,.5473693D-01,.4855317D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PARRST/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4859055D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ C Parameters for a nonlocal potential with non-relativistic deuteron DATA PNR3C1/ 1 .2706937D+00,.2588931D+01,.2120162D+01,.8334762D+00,.2751792D+01 2,.0000000D+00,.2521404D+01,.5473693D-01,.4889314D+01,.44371806D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.8848614D+03,.6465221D+03,.0000000D+00,.6990612D+03,.7600000D+03 6,.1690008D+00,.5831292D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DO 1 I=1,6 DO 1 J=1,5 IF(PHNAME.EQ.'1S0') THEN IF(ISIGN.EQ. 1) PAR(I,J)=APP1S0(J,I) IF(ISIGN.EQ.-1) PAR(I,J)=ANP1S0(J,I) ELSEIF(PHNAME.EQ.'1P1') THEN PAR(I,J)=PAR1P1(J,I) ELSEIF(PHNAME.EQ.'1D2') THEN PAR(I,J)=PAR1D2(J,I) ELSEIF(PHNAME.EQ.'1F3') THEN PAR(I,J)=PAR1F3(J,I) ELSEIF(PHNAME.EQ.'1G4') THEN PAR(I,J)=PAR1G4(J,I) ELSEIF(PHNAME.EQ.'3P0') THEN PAR(I,J)=PAR3P0(J,I) ELSEIF(PHNAME.EQ.'3P1') THEN PAR(I,J)=PAR3P1(J,I) ELSEIF(PHNAME.EQ.'3D2') THEN PAR(I,J)=PAR3D2(J,I) ELSEIF(PHNAME.EQ.'3F3') THEN PAR(I,J)=PAR3F3(J,I) ELSEIF(PHNAME.EQ.'3G4') THEN PAR(I,J)=PAR3G4(J,I) ELSEIF(PHNAME.EQ.'3C1') THEN IF(.NOT.NONREL) PAR(I,J)=PAR3C1(J,I) IF(NONREL) PAR(I,J)=PNR3C1(J,I) ELSEIF(PHNAME.EQ.'3C2') THEN PAR(I,J)=PAR3C2(J,I) ELSEIF(PHNAME.EQ.'3C3') THEN PAR(I,J)=PAR3C3(J,I) ELSEIF(PHNAME.EQ.'3C4') THEN PAR(I,J)=PAR3C4(J,I) ELSE PAR(I,J)=PARRST(J,I) ENDIF 1 CONTINUE RETURN END ************************************************************************ C This subroutine reads the parameters of the 0-350 MeV fitted local C potential (without Q**2), Jan 1993, chi**2/data=1.03 Nijm II C----------------------------------------------------------------------- SUBROUTINE PHSLOC(ISIGN) IMPLICIT REAL*8 (A-H,O-Z) CHARACTER PHNAME*3 LOGICAL NONREL COMMON/PARAMS/PAR(6,5) COMMON/EMANHP/PHNAME COMMON/RELKIN/NONREL DIMENSION APP1S0(5,6),ANP1S0(5,6),PAR1P1(5,6),PAR1D2(5,6), . PAR1F3(5,6),PAR1G4(5,6),PAR3P1(5,6),PAR3D2(5,6), . PAR3F3(5,6),PAR3G4(5,6),PAR3P0(5,6),PAR3C1(5,6), . PAR3C2(5,6),PAR3C3(5,6),PAR3C4(5,6),PARRST(5,6), . PNR3C1(5,6) DATA APP1S0/ 1 .2702427D+00,.6723103D+00,.4006116D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4669822D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.3909783D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4692497D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA ANP1S0/ 1 .2702427D+00,.6723103D+00,.7407998D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4626459D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2196159D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4722381D+03,.6130152D+00,.1021139D+04,.6400000D+03/ C DATA PAR1P1/ C 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 C 2,.2871299D+00,.2849628D+01,.5473693D-01,.5411714D+01,.4437220D+00 C 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 C 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 C 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 C 6,.1690008D+00,.5217623D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1P1/ . .2702427D+00,.6723103D+00,.0000000D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.3430827D+01,.5473693D-01,.2825877D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+03,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1D2/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4138573D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2243917D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4367031D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1F3/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.3983000D+01,.5473693D-01,.5627977D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR1G4/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4859055D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2037620D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3P1/ 1 .2702427D+00,.6723103D+00,.0000000D+00,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4171550D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.3368384D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4530824D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3D2/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.5469270D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.1847244D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3F3/ 1 .2702427D+00,.6723103D+00,.6012926D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.5530460D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3G4/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.3663270D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3P0/ 1 .2702427D+00,.6723103D+00,.2761025D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.3041218D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.8275346D+03,.0000000D+00,.1134832D+04,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C1/ 1 .2702427D+00,.1607944D+01,.1841778D+01,.5469244D+00,.3469472D+01 2,.2871299D+00,.2240543D+01,.5473693D-01,.4035077D+01,.4437213D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.5151821D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.8275346D+03,.0000000D+00,.8044237D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C2/ 1 .2702427D+00,.6723103D+00,.5816373D+01,.8334762D+00,.3957678D+01 2,.2353573D+00,.0000000D+00,.5473693D-01,.4143714D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2781205D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.8275346D+03,.0000000D+00,.6121468D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C3/ 1 .2702427D+00,.6723103D+00,.4050335D+01,.8334762D+00,.4316501D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.5048592D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PAR3C4/ 1 .2702427D+00,.6723103D+00,.7347855D+01,.8334762D+00,.2579081D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.5157279D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DATA PARRST/ 1 .2702427D+00,.6723103D+00,.4728635D+01,.8334762D+00,.2751792D+01 2,.2871299D+00,.2849628D+01,.5473693D-01,.4859055D+01,.4437220D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.2579522D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.0000000D+00,.0000000D+00,.8275346D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ C Parameters for a local potential with non-relativistic deuteron DATA PNR3C1/ 1 .2702427D+00,.1710842D+01,.1781765D+01,.5362515D+00,.3461562D+01 2,.2871299D+00,.2247159D+01,.5473693D-01,.4028595D+01,.44387996D+00 3,.3550000D+00,.1000000D+01,.0000000D+00,.8389363D+00,.5184792D+03 4,-.230000D+02,.3750000D+02,.0000000D+00,.3790000D+02,.0000000D+00 5,.8275346D+03,.0000000D+00,.8044237D+03,.0000000D+00,.7600000D+03 6,.1690008D+00,.4878179D+03,.6130152D+00,.1021139D+04,.6400000D+03/ DO 1 I=1,6 DO 1 J=1,5 IF(PHNAME.EQ.'1S0') THEN IF(ISIGN.EQ. 1) PAR(I,J)=APP1S0(J,I) IF(ISIGN.EQ.-1) PAR(I,J)=ANP1S0(J,I) ELSEIF(PHNAME.EQ.'1P1') THEN PAR(I,J)=PAR1P1(J,I) ELSEIF(PHNAME.EQ.'1D2') THEN PAR(I,J)=PAR1D2(J,I) ELSEIF(PHNAME.EQ.'1F3') THEN PAR(I,J)=PAR1F3(J,I) ELSEIF(PHNAME.EQ.'1G4') THEN PAR(I,J)=PAR1G4(J,I) ELSEIF(PHNAME.EQ.'3P0') THEN PAR(I,J)=PAR3P0(J,I) ELSEIF(PHNAME.EQ.'3P1') THEN PAR(I,J)=PAR3P1(J,I) ELSEIF(PHNAME.EQ.'3D2') THEN PAR(I,J)=PAR3D2(J,I) ELSEIF(PHNAME.EQ.'3F3') THEN PAR(I,J)=PAR3F3(J,I) ELSEIF(PHNAME.EQ.'3G4') THEN PAR(I,J)=PAR3G4(J,I) ELSEIF(PHNAME.EQ.'3C1') THEN IF(.NOT.NONREL) PAR(I,J)=PAR3C1(J,I) IF(NONREL) PAR(I,J)=PNR3C1(J,I) ELSEIF(PHNAME.EQ.'3C2') THEN PAR(I,J)=PAR3C2(J,I) ELSEIF(PHNAME.EQ.'3C3') THEN PAR(I,J)=PAR3C3(J,I) ELSEIF(PHNAME.EQ.'3C4') THEN PAR(I,J)=PAR3C4(J,I) ELSE PAR(I,J)=PARRST(J,I) ENDIF 1 CONTINUE RETURN END ************************************************************************ SUBROUTINE NYMCOP(ISIGN) IMPLICIT REAL*8 (A-H,O-Z) COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2, . ALVE,THV ,GV1, GRO,GOM,GFI, GRO2,GOM2,GFI2, . ALVM, FV1, FRO,FOM,FFI, FRO2,FOM2,FFI2, . ALGS,THGS,GS1, GA0,GEP,GF0, GA02,GEP2,GF02, . GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI, . FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2 COMMON/SCALIN/AMT,AMPV DATA SR3/1.7320508075688772935D0/ ** pseudovector coupling constants PV8 = FPI * (4D0*ALPV-1D0)/SR3 COST = DCOS(THPV) SINT = DSIN(THPV) FETA = COST*PV8 - SINT*PV1 FETP = SINT*PV8 + COST*PV1 FPI2 = FPI*FPI/(AMPV*AMPV) * ISIGN FETA2= FETA*FETA/(AMPV*AMPV) FETP2= FETP*FETP/(AMPV*AMPV) ** vector coupling constants GV8 = GRO * (4D0*ALVE-1D0)/SR3 COST = DCOS(THV) SINT = DSIN(THV) GFI = COST*GV8 - SINT*GV1 GOM = SINT*GV8 + COST*GV1 GRO2 = GRO*GRO * ISIGN GOM2 = GOM*GOM GFI2 = GFI*GFI ** tensor coupling constants COST = DCOS(THV) SINT = DSIN(THV) FV8 = COST*(GFI+FFI) + SINT*(GOM+FOM) - GV8 FV1 =-SINT*(GFI+FFI) + COST*(GOM+FOM) - GV1 ALVM = (SR3*(GV8+FV8)+(GRO+FRO))/(4D0*(GRO+FRO)) FRO2 = FRO*FRO/(4D0*AMT*AMT) * ISIGN FOM2 = FOM*FOM/(4D0*AMT*AMT) FFI2 = FFI*FFI/(4D0*AMT*AMT) GFPRO = (GRO*FRO+FRO*GRO)/(2D0*AMT) * ISIGN GFPOM = (GOM*FOM+FOM*GOM)/(2D0*AMT) GFPFI = (GFI*FFI+FFI*GFI)/(2D0*AMT) GFMRO = (GRO*FRO-FRO*GRO)/(2D0*AMT) * ISIGN GFMOM = (GOM*FOM-FOM*GOM)/(2D0*AMT) GFMFI = (GFI*FFI-FFI*GFI)/(2D0*AMT) ** scalar coupling constants COST = DCOS(THGS) SINT = DSIN(THGS) GS8 = COST*GEP + SINT*GF0 GS1 =-SINT*GEP + COST*GF0 ALGS = (SR3*GS8+GA0)/(4D0*GA0) GA02 = GA0*GA0 * ISIGN GEP2 = GEP*GEP GF02 = GF0*GF0 RETURN END ************************************************************************ SUBROUTINE PSVECT(TYPE,NLOC) C* Partial-wave momentum-space potentials for pseudoscalar exchange C* NLOC=0,1: no nonlocal (q^2+k^2/4) contributions C* 2: nonlocal contributions in spin-spin and tensor C----------------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) CHARACTER TYPE*2 COMMON/POTMOM/VC(-1:1),VS(-1:1),VT(-2:2),VLS(-2:2),VLA(-2:2), . VQ(-3:3),VQP(-2:2) COMMON/JFAC/ JMMM,JMM,JM,J,JP,JPP,JPPP COMMON/QFAC/ QI2,QF2,QIQF,QI2F2,QI2F22,QIPF2 COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0, . AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2, . ALVE,THV ,GV1, GRO,GOM,GFI, GRO2,GOM2,GFI2, . ALVM, FV1, FRO,FOM,FFI, FRO2,FOM2,FFI2, . ALGS,THGS,GS1, GA0,GEP,GF0, GA02,GEP2,GF02, . GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI, . FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2 COMMON/CUTOFF/ALAM,ALAMP,ALAMV,ALAMS COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2, . AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2, . AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2 DIMENSION U(-3:12), R(-3:12), ELN(15), ERN(15) DATA U/16*0D0/, R/16*0D0/ DATA PI/3.14159265358979323846D0/ XCOM=0.5D0*QI2F2/QIQF ALAM2=ALAMP*ALAMP Y=2D0*QIQF/ALAM2 JMAX=J+2 KOM=1 FAC=2D0*PI/QIQF C** Pseudoscalar mesons: pi, eta, eta' DO 1000 IN=1,3 IF(IN.EQ.1) THEN AMES=AMPI FP2 =FPI2 ELSEIF(IN.EQ.2) THEN AMES=AMETA FP2 =FETA2 ELSEIF(IN.EQ.3) THEN AMES=AMETP FP2 =FETP2 ENDIF AMES2=AMES*AMES X=XCOM+0.5D0*AMES2/QIQF CALL SEX(X,Y,AMES/ALAMP,JMAX,KOM,ELN,ERN) XPS =-FAC*FP2*QI2F2/3D0 YPS = FAC*FP2*2D0*QIQF/3D0 XPT =-FAC*FP2 XNLS= FAC*FP2*AMYNI*QI2F22/12D0 YNLS=-FAC*FP2*AMYNI*QI2F2*QIQF/6D0 XNLT= FAC*FP2*AMYNI*QI2F2/4D0 DO 5 K=0,JMAX U(K)=ELN(K+1) R(K)=ERN(K+1) 5 CONTINUE VS(-1) = VS(-1) + (XPS+X*YPS)*U(JM) - YPS*R(JM) VS( 0) = VS( 0) + (XPS+X*YPS)*U(J) - YPS*R(J) VS( 1) = VS( 1) + (XPS+X*YPS)*U(JP) - YPS*R(JP) VT(-2) = VT(-2) + XPT*U(JMM) VT(-1) = VT(-1) + XPT*U(JM) VT( 0) = VT( 0) + XPT*U(J) VT( 1) = VT( 1) + XPT*U(JP) VT( 2) = VT( 2) + XPT*U(JPP) IF(NLOC.EQ.2) THEN VS(-1) = VS(-1) + (XNLS+X*YNLS)*U(JM) - YNLS*R(JM) VS( 0) = VS( 0) + (XNLS+X*YNLS)*U(J) - YNLS*R(J) VS( 1) = VS( 1) + (XNLS+X*YNLS)*U(JP) - YNLS*R(JP) VT(-2) = VT(-2) + XNLT*U(JMM) VT(-1) = VT(-1) + XNLT*U(JM) VT( 0) = VT( 0) + XNLT*U(J) VT( 1) = VT( 1) + XNLT*U(JP) VT( 2) = VT( 2) + XNLT*U(JPP) ENDIF 1000 CONTINUE IF(TYPE.EQ.'PP' .OR. TYPE.EQ.'NN') RETURN C** Charged pseudoscalar: pi+ AMES=AWPIC FP2 =FPIC2 AMES2=AMES*AMES X=XCOM+0.5D0*AMES2/QIQF CALL SEX(X,Y,AMES/ALAMP,JMAX,KOM,ELN,ERN) XPS =-FAC*FP2/ALF*QI2F2/3D0 YPS = FAC*FP2/ALF*2D0*QIQF/3D0 XPT =-FAC*FP2/ALF XNLS= FAC*FP2/ALF*AYPNI2*QI2F22/24D0 YNLS=-FAC*FP2/ALF*AYPNI2*QI2F2*QIQF/12D0 XNLT= FAC*FP2/ALF*AYPNI2*QI2F2/8D0 DO 6 K=0,JMAX U(K)=ELN(K+1) R(K)=ERN(K+1) 6 CONTINUE VS(-1) = VS(-1) + (XPS+X*YPS)*U(JM) - YPS*R(JM) VS( 0) = VS( 0) + (XPS+X*YPS)*U(J) - YPS*R(J) VS( 1) = VS( 1) + (XPS+X*YPS)*U(JP) - YPS*R(JP) VT(-2) = VT(-2) + XPT*U(JMM) VT(-1) = VT(-1) + XPT*U(JM) VT( 0) = VT( 0) + XPT*U(J) VT( 1) = VT( 1) + XPT*U(JP) VT( 2) = VT( 2) + XPT*U(JPP) IF(NLOC.EQ.2) THEN VS(-1) = VS(-1) + (XNLS+X*YNLS)*U(JM) - YNLS*R(JM) VS( 0) = VS( 0) + (XNLS+X*YNLS)*U(J) - YNLS*R(J) VS( 1) = VS( 1) + (XNLS+X*YNLS)*U(JP) - YNLS*R(JP) VT(-2) = VT(-2) + XNLT*U(JMM) VT(-1) = VT(-1) + XNLT*U(JM) VT( 0) = VT( 0) + XNLT*U(J) VT( 1) = VT( 1) + XNLT*U(JP) VT( 2) = VT( 2) + XNLT*U(JPP) ENDIF RETURN END ************************************************************************ SUBROUTINE VECTOR(TYPE,NLOC,NQ12) C* Partial-wave momentum-space potentials for vector exchange C* NQ12=0: inexact Fourier of quadratic spin-orbit Q12 operator C* 1: exact Fourier Q12 from coordinate to momentum C* NLOC=0: no nonlocal (q^2+k^2/4) contributions C* 1: nonlocal contributions in central potential only C* 2: nonlocal contributions also in spin-spin (only for NP) C----------------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) CHARACTER TYPE*2 COMMON/POTMOM/VC(-1:1),VS(-1:1),VT(-2:2),VLS(-2:2),VLA(-2:2), . VQ(-3:3),VQP(-2:2) COMMON/JFAC/ JMMM,JMM,JM,J,JP,JPP,JPPP COMMON/QFAC/ QI2,QF2,QIQF,QI2F2,QI2F22,QIPF2 COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0, . AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2, . ALVE,THV ,GV1, GRO,GOM,GFI, GRO2,GOM2,GFI2, . ALVM, FV1, FRO,FOM,FFI, FRO2,FOM2,FFI2, . ALGS,THGS,GS1, GA0,GEP,GF0, GA02,GEP2,GF02, . GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI, . FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2 COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2, . AVSC1,AVSC2, AEPS,AME1,BEPS,AME2 COMMON/CUTOFF/ALAM,ALAMP,ALAMV,ALAMS COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2, . AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2, . AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2 DIMENSION U(-3:12), R(-3:12), S(-3:12), G(-3:12), ELN(15), ERN(15) DATA U/16*0D0/, R/16*0D0/, S/16*0D0/, G/16*0D0/ DATA PI/3.14159265358979323846D0/ XCOM=0.5D0*QI2F2/QIQF ALAM2=ALAMV*ALAMV Y=2D0*QIQF/ALAM2 JMAX=J+3 KOM=1 FAC=2D0*PI/QIQF C** Vector mesons: broad rho (2 Yukawa's), omega, fi DO 1000 IN=1,4 IF(IN.EQ.1) THEN AMES=AMR1 GV2 =GRO2 *ARO FV2 =FRO2 *ARO GFV =GFPRO*ARO GFM =GFMRO*ARO ELSEIF(IN.EQ.2) THEN AMES=AMR2 GV2 =GRO2 *BRO FV2 =FRO2 *BRO GFV =GFPRO*BRO GFM =GFMRO*BRO ELSEIF(IN.EQ.3) THEN AMES=AMOM GV2 =GOM2 FV2 =FOM2 GFV =GFPOM GFM =GFMOM ELSEIF(IN.EQ.4) THEN AMES=AMFI GV2 =GFI2 FV2 =FFI2 GFV =GFPFI GFM =GFMFI ENDIF AMES2=AMES*AMES X=XCOM+0.5D0*AMES2/QIQF CALL SEX(X,Y,AMES/ALAMV,JMAX+1,KOM,ELN,ERN) XVC = FAC * (GV2*(1D0-AMYNI/2D0/ALF*QI2F2) - GFV/4D0/REDM*QI2F2 . + FV2*AMYNI/4D0*QI2F22) YVC = FAC * (GV2*AMYNI/ALF*QIQF+GFV/2D0/REDM*QIQF . - FV2*AMYNI*QIQF*QI2F2) ZVC = FAC * FV2*AMYNI*QI2*QF2 XVS =-FAC * ((GV2+AMYPN*GFV+4D0*AMYN*FV2)*AMYNI/6D0*QI2F2 . - FV2*AMYNI/12D0*QI2F22) YVS = FAC * ((GV2+AMYPN*GFV+4D0*AMYN*FV2)*AMYNI/3D0*QIQF . - FV2*AMYNI/3D0*QIQF*QI2F2) ZVS = FAC * FV2*AMYNI/3D0*QI2*QF2 XVT = FAC * ((GV2+AMYPN*GFV+4D0*AMYN*FV2)*AMYNI/4D0 . - FV2*AMYNI/8D0*QI2F2) YVT = FAC * FV2*AMYNI/4D0*QIQF XVLS=-FAC * ((GV2*(0.5D0+1D0/ALF)+AMYPN*GFV)*AMYNI . - FV2*(0.5D0+1D0/ALF)*AMYNI*QI2F2) YVLS=-FAC * FV2*(1D0+2D0/ALF)*AMYNI*QIQF XVQ =-FAC * (GV2+4D0*AMYPN*GFV+8D0*AMYPN2*FV2)*AMYNI2/16D0 XVQ2= XVQ * QIQF XVNL= FAC * GV2*(AYPNI2+AMYNI)/4D0*QI2F2 DO 5 K=1,JMAX U(K)=ELN(K+1) R(K)=ERN(K+1) S(K)=(K*ERN(K)+(K+1)*ERN(K+2))/(2*K+1) G(K)=(ELN(K+2)-ELN(K))/(2*K+1) 5 CONTINUE U(0)=ELN(1) R(0)=ERN(1) S(0)=ERN(2) ARG=(QIPF2+AMES2)/ALAM2 G(0)=(ELN(2)-ELN(1))-FDEXP(-QIPF2/ALAM2)*CE1(ARG) VC(-1) = VC(-1) + (XVC+X*YVC+X*X*ZVC)*U(JM) . - (YVC+X*ZVC)*R(JM) - ZVC*S(JM) VC( 0) = VC( 0) + (XVC+X*YVC+X*X*ZVC)*U(J) . - (YVC+X*ZVC)*R(J) - ZVC*S(J) VC( 1) = VC( 1) + (XVC+X*YVC+X*X*ZVC)*U(JP) . - (YVC+X*ZVC)*R(JP) - ZVC*S(JP) VS(-1) = VS(-1) + (XVS+X*YVS+X*X*ZVS)*U(JM) . - (YVS+X*ZVS)*R(JM) - ZVS*S(JM) VS( 0) = VS( 0) + (XVS+X*YVS+X*X*ZVS)*U(J) . - (YVS+X*ZVS)*R(J) - ZVS*S(J) VS( 1) = VS( 1) + (XVS+X*YVS+X*X*ZVS)*U(JP) . - (YVS+X*ZVS)*R(JP) - ZVS*S(JP) VT(-2) = VT(-2) + (XVT+X*YVT)*U(JMM)- YVT*R(JMM) VT(-1) = VT(-1) + (XVT+X*YVT)*U(JM) - YVT*R(JM) VT( 0) = VT( 0) + (XVT+X*YVT)*U(J) - YVT*R(J) VT( 1) = VT( 1) + (XVT+X*YVT)*U(JP) - YVT*R(JP) VT( 2) = VT( 2) + (XVT+X*YVT)*U(JPP)- YVT*R(JPP) VLS(-2)= VLS(-2)+ (XVLS+X*YVLS)*U(JMM)- YVLS*R(JMM) VLS(-1)= VLS(-1)+ (XVLS+X*YVLS)*U(JM) - YVLS*R(JM) VLS( 0)= VLS( 0)+ (XVLS+X*YVLS)*U(J) - YVLS*R(J) VLS( 1)= VLS( 1)+ (XVLS+X*YVLS)*U(JP) - YVLS*R(JP) VLS( 2)= VLS( 2)+ (XVLS+X*YVLS)*U(JPP)- YVLS*R(JPP) VQ(-3) = VQ(-3) + XVQ*U(JMMM) VQ(-2) = VQ(-2) + XVQ*U(JMM) VQ(-1) = VQ(-1) + XVQ*U(JM) VQ( 0) = VQ( 0) + XVQ*U(J) VQ( 1) = VQ( 1) + XVQ*U(JP) VQ( 2) = VQ( 2) + XVQ*U(JPP) VQ( 3) = VQ( 3) + XVQ*U(JPPP) IF(NLOC.NE.0) THEN VC(-1) = VC(-1) + XVNL*U(JM) VC( 0) = VC( 0) + XVNL*U(J) VC( 1) = VC( 1) + XVNL*U(JP) ENDIF IF(NQ12.EQ.1) THEN C** Extra contribution from inverse Fourier transform of Q12 VQP(-2) = VQP(-2) + XVQ2*G(JMM) VQP(-1) = VQP(-1) + XVQ2*G(JM) VQP( 0) = VQP( 0) + XVQ2*G(J) VQP( 1) = VQP( 1) + XVQ2*G(JP) VQP( 2) = VQP( 2) + XVQ2*G(JPP) ENDIF 1000 CONTINUE IF(TYPE.EQ.'PP' .OR. TYPE.EQ.'NN') RETURN C** Charged vector: broad rho+ (2 Yukawa's) DO 2000 IN=1,2 IF(IN.EQ.1) THEN AMES=AWRC1 GV2 =GROC2 *AROC FV2 =FROC2 *AROC GFV =GFPROC*AROC GFM =GFMROC*AROC FFAC=FAC*AVSC1 ELSEIF(IN.EQ.2) THEN AMES=AWRC2 GV2 =GROC2 *BROC FV2 =FROC2 *BROC GFV =GFPROC*BROC GFM =GFMROC*BROC FFAC=FAC*AVSC2 ENDIF AMES2=AMES*AMES X=XCOM+0.5D0*AMES2/QIQF CALL SEX(X,Y,AMES/ALAMV,JMAX+1,KOM,ELN,ERN) XVC = FAC * (GV2*(ALF-AMYNI/2D0*QI2F2) - ALF* . (GFV*AMYPN+FV2*AMYMN2)*AYPNI2/8D0*QI2F2 + FV2*AMYNI/4D0* . (1D0+(AY2MN2**2+2D0*AMYN*AMYMN2)*AMYNI2/16D0)*QI2F22) YVC = FAC * (GV2*AMYNI*QIQF + ALF*(GFV*AMYPN+FV2*AMYMN2)* . AYPNI2/4D0*QIQF - FV2*AMYNI* . (1D0+(AY2MN2**2+2D0*AMYN*AMYMN2)*AMYNI2/16D0)*QI2F2*QIQF) ZVC = FAC * FV2*AMYNI* . (1D0+(AY2MN2**2+2D0*AMYN*AMYMN2)*AMYNI2/16D0)*QI2*QF2 XVS = FAC * (-(GV2+AMYPN*GFV+AMYPN2*FV2)* . (AMYNI/6D0+ALF*AMYMN2*AMYNI2/16D0)*QI2F2 + . FV2/ALF*AYPNI2/24D0*QI2F22) YVS = FAC * ((GV2+AMYPN*GFV+AMYPN2*FV2)* . (AMYNI/3D0+ALF*AMYMN2*AMYNI2/8D0)*QIQF - . FV2/ALF*AYPNI2/6D0*QIQF*QI2F2) ZVS = FAC * FV2/ALF*AYPNI2/6D0*QI2*QF2 XVT = FAC * ((GV2+AMYPN*GFV+AMYPN2*FV2)*AMYNI/4D0 - . FV2/ALF*AYPNI2/16D0*QI2F2) YVT = FAC* FV2/ALF*AYPNI2/8D0*QIQF XVLS= FAC * (-(GV2*(2D0-ALF/2D0)+2D0*AY2PN2/AMYPN*GFV + . 2D0*AMYMN2*FV2)*AMYNI + . FV2*((2D0/ALF-0.5D0)/ALF-AMYMN2*AMYNI/8D0)*AMYNI*QI2F2) YVLS=-FAC* FV2*((4D0/ALF-1D0)/ALF-AMYMN2*AMYNI/4D0)*AMYNI*QIQF XVQ =-FAC * ALF*(GV2+4D0*AMYPN*GFV+8D0*AMYPN2*FV2)*AMYNI2/16D0 XVNLC= FAC * (ALF*GV2*(AYPNI2+AMYNI)/2D0+GFV*AMYMN2/AMYPN*AMYNI+ . FV2*(4D0/ALF-2D0*ALF*AY2PN2*AMYNI))*QI2F2/2D0 XVNLS= FAC * ALF*AMYMN2*AMYNI2/8D0*(GV2+AMYPN*GFV+AMYPN2*FV2)* . QI2F2 C** Second part of vector-vector potential (same as scalar) XVC = XVC + FFAC * ALF*GV2*(1D0+AMYNI/4D0*QI2F2) YVC = YVC - FFAC * ALF*GV2*AMYNI/2D0*QIQF XVLS = XVLS + FFAC * ALF*GV2*AMYNI/2D0 XVQ = XVQ - FFAC * ALF*GV2*AMYNI2/16D0 XVNLC= XVNLC+ FFAC * ALF*GV2*(AYPNI2/4D0-AMYNI)*QI2F2/2D0 XVQ2= XVQ * QIQF DO 6 K=1,JMAX U(K)=ELN(K+1) R(K)=ERN(K+1) S(K)=(K*ERN(K)+(K+1)*ERN(K+2))/(2*K+1) G(K)=(ELN(K+2)-ELN(K))/(2*K+1) 6 CONTINUE U(0)=ELN(1) R(0)=ERN(1) S(0)=ERN(2) ARG=(QIPF2+AMES2)/ALAM2 G(0)=(ELN(2)-ELN(1))-FDEXP(-QIPF2/ALAM2)*CE1(ARG) VC(-1) = VC(-1) + (XVC+X*YVC+X*X*ZVC)*U(JM) . - (YVC+X*ZVC)*R(JM) - ZVC*S(JM) VC( 0) = VC( 0) + (XVC+X*YVC+X*X*ZVC)*U(J) . - (YVC+X*ZVC)*R(J) - ZVC*S(J) VC( 1) = VC( 1) + (XVC+X*YVC+X*X*ZVC)*U(JP) . - (YVC+X*ZVC)*R(JP) - ZVC*S(JP) VS(-1) = VS(-1) + (XVS+X*YVS+X*X*ZVS)*U(JM) . - (YVS+X*ZVS)*R(JM) - ZVS*S(JM) VS( 0) = VS( 0) + (XVS+X*YVS+X*X*ZVS)*U(J) . - (YVS+X*ZVS)*R(J) - ZVS*S(J) VS( 1) = VS( 1) + (XVS+X*YVS+X*X*ZVS)*U(JP) . - (YVS+X*ZVS)*R(JP) - ZVS*S(JP) VT(-2) = VT(-2) + (XVT+X*YVT)*U(JMM)- YVT*R(JMM) VT(-1) = VT(-1) + (XVT+X*YVT)*U(JM) - YVT*R(JM) VT( 0) = VT( 0) + (XVT+X*YVT)*U(J) - YVT*R(J) VT( 1) = VT( 1) + (XVT+X*YVT)*U(JP) - YVT*R(JP) VT( 2) = VT( 2) + (XVT+X*YVT)*U(JPP)- YVT*R(JPP) VLS(-2)= VLS(-2)+ (XVLS+X*YVLS)*U(JMM)- YVLS*R(JMM) VLS(-1)= VLS(-1)+ (XVLS+X*YVLS)*U(JM) - YVLS*R(JM) VLS( 0)= VLS( 0)+ (XVLS+X*YVLS)*U(J) - YVLS*R(J) VLS( 1)= VLS( 1)+ (XVLS+X*YVLS)*U(JP) - YVLS*R(JP) VLS( 2)= VLS( 2)+ (XVLS+X*YVLS)*U(JPP)- YVLS*R(JPP) VQ(-3) = VQ(-3) + XVQ*U(JMMM) VQ(-2) = VQ(-2) + XVQ*U(JMM) VQ(-1) = VQ(-1) + XVQ*U(JM) VQ( 0) = VQ( 0) + XVQ*U(J) VQ( 1) = VQ( 1) + XVQ*U(JP) VQ( 2) = VQ( 2) + XVQ*U(JPP) VQ( 3) = VQ( 3) + XVQ*U(JPPP) IF(NLOC.NE.0) THEN VC(-1) = VC(-1) + XVNLC*U(JM) VC( 0) = VC( 0) + XVNLC*U(J) VC( 1) = VC( 1) + XVNLC*U(JP) ENDIF IF(NLOC.EQ.2) THEN VS(-1) = VS(-1) + XVNLS*U(JM) VS( 0) = VS( 0) + XVNLS*U(J) VS( 1) = VS( 1) + XVNLS*U(JP) ENDIF IF(NQ12.EQ.1) THEN C** Extra contribution from inverse Fourier transform of Q12 VQP(-2) = VQP(-2) + XVQ2*G(JMM) VQP(-1) = VQP(-1) + XVQ2*G(JM) VQP( 0) = VQP( 0) + XVQ2*G(J) VQP( 1) = VQP( 1) + XVQ2*G(JP) VQP( 2) = VQP( 2) + XVQ2*G(JPP) ENDIF 2000 CONTINUE RETURN END ************************************************************************ SUBROUTINE SCALAR(TYPE,NLOC,NQ12) C* Partial-wave momentum-space potentials for scalar exchange C* NQ12=0: inexact Fourier of quadratic spin-orbit Q12 operator C* 1: exact Fourier Q12 from coordinate to momentum C* NLOC=0: no nonlocal (q^2+k^2/4) contributions C* 1: nonlocal contributions in central potential C----------------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) CHARACTER TYPE*2 COMMON/POTMOM/VC(-1:1),VS(-1:1),VT(-2:2),VLS(-2:2),VLA(-2:2), . VQ(-3:3),VQP(-2:2) COMMON/JFAC/ JMMM,JMM,JM,J,JP,JPP,JPPP COMMON/QFAC/ QI2,QF2,QIQF,QI2F2,QI2F22,QIPF2 COMMON/MESONM/AMPI,AMETA,AMETP, AMRO,AMOM,AMFI, AMA0,AMEP,AMF0, . AMPIC,AMROC,AMA0C,AWPIC,AWROC,AWA0C, AVSC COMMON/COPLNG/ALPV,THPV,PV1, FPI,FETA,FETP,FPI2,FETA2,FETP2, . ALVE,THV ,GV1, GRO,GOM,GFI, GRO2,GOM2,GFI2, . ALVM, FV1, FRO,FOM,FFI, FRO2,FOM2,FFI2, . ALGS,THGS,GS1, GA0,GEP,GF0, GA02,GEP2,GF02, . GFPRO,GFPOM,GFPFI, GFMRO,GFMOM,GFMFI, . FPIC2,GROC2,FROC2,GFPROC,GFMROC,GA0C2 COMMON/YUKEFF/ARO,AMR1,AROC,AMRC1,AWRC1,BRO,AMR2,BROC,AMRC2,AWRC2, . AVSC1,AVSC2, AEPS,AME1,BEPS,AME2 COMMON/CUTOFF/ALAM,ALAMP,ALAMV,ALAMS COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2, . AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2, . AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2 DIMENSION U(-3:12), R(-3:12), G(-3:12), ELN(15), ERN(15) DATA U/16*0D0/, R/16*0D0/, G/16*0D0/ DATA PI/3.14159265358979323846D0/ XCOM=0.5D0*QI2F2/QIQF ALAM2=ALAMS*ALAMS Y=2D0*QIQF/ALAM2 JMAX=J+3 KOM=1 FAC=2D0*PI/QIQF C** Scalar mesons: a0, broad epsilon (2 Yukawa's), f0 DO 1000 IN=1,4 IF(IN.EQ.1) THEN AMES=AMA0 GS2 =GA02 ELSEIF(IN.EQ.2) THEN AMES=AME1 GS2 =GEP2*AEPS ELSEIF(IN.EQ.3) THEN AMES=AME2 GS2 =GEP2*BEPS ELSEIF(IN.EQ.4) THEN AMES=AMF0 GS2 =GF02 ENDIF AMES2=AMES*AMES X=XCOM+0.5D0*AMES2/QIQF CALL SEX(X,Y,AMES/ALAMS,JMAX+1,KOM,ELN,ERN) XSC =-FAC * GS2*(1D0+AYPNI2/8D0*QI2F2) YSC = FAC * GS2*AYPNI2/4D0*QIQF XSLS=-FAC * GS2*AYPNI2/4D0 XSQ = FAC * GS2*AMYNI2/16D0 XSNL= FAC * GS2*AMYNI/4D0*QI2F2 XSQ2= XSQ * QIQF DO 5 K=1,JMAX U(K)=ELN(K+1) R(K)=ERN(K+1) G(K)=(ELN(K+2)-ELN(K))/(2*K+1) 5 CONTINUE U(0)=ELN(1) R(0)=ERN(1) ARG=(QIPF2+AMES2)/ALAM2 G(0)=(ELN(2)-ELN(1))-FDEXP(-QIPF2/ALAM2)*CE1(ARG) VC(-1) = VC(-1) + (XSC+X*YSC)*U(JM) - YSC*R(JM) VC( 0) = VC( 0) + (XSC+X*YSC)*U(J) - YSC*R(J) VC( 1) = VC( 1) + (XSC+X*YSC)*U(JP) - YSC*R(JP) VLS(-2)= VLS(-2)+ XSLS*U(JMM) VLS(-1)= VLS(-1)+ XSLS*U(JM) VLS( 0)= VLS( 0)+ XSLS*U(J) VLS( 1)= VLS( 1)+ XSLS*U(JP) VLS( 2)= VLS( 2)+ XSLS*U(JPP) VQ(-3) = VQ(-3) + XSQ*U(JMMM) VQ(-2) = VQ(-2) + XSQ*U(JMM) VQ(-1) = VQ(-1) + XSQ*U(JM) VQ( 0) = VQ( 0) + XSQ*U(J) VQ( 1) = VQ( 1) + XSQ*U(JP) VQ( 2) = VQ( 2) + XSQ*U(JPP) VQ( 3) = VQ( 3) + XSQ*U(JPPP) IF(NLOC.NE.0) THEN VC(-1) = VC(-1) + XSNL*U(JM) VC( 0) = VC( 0) + XSNL*U(J) VC( 1) = VC( 1) + XSNL*U(JP) ENDIF IF(NQ12.EQ.1) THEN C** Extra contribution from inverse Fourier transform of Q12 VQP(-2) = VQP(-2) + XSQ2*G(JMM) VQP(-1) = VQP(-1) + XSQ2*G(JM) VQP( 0) = VQP( 0) + XSQ2*G(J) VQP( 1) = VQP( 1) + XSQ2*G(JP) VQP( 2) = VQP( 2) + XSQ2*G(JPP) ENDIF 1000 CONTINUE IF(TYPE.EQ.'PP' .OR. TYPE.EQ.'NN') RETURN C** Charged scalar: a0+ AMES=AWA0C GS2 =GA0C2 AMES2=AMES*AMES X=XCOM+0.5D0*AMES2/QIQF CALL SEX(X,Y,AMES/ALAMS,JMAX+1,KOM,ELN,ERN) XSC =-FAC * ALF*GS2*(1D0+AMYNI/4D0*QI2F2) YSC = FAC * ALF*GS2*AMYNI/2D0*QIQF XSLS=-FAC * ALF*GS2*AMYNI/2D0 XSQ = FAC * ALF*GS2*AMYNI2/16D0 XSNL=-FAC * ALF*GS2*(AYPNI2/4D0-AMYNI)*QI2F2/2D0 XSQ2= XSQ * QIQF DO 6 K=1,JMAX U(K)=ELN(K+1) R(K)=ERN(K+1) G(K)=(ELN(K+2)-ELN(K))/(2*K+1) 6 CONTINUE U(0)=ELN(1) R(0)=ERN(1) ARG=(QIPF2+AMES2)/ALAM2 G(0)=(ELN(2)-ELN(1))-FDEXP(-QIPF2/ALAM2)*CE1(ARG) VC(-1) = VC(-1) + (XSC+X*YSC)*U(JM) - YSC*R(JM) VC( 0) = VC( 0) + (XSC+X*YSC)*U(J) - YSC*R(J) VC( 1) = VC( 1) + (XSC+X*YSC)*U(JP) - YSC*R(JP) VLS(-2)= VLS(-2)+ XSLS*U(JMM) VLS(-1)= VLS(-1)+ XSLS*U(JM) VLS( 0)= VLS( 0)+ XSLS*U(J) VLS( 1)= VLS( 1)+ XSLS*U(JP) VLS( 2)= VLS( 2)+ XSLS*U(JPP) VQ(-3) = VQ(-3) + XSQ*U(JMMM) VQ(-2) = VQ(-2) + XSQ*U(JMM) VQ(-1) = VQ(-1) + XSQ*U(JM) VQ( 0) = VQ( 0) + XSQ*U(J) VQ( 1) = VQ( 1) + XSQ*U(JP) VQ( 2) = VQ( 2) + XSQ*U(JPP) VQ( 3) = VQ( 3) + XSQ*U(JPPP) IF(NLOC.NE.0) THEN VC(-1) = VC(-1) + XSNL*U(JM) VC( 0) = VC( 0) + XSNL*U(J) VC( 1) = VC( 1) + XSNL*U(JP) ENDIF IF(NQ12.EQ.1) THEN C** Extra contribution from inverse Fourier transform of Q12 VQP(-2) = VQP(-2) + XSQ2*G(JMM) VQP(-1) = VQP(-1) + XSQ2*G(JM) VQP( 0) = VQP( 0) + XSQ2*G(J) VQP( 1) = VQP( 1) + XSQ2*G(JP) VQP( 2) = VQP( 2) + XSQ2*G(JPP) ENDIF RETURN END ************************************************************************ SUBROUTINE DIFRAC(TYPE,NLOC,NQ12,ISIGN,FACISO) C* Partial-wave momentum-space potentials for diffractive part C* NQ12=0: inexact Fourier of quadratic spin-orbit Q12 operator C* 1: exact Fourier Q12 from coordinate to momentum C* NLOC=0: no nonlocal (q^2+k^2/4) contributions C* 1: nonlocal contributions in central potential C----------------------------------------------------------------------- IMPLICIT REAL*8(A-H,O-Z) CHARACTER TYPE*2 COMMON/POTMOM/VC(-1:1),VS(-1:1),VT(-2:2),VLS(-2:2),VLA(-2:2), . VQ(-3:3),VQP(-2:2) COMMON/JFAC/ JMMM,JMM,JM,J,JP,JPP,JPPP COMMON/QFAC/ QI2,QF2,QIQF,QI2F2,QI2F22,QIPF2 COMMON/POMRON/GPOM,FPOM2,AMPOM,AMPOM2,AMPOM4, . GA2D,FA2D2,AMA2D,AMA2D2,AMA2D4 COMMON/AMCOEF/ALF,REDM, AMY,AMY2,AMYI,AMYI2, AMN,AMN2,AMNI,AMNI2, . AMYPN,AMYMN, AMYPN2,AMYMN2, AY2PN2,AY2MN2, . AMYN,AMYNI,AMYNI2, AYPNI2,AYMNI2 DIMENSION R(-3:12), S(-3:12), G(-3:12), ELN(15), ERN(15) DATA R/16*0D0/, S/16*0D0/, G/16*0D0/ DATA PI/3.14159265358979323846D0/ XCOM=0.5D0*QI2F2/QIQF JMAX=J+3 KOM=2 FAC=4D0*PI C** Diffractive contribution: pomeron (f, f', A2), pomeron' DO 1000 IN=1,2 IF(IN.EQ.1) THEN AMES2=AMPOM2 GD2=FPOM2 ELSEIF(IN.EQ.2) THEN AMES2=AMA2D2 GD2=FA2D2*ISIGN ENDIF Y=0.5D0*QIQF/AMES2 X=XCOM CALL SEX(X,Y,0D0,JMAX+1,KOM,ELN,ERN) XDC = FAC * GD2*(1D0+AYPNI2/8D0*QI2F2) YDC =-FAC * GD2*AYPNI2/4D0*QIQF XDLS= FAC * GD2*AYPNI2/4D0 XDQ =-FAC * GD2*AMYNI2/16D0 XDNL=-FAC * GD2*AMYNI/4D0*QI2F2 XDQ2= XDQ * QIQF DO 5 K=1,JMAX R(K)=ERN(K+1) S(K)=(K*ERN(K)+(K+1)*ERN(K+2))/(2*K+1) G(K)=(ERN(K+2)-ERN(K))/(2*K+1) 5 CONTINUE R(0)=ERN(1) S(0)=ERN(2) G(0)=(ERN(2)-ERN(1))-2D0*FDEXP(-QIPF2/4D0/AMES2)*AMES2/QIQF VC(-1) = VC(-1) + XDC*R(JM) + YDC*S(JM) VC( 0) = VC( 0) + XDC*R(J) + YDC*S(J) VC( 1) = VC( 1) + XDC*R(JP) + YDC*S(JP) VLS(-2)= VLS(-2)+ XDLS*R(JMM) VLS(-1)= VLS(-1)+ XDLS*R(JM) VLS( 0)= VLS( 0)+ XDLS*R(J) VLS( 1)= VLS( 1)+ XDLS*R(JP) VLS( 2)= VLS( 2)+ XDLS*R(JPP) VQ(-3) = VQ(-3) + XDQ*R(JMMM) VQ(-2) = VQ(-2) + XDQ*R(JMM) VQ(-1) = VQ(-1) + XDQ*R(JM) VQ( 0) = VQ( 0) + XDQ*R(J) VQ( 1) = VQ( 1) + XDQ*R(JP) VQ( 2) = VQ( 2) + XDQ*R(JPP) VQ( 3) = VQ( 3) + XDQ*R(JPPP) IF(NLOC.NE.0) THEN VC(-1) = VC(-1) + XDNL*R(JM) VC( 0) = VC( 0) + XDNL*R(J) VC( 1) = VC( 1) + XDNL*R(JP) ENDIF IF(NQ12.EQ.1) THEN C** Extra contribution from inverse Fourier transform of Q12 VQP(-2) = VQP(-2) + XDQ2*G(JMM) VQP(-1) = VQP(-1) + XDQ2*G(JM) VQP( 0) = VQP( 0) + XDQ2*G(J) VQP( 1) = VQP( 1) + XDQ2*G(JP) VQP( 2) = VQP( 2) + XDQ2*G(JPP) ENDIF 1000 CONTINUE IF(TYPE.EQ.'PP' .OR. TYPE.EQ.'NN') RETURN C** Charged diffractive contribution: pomeron' AMES2=AMA2D2 GD2=FA2D2*FACISO Y=0.5D0*QIQF/AMES2 X=XCOM CALL SEX(X,Y,0D0,JMAX+1,KOM,ELN,ERN) XDC = FAC * ALF*GD2*(1D0+AMYNI/4D0*QI2F2) YDC =-FAC * ALF*GD2*AMYNI/2D0*QIQF XDLS= FAC * ALF*GD2*AMYNI/2D0 XDQ =-FAC * ALF*GD2*AMYNI2/16D0 XDNL= FAC * ALF*GD2*(AYPNI2/4D0-AMYNI)*QI2F2/2D0 XDQ2= XDQ * QIQF DO 6 K=1,JMAX R(K)=ERN(K+1) S(K)=(K*ERN(K)+(K+1)*ERN(K+2))/(2*K+1) G(K)=(ERN(K+2)-ERN(K))/(2*K+1) 6 CONTINUE R(0)=ERN(1) S(0)=ERN(2) G(0)=(ERN(2)-ERN(1))-2D0*FDEXP(-QIPF2/4D0/AMES2)*AMES2/QIQF VC(-1) = VC(-1) + XDC*R(JM) + YDC*S(JM) VC( 0) = VC( 0) + XDC*R(J) + YDC*S(J) VC( 1) = VC( 1) + XDC*R(JP) + YDC*S(JP) VLS(-2)= VLS(-2)+ XDLS*R(JMM) VLS(-1)= VLS(-1)+ XDLS*R(JM) VLS( 0)= VLS( 0)+ XDLS*R(J) VLS( 1)= VLS( 1)+ XDLS*R(JP) VLS( 2)= VLS( 2)+ XDLS*R(JPP) VQ(-3) = VQ(-3) + XDQ*R(JMMM) VQ(-2) = VQ(-2) + XDQ*R(JMM) VQ(-1) = VQ(-1) + XDQ*R(JM) VQ( 0) = VQ( 0) + XDQ*R(J) VQ( 1) = VQ( 1) + XDQ*R(JP) VQ( 2) = VQ( 2) + XDQ*R(JPP) VQ( 3) = VQ( 3) + XDQ*R(JPPP) IF(NLOC.NE.0) THEN VC(-1) = VC(-1) + XDNL*R(JM) VC( 0) = VC( 0) + XDNL*R(J) VC( 1) = VC( 1) + XDNL*R(JP) ENDIF IF(NQ12.EQ.1) THEN C** Extra contribution from inverse Fourier transform of Q12 VQP(-2) = VQP(-2) + XDQ2*G(JMM) VQP(-1) = VQP(-1) + XDQ2*G(JM) VQP( 0) = VQP( 0) + XDQ2*G(J) VQP( 1) = VQP( 1) + XDQ2*G(JP) VQP( 2) = VQP( 2) + XDQ2*G(JPP) ENDIF RETURN END ************************************************************************ FUNCTION FDEXP(X) IMPLICIT REAL*8(A-Z) IF(X.LE.-100D0) THEN FDEXP=0D0 ELSE FDEXP=DEXP(X) ENDIF RETURN END ************************************************************************ DOUBLE PRECISION FUNCTION CE1(X) C** CE1 calculates the function : DEXP(X) * E1(X) IMPLICIT REAL*8 (A-H,O-Z) DATA EPS/10D-14/ DATA GAM/0.577215664901532861D0/ DATA AP0,AP1,AP2,AP3,AP4,AP5,AP6 / . 0.463996004278035D+01, 0.127788778637147D+03, . 0.735910238555843D+03, 0.139583023127254D+04, . 0.101614779141469D+04, 0.286647946600883D+03, . 0.256489038620717D+02/ DATA AQ1,AQ2,AQ3,AQ4,AQ5,AQ6,AQ7 / . 0.512251050448444D+02, 0.503800829553457D+03, . 0.165169408854742D+04, 0.220150290642078D+04, . 0.127719811988873D+04, 0.312294439564262D+03, . 0.256489625816454D+02/ DATA BP0,BP1,BP2,BP3,BP4,BP5 / . 0.335956527252693D+01, 0.204955591333077D+02, . 0.267757325223533D+02, 0.112883678215773D+02, . 0.164680678114210D+01, 0.655193572680895D-01/ DATA BQ1,BQ2,BQ3,BQ4,BQ5,BQ6 / . 0.143836492361913D+02, 0.400563387674630D+02, . 0.366148021121537D+02, 0.128696120312766D+02, . 0.171232738644327D+01, 0.655193403549186D-01/ DATA CP0,CP1,CP2,CP3,CP4 / . 0.298014627030798D+01, 0.113803314436134D+02, . 0.947288802836929D+01, 0.247747160891423D+01, . 0.188516317695352D+00/ DATA CQ1,CQ2,CQ3,CQ4,CQ5 / . 0.988019055335016D+01, 0.189408176576544D+02, . 0.117618585876339D+02, 0.266598761793551D+01, . 0.188516320637495D+00/ DATA DP0,DP1,DP2,DP3 / . 0.242331331460798D+01, 0.432777141801875D+01, . 0.160959648287707D+01, 0.148720388893508D+00/ DATA DQ1,DQ2,DQ3,DQ4 / . 0.558734308280980D+01, 0.578865453197840D+01, . 0.175831677540018D+01, 0.148720389489176D+00/ DATA EP0,EP1,EP2,EP3 / . 0.226526458912153D+01, 0.332000741007556D+01, . 0.104761178441346D+01, 0.837423061701825D-01/ DATA EQ1,EQ2,EQ3,EQ4 / . 0.478887726713541D+01, 0.428387700117901D+01, . 0.113135408983342D+01, 0.837423061723804D-01/ DATA FP0,FP1,FP2,FP3 / . 0.190053654321203D+01, 0.151285969203750D+01, . 0.205314346964057D+00, 0.264152351883344D-03/ DATA FQ1,FQ2,FQ3,FQ4 / . 0.320887608816311D+01, 0.171790987670629D+01, . 0.205578499347658D+00, 0.264152351839874D-03/ PQX=1D0 IF(X.LT.1D0) THEN E1=-GAM-DLOG(X)+X N=1 IX=1 TERM=X 1 N=N+1 IX=-IX TERM=TERM*X/N IF(TERM.LT.EPS) GOTO 10 E1=E1+IX*TERM/N GOTO 1 10 PPX=FDEXP(X)*E1 ELSEIF(X.LT.3D0) THEN PPX=AP0+X*(AP1+X*(AP2+X*(AP3+X*(AP4+X*(AP5+X*AP6))))) PQX=1D0+X*(AQ1+X*(AQ2+X*(AQ3+X*(AQ4+X*(AQ5+X*(AQ6+X*AQ7)))))) ELSEIF(X.LT.6D0) THEN PPX=BP0+X*(BP1+X*(BP2+X*(BP3+X*(BP4+X*BP5)))) PQX=1D0+X*(BQ1+X*(BQ2+X*(BQ3+X*(BQ4+X*(BQ5+X*BQ6))))) ELSEIF(X.LT.14D0) THEN PPX=CP0+X*(CP1+X*(CP2+X*(CP3+X*CP4))) PQX=1D0+X*(CQ1+X*(CQ2+X*(CQ3+X*(CQ4+X*CQ5)))) ELSEIF(X.LT.25D0) THEN PPX=DP0+X*(DP1+X*(DP2+X*DP3)) PQX=1D0+X*(DQ1+X*(DQ2+X*(DQ3+X*DQ4))) ELSEIF(X.LT.70D0) THEN PPX=EP0+X*(EP1+X*(EP2+X*EP3)) PQX=1D0+X*(EQ1+X*(EQ2+X*(EQ3+X*EQ4))) ELSEIF(X.LT.165D0) THEN PPX=FP0+X*(FP1+X*(FP2+X*FP3)) PQX=1D0+X*(FQ1+X*(FQ2+X*(FQ3+X*FQ4))) ELSEIF(X.GE.165D0) THEN Y=1D0/X N=0 IN=1 K=1 CE1=Y 2 N=N+1 IN=-IN Y=Y/X K=K*N TERM=Y*K IF(TERM/CE1.LT.EPS) GOTO 20 CE1=CE1+IN*TERM GOTO 2 20 PPX=CE1 ENDIF CE1=PPX/PQX RETURN END cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Compute c c eln = 0.5*exp(ames**2) * Integral(-1,+1) (Pn(z)/(x-z)*f(x,y)) dz c ern = 0.5*exp(ames**2) * Integral(-1,+1) (Pn(z)*f(x,y)) dz c c with Pn(z) = Legendre function c f(x,y) = exp(-y*(x-z)) c c Input: x = (p**2 + q**2 + m**2)/(2*p*q) c y = 2*p*q/(cutoff**2) c j = number of partial wave + 1 c kom = 2 : compute only rn ( for "diffrac" ) c : else only ln c ames = mass / lambda c c Updated November 1999 subroutine sex(x,y,ames,j,kom,eln,ern) implicit none real*8 eln(15),ern(15),x,y,ames,xh integer i,max,j,kom,maxl logical first data first /.true./ data maxl /20/ if (first) then call plcfs(maxl) first = .false. endif if (y.le.1.d-3) then call smllex(eln,ern,x,y,ames,j,kom) elseif(x.ge.1.d0) then call sexn(eln,ern,x,y,ames,j,kom) elseif(dabs(x).lt.1.d0) then call wex(eln,ern,x,y,ames,j,kom) elseif(x.lt.-1.d0) then xh=-x call sexn(eln,ern,xh,y,ames,j,kom) max=j+2 do 1 i=1,max eln(i)=(-1.d0)**i*eln(i) ern(i)=(-1.d0)**i*ern(i) 1 continue endif end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine sexn(eln,ern,x,y,ames,j,kom) implicit none real*8 x,y,ames real*8 eln(15),ern(15) integer j,kom real*8 un(30),tn(30),anl,anlh real*8 t,ce1p,ce1,ce1m,eps2,sum,xi,dsum,eps integer i,ii,l1,n1,k1,k2,mmax,max real*8 ep,em,xyh,yh,sh,ch,check,amax,dum,grens1,grens2 real*8 grens,yp,yx,enu,xyhh,xh,xmin,xmed,ym,fdexp common /lptab/ anl(24,24) data xyhh /1.d20/ data xh /1.d20/ data xmin /-160.d0/ data xmed /5.d0/ data eps /1.d-13/ data eps2 /1.d-06/ data yh /1.d20/ save tn ym=y-y*x yp=y+y*x yx=y*x if((ym-ames**2).lt.xmin) goto 500 enu=fdexp(ames*ames) ep=0.d0 if(ym.gt.xmin) ep=fdexp(ym)*enu/2.d0 em=0.d0 if(-yp.gt.xmin) em=fdexp(-yp)*enu/2.d0 c same momenta, only different mesonmass -> calculate only un again xyh=yx-ames*ames if(dabs(xyh-xyhh).lt.1.d-10*xyh.and.y.eq.yh) goto 50 9 xyhh=xyh yh=y sh=(ep-em)/y ch=(ep+em)/y max=j+2 if(x.gt.xmed) max=max+20 amax=max check=0.434*(-amax+(y+amax+0.5)*dlog(1+(amax+1)/y)) if (check.gt.7.d0) goto 10 c calculate tn with rec. relation forwards tn(1)=sh do 1 i=2,max tn(i)=ch-(i-1)*tn(i-1)/y dum=ch ch=sh sh=dum 1 continue goto 51 c calculate tn with rec relation backwards 10 if(y.lt.1.d-04) then mmax=amax+4 else grens1=9.197/y grens2=amax*dlog((amax+1)/2.718/y)/2.718/y+8.829/y grens1=grens(grens1)*2.718*y grens2=grens(grens2)*2.718*y k1=grens1+1 k2=grens2+1 mmax=max0(k1,k2) check=0.434*(y*dlog(y)+y+1-(y+0.5)*dlog(y+1)) c check: on loss of more then 7 significant digits if(check.gt.7.d0) print 5003,check endif c now: rec relation backwards t=0.d0 if(mod(mmax,2).eq.0) mmax=mmax+1 sh=enu*dsinh(y)*fdexp(-yx) ch=ep+em do 11 ii=2,mmax i=mmax-ii+2 t=(sh-y*t)/(i-1) if(i.le.(max+1)) tn(i-1)=t dum=sh sh=ch ch=dum 11 continue goto 51 50 if(x.gt.xmed.and.xh.le.xmed) goto 9 c if mesonmass differs but x has passed critical value more c tn's have to be calculated 51 xh=x c FOR DIFFRACTIVE MESONS CALCULATION OF RN IS ENOUGH if(kom.eq.2) goto 800 if(x.gt.xmed) goto 60 c CALC FOR X < XMED USING REC. RELATION FORWARDS ce1p=ep*ce1(-ym) ce1m=em*ce1(yp) c CHECK ON LOSS OF SIGNIFICANT DIGITS if (ce1p.ne.0.d0.or.ce1m.ne.0.d0) then c if (dabs((ce1p-ce1m)/(ce1p+ce1m)).lt.eps2) print 5001 endif un(1)=ce1p-ce1m max=j+1 do i=2,max un(i)=x*un(i-1)-tn(i-1) enddo goto 600 c CALC FOR X > XMED USING REC. RELATION BACKWARDS 60 max=j+1 c FIRST CALCULATE UN(MAX) WITH EXPANSION IN TN sum=tn(max)/x if(sum.eq.0.d0) goto 62 xi=1/x dsum=0.d0 do i=1,20 sum=sum+dsum xi=xi/x dsum=tn(max+i)*xi if(dabs(dsum/sum).lt.eps) goto 62 enddo c print 5002 62 un(max)=sum do ii=2,max i=max+2-ii un(i-1)=(un(i)+tn(i-1))/x enddo goto 600 c Form linear combinations to obtain ln and rn 500 max=j+1 do i=1,30 un(i)=0.d0 tn(i)=0.d0 enddo do i=1,max eln(i)=0.d0 ern(i)=0.d0 enddo return 600 max=j+1 do l1=1,max eln(l1)=0.d0 ern(l1)=0.d0 do n1=1,l1 anlh = anl(n1,l1) eln(l1)=eln(l1)+anlh*un(n1) ern(l1)=ern(l1)+anlh*tn(n1) enddo enddo return 800 max=j+1 do l1=1,max ern(l1)=0.d0 do n1=1,l1 anlh = anl(n1,l1) ern(l1)=ern(l1)+anlh*tn(n1) enddo enddo return 5001 format(1x,'**** SEXN **** LOSS OFF SIGNIFICANT DIGITS UN(1) ') 5002 format(1x,'**** SEXN **** UN MAX NOT ACCURATE ') 5003 format(1x,'**** SEXN **** TN NOT ACCURATE ; LOSS OF :',d10.3, c ' DECIMAL DIGITS ') end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine smllex(eln,ern,x,y,ames,j,kom) implicit none real*8 eln(15),ern(15),qleg(15) real*8 x,y,ames,sign,term,akfac,ak,ank real*8 yx,enu,qx,qxm,tll,tlp,tlm integer j,kom,n,aln,nn,kk,kmin,kmax,ll,lint,k,nmax,max real*8 anl,fdexp common /lptab/ anl(24,24) data qleg /15*0.d0/ data nmax /13/ max = j+3 yx=y*x enu=fdexp(ames*ames) do 10 k=1,max if (ames.ne.0.d0) then lint = k-1 call qj99(lint,x,qx,qxm) qleg(k) = qx else qleg(k) = 0.d0 endif 10 continue do 1 ll=1,max eln(ll) = 0.d0 ern(ll) = 0.d0 lint = ll-1 c COMPUTATION ELN'S: tll = 2.d0*lint+1.d0 tlp = 1.d0*(lint+1)/tll tlm = 1.d0*lint/tll if (lint.eq.0) then eln(ll) = fdexp(-yx)*(qleg(ll)+y*qleg(ll+1)) endif if (lint.gt.0) then eln(ll) = fdexp(-yx)*(qleg(ll)+y*( . tlp*qleg(ll+1)+tlm*qleg(ll-1) ) ) endif c COMPUTATION ERN'S: kmax = 20 kmin = ll akfac = 1.d0 do 21 kk=1,kmin if(kk.ne.1.and.kk.ne.kmin) akfac=(kk-1)*akfac 21 continue do 2 kk=kmin,kmax k = kk-1 if(mod(k,2).eq.0) sign = 1.d0 if(mod(k,2).eq.1) sign =-1.d0 if(k.ne.0) akfac = k*akfac term = 0.d0 do 22 nn=1,nmax if (sign.ne.-1.d0) then aln = anl(nn,ll) if ( aln.ne.0.d0 ) then n = nn - 1 ak = 1.d0*(n+k+1) ank = 1.d0/ak term = term + aln*ank endif endif sign = -sign 22 continue ern(ll) = ern(ll)+fdexp(-yx)*y**k*term/akfac 2 continue eln(ll) = enu*eln(ll) ern(ll) = enu*ern(ll) 1 continue end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c COMPUTATION LEGENDRE FUNCTIONS OF THE SECOND KIND, ARGUMENT X.GT.1 c QXM = FUNCTION QX(JJ-1) , INDEX = JJ c QX = FUNCTION QX(JJ) , INDEX = JJ+1 c X IS INPUT , QX IS OUTPUT subroutine qj99(jj,x,qx,qxm) implicit none real*8 x,qx,qxm,qm2,xj,hyprg99,en,wun,two,ax,eps real*8 tx,cx,bx,zx,z,term0,term1,tj,ak,xi2,fact integer jj,i,max,n,j,nmax data wun /1.d0/ data two /2.d0/ data ax /0.5d0/ data nmax /400/ data eps /1.d-40/ j=jj+1 tj=j if(dabs(x).lt.1.d+04) then qx=dlog((x+wun)/(x-wun))/two if(j.eq.1) return if(j.eq.2) goto 1 endif if(dabs(x).gt.1.d+04) then ak = 1.d0 xi2 = 1.d0/x/x fact = 1.d0/x if(j.eq.1) qx = fact/ak if(j.eq.2) qxm= fact/ak if(j.eq.2) qx = fact/x/(ak+2.d0) 20 ak = ak + 2.d0 fact = fact*xi2 term0 = fact/ak term1 = fact/x/(ak+2.d0) if(dabs(term0).lt.eps) goto 21 if(j.eq.1) qx = qx + term0 if(j.eq.2) qxm= qxm + term0 if(j.eq.2) qx = qx + term1 goto 20 21 if(j.le.2) return endif if(dabs(x).gt.1) then z=wun/(x+dsqrt(x*x-wun)) zx=z*z j=j-1 tj=tj-wun n=0 5 bx=tj cx=bx+ax tx=wun en=0.0d0 max=j-1 do 4 i=1,max en=en+wun tx=tx*en/(two*en+wun) 4 continue qx=tx*(two*z)**j*hyprg99(ax,bx,cx,zx,nmax) if(n.eq.1) return qxm=qx j=j+1 tj=tj+wun n=1 goto 5 endif 1 qxm=qx qx=x*qx-wun if(j.eq.2) return xj=wun do 10 i=3,j qm2=qxm xj=xj+wun qxm=qx qx=((two*xj-wun)*x*qx-(xj-wun)*qm2)/xj 10 continue return end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function hyprg99(a,b,c,z,nmax) implicit none real*8 hyprg99,a,b,c,z,az,bz,cz,dz,s,t,tst,wun,eps integer nmax,i data wun /1.d0/ data eps /1.d-07/ az=a bz=b cz=c dz=wun t=a*b*z/c s=t+wun do 1 i=2,nmax az=az+wun bz=bz+wun cz=cz+wun dz=dz+wun t= t*az*bz/cz/dz*z tst=dabs(t/s) if (tst.le.eps) goto 2 s=s+t 1 continue print 1000,tst 1000 format(///'MAXIMAL NUMBER OF ITERATIONS IN HYPGR , SO STOP ., TST=',d20.5///) 2 hyprg99=s end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine wex(eln,ern,x,y,ames,j,kom) implicit none real*8 eln,ern,un,tn,wn,xx,wx real*8 xyhh,xh,xmin,xmed,yh,anl,fdexp integer nogp integer i,j,k1,k2,max,amax,mmax,ii,kom,ir,lmax,l1 real*8 grens1,grens2,ch,sh,em,xyh,dum,check,t,rint real*8 grens,ym,y,x,yp,yx,ames,enu,ep integer n1,ix logical first dimension eln(15),ern(15) dimension un(30),tn(30),wn(30) dimension xx(20),wx(20) common /lptab/ anl(24,24) data first /.true./ data nogp /20/ data xyhh /1.d20/ data xh /1.d20/ data xmin /-160.d0/ data xmed /5.d0/ data yh /1.d20/ save xx,wx save tn max=j+1 do i=1,max eln(i)=0.d0 ern(i)=0.d0 enddo if (first) then call gauss(-1.d0,1.d0,xx,wx,nogp) first = .false. endif c Calculate basic quantities ym=y-y*x yp=y+y*x yx=y*x if((ym-ames**2).lt.xmin) return enu=fdexp(ames**2) ep=0.d0 if(ym.gt.xmin) ep=fdexp(ym)*enu/2.d0 em=0.d0 if(-yp.gt.xmin) em=fdexp(-yp)*enu/2.d0 c Computation of TN functions: c If same momenta, only diff. mesonmass calculate only un again xyh=yx-ames**2 if(dabs(xyh-xyhh).lt.1.d-10*xyh.and.y.eq.yh) goto 50 9 xyhh=xyh yh=y sh=(ep-em)/y ch=(ep+em)/y max=j+2 if(x.gt.xmed) max=max+20 amax=max check=0.434*(-amax+(y+amax+0.5)*dlog(1+(amax+1)/y)) if (check.le.7.d0) then c Calculate TN with rec. relation forwards tn(1)=sh do i=2,max tn(i)=ch-(i-1)*tn(i-1)/y dum=ch ch=sh sh=dum enddo else c Calculate TN with rec. relation backwards c First: calculate starting point 10 if(y.le.1.d-04) then mmax=amax+4 else grens1=9.197/y grens2=amax*dlog((amax+1)/2.718/y)/2.718/y+8.829/y grens1=grens(grens1)*2.718*y grens2=grens(grens2)*2.718*y k1=grens1+1 k2=grens2+1 mmax=max0(k1,k2) check=0.434*(y*dlog(y)+y+1-(y+0.5)*dlog(y+1)) c Check: on loss of more then 7 significant digits c if(check.gt.7.d0) print 5003,check endif c Now: rec. relation backwards t=0.d0 if(mod(mmax,2).eq.0) mmax=mmax+1 sh=enu*dsinh(y)*fdexp(-yx) ch=ep+em do ii=2,mmax i=mmax-ii+2 t=(sh-y*t)/(i-1) if(i.le.(max+1)) tn(i-1)=t dum=sh sh=ch ch=dum enddo endif goto 51 50 if(x.gt.xmed.and.xh.le.xmed) goto 9 51 xh=x if (kom.eq.1) then c Computation WN functions lmax = j + 1 wn(1) = 0.d0 if (j.ne.1) then do n1 = 1,lmax wn(n1+1)=0.d0 do ir=1,n1 wn(n1+1) = wn(n1+1) + x**(n1-ir)*tn(ir) enddo enddo endif c Computation rest integral: c ENU * DEXP(-Y*X) * 0.5 * INT(-1,+1) (DEXP(Y*X)- DEXP(Y*Z))/(X-Z) DZ rint=0.d0 do ix=1,nogp rint = rint + 0.5d0*wx(ix)* . (1.d0-fdexp(y*(xx(ix)-x)))/(x-xx(ix)) enddo c Computation UN functions do n1=0,lmax un(n1+1) = - wn(n1+1) un(n1+1) = un(n1+1) + 0.5d0*enu*x**n1* . (dlog(dabs((x+1.d0)/(x-1.d0)))-rint) enddo c Computation ELN and ERN functions max=j+1 do l1=1,max do n1=1,l1 eln(l1)=eln(l1)+anl(n1,l1)*un(n1) ern(l1)=ern(l1)+anl(n1,l1)*tn(n1) enddo enddo elseif (kom.eq.2) then c Computation ERN functions max=j+1 do l1=1,max do n1=1,l1 ern(l1)=ern(l1)+anl(n1,l1)*tn(n1) enddo enddo endif 5001 format(1x,'**** WEX **** LOSS OFF SIGNIFICANT DIGITS UN(1) ') 5003 format(1x,'**** WEX **** TN NOT ACCURATE ; LOSS OF :',d10.3, . ' DECIMAL DIGITS ') end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Calculates abscissas and weigths for Gauss-Legendre c integration (n-points Gauss-Legendre quadrature). c Input: x1,x2: begin and end of integral c n: number of integration points c Output: x: abscissas c w: weights c Literature: * Numerical Recipes in Fortran, paragraph 4.5 c * Abramowitz and Stegun, 25.4.29 and 25.4.30 c Only for integrals with finite integration range; c For integrals to infinity; perform a mapping first subroutine gauss(x1,x2,x,w,n) implicit none integer i,j,m,n real*8 x1,x2,x(n),w(n) real*8 eps,xm,xl,z,z1,p1,p2,p3,pp,pi data eps /3.0d-14/ data pi /3.1415926535897932384626433832795028d0/ m=(n+1)/2 xm=0.5d0*(x2+x1) xl=0.5d0*(x2-x1) do i=1,m z=dcos(pi*(i-0.25d0)/(n+0.5d0)) 1 continue p1=1.d0 p2=0.d0 do j=1,n p3=p2 p2=p1 p1=((2.d0*j-1.d0)*z*p2-(j-1.d0)*p3)/j enddo pp=n*(z*p1-p2)/(z*z-1.d0) z1=z z=z1-p1/pp if(dabs(z-z1).le.eps) goto 2 goto 1 2 continue x(i)=xm-xl*z x(n+1-i)=xm+xl*z w(i)=2.d0*xl/((1.d0-z*z)*pp*pp) w(n+1-i)=w(i) enddo end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c P(l) = Sum A_n(l) * z**n c Recurrence relation P(l): P(l+1) = ( (2l+1)*z*P(l) - l*P(l-1) ) / (l+1) c ==> c Recurrence relation A_n(l): c A_n(l+1) = ( (2l+1)*A_{n-1}(l) - l*A_n(l-1) ) / (l+1) cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine plcfs(lmax) implicit none integer lmax integer n1,ll,nn,l,l1,n real*8 anl,tll,tlp,tl,bb,cc common /lptab/ anl(24,24) n1=lmax+1 do 1 nn=1,n1 do 2 ll=1,n1 anl(nn,ll) = 0.d0 2 continue 1 continue anl(1,1) = 1.d0 anl(2,2) = 1.d0 do 3 ll=2,lmax l = ll-1 tll = (2*l+1)*1.d0 tlp = (l+1)*1.d0 tl = l*1.d0 l1 = ll+1 do 4 nn=1,l1 n = nn-1 bb = 0.d0 if((n-1).ge.0) then if(l.ge.(n-1)) bb=anl(n,l+1) endif cc = 0.d0 if(l.ge.1) then if((n+1).le.l) cc=anl(n+1,l) endif anl(n+1,ll+1) = ( tll*bb - tl*cc )/tlp 4 continue 3 continue end cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function grens(x) implicit none real*8 grens,x,a0,a1,a2,a3 data a0 /-0.589654d0/ data a1 /-0.0595734d0/ data a2 /0.649227d0/ data a3 /0.1809910d0/ if(x.gt.1.2d6) then grens=x**0.847 elseif(x.gt.9d4) then grens=x**0.825 elseif(x.gt.6900d0) then grens=x**0.806 elseif(x.gt.460d0) then grens=x**0.781 elseif(x.gt.23d0) then grens=x**0.751 else grens=(x*a3-a1+dsqrt((a1-x*a3)**2-4d0*a2*(a0-x)))/2d0/a2 endif end