PROGRAM B8MATRIX C******************************************************************** C Calculates electromagnetic matrix elements of Boron-8 C Ground-state to continuum transitions C Written by C.A. Bertulani -Last revised: 12/12/95 C******************************************************************** C IMPLICIT REAL*8(A-H,O-Z) COMMON/SNBAR/ SNBAR,PSHIFT COMMON/DER/AM2,AUX1 COMMON/INTE/DR,RI COMMON/INTE2/NAX COMMON/ET/ETA,XK COMMON/POT/V0R,R0R,AR,FS0,RS0,AM,A13,ECM COMMON/SPIN/L,J,ALL,ALS COMMON/COU/Z1,Z2,A1,A2,RC,RC2,AUX6 COMMON/AUX/AUX4,AUX5 COMPLEX*16 PSI(1600),SUM,E1S,E1D3,E1D5,E2P,E2P3,E2F5,E2F7, & E2PNR,M1P1,M1P3,PSIR1(1600),PSIR3(1600) REAL*8 PSIGS(1600),RR(1600),GINT(1600) c OPEN(11,FILE='GS.IN',STATUS='OLD') OPEN(12,FILE='DBEB8.OUT',STATUS='UNKNOWN') OPEN(13,FILE='CAPB8.OUT',STATUS='UNKNOWN') OPEN(14,FILE='SB8.OUT',STATUS='UNKNOWN') c Potential parameters: V0R,R0R,AR - Real part c FS0,RS0 - Spin-orbit c RC0 - Coulomb radius C c potential parameters AR=0.52 R0R=0.8208 FS0=0.351 RS0=1.25 RC0=R0R c Depth for non-resonant transitions V0RN=-42.48 c Depth for 1+ resonance C V0R1=V0RN c Depth for 3+ resonance V0R3=-36.05 c c System (Z1,A1,Z2,A2), Z1=1. A1=1. Z2=4. A2=7. E2=1.44 c energy range and step (EI,EF,DE) EI=0.0005 EF=3. DE=0.03 c Radial step (DR) and maximum number of steps (NAX) DR=0.14 NAX=1599 IF(MOD(NAX,2).EQ.0) NAX=NAX+1 c enter ground-state wavefunction DO 200 K=1,NAX READ(11,*)RR(K),PSIGS(K) 200 CONTINUE c general parameters PI=3.141519 PI2=2.*PI AM0=A1*A2/(A1+A2) A13=A1**.333333+A2**.333333 RC=RC0*A13 RC2=RC*2. AM=931.5*AM0 AM2=2.*AM H2=197.33**2 AUX4=R0R*A13 AUX6=1.44*Z1*Z2 EEFF=(Z2/A2-Z1/A1)*AM0 EEFF2=(Z2/A2**2+Z1/A1**2)*AM0**2 EEFFM=(Z2/A2+Z1/A1)*AM0 GPROT=2.79 H2MC=SQRT(H2)/2./938. AMH2=AM/H2 c repeated indices for the clebsh and racah coefficients AL0=1. AJ0=1.5 AJ0M=0.5 ALAM1=1. ALAM2=2. ALAMM=0. AJM=0.5 AL0HAT=SQRT(2.*AL0+1.) AJ0HAT=SQRT(2.*AJ0+1.) AJ0J0=SQRT(AJ0**2+AJ0) c common factors for strength functions and cross sections SQ1=SQRT((2.*ALAM1+1)*(2.*AJ0+1)/4./PI) SQ2=SQRT((2.*ALAM2+1)*(2.*AJ0+1)/4./PI) FACE1=3. FACE2=15. FPHOT1=8.*PI**3*(ALAM1+1.)/ALAM1/FACE1**2/SQRT(H2)*E2 FPHOT2=8.*PI**3*(ALAM2+1.)/ALAM2/FACE2**2/SQRT(H2)**3*E2 FM1=H2MC*EEFFM*2.*AJ0J0/AL0HAT**2*AL0 IDUM=NINT(AL0+0.5-AJ0) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(l0+0.5-j0) FM2=SIGN*GPROT*H2MC/AL0HAT**2*AJ0J0 c spins of 8B, 7Be, and proton AJCORE=1.5 AJPART=0.5 AJB8=2. FCAP=1./H2*2.*(2.*AJB8+1.)/(2.*AJPART+1.)/(2.*AJCORE+1.) c coefficients used in Simpson's integrations GINT(1)=DR/3. GINT(NAX)=DR/3. RR(1)=DR RR(NAX)=NAX*DR G=4. DO 201 K=2,NAX-1 RR(K)=K*DR GINT(K)=G*DR/3. G=6.-G 201 CONTINUE C C******************** Energy mesh ********************** C IE=0. DO 2000 ECM=EI,EF,DE IE=IE+1 AUX1=AM2*ECM/H2 XK=DSQRT(AUX1) AUX1=AUX1/ECM ETA=1.44*Z1*Z2*AM/XK/H2 C C matrix elements for s,p,d, and f waves C C******************************************************************** C ****** E1 E1 E1 E1 E1 E1 E1 E1 ******** C******************************************************************** C ****** s and d waves (E1 transitions) ****** V0R=V0RN C C ***** s wave C L=0 J=1 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=0, j=1/2, continuum (Simpson) SUM=0. DO 3000 K=1,NAX SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K) 3000 CONTINUE c C Single particle reduced matrix elements (E1: g.s. --> s-wave) c AL=0. AJ=0.5 FACT=EEFF*CLEBSH(AJ0,ALAM1,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM1) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E1S=((1.+SIGN)/2.)*SQ1*FACT*SUM IDUM=NINT(AL0+AL+AJ0-AJ) ! (-1)**(l0+l+j0-j) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) E1S=SIGN*E1S C C ***** d wave (j=3/2) C L=2 J=3 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=2, j=3/2, continuum (Simpson) SUM=0. DO 3001 K=1,NAX SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K) 3001 CONTINUE c C Single particle reduced matrix elements (E1: g.s. --> d-wave, j=3/2) c AL=2. AJ=1.5 FACT=EEFF*CLEBSH(AJ0,ALAM1,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM1) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E1D3=(1.+SIGN)/2.*SQ1*FACT*SUM IDUM=NINT(AL0+AL+AJ0-AJ) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l0+l+j0-j) E1D3=SIGN*E1D3 C C ***** d wave (j=5/2) C L=2 J=5 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=2, j=5/2, continuum (Simpson) SUM=0. DO 3002 K=1,NAX SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K) 3002 CONTINUE c C Single particle reduced matrix elements (E1: g.s. --> d-wave, j=5/2) c AL=2. AJ=2.5 FACT2=EEFF*CLEBSH(AJ0,ALAM1,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM1) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E1D5=(1.+SIGN)/2.*SQ1*FACT2*SUM IDUM=NINT(AL0+AL+AJ0-AJ) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l0+l+j0-j) E1D5=SIGN*E1D5 C C E1 Strength function C DBE1S=AMH2*ABS(E1S/AJ0HAT)**2/XK DBE1D3=AMH2*ABS(E1D3/AJ0HAT)**2/XK DBE1D5=AMH2*ABS(E1D5/AJ0HAT)**2/XK DBE1D=DBE1D3+DBE1D5 DBE1=DBE1S+DBE1D C C Photo absorption cross section PHOTE1=FPHOT1*(ECM+0.14)*DBE1 PHTE1S=FPHOT1*(ECM+0.14)*DBE1S PHTE1D=FPHOT1*(ECM+0.14)*DBE1D C Radiative capture cross section CAPE1=(ECM+0.14)**2*FCAP/XK**2*PHOTE1 CAPE1S=(ECM+0.14)**2*FCAP/XK**2*PHTE1S CAPE1D=(ECM+0.14)**2*FCAP/XK**2*PHTE1D C S-factor SE1=ECM*CAPE1*EXP(PI2*ETA) SE1S=ECM*CAPE1S*EXP(PI2*ETA) SE1D=ECM*CAPE1D*EXP(PI2*ETA) C C******************************************************************** C ****** E2 E2 E2 E2 E2 E2 E2 E2 ******** C******************************************************************** C ****** p and f waves (E2 transitions) ****** C C ***** p(3/2) wave (resonant, at 0.63 keV) C V0R=V0RN C L=1 J=3 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=1, j=3/2 , J=1+, resonance (Simpson) SUM=0. DO 3003 K=1,NAX PSIR1(K)=PSI(K) SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K)**2 3003 CONTINUE c C Single particle reduced matrix elements (E2: g.s. --> p-wave, at 0.63 MeV) c AL=1. AJ=1.5 FACT2=EEFF2*CLEBSH(AJ0,ALAM2,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E2P=((1.+SIGN)/2.)*SQ2*FACT2*SUM IDUM=NINT(AL0+AL+AJ0-AJ) ! (-1)**(l0+l+j0-j) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) E2P=SIGN*E2P C C ***** p(3/2) wave (resonant, at 2.17 MeV) C V0R=V0R3 C L=1 J=3 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=1, j=3/2, J=3+ continuum (Simpson) SUM=0. DO 3004 K=1,NAX PSIR3(K)=PSI(K) SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K)**2 3004 CONTINUE c C Single particle reduced matrix elements (E2: g.s. --> p-wave, at 2.17 MeV) c AL=1. AJ=1.5 FACT2=EEFF2*CLEBSH(AJ0,ALAM2,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E2P3=((1.+SIGN)/2.)*SQ2*FACT2*SUM IDUM=NINT(AL0+AL+AJ0-AJ) ! (-1)**(l0+l+j0-j) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) E2P3=SIGN*E2P3 C C ***** p(1/2) wave (non-resonant) C V0R=V0RN C L=1 J=1 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=1, j=3/2, J=3+ continuum (Simpson) SUM=0. DO 3014 K=1,NAX PSIR3(K)=PSI(K) SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K)**2 3014 CONTINUE c C Single particle reduced matrix elements (E2: g.s. --> (p1/2)-wave) c AL=1. AJ=0.5 FACT2=EEFF2*CLEBSH(AJ0,ALAM2,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E2PNR=((1.+SIGN)/2.)*SQ2*FACT2*SUM IDUM=NINT(AL0+AL+AJ0-AJ) ! (-1)**(l0+l+j0-j) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) E2PNR=SIGN*E2PNR C C ***** f wave (j=5/2) C V0R=V0RN L=3 J=5 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=3, j=5/2, continuum (Simpson) SUM=0. DO 3005 K=1,NAX SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K)**2 3005 CONTINUE c C Single particle reduced matrix elements (E2: g.s. --> f-wave, j=5/2) c AL=3. AJ=2.5 FACT2=EEFF2*CLEBSH(AJ0,ALAM2,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E2F5=(1.+SIGN)/2.*SQ2*FACT2*SUM IDUM=NINT(AL0+AL+AJ0-AJ) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l0+l+j0-j) E2F5=SIGN*E2F5 C C ***** f wave (j=7/2) C L=3 J=7 ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) C radial integral for g.s. --> l=3, j=7/2, continuum (Simpson) SUM=0. DO 3006 K=1,NAX SUM=SUM+GINT(K)*PSI(K)*PSIGS(K)*RR(K)**2 3006 CONTINUE c C Single particle reduced matrix elements (E2: g.s. --> f-wave, j=7/2) c AL=3. AJ=3.5 FACT2=EEFF2*CLEBSH(AJ0,ALAM2,AJ,AJ0M,ALAMM,AJM) IDUM=NINT(AL0+AL+ALAM2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l+l0+lambda) E2F7=(1.+SIGN)/2.*SQ2*FACT2*SUM IDUM=NINT(AL0+AL+AJ0-AJ) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(l0+l+j0-j) E2F7=SIGN*E2F7 C C E2 Strength function C C p(3/2) waves at 0.63 keV and 2.17 MeV SUM=0. DO 3500 J=1,3,2 AJ1=1.5 AJ2=DFLOAT(J) AJ3=1.5 AL1=2. AL2=1.5 AL3=2. IDUM=NINT(AJ1+AJ2+AL1+AL2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(j1+j2+l1+l2) SIXJ=SIGN*RACAHR(AJ1,AJ2,AL2,AL1,AJ3,AL3) IF (J.EQ.1) THEN SUM=SUM+(2.*AJ2+1.)*SIXJ**2*ABS(E2P)**2 ELSE SUM=SUM+(2.*AJ2+1.)*SIXJ**2*ABS(E2P3)**2 ENDIF 3500 CONTINUE DBE2P=AMH2*SUM/XK C C add non-resonant p1/2 continuum C DBE2P=DBE2P+AMH2*ABS(E2PNR/AJ0HAT)**2/XK C C f(5/2) and f(7/2) waves C DBE2F5=AMH2*ABS(E2F5/AJ0HAT)**2/XK DBE2F7=AMH2*ABS(E2F7/AJ0HAT)**2/XK DBE2F=DBE2F5+DBE2F7 DBE2=DBE2P+DBE2F C C C Photo absorption cross section PHOTE2=FPHOT2*(ECM+0.14)**3*DBE2 C Radiative capture cross section CAPE2=(ECM+0.14)**2*FCAP/XK**2*PHOTE2 C S-factor SE2=ECM*CAPE2*EXP(PI2*ETA) c C C M1 matrix elements for p waves C C******************************************************************** C ****** M1 M1 M1 M1 M1 M1 M1 M1 ******** C******************************************************************** C ****** p waveS (M1 transitions) ****** C C ***** p(3/2) wave (resonant, at 0.63 keV) C C wavefunction was calculated and stored above C C radial integral for g.s. --> l=1, j=3/2 , J=1+, continuum (Simpson) SUM=0. DO 3007 K=1,NAX SUM=SUM+GINT(K)*PSIR1(K)*PSIGS(K) 3007 CONTINUE c C Single particle reduced matrix elements (M1: g.s. --> p-wave, at 0.63 MeV) c AL=1. AJ=1.5 AJ1=AJ AJ2=1. AJ3=1.5 AL1=2. AL2=AJ0 AL3=1. IDUM=NINT(AJ1+AJ2+AL1+AL2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(j1+j2+l1+l2) SIXJ=SIGN*RACAHR(AJ1,AJ2,AL2,AL1,AJ3,AL3) FACT1=SQRT((2.*AJ2+1.)*(2.*AL1+1))*SIXJ IDUM=NINT(AJ1+AJ3+AL1+1.) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(j1+j3+l1+1) FACT1=SIGN*FACT1 FACT1=FACT1*(FM1+FM2) C AJ1=1.5 AJ2=1. AJ3=AJ0 AL1=2. AL2=1.5 AL3=1. IDUM=NINT(AJ1+AJ2+AL1+AL2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(j1+j2+l1+l2) SIXJ=SIGN*RACAHR(AJ1,AJ2,AL2,AL1,AJ3,AL3) IDUM=NINT(AJ1+AJ3+AJ2+1.) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(j1+j2+j3+1) FACT2=SQRT((2.*AL1+1.)*(2.*AJ2+1.)*(2*AL2+1.)*(AJ1**2+AJ1))* & SIXJ*(-1.7*H2MC)*SIGN M1P1=(FACT1+FACT2)*SUM C C ***** p(3/2) wave (resonant, at 2.17 MeV) C C wavefunction was calculated and stored above C C radial integral for g.s. --> l=1, j=3/2, J=3+ continuum (Simpson) SUM=0. DO 3008 K=1,NAX SUM=SUM+GINT(K)*PSIR3(K)*PSIGS(K) 3008 CONTINUE c C Single particle reduced matrix elements (M1: g.s. --> p-wave, at 2.17 MeV) c AL=1. AJ=1.5 AJ1=AJ AJ2=3. AJ3=1.5 AL1=2. AL2=AJ0 AL3=1. IDUM=NINT(AJ1+AJ2+AL1+AL2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(j1+j2+l1+l2) SIXJ=SIGN*RACAHR(AJ1,AJ2,AL2,AL1,AJ3,AL3) FACT1=SQRT((2.*AJ2+1.)*(2.*AL1+1))*SIXJ IDUM=NINT(AJ1+AJ3+AL1+1.) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(j1+j3+l1+1) FACT1=SIGN*FACT1 IDUM=NINT(AL0+0.5-AJ) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(l0+0.5-j) FACT1=FACT1*(FM1+FM2) C AJ1=1.5 AJ2=3. AJ3=AJ0 AL1=2. AL2=1.5 AL3=1. IDUM=NINT(AJ1+AJ2+AL1+AL2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) !(-1)**(j1+j2+l1+l2) SIXJ=SIGN*RACAHR(AJ1,AJ2,AL2,AL1,AJ3,AL3) IDUM=NINT(AJ1+AJ3+AJ2) SIGN=DFLOAT(1-2*IABS(IDUM-2*(IDUM/2))) ! (-1)**(j1+j2+j3) FACT2=SQRT((2.*AL1+1.)*(2.*AJ2+1.)*(2*AL2+1.)*(AJ1**2+AJ1))* & SIXJ*(-1.7*H2MC)*SIGN M1P3=(FACT1+FACT2)*SUM C C C p(3/2),(J=1+) and p(3/2),(J=3+) resonances C DBM1P1=AMH2*ABS(M1P1/AJ0HAT)**2/XK DBM1P3=AMH2*ABS(M1P3/AJ0HAT)**2/XK DBM1=DBM1P1+DBM1P3 C C C Photo absorption cross section PHOTM1=FPHOT1*(ECM+0.14)*DBM1 C Radiative capture cross section CAPM1=(ECM+0.14)**2*FCAP/XK**2*PHOTM1 C S-factor SM1=ECM*CAPM1*EXP(PI2*ETA) C WRITE(6,40)ECM,SE1*1.D4,SE2*1.D4,SM1*1.D4 WRITE(12,40)ECM,DBE1S,DBE1,DBE2F,DBE2 WRITE(13,40)ECM,CAPE1S*10,CAPE1D*10,CAPE1*10,CAPE2*10,CAPM1*10 WRITE(14,40)ECM,SE1S*1.D4,SE1D*1.D4,SE1*1.D4,SE2*1.D4,SM1*1.D4 40 FORMAT(F8.3,1X,D12.4,1X,D12.4,1X,D12.4,1X,D12.4,1X,D12.4) C40 FORMAT(' / E, Strength, S-E1 (eV.b) / =', C & F8.3,1X,D12.4,1X,D12.4,1X,D12.4,1X,D12.4) 2000 CONTINUE C C******************** End energy mesh ********************** C PRINT*,'Output:' PRINT*,' DBEB8.OUT - Strength function for E1 (e^2 fm^2 MeV^-1)' PRINT*,' and E2 (e^2 fm^4 MeV^-1):' PRINT*,' s (2nd column), and s+d waves (third)' PRINT*,' f (4nd column), and f+p waves (fifth)' PRINT*,' CAPB8.OUT - Capt. cross sections (mb): E1S (2nd col.),' PRINT*,' E1D (3rd), E1 (4th), E2 (5th) and M1 (6th)' PRINT*,' SB8.OUT - S-factors (eV.b) for E1S (2nd), E1D (3rd)' PRINT*,' E1D (3rd), E1 (4th), E2 (5th) and M1 (6th)' STOP END C C******************************************************************** c Real part of the potential + spin-orbit c FUNCTION VNR(R) IMPLICIT REAL*8(A-H,O-Z) COMMON/AUX/AUX4,AUX5 COMMON/POT/V0R,R0R,AR,FS0,RS0,AM,A13,ECM COMMON/SPIN/L,J,ALL,ALS TEMP=(R-AUX4)/AR VNR=0. IF(TEMP.LT.30)CR=DEXP(TEMP) VNR=V0R/(1.+CR) c add spin-orbit term DVNR=CR/(1.+CR)**2/AR IF(R.GT.1.E-30) VS0=V0R*FS0*ALS*RS0/2./R*DVNR VNR=VNR+VS0 RETURN END C C******************************************************************** c Coulomb potential c FUNCTION VCOU(R) IMPLICIT REAL*8(A-H,O-Z) COMMON/COU/Z1,Z2,A1,A2,RC,RC2,AUX6 IF(R.GT.RC)GO TO 1 VCOU=AUX6*(3.-(R/RC)**2)/RC2 RETURN 1 VCOU=AUX6/R RETURN END C C******************************************************************** c Local momentum function c FUNCTION DERIV(R) IMPLICIT REAL*8(A-H,O-Z) REAL*8 DERIV COMMON/DER/AM2,AUX1 COMMON/COU/Z1,Z2,A1,A2,RC,RC2,AUX6 COMMON/POT/V0R,R0R,AR,FS0,RS0,AM,A13,ECM COMMON/SPIN/L,J,ALL,ALS DERIV=ALL/(R*R)-AUX1*(ECM-VCOU(R)-VNR(R)) RETURN END C C******************************************************************** c Solves Schroedinger equation using a Runge-Kutta method c SUBROUTINE WAVE(PSI) IMPLICIT REAL*8(A-H,O-Z) COMMON/SNBAR/SNBAR,PSHIFT COMMON/COU/Z1,Z2,A1,A2,RC,RC2,AUX6 COMMON/INTE/DR,RI COMMON/INTE2/NAX COMMON/ET/ETA,XK COMMON/POT/V0R,R0R,AR,FS0,RS0,AM,A13,ECM COMMON/SPIN/L,J,ALL,ALS DIMENSION F(71),FP(71),G(71),GP(71),CNN(1600) COMPLEX*16 HP,HM,DHP,DHM,DLG,ANORM,SNUC COMPLEX*16 PSI(1600) COMPLEX*16 YY(2) C DR1=DR/10. N=NAX+1 PSI(1)=DCMPLX(0.D0,0.D0) YY(1)=DCMPLX(0.D0,0.D0) YY(2)=DCMPLX(1.D0,0.D0) R=1.D-10 CNN(1)=0.D0 c c begin integration up to r=RMAX c DO 1000 I=2,N CNN(I)=CNN(I-1) CALL RK4(YY,2,R,DR) CHN=ABS(YY(1)) c Renormalization c IF(CHN.GT.1.D30)THEN CNN(I)=CNN(I)+DLOG(CHN) YY(1)=YY(1)/CHN YY(2)=YY(2)/CHN END IF c PSI(I)=YY(1) 1000 CONTINUE c c End integration and match with Coulomb waves at r=RMAX c RMAX=DR*DBLE(NAX) LP=L+1 ROM=RMAX*XK ACCUR=1.D0 STEP=100.D0 CALL RCWFN(ROM,ETA,L,L,F,FP,G,GP,ACCUR,STEP) DEVI=R-RMAX IF (DABS(DEVI).GT.DR1) WRITE(6,707)R,RMAX,DEVI 707 FORMAT(/,' DEU ZEBRA!! RMATCH (=',f8.4,' fm) DIFERENTE 1 DE RMAX (=',f8.4,' fm) !! Delta R=',g10.2//) c c S-matrix c HP=DCMPLX(G(LP),F(LP)) HM=DCONJG(HP) DHP=DCMPLX(GP(LP),FP(LP)) DHM=DCONJG(DHP) DLG=YY(2)/YY(1) SNUC=(XK*DHM-DLG*HM)/(XK*DHP-DLG*HP) c normalize wavefunction: sin(kr+ p.shift) as r --> infty c ANORM=.5*(HM-SNUC*HP)*DCMPLX(CEXP(CMPLX(0.,-1.57)))/YY(1) c normalize wavefunction: (2/pi)**{1/2}*sin(kr+ p.shift) as r --> infty ANORM=0.5*(HM-SNUC*HP)*DCMPLX(CEXP(CMPLX(0.,-1.57)))/YY(1) ANORM=SQRT(2./3.141519)*ANORM DO 2000 I=2,N RENO=EXP(CNN(N)-CNN(I)) PSI(I)=PSI(I)*ANORM/RENO 2000 CONTINUE c c phase shift c SNBAR=CDABS(SNUC) PSHIFT=DATAN2(DIMAG(SNUC),DBLE(REAL(SNUC)))/2. RETURN END C C******************************************************************** C Wavefunction (Y1) and its derivative at X (Y2) c SUBROUTINE DERIVS(X,Y,DYDX) IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 Y(2),DYDX(2) DYDX(1)=Y(2) DYDX(2)=Y(1)*DERIV(X) RETURN END C C******************************************************************** c Runge-Kutta: iterates to get Y(X_(i+1)) from Y(X_i) c SUBROUTINE RK4(Y,N,X,H) PARAMETER (NMAX=10) IMPLICIT REAL*8(A-H,O-Z) COMPLEX*16 Y,DYDX,YT,DYT,DYM DIMENSION Y(N),DYDX(NMAX),YT(NMAX),DYT(NMAX),DYM(NMAX) HH=H*.5 H6=H/6. XH=X+HH c CALL DERIVS(X,Y,DYDX) DO 11 I=1,N YT(I)=Y(I)+HH*DYDX(I) 11 CONTINUE CALL DERIVS(XH,YT,DYT) DO 12 I=1,N YT(I)=Y(I)+HH*DYT(I) 12 CONTINUE CALL DERIVS(XH,YT,DYM) DO 13 I=1,N YT(I)=Y(I)+H*DYM(I) DYM(I)=DYT(I)+DYM(I) 13 CONTINUE CALL DERIVS(X+H,YT,DYT) DO 14 I=1,N Y(I)=Y(I)+H6*(DYDX(I)+DYT(I)+2.*DYM(I)) 14 CONTINUE X=X+H RETURN END C C******************************************************************** c Coulomb waves SUBROUTINE RCWFN(RHO,ETA,MINL,MAXL,FC,FCP,GC,GCP,ACCUR,STEP) IMPLICIT REAL*8(A-H,O-Z) REAL*4 K,K1,K2,K3,K4,M1,M2,M3,M4 DIMENSION FC(71),FCP(71),GC(71),GCP(71) C*** COULOMB WAVEFUNCTIONS CALCULATED AT R=RHO BY THE C*** CONTINUED FRACTION METHOD OF STEED MINL,MAXL ARE ACTUAL L-VALUES C*** SEE BARNETT FENG STEED AND GOLDFARB COMPUTER PHYSICS COMMUN 1974 PACE=STEP ACC=ACCUR IF(PACE.LT.100.0)PACE=100.0 IF(ACC.LT.1.0E-15.OR.ACC.GT.1.0E-6)ACC=1.0E-6 R=RHO KTR=1 LMAX=MAXL LMIN1=MINL+1 XLL1=DBLE(MINL*LMIN1) ETA2=ETA*ETA TURN=ETA+DSQRT(ETA2+XLL1) IF(R.LT.TURN.AND.DABS(ETA).GE.1.0E-6)KTR=-1 KTRP=KTR GO TO 2 1 R=TURN TF=F TFP=FP LMAX=MINL KTRP=1 2 ETAR=ETA*R RHO2=R*R PL=DBLE(LMAX+1) PMX=PL+0.5 C*** CONTINUED FRACTION FOR FP(MAXL)/F(MAXL) XL IS F XLPRIME IS FP FP=ETA/PL+PL/R DK=ETAR*2.0 DEL=0.0 D=0.0 F=1.0 K=(PL*PL-PL+ETAR)*(2.0*PL-1.0) IF(PL*PL+PL+ETAR.NE.0.0)GO TO 3 R=R+1.0E-6 GO TO 2 3 H=(PL*PL+ETA2)*(1.0-PL*PL)*RHO2 K=K+DK+PL*PL*6.0 D=1.0/(D*H+K) DEL=DEL*(D*K-1.0) IF(PL.LT.PMX)DEL=-R*(PL*PL+ETA2)*(PL+1.0)*D/PL PL=PL+1.0 FP=FP+DEL IF(D.LT.0.0)F=-F IF(PL.GT.20000.)GO TO 11 IF(DABS(DEL/FP).GE.ACC)GO TO 3 FP=F*FP IF(LMAX.EQ.MINL)GO TO 5 FC(LMAX+1)=F FCP(LMAX+1)=FP C*** DOWNWARD RECURSION TO MINL FOR F AND FP,ARRAYS GC,GCP ARE STORAGE L=LMAX DO 4 LP=LMIN1,LMAX PL=DBLE(L) GC(L+1)=ETA/PL+PL/R GCP(L+1)=DSQRT(ETA2+PL*PL)/PL FC(L)=(GC(L+1)*FC(L+1)+FCP(L+1))/GCP(L+1) FCP(L)=GC(L+1)*FC(L)-GCP(L+1)*FC(L+1) 4 L=L-1 F=FC(LMIN1) FP=FCP(LMIN1) 5 IF(KTRP.EQ.-1)GO TO 1 C*** REPEAT FOR R=TURN IF RHO LT TURN C*** NOW OBTAIN P+I.Q FOR MINL FROM CONTINUED FRACTION (32) C*** REAL ARITHMETIC TO FACILITATE CONVERSION TO IBM USING REAL*8 P=0.0 Q=R-ETA PL=0.0 AR=-(ETA2+XLL1) AI=ETA BR=2.0*Q BI=2.0 WI=2.0*ETA DR=BR/(BR*BR+BI*BI) DI=-BI/(BR*BR+BI*BI) DP=-(AR*DI+AI*DR) DQ=(AR*DR-AI*DI) 6 P=P+DP Q=Q+DQ PL=PL+2.0 AR=AR+PL AI=AI+WI BI=BI+2.0 D=AR*DR-AI*DI+BR DI=AI*DR+AR*DI+BI T=1.0/(D*D+DI*DI) DR=T*D DI=-T*DI H=BR*DR-BI*DI-1.0 K=BI*DR+BR*DI T=DP*H-DQ*K DQ=DP*K+DQ*H DP=T IF(PL.GT.46000.)GO TO 11 IF(DABS(DP)+DABS(DQ).GE.(DABS(P)+DABS(Q))*ACC)GO TO 6 P=P/R Q=Q/R C*** SOLVE FOR FP,G,GP AND NORMALISE F AT L=MINL G=(FP-F*P)/Q GP =P*G-Q*F W=1.0/DSQRT(FP*G-F*GP) G=W*G GP=W*GP IF(KTR.EQ.1)GO TO 8 F=TF FP=TFP LMAX=MAXL C*** RUNGE-KUTTA INTEGRATION OF G(MINL) AND GP(MINL) INWARDS FROM TURN C*** SEE FOX AND MAYERS 1968 PG 202 IF(RHO.LT.0.2*TURN)PACE=999.0 R3=1.0/3.0 H=(RHO-TURN)/(PACE+1.0) H2=0.5*H I2=(PACE+0.001) ETAH=ETA*H H2LL=H2*XLL1 S=(ETAH+H2LL/R)/R-H2 7 RH2=R+H2 T=(ETAH+H2LL/RH2)/RH2-H2 K1=H2*GP M1=S*G K2=H2*(GP+M1) M2=T*(G+K1) K3=H*(GP+M2) M3=T*(G+K2) M3=M3+M3 K4=H2*(GP+M3) RH=R+H S=(ETAH+H2LL/RH)/RH-H2 M4=S*(G+K3) G=G+(K1+K2+K2+K3+K4)*R3 GP=GP+(M1+M2+M2+M3+M4)*R3 R=RH I2=I2-1 IF(I2.GE.0)GO TO 7 W=1.0/(FP*G-F*GP) C*** UPWARD RECURSION FROM GC(MINL) AND GCP(MINL),STORED VALUES ARE R,S C*** RENORMALISE FC,FCP FOR EACH L-VALUE 8 GC(LMIN1)=G GCP(LMIN1)=GP IF(LMAX.EQ.MINL)GO TO10 DO 9 L=LMIN1,LMAX T=GC(L+1) GC(L+1)=(GC(L)*GC(L+1)-GCP(L))/GCP(L+1) GCP(L+1)=GC(L)*GCP(L+1)-GC(L+1)*T FC(L+1)=W*FC(L+1) 9 FCP(L+1)=W*FCP(L+1) FC(LMIN1)=FC(LMIN1)*W FCP(LMIN1)=FCP(LMIN1)*W RETURN 10 FC(LMIN1)=W*F FCP(LMIN1)=W*FP RETURN 11 W=0.0 G=0.0 GP=0.0 GO TO 8 END C C******************************************************************** C C******************************************************************** FUNCTION RACAH(JAD,JBD,JCD,JDD,JED,JFD) IMPLICIT REAL*8(A-H,O-Z) C C calculates racah coefficients C C for six-j coefficients, use relation C sixj(j1,j2,j3,l1,l2,l3)=(-1)**(j1+j2+l1+l2)*W(j1,j2,l2,l1,j3,l3) C C ENTRIES : RACAH , RACAHI , RACAHR C RACAH : INTEGER ARGUMENTS = 2*J C RACAHI : ARGUMENTS = TRUE INTEGER VALUE C RACAHR : ARGUMENTS = TRUE REAL VALUE C EXTERNAL : THFACI , generates factorial table C DIMENSION I(16) LOGICAL FIRST C REAL*8 G,S PARAMETER (LFACTC=200) COMMON / LOGFAC1 / LFACT COMMON / LOGFAC2 / G(LFACTC) EQUIVALENCE(I(1),I1),(I(2),I2),(I(3),I3),(I(4),I4),(I(5),I5), 1 (I(6),I6),(I(7),I7),(I(8),I8),(I(9),I9),(I(10),I10),(I(11),I11), 2 (I(12),I12),(I(13),I13),(I(14),I14),(I(15),I15),(I(16),I16) DATA FIRST /.TRUE./ C MAKE USEFULL COMBINATIONS K=JAD+JBD-JED+2 I1=K/2 IF((2*I1).NE.K) GOTO 300 K=JCD+JDD-JED+2 I4=K/2 IF((2*I4).NE.K) GOTO 300 K=JAD+JCD-JFD+2 I7=K/2 IF((2*I7).NE.K) GOTO 300 K=JBD+JDD-JFD+2 I10=K/2 IF((2*I10).NE.K) GOTO 300 I13=I1+JED I14=I4+JED I15=I7+JFD I16=I10+JFD I2=I13-JAD I3=I13-JBD I5=I14-JCD I6=I14-JDD I8=I15-JAD I9=I15-JCD I11=I16-JBD I12=I16-JDD C CHECK TRIANGULAR INEQUALITIES,FIND NO. OF TERMS IN SUM N=MIN(I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12)-1 IF(N) 300,2,2 C FIND MINIMUM VALUE OF SUMMATION INDEX 2 IL=MAX(I13,I14,I15,I16) IF(MIN(JAD,JBD,JCD,JDD,JED,JFD)) 300,20,1 C .............. ENTRY RACAHI(JA1,JB1,JC1,JD1,JE1,JF1) C MAKE USEFULL COMBINATIONS I13=JA1+JB1+JE1+1 I14=JC1+JD1+JE1+1 I15=JA1+JC1+JF1+1 I16=JB1+JD1+JF1+1 I1=I13-JE1*2 I2=I13-JA1*2 I3=I13-JB1*2 I4=I14-JE1*2 I5=I14-JC1*2 I6=I14-JD1*2 I7=I15-JF1*2 I8=I15-JA1*2 I9=I15-JC1*2 I10=I16-JF1*2 I11=I16-JB1*2 I12=I16-JD1*2 C CHECK TRIANGULAR INEQUALITIES,FIND NO. OF TERMS IN SUM N=MIN(I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12)-1 IF(N) 300,4,4 C FIND MINIMUM VALUE OF SUMMATION INDEX 4 IL=MAX(I13,I14,I15,I16) LMIN=MIN(JA1,JB1,JC1,JD1,JE1,JF1) IF(LMIN)300,20,1 C ............ ENTRY RACAHR(A,B,C,D,E,F) C CONVERT ARGUMENTS TO INTEGER JA=NINT(2.*A) JB=NINT(2.*B) JC=NINT(2.*C) JD=NINT(2.*D) JE=NINT(2.*E) JF=NINT(2.*F) C MAKE USEFULL COMBINATIONS K=JA+JB-JE+2 I1=K/2 IF((2*I1-K).NE.0) GOTO 300 K=JC+JD-JE+2 I4=K/2 IF((2*I4-K).NE.0) GOTO 300 K=JA+JC-JF+2 I7=K/2 IF((2*I7-K).NE.0) GOTO 300 K=JB+JD-JF+2 I10=K/2 IF((2*I10-K).NE.0) GOTO 300 I13=I1+JE I14=I4+JE I15=I7+JF I16=I10+JF I2=I13-JA I3=I13-JB I5=I14-JC I6=I14-JD I8=I15-JA I9=I15-JC I11=I16-JB I12=I16-JD C CHECK TRIANGULAR INEQUALITIES,FIND NO. OF TERMS IN SUM N=MIN(I1,I2,I3,I4,I5,I6,I7,I8,I9,I10,I11,I12)-1 IF(N) 300,3,3 C FIND MINIMUM VALUE OF SUMMATION INDEX 3 IL=MAX(I13,I14,I15,I16) LMIN=MIN(JA,JB,JC,JD,JE,JF) IF(LMIN)300,20,1 C ------------ 1 IF(FIRST) GOTO 8 9 IF(IL.GE.LFACT) STOP 'RACAH: LENGTH FACTORIAL TABLE INSUFFICIENT' J1=IL-I13+1 J2=IL-I14+1 J3=IL-I15+1 J4=IL-I16+1 J5=I13+I4-IL J6=I15+I5-IL J7=I16+I6-IL PH=1. IF(2*(J5/2).EQ.J5) PH=-1. H=PH*EXP ((G(I1)+G(I2)+G(I3)-G(I13+1)+G(I4)+G(I5)+G(I6)- 1G(I14+1)+G(I7)+G(I8)+G(I9)-G(I15+1)+G(I10)+G(I11)+G(I12)-G(I16+1)) 2*.5+G(IL+1)-G(J1)-G(J2)-G(J3)-G(J4)-G(J5)-G(J6)-G(J7)) IF(N)300,110,120 C 110 RACAH=H RACAHI=H RACAHR=H RETURN C 120 S=1. K=N-1 KL=IL+1 J5=J5-1 J6=J6-1 J7=J7-1 DO 130 J=1,N C! K=N-J S=1.-((KL+K)*(J5-K)*(J6-K)*(J7-K))*S/((J1+K)*(J2+K)*(J3+K)*(J4+K)) K=K-1 130 CONTINUE RACAH=H*S RACAHI=H*S RACAHR=H*S RETURN C C ONE OF THE ARGUMENTS =0 20 IAD=IL IBD=IL DO 21 J=13,16 IF(IAD.LT.I(J)) GOTO 22 IF(IAD.LT.IBD) IBD=IAD IAD=I(J) GOTO 21 22 IF(IBD.GT.I(J)) IBD=I(J) 21 CONTINUE J5=I13+I4-IL PH=1. IF(2*(J5/2).EQ.J5) PH=-1. RACAH=PH/SQRT(FLOAT(IAD*IBD)) RACAHI=RACAH RACAHR=RACAH RETURN C C IMPOSSIBLE COMBINATION OF ARGUMENTS 300 RACAH=0. RACAHI=0. RACAHR=0. RETURN 8 CALL THFACI FIRST=.FALSE. GOTO 9 END c C******************************************************************** C SUBROUTINE THFACI IMPLICIT REAL*8(A-H,O-Z) C ... Set up log of factorial table. Used in calculation of Racah coefficients PARAMETER (LFACTC=200) COMMON / LOGFAC1 / LFACT COMMON / LOGFAC2 / FACLOG(LFACTC) C REAL*8 FACLOG,FN LFACT=LFACTC FACLOG(1)=0.0 FACLOG(2)=0.0 FN=1.0 DO 10 I=3,LFACTC FN=FN+1.0 FACLOG(I)=FACLOG(I-1)+LOG(FN) 10 CONTINUE RETURN END c c*********************************************************************** c Returns Wigner 3j-coefficients (for Clebsh coeficients change c name to clebsh and delete line b3=-b3) c a1,a2,a3 are j1,j2,j3 and b1,b2,b3 are m1,m2,m3. c array f(n+1)=n!/a**n, where a is an arbitrary number c that cancels in the calculation of clebsh; its purpose c is to prevent under/overflows in array f c c*********************************************************************** c c real*8 function wigner(a1,a2,a3,b1,b2,b3) real*8 function clebsh(a1,a2,a3,b1,b2,b3) c function clebsh(a1,a2,a3,b1,b2,b3) c implicit real*8 (a-h,o-z) dimension f(85) save f data ifirst /1/ c if (ifirst .eq. 1) then f(1)=1.0d0 do 10 i=1,84 10 f(i+1)=f(i)*dfloat(i)/15.d0 ifirst=0 endif c c delete this line to run this program for clebsh coefficients c b3=-b3 c clebsh=0.d0 c wigner=0.d0 if(abs(b1)-a1 .gt. 0.01d0) goto 100 if(abs(b2)-a2 .gt. 0.01d0) goto 100 if(abs(b3)-a3 .gt. 0.01d0) goto 100 if(abs(b1+b2-b3).gt.0.01d0) goto 100 if(a3-a1-a2 .gt. 0.01d0) goto 100 if(abs(a1-a2)-a3 .gt. 0.01d0) goto 100 l1=a1+a2-a3+1.01d0 l2=a1+b1+1.01d0 l3=a2+b2+1.01d0 l4=a2-b2+1.01d0 l5=a3+a2-a1+1.01d0 l6=a1-b1+1.01d0 l7=a3+b3+1.01d0 l8=a3-b3+1.01d0 l9=a3+a1-a2+1.01d0 l10=a3+a1+a2+1.01d0 m1=l1 m2=l6 m3=l3 m4=(l7+l8-l3-l4+l2-l6)/2+1 m5=(l7+l8-l2-l6-l3+l4)/2+1 m6=1 zw=sqrt(f(l1)*f(l2)*f(l6)) tw=sqrt(f(l9)/f(l10)*(2.0d0*a3+1.d0)/float(l10)) sw=tw*sqrt(f(l3)*f(l7)*f(l8)*f(l4)*f(l5)) min=max0(-m4,-m5,-m6)+2 max=min0(m1,m2,m3) fr=-(-1.)**min if(max.lt.min) return do m=min,max mw=m-1 clebsh=clebsh+sw*fr/f(m4+mw)/(f(m2-mw)/zw)/f(m3-mw)/ 1 f(m1-mw)/f(m5+mw)/f(m6+mw) fr=-fr enddo c here returns wigner coeeficients c wigner=(-1)**nint(a1-a2+b3)/sqrt(2.d0*a3+1.d0)*clebsh 100 return end