PtolemyGUI/dwuck4/ADWCK4.FOR
2025-02-12 23:34:08 -05:00

1076 lines
30 KiB
Fortran

c$debug
c***********************************************************************
SUBROUTINE ADWUCK4
c
c Control program for first portion of DWUCK4 which reads in data
c and calculates the form factors
c***********************************************************************
c
parameter(ispc0 = 4010, ispc1 = 8000, ispc2 = 8000)
IMPLICIT REAL*8(A-H,O-Z)
COMMON ALPHA(15),IDAT(6),ICON(20),ANGLE(5),HBARC,AMU,AMASS,CHSQ,AA
1,DRF,Z(3),ZA(3),FM(3),FMA(3),RC(3),AC(3),PNLOC(3),FS(3),ECM(3)
2,FK(3),FK2(3),ETA(3),DR(3),FMU(3),FN,FL,FJ2,FSS,VCE,FNRNG,RSIG(2)
3,K,KZ,LPLUS,LPL2,IS(3),NS(3),NLTR,LTRT(8),JTRT(8),ISTRT(8),IBF(8)
4,KC,IBUFF,IWORD,ILINE,JLINE
Common/array0/space0(ispc0)
Common/array1/space1(ispc1)
Common/array2/space2(ispc2)
DIMENSION U(800,2), V(800,2), UB(800), FF(10), VB(800)
EQUIVALENCE (space0( 1), U ), (space0(1601), V)
1 ,(space0(3201), UB), (space0(4001), FF)
2 ,(space1( 1), VB)
C
IF(ICON(1).EQ.0) then
c
c Use standard angle data
c
II=(ANGLE(1)-ANGLE(2))/ANGLE(3)+1.0
IF(II.LT.0) II=0
ANGLE(1)=II
else
C
C READ CARD SET 2, Angle data
c
READ (5,9001)ANGLE
endif
WRITE(6,9010)ANGLE
C
C
C READ CARD SET 3, Angular momentum parameters
C
c L Maximum L for partial waves
c Nltr Number of angular momentum transfers
c Ltrt Orbital angular momentum transfer
c Jtrt Total angular momentum transfer
c
READ (5,9002)L,NLTR,(LTRT(I),I=1,NLTR),(JTRT(I),I=1,NLTR)
WRITE(6,9008)L ,(LTRT(I),I=1,NLTR)
WRITE(6,9009)NLTR,(JTRT(I),I=1,NLTR)
LC=0
IF(L.LT.0) LC=IABS(L)
L=IABS(L)
C
C
C READ CARD SET 4
c
c Drf Radial increment
c Rz Lower radial cutoff
c Rmaz Upper Radial cutoff
c Vce Coulomb excitation scale factor for inelastic scattering reactions
c Fnrng Finite range parameter for stripping reactions
c Amass Mass for scaling radial increment, if zero, then scaling automatic
c
READ (5,9001)DRF,RZ,RMAX,VCE,FNRNG,AMASS
IF(ABS(RMAX).LT.ABS(RZ)) THEN
temp=RZ
RZ=RMAX
RMAX=temp
ENDIF
WRITE(6,9011)DRF,RZ,RMAX,VCE,FNRNG
IF(AMASS.NE.0.0) WRITE(6,9012)AMASS
C default value for drf
IF(DRF.EQ.0.0) DRF=0.1
KZ=ABS(RZ )/DRF+1.0E-08
K =ABS(RMAX)/DRF+1.0E-08
KC=2*((K+1)/2)
LPLUS=L+1
IND=0
INCR=0
IBF(5)=0
IBF(6)=0
WRITE(6,'(a)')' PARTICLE DATA '
C
C READ IN DISTORTED WAVE INFORMATION
C
c ibf(4) is flag to shut off spins if no spin orbit potential
c
IBF(4)=0
c Initial wave
CALL FORMF(U(1,1),V(1,1),1,0,0,IVB)
IF(IBF(4).NE.0) IBF(5)=1
IBF(4)=0
c Final wave
CALL FORMF(U(1,2),V(1,2),2,0,0,IVB)
IF(IBF(4).NE.0) IBF(6)=1
C
C FORM FACTOR LOOP, up to 8 form factors may be used
C
ILINE=0
DO 35 II=1,NLTR
IBF(4)=0
WRITE(6,9999)ALPHA,(IDAT(I),I=1,3)
IF(ICON(3).NE.0.or.II.EQ.1) then
DO 10 I=1,800
UB(I)=0.0
10 CONTINUE
ivb=0
do 11 I=1,10
FF(I) = 0.0
11 continue
c IF
IF(ICON(2).EQ.0.or.ICON(2).EQ.3) then
C
C Here for collective or single particle form factors
C
WRITE(6,'(a)')'0FORM FACTOR DATA '
IF(RMAX.LT.0.0) K=KC
C
CALL FORMF(UB,VB,3,LTRT(II),0,IVB)
if(FK(3) .ne. 0.0) then
FF( 7)=eta(3)*fk(3)
FF( 8)=float(IVB)
FF( 9)=float(kc)
else
FK(3) = sqrt(abs(FF(5)))
IVB = FF(8)
endif
C
C HERE FOR 2 NUCLEON FORM FACTORS
C
c ELSEIF
ELSEIF(ICON(2).EQ.2.OR.ICON(2).EQ.3) THEN
CALL CATHEN(LTRT(II),JTRT(II),space2( 1),space2(1601),UB)
ELSE
write(6,'(a)')'0Warning, no form factor computed'
endif
endif
c
C L.S FORM FACTOR KLUDGE
C Move Spin-Orbit form factor to UB
IF(IVB.EQ.1) THEN
DO 14 M=1,800
UB(M)=VB(M)
14 CONTINUE
FF( 8)=float(ivb)
ENDIF
ISTRT(II)=IS(3)
IF(ICON(4).EQ.0.or.icon(4).eq.3) then
IWRITE=2
call ffprint(iwrite,k,alpha,ub(1),dr(3)
1 ,ltrt(ii),istrt(ii),jtrt(ii),idat)
endif
C
IF(ICON(3).NE.0.or.II.EQ.1) then
C Apply non-local correction to form factor
CALL FNLOC(U,V,UB,VB,LTRT(II),Z,FNRNG,K)
endif
IF(VCE.NE.0.0) then
c Store end point for Vcx (= Rmax*FF) for integration into complex R plane
ivb = 4
FF( 1)=3.0*VCE*CHSQ*Z(3)*ZA(3)
1 *(RC(3)/(FLOAT(KC)*DR(3)))**LTRT(II)/FLOAT(2*LTRT(II)+1)
FF( 2)=0.0
FF( 3)=0.0
FF( 4)=1.0
FF( 5)=0.0
FF( 6)=0.0
FF( 7)=0.0
FF( 8)=float(IVB)
if(ltrt(ii).eq.0) ff( 1)=ff( 1)/3.0
ENDIF
IF(IVB.EQ.2.OR.IVB.EQ.4) THEN
C
C SET ERROR TRAPS FOR UNBOUND STRIPPING AND COULOMB EXCITATION CASES
C
ILINE=-1
C RESTRICT MAX L TRANSFER IN UNBOUND STRIPPING - STORAGE LIMITED
IF(LTRT(II).GT.7) THEN
WRITE(6,9903)LTRT(II)
IBF(3)=1
ENDIF
IF(FM(1).LT.FM(2)) THEN
C SET ERROR FLAG
IBF(3)=1
WRITE(6,'(a)')'0ERROR * UNBOUND STATE FOR PICK UP NOT ALLOWED'
ENDIF
C
IF(IVB.EQ.2) THEN
T1=(SQRT(FLOAT(LTRT(II)*(LTRT(II)+1))+ETA(3)**2)+ETA(3))/FK(3)
IF(ABS(RMAX).LT.T1) WRITE(6,9904)T1,RMAX
ENDIF
L = FK(1)*FLOAT(KC)*DR(1) + 12.5
IF(L.LT.LPLUS) WRITE(6,9905)
C
ENDIF
C
C STORE FF PARAMETERS FOR USE IN UNBND STRIP. AND Coulomb excitation
C
c write UB and FF onto disk
WRITE(4) UB,FF
IS(3)=ISTRT(II)
IF(IS(3).EQ.0) JTRT(II)=LTRT(II)+LTRT(II)
IND=MAX0(LTRT(II)+1,IND)
INCR=MAX0(INCR,JTRT(II))
C
IF(ICON(4).EQ.3) then
IWRITE=1
call ffprint(iwrite,k,alpha,ub(1),dr(3)
1 ,ltrt(ii),istrt(ii),jtrt(ii),idat)
endif
C
C CHECK ON SPINS AND STATISTICS
C
LTR=LTRT(II)
JTR=JTRT(II)
IF(PHASEF(IS(1)+IS(2)+IS(3)).LT.0.0) GO TO 33
IF(PHASEF( JTR+IS(3)).LT.0.0) GO TO 33
GO TO 35
33 CONTINUE
C
C SET ERROR FLAG
C
IBF(3)=1
WRITE(6,9906)IS,LTR,JTR
35 CONTINUE
c
IF(NLTR.GT.8) IBF(3)=1
C
C CHECK ON COMPATIBILITY OF LMAX AND STORAGE
C
C # Partial waves # Radial matrix elements
LPLUS=MIN0(LPLUS,400/(NS(1)+NS(2)),4000/(NS(1)*NS(2)*(IND+1)))
L=LPLUS-1
LPL2=LPLUS*2
IF(RMAX.LT.0.0) K=KC
IF(ILINE.NE.0) GO TO 40
K=2*(K/2)
L=FK(1)*float(K)*DR(1)+4.0
L=MIN0(LPLUS-1,L)
IF(LC.NE.0) L=MIN0(L,LC)
LPLUS=L+1
LPL2=LPLUS*2
40 WRITE(6,9506)L,K,NLTR
INC=(INCR+IS(1)+IS(2))/2+1
c # Plm's # Fll's
IBUF=MAX0(IBUF, INC*LPLUS+1, LPL2*IND+1)
IBF(1)=IBUF
C
IBF(7)=IS(1)
IF(IBF(5).NE.0) GO TO 51
IS(1)=0
FS(1)=0.0
NS(1)=1
51 CONTINUE
IBF(8)=IS(2)
IF(IBF(6).NE.0) GO TO 52
IS(2)=0
FS(2)=0.0
NS(2)=1
52 CONTINUE
IF(ICON(15).NE.0) THEN
C
C PRINT OUT K(R)**2 FOR DISTORTED WAVES
C
WRITE(6,9999)ALPHA,(IDAT(I),I=1,3)
DO 55 I=1,2
DO 55 N=1,2
IF(N.EQ.1) THEN
WRITE(6,9057)I
ELSE
WRITE(6,9058)I
ENDIF
IND=0
R=DR(I)
DO 54 M=1,K,5
MK=M+M-1
MK4=MIN0(MK +9,K+K)
IF(N.EQ.1) THEN
WRITE(6,9052)R,(U(J,I),J=MK,MK4)
ELSE
WRITE(6,9052)R,(V(J,I),J=MK,MK4)
ENDIF
R=R+5.0*DR(I)
54 CONTINUE
55 CONTINUE
ENDIF
56 CONTINUE
RETURN
C
9001 FORMAT(10F8.4)
9002 FORMAT(18I3)
9008 FORMAT(18H0CARD SET 3 DATA ,9H LMAX =,I4,14H LTR=,8I4)
9009 FORMAT(18H ,9H NLTR =,I4,14H 2*JTR=,8I4)
9010 FORMAT(18H0ANGLE DATA ,9H THETN=,F9.4,9H THET1=,F9.4
1, 9H DTHET=,F9.4,9H A-ANG=,F9.4,9H B-ANG=,F9.4)
9011 FORMAT(18H0CARD SET 4 DATA ,9H DRF =,F9.4,9H RZ =,F9.4
1, 9H RMAX =,F9.4,9H VCE =,F9.4,9H FNRNG=,F9.4)
9012 FORMAT(18X,9H AMASS=,F9.4)
9052 FORMAT(1H ,F6.2,1P10E12.4)
9057 FORMAT(38H0Central K(R)**2 FOR DISTORTED WAVE,I2)
9058 FORMAT(38H0Spin orbit K(R)**2 FOR DISTORTED WAVE,I2)
9500 FORMAT(1P5E16.7)
9506 FORMAT(1H0,17X,9H LMAX =,I4,14H NSTEP=,I4
1, 14H NLTR =,I4)
c 9803 FORMAT(1H ,16I2,4X,15A4,I2,2(1H/,I2.2),I4,2(1H.,I2.2))
9903 FORMAT('0ERROR ** LTR=',I3,' is too large' )
9904 FORMAT('0***** Warning ***** CLASSICAL TURNING POINT FOR BOUND'
1,'STATE =',F8.4,' IS GREATER THAN RMAX =',F8.4)
9905 FORMAT('0***** Warning ***** LMAX too large for reliable '
1,'calculation of unbound stripping or coulomb excitation cases.'
2,' Increase RMAX ')
9906 FORMAT(1H0,28HSPIN STATISTICS NOT CORRECT ,7H 2*IS1=,I3,7H 2*IS2=
1,I3,7H 2*IS3=,I3,7H LTR=,I3,7H 2*JTR=,I3)
9999 FORMAT(1H1,15A4,I4,2(1H/,I2),I4,2(1H.,I2))
END
c***********************************************************************
subroutine ffprint(ii,k,alpha,ub,drf,ltrt,istrt,jtrt,idat)
C
C Prints form factor
c***********************************************************************
c
IMPLICIT REAL*8(A-H,O-Z)
dimension ub(800),idat(6),alpha(15)
c
if(ii.eq.1) then
WRITE(6,9999)ALPHA,(IDAT(I),I=1,3)
endif
R=DRF
WRITE(6,9100)LTRT,ISTRT,JTRT
WRITE(6,9101)
DO 3006 M=1,K,5
MK=M+M-1
MK4=MIN0(MK+9,K+K)
WRITE(6,9052)R ,(UB(N),N=MK,MK4)
R=R+5.0*DRF
3006 CONTINUE
return
9052 FORMAT(1H ,F6.2,1P10E12.4)
9100 FORMAT(12H0FORM FACTOR,6X,9H LTR =,I4,5X,9H 2*STR=,I4,5X
1, 9H 2*JTR=,I4)
9101 FORMAT(55H0 R RL,R IM,R RL,R+DR*1 IM,R+DR*1
1, 48H RL,R+DR*2 IM,R+DR*2 RL,R+DR*3 IM,R+DR*3
2, 24H RL,R+DR*4 IM,R+DR*4 )
9999 FORMAT(1H1,15A4,I4,2(1H/,I2.2),I4,2(1H.,I2.2))
end
c***********************************************************************
SUBROUTINE FORMF(U,V,N,LAM,IK,IVB)
c
c Processes kinematic input for waves and form factors and
c calculates potentials or form factors
c***********************************************************************
c
IMPLICIT REAL*8(A-H,O-Z)
C
COMMON ALPHA(15),IDAT(6),ICON(20),ANGLE(5),HBARC,AMU,AMASS,CHSQ,AA
1,DRF,Z(3),ZA(3),FM(3),FMA(3),RC(3),AC(3),PNLOC(3),FS(3),ECM(3)
2,FK(3),FK2(3),ETA(3),DR(3),FMU(3),FN,FL,FJ2,FSS,VCE,FNRNG,RSIG(2)
3,K,KZ,LPLUS,LPL2,IS(3),NS(3),NLTR,LTRT(8),JTRT(8),ISTRT(8),IBF(8)
4,KC
COMMON/POTTER/DRX,AFACT(2),VFACT,SFACT,E,RM,G(4),ETAX,FKX,ETAKX
1 ,RCX,HBARC2,ABETA(3),FLDF(3)
2 ,NX,LAMX,KMXX,KX,IBX,LPLUSX,ICON4,NSPC,IDIRAC,ICHK
DIMENSION U(810),V(800),X(4)
EQUIVALENCE (X,FN)
C
C READ IN POTENTIAL CARDS FOR CARD SETS 5,6,OR 7
C
ETA6=60.
KMAX=400
ICON4=ICON(4)
NSPC=N
IVB=0
KB=K
READ (5,9000)E,FM(N),Z(N),FMA(N),ZA(N),RY,AC(N),PNLOC(N),FS(N),QCD
E=E+QCD
IF(N.EQ.2) THEN
Q=E
E=(ECM(1)+Q)*(FM(2)+FMA(2))/FMA(2)
ENDIF
IS(N)=FS(N) ! 2*spin, FS(N) will divided by 2 later
NS(N)=IS(N)+1 ! 2*spin+1, number of m-state
IF(AMASS.EQ.0.0) AMASS=FMA(1)
IF(IK.EQ.0) DR(N)=DRF*AMASS/FMA(N)
KMXX=KMAX
DRX=DR(N)
AFACT(1)=FMA(N)**.333333333
AFACT(2)=FM (N)**.333333333
RC(N)=ABS(RY)*AFACT(1)
IF(RY.LT.0.0) RC(N)=RC(N)+ABS(RY)*AFACT(2)
DO 12 M=1,KMAX
MK=M+M-1
U(MK )=0.0
U(MK+1)=0.0
V(MK )=0.0
V(MK+1)=0.0
12 CONTINUE
RM=0.0
IF(E.EQ.0.0) GO TO 66
C
C ICON(10).NE.0 GIVES RELATIVISTIC KINEMATICS
C
IF(ICON(10).NE.0 .or. (fm(n).eq.0.0 .and. is(n).eq.2)) THEN
IF(N.NE.2) GO TO 26
IF(QCD.GT.0.0) GO TO 26
E=E+(ECM(1)+Q)**2/(2.0*FMA(2)*AMU)
26 CONTINUE
FM1=FM(N)*AMU
FM2=FMA(N)*AMU
FMT=FM1+FM2
C KLUDGE FAKE KE-LAB FOR BOUND STATES
IF(N.GE.3) E=((E+FMT)**2-FMT**2)/(2.0*FM2)
WLAB=E+FMT
WCM =SQRT(2.0*E*FM2+FMT**2)
GAMM=WLAB/WCM
W1=GAMM*(FMT*FM1+FM2*E)/WLAB
W2=GAMM* FM2
ECM(N)=WCM-FMT
IF(ICON(10).EQ.1) THEN
C OLD KINEMATICS SELECTED BY ICON(10) = 1
FMU(N)=W1*W2/(W1+W2)
VFACT=2.0*FMU(N)/HBARC**2
SFACT=2.0*fm(n)*fma(n)/(fm(n)+fma(n))/hbarc**2
EFACT=VFACT
ELSE
C NEW KINEMATICS SELECTED BY ICON(10) > 1
FMU(N)=W1
VFACT= 2.0*W1 /HBARC**2
SFACT= 2.0*FM1/HBARC**2
EFACT=VFACT
ENDIF
FMU(N)=FMU(N)/AMU
FK2(N)=(W1**2-FM1**2)/HBARC**2
ELSE
C NON RELATIVISTIC KINEMATICS
FMU(N)=FM(N)*(FMA(N)/(FM(N)+FMA(N)))
IF(N.LT.3) ECM(N)=E*(FMA(N)/(FM(N)+FMA(N)))
IF(N.GE.3) ECM(N)=E
VFACT=2.0*FMU(N)*AMU/HBARC**2
SFACT=VFACT
EFACT=VFACT
FK2(N)=SFACT*ECM(N)
ENDIF
FK(N)=SQRT(ABS(FK2(N)))
ETAK=CHSQ*Z(N)*ZA(N)*EFACT
ETA(N)=ETAK*0.5/FK(N)
HBARC2=HBARC**2
C
C ADD COULOMB AND KINETIC ENERGIES TO U
C
RCX=RC(N)
IF(RCX.EQ.0.0) RCX=DR(N)
R=0.0
FCOU=0.5*ETAK/RCX
DO 42 M=1,kmax
MK=M+M-1
R=R+DR(N)
IF(R.GT.RCX) GO TO 40
FC=FCOU*(3.0-(R/RCX)**2)
GO TO 41
40 CONTINUE
FC=ETAK/R
41 CONTINUE
IF(N.NE.3) U(MK )=U(MK )+FK2(N)-FC
IF(N.EQ.3) U(MK+1)=U(MK+1)+FK2(N)-FC
42 CONTINUE
GO TO 67
66 CONTINUE
FK(N)=0.0
ETA(N)=0.0
FK2(N)=0.0
ECM(N)=0.0
FMU(N)=FM(N)*FMA(N)/(FM(N)+FMA(N))
VFACT=2.0*FMU(N)*AMU/HBARC**2
SFACT=VFACT
ETAK=CHSQ*Z(N)*ZA(N)*VFACT
67 CONTINUE
IF(N.LT.3) GO TO 68
IF(ICON(4).EQ.2) GO TO 69
68 CONTINUE
IF(N.GT.2) THEN
Q=ECM(N)
ELSE
Q=ECM(N)-ECM(1)
ENDIF
WRITE(6,9010)N
WRITE(6,9503)E,RY,AC(N),FS(N)
WRITE(6,9504)FM(N),FMA(N),Q
WRITE(6,9505)Z(N),ZA(N),PNLOC(N)
WRITE(6,9500)
RHO=FK(N)*RC(N)
WRITE(6,9506)ECM(N),RC(N),RHO
WRITE(6,9507)FK(N),ETA(N),DR(N)
WRITE(6,9008)
69 CONTINUE
FS(N)=FS(N)/2.0
C Set up variables for potential routine
IBX=0
ETAX=ETA(N)
FKX=FK(N)
ETAKX=ETAK
RCX=RC(N)
LAMX=LAM
NX=N
LPLUSX=LPLUS
KX=K
ICHK=0
IDIRAC=0
C
CALL POTS(U,V)
C
DR(N)=DRX
LPLUS=LPLUSX
IBF(4)=IBX
C Set nonlocality for Dirac-Darwin term
IF(N.LE.2.AND.IDIRAC.NE.0.AND.PNLOC(N).EQ.0.0) PNLOC(N)=-1.0
K=MIN0(MAX0(K,KX),KMAX)
IF(N.LE.2) GO TO 3000
IF(E.EQ.0.0) THEN
IF(IBF(4).NE.0) u(808) = 1
IF(IVB.EQ.4) THEN
U(801)=U(2*KC-1)*FLOAT(KC)*DR(3)
U(802)=0.0
U(803)=0.0
U(804)=1.0
ENDIF
GO TO 3000
ENDIF
C
C SINGLE PARTICLE ORBITAL
C
C
C READ IN QUANTUM NUMBERS FOR SINGLE PARTICLE ORBITAL
C
2000 CONTINUE
C Set flags for unbound stripping case
IF(E.GT.0.0) THEN
IVB=2
K=KC
ENDIF
READ(5,9000)G,VTRIAL,FISW
FN =G(1)
FL =G(2)
FJ2=G(3)
FSS=G(4)
ISW=FISW
WRITE(6,9500)
FJ0=FJ2/2.0
FS0=FSS/2.0
IF(VTRIAL.EQ.0.0) VTRIAL=ETA6
WRITE(6,9508)G,FISW
FACT=(FJ0**2+FJ0-FL**2-FL-FS0**2-FS0)*0.5
c
DO 2028 M=1,kmax
MK=M+M-1
V(MK )=U(MK )+V(MK )*FACT
V(MK+1)=U(MK+1)+V(MK+1)*FACT
2028 CONTINUE
WRITE(6,9500)
c
CALL BIND(V, U ,DR(3),RM,FN,FL,K,FK(3),ETA(3),VTRIAL,ECM(3)
1,FK2(3),ISW,IBF(3),U(801))
C
IBF(2)=RM/DR(3)
DO 2050 M=1,K
MK=M+M-1
V(M)=VTRIAL*V(MK )+V(MK+1)
2050 CONTINUE
Anorm=1.0
FACT=PNLOC(3)**2/8.0
IF(FACT.ne.0.0) then
C
C NON-LOCAL CORRECTION FOR SINGLE PARTICLE FUNCTION
C
SUM=0.0
R=0.0
DO 2075 M=1,K
MK=M+M-1
R=R+DR(3)
U(M)=U(M)*EXP(FACT*(FK2(3)-V(M)))
SUM=SUM+(U(M)*R)**2
2075 CONTINUE
C
C DO NOT RENORMALIZE FOR POSITIVE ENERGY
C
IF(FK2(3).LT.0.0) then
Anorm=1.0/SQRT(SUM*DR(3))
ELSE
Anorm=1.0
ENDIF
ENDIF
DO 2100 M=kmax,1,-1
MK=M+M-1
IF(M.GT.K) then
V(M )=0.0
U(MK )=0.0
ELSE
U(MK )=U(M)*Anorm
ENDIF
U(MK+1)=0.0
2100 CONTINUE
c
3000 CONTINUE
c Store quantum numbers for transfer back
DO 3020 M=1,4
X(M)=G(M)
3020 CONTINUE
RETURN
9000 FORMAT(10F8.4)
9008 FORMAT(21H0POTENTIAL PARAMETERS )
9010 FORMAT( 9H0PARTICLE,I2,115(1H*))
9500 FORMAT(1H ,3A6,5(3X,A6,F9.4))
9503 FORMAT(18H INPUT DATA ,9H ELAB =,F9.4,9H RC0 =,F9.4
1 ,9H AC =,F9.4,9H 2*STR=,F9.4)
9504 FORMAT(18X ,9H MASSP=,F9.4,9H MASST=,F9.4,9H Q =,F9.4)
9505 FORMAT(18X ,9H ZP =,F9.4,9H ZT =,F9.4,9H PNLOC=,F9.4)
9506 FORMAT(18H DERIVED DATA ,9H ECM =,F9.4,9H RC =,F9.4
1 ,9H RHO =,F9.4)
9507 FORMAT(18X ,9H K =,F9.4,9H ETA =,F9.4,9H DR =,F9.4)
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)
9511 FORMAT(18X ,9H DAMP =,F9.4)
END
c***********************************************************************
SUBROUTINE FNLOC(U,V,W,VB,LTR,BKIN,FNRNG,KT)
c
c Calculates non-local and finite range correction for distorted waves
c
c U = Optical potentials for distorted waves
c V = Spin orbit potentials for distorted waves
c W = Form factor
c VB = Bound State potential
c LTR = Angular momentum transfer ( of bound state)
c BKIN = Kinematic quatities stored in blank common
c FNRNG = Finite range parameter
c KT = Number of points
c***********************************************************************
c
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION U(800,2),V(800,2),W(800),VB(800),BKIN(3,14)
1,RT(8,7),FM(3),FMA(3),RC(3),PNLOC(3),FMU(3),FK(3),FK2(3)
1 ,ETA(3),DR(3)
DATA RT/ 0.0 , 4.4934, 7.7253, 10.9041, 14.0662, 17.2208,
1 20.3713, 23.5195, 2.0816, 5.9404, 9.2058, 12.4044, 15.5792,
2 18.7426, 21.8997, 25.0528, 3.3421, 7.2899, 10.6139, 13.8461,
3 17.0429, 20.2219, 23.3905, 26.5526, 4.5141, 8.5838, 11.9727,
4 15.2445, 18.4681, 21.6666, 24.8501, 28.0239, 5.6467, 9.8404,
5 13.2956, 16.6093, 19.8624, 23.0828, 26.2833, 29.4706, 6.7565,
6 11.0702, 14.5906, 17.9472, 21.2311, 24.4748, 27.6937, 30.8960,
7 7.8511, 12.2793, 15.8632, 19.2627, 22.5781, 25.8461, 29.0843,
8 32.3025/
C
C FINITE RANGE CORRECTION
C
C
C POSITIVE FR PARAMETER CALCULATES A HULTHEN TYPE FR CORRECTION
C NEGATIVE FR PARAMETER CALCULATES A EXPONENTIAL TYPE FR CORRECTION
C
K=KT
DO 20 I=1,3
FM(I) =BKIN(I,3)
FMA(I) =BKIN(I,4)
RC(I) =BKIN(I,5)
PNLOC(I)=BKIN(I,7)
FK(I) =BKIN(I,10)
FK2(I) =BKIN(I,11)
ETA(I) =BKIN(I,12)
DR(I) =BKIN(I,13)
FMU(I) =BKIN(I,14)
20 CONTINUE
IF(FNRNG.EQ.0.0) GO TO 398
IF(FM(1)-FM(2))143,398,142
142 I1=1
I2=2
GO TO 144
143 I1=2
I2=1
C
C I1 IS FOR LARGER MASS , I2 IS FOR SMALLER MASS
C
144 CONTINUE
FM1=FMU(I1)
FM2=FMU(I2)
TEMP=FM(I2)*FNRNG**2/FM(I1)
FACT=TEMP*FMA(I1)/FMA(I2)
FMX=FM(I1)-FM(I2)
TEMP=TEMP*FMX
IF(FMX.LT.1.9) GO TO 350
C
C FORM FACTOR PART BY ROST-KUNZ METHOD
C THIS OPTION CALCULATES AN EXPONENTIAL TYPE FR CORRECTION
C FOR MULTI -NUCLEON TRANSFER REACTIONS
C
DR2=DR(3)**2
FLF=LTR*(LTR+1)
R3=DR(3)*FLOAT(K+1)
D3=0.0
A3=0.0
KM=K+K+1
R2=R3-DR(3)
D2=-FLF/R3
A2=W (KM-2)*R2
DO 220 M=1,K
KM=KM-2
R1=R2-DR(3)
A1=W (KM-2)*R1
D1=-FLF/R2**2
IF(A3.EQ.0.0) GO TO 215
D1=D1+(A3-2.0*A2+A1)/(A2*DR2)
W (KM+2)=W (KM+2)*EXP(D2*FACT)
DX2=(D3*A3-2.0*D2*A2+D1*A1)/(R2*DR2)-A2*D2**2/R2
W (KM+2)=W (KM+2)+0.5*DX2*FACT**2
215 CONTINUE
IF(ABS(W (KM-2)).LT.ABS(W (KM ))) GO TO 225
R3=R2
R2=R1
D3=D2
D2=D1
A3=A2
A2=A1
220 CONTINUE
225 CONTINUE
KM=KM/2+1
DO 227 M=1,KM
MK=M+M-1
VB(M)=W(MK)
W(MK)=0.0
227 CONTINUE
RKM=FLOAT(KM)*DR(3)
LL=MIN0(LTR,6)
FLF=LL*(LL+1)
DR12=DR(3)**2/12.0
DO 250 N=1,8
FKX=(RT(N,LL+1)/RKM)**2
A1=0.0
D1=0.0
R=DR(3)
A2=DR(3)
D2=(1.0+DR12*(FKX-FLF/R**2))*A2
FNORM=0.0
AN=0.0
DO 230 M=1,KM
AN=AN+VB(M)*A2*R
FNORM=FNORM+A2**2
VB(M+400)=A2/R
R=R+DR(3)
D3=1.0+DR12*(FKX-FLF/R**2)
A3=(12.0*A2-10.0*D2-D1)/D3
D1=D2
D2=D3*A3
A1=A2
A2=A3
230 CONTINUE
AN=AN-0.5*VB(KM)*A1*(R-DR(3))
FNORM=FNORM-0.5*A1**2
FNORM=AN*EXP(-FACT*FKX)/FNORM
DO 240 M=1,KM
MK=M+M-1
W (MK)=W (MK)+VB(M+400)*FNORM
240 CONTINUE
250 CONTINUE
340 CONTINUE
C
C DISTORTED WAVE PART, GAUSSIAN FACTOR
C
DO 345 M=1,K
MK=M+M-1
CTEMP1=TEMP*(U(MK ,I2)/FM2-U(MK ,I1)/FM1)
CTEMP2=TEMP*(U(MK+1,I2)/FM2-U(MK+1,I1)/FM1)
FACT=EXP(-CTEMP1)
UT1= FACT*COS(CTEMP2)
UT2=-FACT*SIN(CTEMP2)
UF1 =UT1*W(MK )-UT2*W(MK+1)
UF2 =UT1*W(MK+1)+UT2*W(MK )
W(MK )=UF1
W(MK+1)=UF2
345 CONTINUE
348 CONTINUE
GO TO 398
350 CONTINUE
IF(FNRNG.GT.0.0) GO TO 360
C
C SINGLE PARTICLE PART, GAUSSIAN FACTOR
C
DO 355 M=1,K
MK=M+M-1
FACT=EXP(-TEMP*VB(M)/FMU(3))
W(MK )=W(MK )*FACT
W(MK+1)=W(MK+1)*FACT
355 CONTINUE
GO TO 340
360 CONTINUE
C
C HULTHEN TYPE
C
DO 365 M=1,K
MK=M+M-1
CTEMP1=VB(M)/FMU(3)+U(MK ,I2)/FM2-U(MK ,I1)/FM1
CTEMP2= U(MK+1,I2)/FM2-U(MK+1,I1)/FM1
CTEMP1=1.0+TEMP*CTEMP1
CTEMP2= TEMP*CTEMP2
DET=CTEMP1**2+CTEMP2**2
UT1=(CTEMP1*W(MK )+CTEMP2*W(MK+1))/DET
UT2=(CTEMP1*W(MK+1)-CTEMP2*W(MK ))/DET
W(MK )=UT1
W(MK+1)=UT2
365 CONTINUE
398 CONTINUE
C
C NON LOCALITY CORRECTION FACTOR FOR DISTORTED WAVES
C
DO 410 I=1,2
CALL FNLOC5(U(1,I),V(1,I),W(1),PNLOC(I),FK2(I),ETA(I)
1, FK(I),RC(I),DR(I),K)
410 CONTINUE
RETURN
END
c***********************************************************************
SUBROUTINE CATHEN(LOPT,JOPT,UB,VB,SI)
c
c Two particle form factor subroutine, calculates microscopic
c inelastic or two-nucleon transfer form factors
c***********************************************************************
c
IMPLICIT REAL*8(A-H,O-Z)
COMMON ALPHA(15),IDAT(6),ICON(20),ANGLE(5),HBARC,AMU,AMASS,CHSQ,AA
1,DRF,Z(3),ZA(3),FM(3),FMA(3),RC(3),AC(3),PNLOC(3),FS(3),ECM(3)
2,FK(3),FK2(3),ETA(3),DR(3),FMU(3),FN,FL,FJ2,FSS,VCE,FNRNG,RSIG(2)
3,K,KZ,LPLUS,LPL2,IS(3),NS(3),NLTR,LTRT(8),JTRT(8),ISTRT(8),IBF(8)
4,KC
DIMENSION UB(800,2),VB(801),SI(800),C(10),QNUM(4,2),G(4),IQN(3)
1 ,JI(2),JJ(2),FJT(2),CON(8)
EQUIVALENCE (G,FN),(C(1),CNTROL),(C(2),QCODE),(C(3),FMUV)
1,(C(4),OPT)
DATA SQR4PI,SQRTEN,PI/3.54490780,3.162277660,3.141592/
C
IQFLG=0
KMFLG=0
IQN(1)=LOPT
IQN(3)=JOPT
KT=399
IBF(2)=1
FK2(3)=0.0
100 CONTINUE
READ(5,9001)C
CNTRL=ABS(CNTROL)
IF(CNTRL.EQ.0.0) GO TO 3010
IF(QCODE)2020,3000,1100
1100 CONTINUE
icode=abs(qcode)
C
C LOOP FOR TWO ORBITALS
C
ML=1
DO 2015 I=1,2
FJT(1)=C(5)
FJT(2)=C(I+5)
IF(I.EQ.2.AND.CNTRL.EQ.1.0) GO TO 2005
IK=0
IF(icode.EQ.7.OR.icode.EQ.8) ik=i-1
G(1)=0.0
G(2)=0.0
G(3)=0.0
G(4)=0.0
C
CALL FORMF(UB(1,2),VB,3,LOPT,IK,IVB)
IVB=0
C
DRX=DR(3)
IF(KMFLG.EQ.0) KM=IBF(2)
KMFLG=1
KT=MIN0(KT,K)
2005 CONTINUE
IF(FJT(1).EQ.0.0) FJT(2)=FJ2
C
C PRINT OUT SINGLE PARTICLE INFORMATION
C
if(icode .ne. 5) then
WRITE(6,9501)
WRITE(6,9100)FJ2,FJT
endif
IF(ICON(4).ne.2 .and. (CNTRL.EQ.2.0 .or. I.EQ.1)) then
R=DRX
K2=K+K
DO 2006 M=1,K2,20
MK4=MIN0(M+19,K2)
WRITE(6,9052)R,(UB(N,2),N=M,MK4,2)
R=R+10.0*DRX
2006 CONTINUE
endif
DO 2008 M=1,400
MK=M+M-1
UB(ML,1)=UB(MK,2)
ML=ML+1
2008 CONTINUE
DO 2009 M=1,4
QNUM(M,I)=G(M)
2009 CONTINUE
JI(I)=FJT(1)
JJ(I)=FJT(2)
2015 CONTINUE
C
2020 ICODE=ABS(QCODE)
KT=MIN0(KT,K)
LTR1=IQN(1)+1
IQN(2)=IS(3)
IF(IS(3).EQ.0) IQN(3)=IQN(1)+IQN(1)
WRITE(6,9051)QCODE,FMUV,OPT
C ENTER TIME REVERSAL PHASE
LVR=QNUM(2,1)+QNUM(2,2)
OPT=OPT*PHASEF((LVR+IQN(1))/2)
IFLAG=1
IF(ICODE.LE.10) GO TO 2025
IFLAG=2
ICODE=ICODE-10
2025 CONTINUE
C
GO TO (2100,2100,2300,2990,2500,2100,2700,2700),ICODE
C
C ICODE=1 YUKAWA POTENTIAL
C ICODE=2 COULOMB POTENTIAL
C ICODE=3 OPEP TENSOR POTENTIAL
C ICODE=5 TWO PARTICLE TRANSFER
C ICODE=6 ZERO RANGE KNOCK OUT
C ICODE=7 NO RECOIL FORM FACTOR
C
2100 LVR=0
TOPT=1.0
GO TO 2350
2300 ICODE=1
LVR=2
TOPT=-SQRTEN
2350 MNNL=IABS(LTR1-LVR-1)+1
MXXL=LTR1+LVR
LTR2=LTR1+LTR1-2
LVR2=LVR+LVR
LLX=QNUM(2,1)+QNUM(2,1)
LPX=QNUM(2,2)+QNUM(2,2)
JLX=QNUM(3,1)
JPX=QNUM(3,2)
IS1=QNUM(4,1)
IS2=QNUM(4,2)
OPT=OPT*PHASEF(LPX/2)
DRX=DR(1)
DO 2390 LAM =MNNL,MXXL,2
LAM2 =LAM +LAM -2
VOPT=TOPT*OPT*SQR4PI*SQRT(FLOAT(LAM2 +1)*FLOAT(IQN(2)+1))
1 *VCC(LAM2 ,LVR2,LTR2 ,0,0)
2 *RACAH(IQN(2),IQN(2),LAM2 ,LTR2 ,LVR2,IQN(3))
RME=0.0
IF(IQN(2).EQ.0) RME=1.0
IF(IQN(2).EQ.2) RME=SQRT(FLOAT(IS(1)*(IS(1)+2))
1*FLOAT(IS1*(IS1+2)))
RME=RME*SQRT(FLOAT(IS1+1))*PHASEF((IS(1)-IS(2))/2)
1 *SQRT(FLOAT((LLX+1)*(LAM2+1)))*VCC(LLX,LAM2,LPX,0,0)
2 *SQRT(FLOAT(JLX+1)*FLOAT(JPX+1)*FLOAT(IQN(3)+1))
rme =rme*WINEJ(LPX,LLX,LAM2 ,IS2,IS1,IQN(2),JPX,JLX,IQN(3))
4 *PHASEF((JI(1)+IQN(3)-JLX-JJ(2))/2)
5 *RACAH(JPX,JJ(2),JLX,JJ(1),JI(1),IQN(3))
6 *SQRT(FLOAT(JJ(1)+1)*FLOAT(JJ(2)+1))
VOPT=VOPT*RME*SQRT(FLOAT(LTR2+1)/FLOAT((JJ(1)+1)*(IQN(3)+1)))
1 *(DR(3)/DR(1))**3
SL=0.0
IF(VOPT.EQ.0.0) GO TO 2390
IF(ICODE.EQ.6 ) GO TO 2365
C
KMAX=400
CALL SLATR (KT,KMAX,DRX,VB( 1),LTR1,FMUV,ICODE)
CALL SLATR (KT,KMAX,DRX,VB(801),LAM ,FMUV,ICODE)
C
C
CALL RADIN
1(KT,KMAX,DRX,FMUV,VB,UB(1,1),UB(401,1),SL,VOPT,SI(IFLAG),KM,SK)
C
GO TO 2380
2365 CONTINUE
SCALE=VOPT/SQR4PI**2
R=0.0
MK=IFLAG
DO 2370 M=1,KT
R=R+DRX
TEMP=UB(M ,1)*UB(M+400,1)*SCALE
SI(MK )=SI(MK )+TEMP
SL =SL +TEMP*R**2
IF(M.EQ.KM) SK=TEMP
MK=MK+2
2370 CONTINUE
SL=SL*DRX
2380 CONTINUE
2385 I=LAM -1
WRITE(6,9102)SL,IQN(1),I,RME
WRITE(6,9002)KM,SK
2390 CONTINUE
GO TO 2990
C
C HERE FOR TWO NUCLEON TRANSFER
C
2500 CONTINUE
c
c c1 = R1 scale
c c2 = R2 scale
c c3 = r1 scale
c c4 = r2 scale
c c5 = r0, integration scale length for relative coordinate
c c6 = Pauli flag
c c7 = order of gaussian integration
c c8 = number of integration points
c
IF(FMUV.EQ.0.0) FMUV=1.7
CON(1)= 1.0
CON(2)= 1.0
CON(3)= 0.5
CON(4)=-0.5
CON(5)= 2.0*fmuv
con(6)= 0.0
con(7)= 0.5
con(8)= 0.0
OPT=OPT*16./PI
CALL DSTRIP
1(IQN,DRX,KT,UB(1,1),UB(401,1),SI(IFLAG),QNUM,OPT,KM,SL,CON)
WRITE(6,9002)KM,SL
IQFLG=1
GO TO 2990
2700 CONTINUE
T1= MIN (FM(1),FM(2))
T2= MAX (FM(1),FM(2))
CON(1)= 1.0
CON(2)= 0.0
CON(3)= T1/T2
CON(4)= 1.0
CON(5)= 1.0
CON(7)= 0.0
CON(8)= 0.0
IF(ICODE.EQ.8) GO TO 2706
T3=HBARC**2/(2.0*AMU*FMU(3))
R=0.0
DO 2705 M=1,KT
R=R+DRX
UB(M+400,1)=UB(M+400,1)*T3* (FK2(3)-VB(M))
2705 CONTINUE
2706 CONTINUE
M=(KT+KM)/2+400
T1= LOG(UB(M-1,1)*FLOAT(M-401)/(UB(M ,1)*FLOAT(M-400)))
FMUV = DRX/T1
CON(5)=fmuv
DO 2710 M=1,KT
UB(M+400,1)=UB(M+400,1)*EXP(T1*FLOAT(M))*DRX*FLOAT(M)/100.
2710 CONTINUE
OPT=OPT*FMUV*FMUV*SQR4PI/2.0
CALL DSTRIP
1(IQN,DRX,KT,UB(1,1),UB(401,1),SI(IFLAG),QNUM,OPT,KM,SL,CON)
WRITE(6,9002)KM,SL
IQFLG=1
GO TO 2990
2990 CONTINUE
3000 CONTINUE
IF(CNTROL.GT.0.0) GO TO 100
IF(IQFLG .EQ.1) KT=FLOAT(KT)-1.50/DRX
K=KT
3010 RETURN
9001 FORMAT(10F8.4)
9002 FORMAT(15H FORM FACTOR,M=,I3,1PE18.6)
9051 FORMAT(18H0PARAMETERS ,9H QCODE=,F9.4,9H RANGE=,F9.4
1 ,9H VZERO=,F9.4)
9052 FORMAT(1H ,F6.2,1P10E12.4)
9100 FORMAT(27H0SINGLE PARTICLE FUNCTIONS ,8HCOUPLING ,6H 2*J1=,F3.0
1 ,9H 2*J2=,F3.0,9H 2*JI=,F3.0)
9102 FORMAT(19H VOLUME INTEGRAL = ,F10.4,8H LTR=,I3,8H LAM=,I3
1,8H RME= ,F8.4)
9501 FORMAT(1H ,18HFORM FACTOR DATA )
END