c***********************************************************************
SUBROUTINE POTS(U,V)
c
c Calculates the potentials or form factors
c***********************************************************************
c
c IMPLICIT REAL*8(A-H,O-Z)
COMMON/POTTER/DRX,AFACT(2),VFACT,SFACT,ENG,RM,G(4),ETAX,FKX,ETAKX
1 ,RCX,HBARC2,ABETA(3),FLDF(3)
2 ,NX,LAMX,KMXX,KX,IBX,LPLUSX,ICON4,NSPC,IDIRAC,ICHK
CHARACTER*18 B(26)
DIMENSION UT(5),CN(16),CP(16),YLAM(16),PLM(8)
1 ,XG(8),WG(8),U(800),V(800),LDFRM(3)
C DIMENSION UD(800)
EQUIVALENCE (YLAM(1),PLM(1))
C
c
DATA b / ' NX=0 No potential', ' NX=1 VOLUME W-S '
1, ' NX=2 SURFAC, W-S ', ' NX=3 2ND DERIV ', ' NX=4 L.S VOLUME '
2, ' NX=5 L.S SURFACE ', ' NX=6 VOL*R**POWR ', ' NX=7 SUR*R**POWR '
3, ' NX=8 EXTERN FORMF', ' NX=9 HARMONIC OSC', ' NX=10 GAUSSIAN '
4, ' NX=11 DEFORM VOL ', ' NX=12 DEFORM SURF', ' NX=13 Hulthen '
5, ' NX=14 Yukawa Lam ', ' NX=15 Yukawa L*S ', ' NX=16 NO OPTION '
6, ' NX=17 NO OPTION ', ' NX=18 NO OPTION ', ' NX=19 NO OPTION '
7, ' NX=20 VECTOR ', ' NX=21 SCALAR ', ' NX=22 Not used '
8, ' NX=23 Not used ', ' NX=24 Not used ', ' NX=25 Not used '
9/
C
DATA NG,NGX/8,0/
C
ETA4 = 6.0
ETA5 =10.0
SQRPI= 1.772453851
PI=3.141592654
C
IVFLAG=0
ISFLAG=0
FACT=VFACT
70 CONTINUE
C
C READ IN CARD SET 5,6,OR 7 POTENTIAL CARDS
C
READ (5,9000)FZ,VR,RY,AR,VSOR,VI,RZ,AI,VSOI,PWR
C
NZ=ABS(FZ)
RR=ABS(RY)*AFACT(1)
IF(RY.LT.0.0) RR=RR+ABS(RY)*AFACT(2)
RI=ABS(RZ)*AFACT(1)
IF(RZ.LT.0.0) RI=RI+ABS(RZ)*AFACT(2)
IF(ICON4.NE.2.OR.NSPC.LT.3) THEN
WRITE(6,9509)B(NZ+1),VR,RY,AR,RR,VSOR
WRITE(6,9510) VI,RZ,AI,RI,VSOI,PWR
ENDIF
C
KFLAG=0
IF(NX.LT.3) THEN
VR=VR*FACT
VI=VI*FACT
KT=FKX * MAX (RR,RI)+ETA5
LPLUSX=MAX0(LPLUSX,KT)
KT=(2.3*ETA4* MAX (AR,AI)+ MAX (RR,RI))/DRX
ELSE
IF(ENG.EQ.0.0) THEN
KT=(2.3*ETA4* MAX (AR,AI)+ MAX (RR,RI))/DRX
ELSE
RM= MAX (RM,RR)
RM= MAX (RM,RI)
IF(RM.EQ.0.0) RM=1.0
VR=VR*FACT
VI=VI*FACT
KT=(2.3*ETA4/SQRT(FKX**2+2.0*ETAX*FKX/RM))/DRX
ENDIF
ENDIF
KX=MIN0(MAX0(KX,KT),KMXX)
83 CONTINUE
IF(AR.EQ.0.0) THEN
F1=0.0
F2=0.0
ELSE
F2=EXP(-DRX /AR)
F1=EXP( RR/AR)
ENDIF
IF(AI.EQ.0.0) THEN
F3=0.0
F4=0.0
ELSE
F4=EXP(-DRX /AI)
F3=EXP( RI/AI)
ENDIF
C
IF(NX.GE.3.AND.ENG.EQ.0.0.AND.NZ.LE.5) THEN
IF(AR.NE.0.0) VR=VR*(RR/AR)**(PWR+1)
IF(AI.NE.0.0) VI=VI*(RI/AI)**(PWR+1)
ENDIF
IF(NZ.EQ.0) GO TO 6050
GO TO ( 100, 200, 300, 400, 500, 600, 700, 800, 900,1000),NZ
GO TO (1100,1200,1300,1400,1500,1600,1700,1800,1900,2000),NZ-10
GO TO (2100),NZ-20
write(6,'(''0Invalid potential option '',i3)')nz
C
C VOLUME WOODS SAXON
C
100 CONTINUE
DO 160 M=1,KX
MK=M+M-1
F1=F1*F2
F3=F3*F4
U(MK )=U(MK )-VR*F1/(1.0+F1)
U(MK+1)=U(MK+1)-VI*F3/(1.0+F3)
160 CONTINUE
GO TO 6000
C
C 1ST DERIVATIVE WOODS SAXON
C
200 CONTINUE
DO 260 M=1,KX
MK=M+M-1
F1=F1*F2
F3=F3*F4
U(MK )=U(MK )+VR*F1/(1.0+F1)**2
U(MK+1)=U(MK+1)+VI*F3/(1.0+F3)**2
260 CONTINUE
GO TO 6000
C
C 2ND DERIVATIVE WOODS SAXON
C
300 CONTINUE
DO 360 M=1,KX
MK=M+M-1
F1=F1*F2
F3=F3*F4
U(MK )=U(MK )-VR*F1*(1.0-F1)/(1.0+F1)**3
U(MK+1)=U(MK+1)-VI*F3*(1.0-F3)/(1.0+F3)**3
360 CONTINUE
GO TO 6000
C
C L.S VOLUME WOODS SAXON
C
400 CONTINUE
IBX=1
IF(AR.NE.0.0) VR=VR/AR
IF(AI.NE.0.0) VI=VI/AI
R=0.0
DO 460 M=1,KX
R=R+DRX
MK=M+M-1
F1=F1*F2
F3=F3*F4
V(MK )=V(MK )-VR*F1/(R*(1.0+F1)**2)
V(MK+1)=V(MK+1)-VI*F3/(R*(1.0+F3)**2)
460 CONTINUE
GO TO 6000
C
C L.S 1ST DERIVATIVE WOODS SAXON
C
500 CONTINUE
IBX=1
IF(AR.NE.0.0) VR=VR/AR
IF(AI.NE.0.0) VI=VI/AI
R=0.0
DO 560 M=1,KX
R=R+DRX
MK=M+M-1
F1=F1*F2
F3=F3*F4
V(MK )=V(MK )+VR*F1*(1.0-F1)/(R*(1.0+F1)**3)
V(MK+1)=V(MK+1)+VI*F3*(1.0-F3)/(R*(1.0+F3)**3)
560 CONTINUE
GO TO 6000
C
C WOOD SAXON*R**PWR
C
600 CONTINUE
R=0.0
DO 660 M=1,KX
R=R+DRX
MK=M+M-1
F1=F1*F2
F3=F3*F4
RPWR=R**PWR
U(MK )=U(MK )-VR*F1*RPWR/(1.0+F1)
U(MK+1)=U(MK+1)-VI*F3*RPWR/(1.0+F3)
660 CONTINUE
GO TO 6000
C
C 1ST DERIVATIVE WOOD SAXON*R**PWR
C
700 CONTINUE
R=0.0
DO 760 M=1,KX
MK=M+M-1
R=R+DRX
F1=F1*F2
F3=F3*F4
RPWR=R**PWR
U(MK )=U(MK )+VR*F1*RPWR/(1.0+F1)**2
U(MK+1)=U(MK+1)+VI*F3*RPWR/(1.0+F3)**2
760 CONTINUE
GO TO 6000
C
C EXTERNAL FORM FACTOR
C
800 CONTINUE
IF(NX.GE.3) THEN
READ (5,9000)G
WRITE(6,9508)G
ENDIF
READ(5,9000)F1,F2
C
C F2 = 0 REAL CENTRAL
C F2 = 1 IMAG CENTRAL
C F2 = 2 REAL SPIN ORBIT
C F2 = 3 IMAG SPIN ORBIT
C
IF(F2.EQ.0.0.OR.F2.EQ.2) THEN
F3=VR
MK=1
ELSE
F3=VI
MK=2
ENDIF
IF(F3.EQ.0.0) F3=1.0
KT=F1
DO 820 M=1,KT,5
READ(5,9100)UT
IF(M.GT.KMXX) GO TO 820
DO 810 I=1,5
IF(F2.LT.2.0) THEN
U(MK )=U(MK )+UT(I)*F3
ELSE
V(MK )=V(MK )+UT(I)*F3
ENDIF
MK=MK+2
810 CONTINUE
820 CONTINUE
C
KX=min0(KT,kmxx)
GO TO 6000
C
C HARMONIC OSCILLATOR, NORMALIZED N*EXP(-(R/RY)**2/2)
C
900 CONTINUE
READ (5,9000)G
WRITE(6,9508)G
F1=1.0/RY**2
F2=F1/RY
F3=0.5
F4=SQRPI*0.5
L=G(2)
IF(L.NE.0) THEN
DO 930 I=1,L
F3=F3+1.0
F4=F4*F3
F2=F2*F1
930 CONTINUE
ENDIF
NN=G(1)
T1=1.0
T2=F4
C LAGUERRE POLYNOMIAL COEFFICIENTS = (Abramowitz-Stegun)*(-1)**nn
CN(1)=(-1.0)**NN
IF(NN.NE.0) THEN
DO 940 I=1,NN
F3=F3+1.0
T1=T1*FLOAT(I)
T2=T2*F3
CN(I+1)=-CN(I)*F1*FLOAT(NN+1-I)/(FLOAT(I)*F3)
940 CONTINUE
ENDIF
ANORM=SQRT(2.0*F2*T1/T2)*T2/(T1*F4)
IF(VR.NE.0.0) ANORM=VR*ANORM
KT=10.0*RY/DRX
KT=MIN0(KT,KMXX)
R=0.0
F1=F1/2.0
DO 960 M=1,KT
MK=M+M-1
R=R+DRX
R2=R**2
F2=CN(1)
IF(NN.EQ.0) GO TO 951
F3=1.0
DO 950 I=1,NN
F3=F3*R2
F2=F2+CN(I+1)*F3
950 CONTINUE
951 CONTINUE
U(MK )=U(MK )+F2*ANORM*EXP(-F1*R2)*R**L
960 CONTINUE
GO TO 6000
C
C EXP(-(R/R0)**2)*R**POWR
C
1000 CONTINUE
IF(VR.NE.0.0) THEN
R=0.0
DO 1060 M=1,KX
MK=M+M-1
R=R+DRX
U(MK )=U(MK )-VR*EXP(-(R/RY)**2)*R**PWR
1060 CONTINUE
endif
IF(VI.NE.0.0) THEN
R=0.0
DO 1070 M=1,KX
MK=M+M-1
R=R+DRX
U(MK+1)=U(MK+1)-VI*EXP(-(R/RZ)**2)*R**PWR
1070 CONTINUE
endif
GO TO 6000
C
C DEFORMED VOLUME OR SURFACE BY YLM EXPANSION
C
1100 CONTINUE
IF(NGX.NE.NG) THEN
CALL LEGAUS(2*NG,XG,WG)
NGX=NG
ENDIF
T2=(-1.0)**LAMX
IF(ICHK.EQ.0) THEN
READ (5,9000) (ABETA(J),FLDF(J) ,J=1,3)
WRITE(6,9000)
WRITE(6,9512) (ABETA(J),FLDF(J) ,J=1,3)
ENDIF
LMAX=LAMX+1
DO 1101 J=1,3
LDFRM(J)=FLDF(J)
LMAX=MAX0(LMAX,LDFRM(J)+1)
1101 CONTINUE
T2=(-1.0)**LAMX
C
DO 1108 I=1,NG
CN(I )=0.0
CN(I+NG)=0.0
P2=0.0
P3=1.0
C
DO 1106 M=1,LMAX
L=M-1
FL=L-1
IF(L.EQ.0) GO TO 1102
P3=((2.0*FL+1.0)*XG(I)*P2-FL*P1)/(FL+1.0)
1102 CONTINUE
DO 1103 J=1,3
IF(ABETA(J).EQ.0.0) GO TO 1103
IF(L.NE.LDFRM(J)) GO TO 1103
FACTOR=P3*ABETA(J)*SQRT(FLOAT(L+L+1)/(4.0*PI))
CN(I )=CN(I )+FACTOR
CN(I+NG)=CN(I+8)+FACTOR*(-1.0)**LDFRM(J)
1103 CONTINUE
IF(L.NE.LAMX) GO TO 1104
YLAM(I )= P3*WG(I)*SQRT(FLOAT(L+L+1)*PI)
YLAM(I+NG)=YLAM(I )*T2
1104 CONTINUE
P1=P2
P2=P3
1106 CONTINUE
1108 CONTINUE
DO 1170 I=1,2
IF(I.EQ.1) THEN
IF(VR.EQ.0.0) GO TO 1170
VX=VR
RX=RR
AX=AR
F1=1.0
IFL=-1
ELSE
IF(VI.EQ.0.0) GO TO 1170
VX=VI
RX=RI
AX=AI
F1=1.0
F2=F4
IFL=0
ENDIF
DO 1135 J=1,16
CP(J)=EXP((1.0+CN(J))*RX/AX)
1135 CONTINUE
IF(LAMX.EQ.0) VX=VX/(SQRPI*2.0)
J=NZ-10
IF(J.EQ.2) GO TO 1151
DO 1150 M=1,KX
MK=M+M+IFL
VTEMP=0.0
F1=F1*F2
DO 1145 J=1,16
F3=F1*CP(J)
VTEMP=VTEMP-YLAM(J)*VX*F3/(1.0+F3)
1145 CONTINUE
U(MK )=U(MK )+VTEMP
1150 CONTINUE
GO TO 1170
1151 CONTINUE
DO 1160 M=1,KX
MK=M+M+IFL
VTEMP=0.0
F1=F1*F2
DO 1155 J=1,16
F3=F1*CP(J)
VTEMP=VTEMP+YLAM(J)*VX*F3/(1.0+F3)**2
1155 CONTINUE
U(MK )=U(MK )+VTEMP
1160 CONTINUE
1170 CONTINUE
GO TO 6000
1200 CONTINUE
GO TO 1100
C
C HULTHEN WAVE FUNCTION= NORM*(EXP(-R/RY)-EXP(-R/RZ))/R
C
1300 CONTINUE
READ (5,9000)G
WRITE(6,9508)G
T1=1.0/RY
T2=1.0/RZ
T3=T2**2-T1**2
T4=SQRT(2.0*T1*(T1+T2)*T2)/(T2-T1)
IF(VR.NE.0.0) KT=16.0* MIN (RY,RZ)/DRX
IF(VR.EQ.0.0) KT=16.0* MAX (RY,RZ)/DRX
KX=MIN0(KX,KMXX)
F1=1.0
F2=EXP(-DRX *T1)
F3=1.0
F4=EXP(-DRX *T2)
R=0.0
DO 1260 M=1,KX
MK=M+M-1
R=R+DRX
F1=F1*F2
F3=F3*F4
V(M)=T3*F3/(F1-F3)
TEMP=1.0
IF(VR.NE.0.0) TEMP=V(M)/FACT
U(MK )=U(MK )+TEMP*T4*(F1-F3)/R
1260 CONTINUE
GO TO 6000
c
c Yukawa L=LAM
c
c v = V-0*exp(-r/rx)/(r/rx) ay < r
c = Wood-Saxon r < ay
c
1400 CONTINUE
lam=lamx
if(ry.eq.0.0) ry=1.0
if(rz.eq.0.0) rz=1.0
if(ar.eq.0.0) ar=drx
if(ai.eq.0.0) ai=drx
f1=1.0
f2=exp(-drx/ry)
xr=ar/ry
t1=1.0+1.0/xr
t2=1.0
c Recur Hankel functions
do 1420 i=0,lam
t0=t1
t1=t2
t2=float(2*i-1)*t1/xr+t0
1420 continue
vzr =2.0*vr*t2*exp(-xr)/xr
betar=2.0*(float(lam+1)*t2/xr+t1)/t2
f3=1.0
f4=exp(-drx/rz)
xi=ai/rz
t1=1.0+1.0/xi
t2=1.0
c Recur Hankel functions
do 1430 i=0,lam
t0=t1
t1=t2
t2=float(2*i-1)*t1/xi+t0
1430 continue
vzi =2.0*vi*t2*exp(-xi)/xi
betai=2.0*(float(lam+1)*t2/xi+t1)/t2
c
if(vr.ne.0.0) then
r =0.0
do 1460 m=1,k
r=r+drx
mk=m+m-1
f1=f1*f2
f3=f3*f4
if(r.lt.ar) then
u(mk )=u(mk )-vzr/(1.0+exp(betar*(r/ry-xr)))
else
xr=r/ry
t1=1.0+1.0/xr
t2=1.0
c Recur Hankel functions
do 1450 i=0,lam
t0=t1
t1=t2
t2=float(2*i-1)*t1/xr+t0
1450 continue
u(mk )=u(mk )-vr*t2*f1/xr
endif
1460 continue
endif
if(vi.ne.0.0) then
r =0.0
do 1480 m=1,k
r=r+drx
mk=m+m-1
if(r.lt.ai) then
u(mk+1)=u(mk+1)-vzi/(1.0+exp(betai*(r/rz-xi)))
else
xi=r/rz
t1=1.0+1.0/xi
t2=1.0
c Recur Hankel functions
do 1470 i=0,lam
t0=t1
t1=t2
t2=float(2*i-1)*t1/xi+t0
1470 continue
u(mk+1)=u(mk+1)-vi*t2*f3/xi
endif
1480 continue
endif
go to 6000
c
c Yukawa L*S
c
1500 CONTINUE
r =0.0
if(ry.eq.0.0) ry=1.0
if(rz.eq.0.0) rz=1.0
if(ar.eq.0.0) ar=1.0
if(ai.eq.0.0) ai=1.0
f1=1.0
f2=exp(-drx/ry)
xr=ar/ry
betar=2.0*(1.0+3.0/xr+3.0/xr**2)/(1.0+1.0/xr)
vzr =2.0*vr*exp(-xr)*(1.0+1.0/xr)/xr**2
f3=1.0
f4=exp(-drx/rz)
xi=ai/rz
betai=2.0*(1.0+3.0/xi+3.0/xi**2)/(1.0+1.0/xi)
vzi =2.0*vi*exp(-xi)*(1.0+1.0/xi)/xi**2
do 1560 m=1,k
r=r+drx
mk=m+m-1
f1=f1*f2
f3=f3*f4
if(r.lt.ar) then
u(mk )=u(mk )-vzr/(1.0+exp(betar*(r/ry-xr)))
else
u(mk )=u(mk )-vr *f1*(1.0+1.0*ry/r)/(r/ry)**2
endif
if(r.lt.ai) then
u(mk+1)=u(mk+1)-vzi/(1.0+exp(betai*(r/rz-xi)))
else
u(mk+1)=u(mk+1)-vi *f3*(1.0+1.0*rz/r)/(r/rz)**2
endif
1560 continue
go to 6000
1600 CONTINUE
1700 CONTINUE
go to 6000
1800 CONTINUE
DO 1860 M=1,KX
MK=M+M-1
F1=F1*F2
F3=F3*F4
V(MK )=V(MK )-VR*F1/(1.0+F1)
V(MK+1)=V(MK+1)-VI*F3/(1.0+F3)
1860 CONTINUE
GO TO 6000
1900 CONTINUE
c
c Coulomb excitation for a deformed uniform charge distribution
c
IF(ICHK.EQ.0) THEN
READ (5,9000) (ABETA(J),FLDF(J) ,J=1,3)
WRITE(6,9000)
WRITE(6,9512) (ABETA(J),FLDF(J) ,J=1,3)
ENDIF
if(vr.eq.0.0) vr=1.0
do 1990 i=1,3
beta=abeta(i)
if(beta.ne.0.0) then
c set flag for unbound stripping evaluation of coulomb ex.
ibx=4
flam1=beta*sqrt(float(2*lamx+1)/(4.0*pi))
flam2=flam1*phasef(lamx)
c
if(ngx.ne.ng) then
ngx=ng
call legaus(2*ng,xg,wg)
endif
c
an=0.0
do 1920 k=1,ng
p1=0.0
p2=1.0
if(lamx.gt.0) then
do 1910 j=1,lamx
p3=(float(2*j-1)*xg(k)*p2 - float(j-1)*p1)/float(j)
p1=p2
p2=p3
1910 continue
endif
plm(k)=p2
c
c calculate normalization
c
r1=rcx*(1.0+flam1*p2)
r2=rcx*(1.0+flam2*p2)
an=an+(r1**3+r2**3)*wg(k)
1920 continue
c
an =an/3.0
c
r =0.0
c
do 1940 m=1,kx
mk=m+m-1
r=r+drx
sum=0.0
do 1930 k=1,ng
r1=rcx*(1.0+flam1*plm(k))
r2=rcx*(1.0+flam2*plm(k))
c
if(r.gt.r1) then
s1=r1**(lamx+3)/(r**(lamx+1)*float(lamx+3))
else
if(lamx.eq.2) then
s1=( log(r1/r)+1.0/float(lamx+3))*r**2
else
s1=(r**lamx/r1**(lamx-2)
1 -float(2*lamx+1)*r**2/float(lamx+3))/float(2-lamx)
endif
endif
c
if(r.gt.r2) then
s2=r2**(lamx+3)/(r**(lamx+1)*float(lamx+3))
else
if(lamx.eq.2) then
s2=( log(r2/r)+1.0/float(lamx+3))*r**2
else
s2=(r**lamx/r2**(lamx-2)
1 -float(2*lamx+1)*r**2/float(lamx+3))/float(2-lamx)
endif
endif
c
sum=sum+(s1+s2*phasef(lamx))*wg(k)*plm(k)
1930 continue
c
sum=sum/an
if(beta.ne.0.0.and.lamx.ne.0) sum=sum/flam1
u(mk )=u(mk )+sum*(vr*etakx/fact)
1940 continue
endif
1990 continue
GO TO 6000
C
C VECTOR POTENTIAL
C VSOR, VSOI ARE THE THIRD PARAMETERS
C IN THE 3 PARAMETER FERMI MODEL
C [1 + VSO?*(R/R?)**2]
C
2000 CONTINUE
IVFLAG=-1
R=0.0
DO 2060 M=1,KX
MK=M+M-1
R=R+DRX
F1=F1*F2
F3=F3*F4
U(MK )=U(MK )-VR*F1*(1.0+VSOR*(R/RR)**2)/(1.0+F1)
U(MK+1)=U(MK+1)-VI*F3*(1.0+VSOI*(R/RI)**2)/(1.0+F3)
2060 CONTINUE
GO TO 6020
C
C SCALAR POTENTIAL
C VSOR, VSOI ARE THE FERMI THIRD PARAMETERS
C IN THE 3 PARAMETER FERMI MODEL
C [1 + VSO?*(R/R?)**2]
C
2100 CONTINUE
ISFLAG=-1
VR=VR*(SFACT/VFACT)
VI=VI*(SFACT/VFACT)
R=0.0
DO 2160 M=1,KX
MK=M+M-1
R=R+DRX
F1=F1*F2
F3=F3*F4
V(MK )=V(MK )-VR*F1*(1.0+VSOR*(R/RR)**2)/(1.0+F1)
V(MK+1)=V(MK+1)-VI*F3*(1.0+VSOI*(R/RI)**2)/(1.0+F3)
2160 CONTINUE
GO TO 6020
C
C END OF POTENTIALS
C
6000 CONTINUE
IDIRAC=1
6020 CONTINUE
IF(KFLAG.NE.0.OR.NZ.GT.5) GO TO 6050
IF(ABS(VSOR)+ABS(VSOI).EQ.0.0) GO TO 6050
NZ=NZ+3
KFLAG=1
VR=VR*VSOR/45.2
VI=VI*VSOI/45.2
GO TO 83
6050 CONTINUE
IF(FZ.GT.0.0) GO TO 70
C
C PROCESS DIRAC POTENTIALS
C
C ENTRY WITH U -> K**2 - VFACT*V VFACT = 2.0*W1 /HBARC**2
C V -> - SFACT*S SFACT = 2.0*FM1/HBARC**2
C
IF(IVFLAG.NE.0.AND.ISFLAG.NE.0) THEN
IF(IDIRAC.EQ.1) THEN
WRITE(6,9515)
ENDIF
IDIRAC=-1
C
C KT2=KX+KX
C WRITE(20,7777)' ENTRY POTENTIALS'
C WRITE(20,7778)(U(I),I=1,KT2)
C WRITE(20,7777)
C WRITE(20,7778)(V(I),I=1,KT2)
W1M1=(VFACT+SFACT)*0.5*HBARC2
DO 6100 M=1,Kx
MK=M+M-1
VVR=(U(MK )-FKX**2)/VFACT
VVI= U(MK+1) /VFACT
VSR= V(MK ) /SFACT
VSI= V(MK+1) /SFACT
U(MK )=U(MK )+V(MK ) +(VVR**2-VVI**2 - VSR**2+VSI**2)/HBARC2
U(MK+1)=U(MK+1)+V(MK+1) +(VVR*VVI - VSR*VSI)*2.0/HBARC2
T1 = W1M1 + VVR-VSR
T2 = VVI-VSI
V(MK )=0.5* LOG (T1**2 + T2**2)
V(MK+1)=ATAN2(T2,T1)
6100 CONTINUE
C WRITE(20,7777)' SECOND POTENTIALS'
C WRITE(20,7778)(U(I),I=1,KT2)
C WRITE(20,7777)
C WRITE(20,7778)(V(I),I=1,KT2)
C
R=FLOAT(KX+1)*DRX
MK=KX+KX-1
D2=V(MK )
D1=D2
A2=V(MK+1)
A1=A2
IBX=1
DO 6150 M=2,KX
MK=MK-2
R=R-DRX
D3=D2
D2=D1
D1=V(MK )
A3=A2
A2=A1
A1=V(MK+1)
C FIRST DERIVATIVE TERMS
DPR=(D3-D1)/(2.0*DRX)
DPI=(A3-A1)/(2.0*DRX)
V(MK+2)=2.0*DPR/R
V(MK+3)=2.0*DPI/R
C SECOND DERIVATIVE TERMS
DPPR=(D3+D1-2.0*D2)/DRX**2
DPPI=(A3+A1-2.0*A2)/DRX**2
UDR =0.5*DPPR-0.25*(DPR**2-DPI**2)+DPR/R
UDI =0.5*DPPI-0.25*(2.0*DPR*DPI )+DPI/R
C UD(MK+2)=UDR
C UD(MK+3)=UDI
U(MK+2)=U(MK+2)+UDR
U(MK+3)=U(MK+3)+UDI
6150 CONTINUE
V(1 )=V(3 )*2.0
V(2 )=V(4 )*2.0
U(1 )=U(1 )+UDR
U(2 )=U(2 )+UDI
C WRITE(20,7777)'TERTIARY POTENTIALS'
C WRITE(20,7777)'CENTRAL POTENTIAL'
C WRITE(20,7778)(U(I),I=1,KT2)
C WRITE(20,7777)'SPIN ORBIT POTENTIAL'
C WRITE(20,7778)(V(I),I=1,KT2)
C WRITE(20,7777)'UD - DARWIN TERM'
C WRITE(20,7778)(UD(I),I=1,KT2)
C 7777 FORMAT(A30)
C 7778 FORMAT(' ',1P10E12.4)
ENDIF
C
IF(IDIRAC.EQ.1) IDIRAC=0
RETURN
C
9000 FORMAT(10F8.4)
9100 FORMAT(5E16.7)
9508 FORMAT(18X,9H NODES=,F9.4,9H L =,F9.4,9H 2*J =,F9.4
1 ,9H 2*S =,F9.4,9H FISW =,F9.4)
9509 FORMAT(A18,9H V RL =,F9.4,9H R0RL =,F9.4,9H A RL =,F9.4
1 ,9H R RL =,F9.4,9H VSOR =,F9.4)
9510 FORMAT(18X,9H V IM =,F9.4,9H R0IM =,F9.4,9H A IM =,F9.4
1 ,9H R IM =,F9.4,9H VSOI =,F9.4,9H POWR =,F9.4)
9512 FORMAT(18X,9H BETA1=,F9.4,9H LDFR1=,F9.4,9H BETA2=,F9.4
1 ,9H LDFR2=,F9.4,9H BETA3=,F9.4,9H LDFR3=,F9.4)
9515 FORMAT('0',20('*'),' WARNING, Mixing of Dirac and non-Dirac'
1, ' potentials may be hazardous to your calculation'
2, 20('*'))
END