C C.A. Bertulani - version 10/Oct/95 C IMPLICIT REAL*8(A-H,O-Z) REAL*8 V(251),VEE(63001),VCNL(251),HB(63001),FNL(12,251),FH(12) 1 , ENN(1),VR,AR,hb1(63001),indx(300) COMMON /WRIT/ IWG /QQXENN/ IQN COMMON/COUL/GAMMA,CONV,BETA C ********** OPEN THE DATA FILE ******* open (unit=5,file='ctage2.dat',status='old') C EQUIVALENCE(HB,FNL) IQN=0 PI=3.141592653589793D0 C C ********* READ Z1,Z2,A1,A2 *********** 100 READ(5,6) NZ1,NZ2,NA1,NA2 WRITE(6,2) NZ1,NZ2,NA1,NA2 C ****Calculate reduced mass, Coulomb parameter, hbar^2/2 amue, Z1*Z2*e^2 AMUE=DFLOAT(NA1*NA2)/DFLOAT(NA1+NA2) GAMMA=DFLOAT(NZ1*NZ2)*AMUE CONV=20.735D0/AMUE CV=DFLOAT(NZ1*NZ2)*1.4409D0 C C ********** Read parameters C N1: number of meshpoints C IWG: dummy C hstep: mesh size C homega: oscillator parameter (in fm) of internal cluster functions 200 READ(5,3) N1,IWG,HSTEP,HOMEGA C ********** parameters ******** C mstore : dimension of matrix C rest are various combinations of internal cluster parameter BZ=HOMEGA HOMEGA=41.47D0/BZ/BZ MSTORE=N1+1 BETA=BZ/DSQRT(AMUE) BB=BETA*BETA B4=BB*BB HOM=2.D0/BB WRITE(6,4) N1,HSTEP,HOMEGA,BETA C C LOCAL POTENTIALS C Nuclear Potential is assumed to be Gaussian, parameters Vr,Ar C Coulomb potential of homogeneous sphere with radius Rc 300 READ(5,1) VR READ(5,1) AR WRITE(6,92) VR,AR RC=3.55D0 RC2=RC*RC DO 10 I=1,MSTORE X=DFLOAT(I)*HSTEP X2=X*X C C ******** Coulomb potential VCNL(I)=CV/X IF(X.LE.RC) VCNL(I)=CV*0.5D0*(3.D0-X2/RC2)/RC C ******** nuclear potential DEX=DEXP(-X2/(AR*AR)) V(I)=VR*DEX 10 CONTINUE C C C **** print out potential at certain mesh points MPRINT=MAX0(MSTORE/20,1) WRITE(6,8) ((I,V(I),VCNL(I)),I=1,MSTORE,MPRINT) C **** Multiply potential with 2 mu/hbar^2 DO 11 I=1,MSTORE 11 V(I)=(V(I)+VCNL(I))/CONV C C C NON-LOCAL POTENTIAL (from Pauli operator) C *** Read in parameters C L : angular momentum of SCATTERING STATE C NlQ : number of energies C N0 : minimum oscillator energy needed in relative wave function C N0=6 for alpha+12C for even L > 0. C NR : number of Pauli forbidden states in scattering wave fct. C 400 READ(5,6) L,NLQ,N0 NR=(N0+1-L)/2 WRITE(6,91) L,NLQ,N0,NR NIJ=N1*N1 IF(NR.GT.0) GOTO 12 DO 13 IJ=1,NIJ 13 VEE(IJ)=0.D0 GOTO 14 C C RAHAOS calculates the eigenfunctions of the norm kernel C these are radial harmonic oscillator functions C 12 CONTINUE DO 121 I=1,N1 X=DFLOAT(I)*HSTEP CALL RAHAOS(NR,L,BETA,X,FH) DO 121 IL=1,NR 121 FNL(IL,I)=FH(IL) C DO 15 J=1,N1 X=DFLOAT(J)*HSTEP VD=1.5D0*HOM-X*X/B4+V(J) DO 15 I=1,N1 IJ=(J-1)*N1+I VEE(IJ)=0.D0 DO 150 IL=1,NR NT=2*(IL-1) IF (DABS(FNL(IL,I)).LT.1.D-15) FNL(IL,I)=0.D0 IF (DABS(FNL(IL,J)).LT.1.D-15) FNL(IL,J)=0.D0 VEE(IJ)=VEE(IJ)-FNL(IL,I)*FNL(IL,J)*(VD+HOM*DFLOAT(NT)) 150 CONTINUE 15 CONTINUE C C C SOLVE SCHRODINGER EQUATION FOR VARIOUS VALUES OF Q C Q is wave number C 14 I0=0 CALL SOLVM(L,NLQ,N1,HSTEP,V,VEE,MSTORE,HB,hb1,I0,indx) READ(5,5) I GOTO(100,200,300,400,500),I 500 close (unit=5) STOP 1 FORMAT(2D10.4) 2 FORMAT('CHARGE AND MASS NUMBERS : Z1 =',I4,' Z2 =',I4,' 1A1 =',I4,' A2 =',I4) 3 FORMAT(2I5,2D10.4) 4 FORMAT(I5,' STEPS SEPARATED BY',+1PD10.3,10X,'OSCILLATOR PARA 1METER =',D18.10,' MEV =',D22.14,' FM') 5 FORMAT(I1) 6 FORMAT(16I5) 8 FORMAT('0 I LOCAL POTENTIAL COULOMB POTENTIAL'// 1 2(' ',I3,2D24.12/)) 9 FORMAT('0',+1PD10.3,11D11.3/3(12D11.3/)) 91 FORMAT('L =',I4/'NO. OF Q-VALUES',I5,10X,'REDUNDANCY LIMIT N0 1 =',I4,5X,'NO. OF REDUNDANT STATES :',I4/) 92 FORMAT('POTENTIAL PARAMETERS :'/' STRENGTHS',8D15.7/' WIDTHS ', 1 8D15.7) END C************************** C This Subroutine solves the non-local scattering problem C using the Numerov method SUBROUTINE SOLVM(L,NQ,NSTORE,HSTEP,V,VEE,MSTORE,HB,hb1,ISUP,indx) IMPLICIT REAL*8(A-H,O-Z) DIMENSION V(MSTORE),VEE(NSTORE,NSTORE),indx(mstore), 1 HB(MSTORE,MSTORE),HILF(321),G(3000),G0(650),hb1(mstore,mstore) DIMENSION F(31),GG(31),GP(31),FP(31),IEXP(31),SIGMA(31) DIMENSION EIG00(50),EIG22(50),EIG20(50) DIMENSION FH(50),GG0(30,300),GG2(30,300) DIMENSION SUM0(30),SUM2(30) COMMON /WRIT/ IWG /QQXENN/ IQN COMMON/COUL/GAMMA,CONV,BETA C *** Various parameters EPS=0.000000001 PI=4.D0*DATAN(1.D0) DFLL=DFLOAT(L*L+L) HST2=HSTEP*HSTEP ONE12=0.08333333333333333D0*HST2 TEN12=HST2-2.D0*ONE12 KSTORE=NSTORE-1 C Read in relative wave function of O16 ground state C Here the stepsize has been set to 0.04 fm C G0(ix) = rel. wave fct. at ix*0.04 fm DO 7011 IX=1,600 READ (5,*) IZ,G0(IX) 7011 CONTINUE C Read in the decomposition of G0 in terms of radial harm. oscillator C functions; this is needed for treatment of Pauli principle C iz : index of harm. oscillator function C SUM0: overlap of G0 with harm. oscillator fct. C SQ : probability, sq=sum0**2 DO 7012 IX=1,30 READ (5,*) IZ,SUM0(IX),SQ 7012 CONTINUE C Normalization of ground state; C SG0 = int dr G0(r)**2 SG0=0.D0 DO 7016 IX=1,600 SG0=SG0+G0(IX)*G0(IX)*0.04D0 7016 CONTINUE C IEND : number of mesh points for scattering state C > N1 (filled up by asymptotic form) IEND=200 C START ENERGY LOOP FOR SCATTERING STATES C EQ : energy in MeV C Q : wave number C eta : Sommerfeld parameter C XN,XN1: radii for boundary condition C G0(XN) = FL (Q*XN) + TAN (DELTA) * GL (Q*XN) C Same for XN1 = XN + Hstep DO 10 IQ=1,NQ READ(5,7) EQ ENERGY=EQ EQ=EQ/CONV Q=DSQRT(EQ) ETA=GAMMA*1.4409D0/(41.47D0*Q) IEXP(L+1)=0 XN=DFLOAT(NSTORE)*HSTEP*Q C DFCOUL calculates reg. (F) and irreg. (GG) Coulomb functions C upto angular momentum L C FP,GP are the derivatives with respect to Q*x C Sigma is Coulomb phase shift C Iexp error parameter CALL RCWFN(XN,ETA,0,L,F,FP,GG,GP,IEXP,500.D0) C CALL DFCOUL(L,ETA,XN,F,FP,GG,GP,IEXP,SIGMA) FN=F(L+1) GN=GG(L+1) XN1=XN+HSTEP*Q CALL RCWFN(XN,ETA,0,L,F,FP,GG,GP,IEXP,500.D0) C CALL DFCOUL(L,ETA,XN1,F,FP,GG,GP,IEXP,SIGMA) FN1=F(L+1) GN1=GG(L+1) C SOLVE PROBLEM DO 11 J=1,NSTORE DO 15 I=1,NSTORE VIJ=VEE(I,J) 15 HILF(I)=VIJ*HSTEP HB(1,J)=-HILF(1)*TEN12-HILF(2)*ONE12 DO 16 I=2,KSTORE 16 HB(I,J)=-(HILF(I-1)+HILF(I+1))*ONE12 - HILF(I)*TEN12 HB(NSTORE,J)=-HILF(NSTORE-1)*ONE12-HILF(NSTORE)*TEN12 X=DFLOAT(J+ISUP)*HSTEP HB(J,J)=HB(J,J)+(EQ-V(J+ISUP)-DFLL/(X*X))*TEN12-2.D0 11 CONTINUE DO 12 I=1,MSTORE X=DFLOAT(I+ISUP)*HSTEP 12 HILF(I)=EQ-V(I+ISUP)-DFLL/(X*X) DO 13 I=2,NSTORE HB(I-1,MSTORE)=0.D0 HB(MSTORE,I-1)=0.D0 HB(I-1,I)=HB(I-1,I)+1.D0+ONE12*HILF(I) HB(I,I-1)=HB(I,I-1)+1.D0+ONE12*HILF(I-1) 13 CONTINUE HB(NSTORE,MSTORE)=GN1 +GN1*ONE12*HILF(MSTORE) HB(MSTORE,NSTORE)=-1.D0 HB(MSTORE,MSTORE)=GN NST2=MSTORE*MSTORE do i=1,mstore do j=1,mstore hb1(i,j)=hb(i,j) enddo enddo C Invert hamilton matrix call dminv(mstore,hb1,hb,indx) RSNST=-FN1-FN1*ONE12*HILF(MSTORE) RSMST=-FN NODE=-1 DS=0.D0 DO 14 I=1,NSTORE G(I)=HB(I,NSTORE)*RSNST+HB(I,MSTORE)*RSMST C Calculate nodes in scattering state IF(DSIGN(1.D0,G(I)).NE.DS) NODE=NODE+1 DS=DSIGN(1.D0,G(I)) 14 CONTINUE C Calculate BL (= Tan (DELTA)) BL=HB(MSTORE,NSTORE)*RSNST+HB(MSTORE,MSTORE)*RSMST IF(IWG.EQ.1) WRITE(6,2) DETH,(G(I),I=1,NSTORE) BL=HB(MSTORE,NSTORE)*RSNST+HB(MSTORE,MSTORE)*RSMST C BUILT up SCATTERING STATE IN ASYMPTOTIC REGION C G = FL(q*x) + Tan(DELTA) * GL(q*x) DO 8734 IX=NSTORE,IEND X=DFLOAT(IX)*Q*HSTEP CALL RCWFN(XN,ETA,0,L,F,FP,GG,GP,IEXP,500.D0) C CALL DFCOUL(L,ETA,X,F,FP,GG,GP,IEXP,SIGMA) G(IX)=F(L+1)+BL*GG(L+1) 8734 CONTINUE C Calculate phase shift Delta (in degrees and Pi's) DEL=DATAN(BL) DELPI=DEL/3.141592653589793D0 C Write ENERGY AND PHASE SHIFT ON FILE write (66,*) energy,del C WRITE to main output WRITE(6,3) L,ENERGY,DEL,DELPI,NODE C *** Define eigenvalues of many-body alpha+C12 norm kernel C eig00 : < L=0 | L=0 > C eig22 : < L=2 | L=2 > C eig20 : < L=2 | L=0 > DO 100 I=1,50 EIG00(I)=1.D0 EIG22(I)=1.D0 EIG20(I)=0.D0 100 CONTINUE EIG00(1)=0.D0 EIG00(1)=0.D0 EIG00(3)=0.296296D0 EIG00(4)=0.316049D0 EIG00(5)=0.5614D0 EIG00(6)=0.740537D0 EIG00(7)=0.856463D0 EIG00(8)=0.924746D0 EIG00(9)=0.961918D0 EIG00(10)=0.9811456D0 EIG00(11)=0.990794D0 EIG00(12)=0.995547D0 EIG22(1)=0.D0 EIG22(2)=0.D0 EIG22(3)=0.23045D0 EIG22(4)=0.5421937D0 EIG22(5)=0.736879D0 EIG22(6)=0.85584D0 EIG22(7)=0.924646D0 EIG22(8)=0.961903D0 EIG22(9)=0.9811434D0 EIG22(10)=0.9908D0 EIG20(1)=0.D0 EIG20(2)=0.D0 EIG20(3)=0.66254D0 EIG20(4)=0.470266D0 EIG20(5)=0.24522D0 EIG20(6)=0.135905D0 EIG20(7)=0.07536D0 EIG20(8)=0.040143D0 EIG20(9)=0.020606D0 EIG20(10)=0.010311D0 SUM00=0.D0 DO 500 I=1,30 SUM00=SUM00+SUM0(I)**2 500 CONTINUE C Calculate eigenfunctions of norm kernel C for ground state (L=0) : GG0 C for scattering (L=2) : GG2 C stepsize hstep, number of meshpoints iend AMUE=48.D0/16.D0 BETA=1.7D0/DSQRT(AMUE) DO 21 I=1,IEND X=DFLOAT(I)*HSTEP CALL RAHAOS(20,0,BETA,X,FH) DO 22 N=1,20 GG0(N,I)=FH(N) 22 CONTINUE CALL RAHAOS(20,2,BETA,X,FH) DO 23 N=1,20 GG2(N,I)=FH(N) 23 CONTINUE 21 CONTINUE C *** DECOMPOSE scattering state in norm eigen functions C *** Only leading terms are required DO 7650 IX=1,20 SUM2(IX)=0.D0 DO 7651 IXX=1,IEND SUM2(IX)=SUM2(IX)+GG2(IX,IXX)*G(IXX)*HSTEP 7651 CONTINUE 7650 CONTINUE C Calculate DIRECT PART of E2 transition C integral dx x^2 G(x) G0(x) C note different step size for G and G0. SUM=0.D0 DO 7013 IX=1,IEND X=DFLOAT(IX)*HSTEP SUM=SUM+G0(3*IX)*G(IX)*X*X 7013 CONTINUE C Calculate many-body contributions FACTOR=DSQRT(45.D0/16.D0/PI) Q2=SUM*FACTOR*HSTEP/DSQRT(SG0) C N0: Sum over oscillator states C three possibilities for rank 2 operator C n2 = n the ground and scattering states C n2 = n-1 have different L and number C n2 = n-2 of Pauli forbidden states C DO 200 N0=3,15 DO 200 IX=1,3 N2=N0+1-IX C Only N2>2 is allowed by Pauli principle IF (N2.LT.3) GOTO 200 SUMMY=0.D0 C many-body contributions C do relative part first C Integral dx x^2 GG0(x) GG2(x) DO 300 I=1,IEND X=DFLOAT(I)*HSTEP QQX=X*X*GG0(N0,I)*GG2(N2,I) SUMMY=SUMMY+QQX*HSTEP 300 CONTINUE AMUENL=DSQRT(EIG00(N0)/EIG22(N2)) IF (N2+1.LT.N0) AMUENL=1.D0/AMUENL FACTOR=DSQRT(45.D0/16.D0/PI) QQ2=SUMMY*(AMUENL-1.D0)*FACTOR C Contribution from internal C12 cluster DSQ=DSQRT(EIG00(N0)*EIG22(N2)*5.D0) IF (N2+1.EQ.N0) QQ2=QQ2-6.07D0*EIG20(N0)/DSQ Q2=Q2+QQ2*SUM0(N0)*SUM2(N2)/DSQRT(SUM00) 200 CONTINUE C BE2 is the many-body transition BE2=Q2*Q2/5.D0 C Include energy factor for E2 EGAM=ENERGY+7.162D0 EGA=EGAM**5 C Transform into transition probabilty PROB=1.223D9*EGA C Normalize the scattering wave function to unit flux C consider that asymptotic is actually C G = COS (DELTA) * FL + SIN (DELTA) * GL QE2=BE2*20.D0*PI*(DCOS(DEL)/Q)**2 VREL=6.D23*DSQRT(ENERGY/(3.D0*(939.507D0+938.213D0))) C Cross section in barn CROSS=PROB*QE2/VREL/100.D0 WRITE (6,432) ENERGY,CROSS C Define astrophysical S-factor (in MeV b) SF=DLOG(ENERGY)+DLOG(CROSS)+650.51D0/DSQRT(1.D3*ENERGY) ASTRO=DEXP(SF) WRITE (6,434) ASTRO write (6,*) C Write S-factor to file write (77,*) energy,astro 10 CONTINUE RETURN 1 FORMAT(4D20.12) 2 FORMAT('0DETERMINANT :',+1PD12.4/8(12D11.3/)) 3 FORMAT(' L =',I3,5X,'E =',D20.12,5X,'PHASE =',D20.12,' =',D20.12, 1' * PI',5X,I3,' NODES') 7 FORMAT(1D20.4) 432 FORMAT(1X,' ENERGY = ',D20.12,' CROSS-SECTION IN BARN =',D20.12) 434 FORMAT(1X,' ASTROPHYSICAL S-FACTOR = ',D20.12) END C ************ C Subroutine for matrix inversion subroutine dminv(np,a,y,indx) implicit real*8 (a-h,o-z) dimension a(np,np),y(np,np),indx(np) n=np do 12 i=1,n do 11 j=1,n y(i,j)=0. 11 continue y(i,i)=1. 12 continue call ludcmp(a,n,np,indx,d) do 13 j=1,n call lubksb(a,n,np,indx,y(1,j)) 13 continue return end SUBROUTINE LUBKSB(A,N,NP,INDX,B) implicit real*8 (a-h,o-z) DIMENSION A(NP,NP),INDX(N),B(N) II=0 DO 12 I=1,N LL=INDX(I) SUM=B(LL) B(LL)=B(I) IF (II.NE.0)THEN DO 11 J=II,I-1 SUM=SUM-A(I,J)*B(J) 11 CONTINUE ELSE IF (SUM.NE.0.) THEN II=I ENDIF B(I)=SUM 12 CONTINUE DO 14 I=N,1,-1 SUM=B(I) IF(I.LT.N)THEN DO 13 J=I+1,N SUM=SUM-A(I,J)*B(J) 13 CONTINUE ENDIF B(I)=SUM/A(I,I) 14 CONTINUE RETURN END SUBROUTINE LUDCMP(A,N,NP,INDX,D) PARAMETER (NMAX=100,TINY=1.0E-20) implicit real*8 (a-h,o-z) DIMENSION A(NP,NP),INDX(N),VV(NMAX) D=1. DO 12 I=1,N AAMAX=0. DO 11 J=1,N IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J)) 11 CONTINUE IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.' VV(I)=1./AAMAX 12 CONTINUE DO 19 J=1,N IF (J.GT.1) THEN DO 14 I=1,J-1 SUM=A(I,J) IF (I.GT.1)THEN DO 13 K=1,I-1 SUM=SUM-A(I,K)*A(K,J) 13 CONTINUE A(I,J)=SUM ENDIF 14 CONTINUE ENDIF AAMAX=0. DO 16 I=J,N SUM=A(I,J) IF (J.GT.1)THEN DO 15 K=1,J-1 SUM=SUM-A(I,K)*A(K,J) 15 CONTINUE A(I,J)=SUM ENDIF DUM=VV(I)*ABS(SUM) IF (DUM.GE.AAMAX) THEN IMAX=I AAMAX=DUM ENDIF 16 CONTINUE IF (J.NE.IMAX)THEN DO 17 K=1,N DUM=A(IMAX,K) A(IMAX,K)=A(J,K) A(J,K)=DUM 17 CONTINUE D=-D VV(IMAX)=VV(J) ENDIF INDX(J)=IMAX IF(J.NE.N)THEN IF(A(J,J).EQ.0.)A(J,J)=TINY DUM=1./A(J,J) DO 18 I=J+1,N A(I,J)=A(I,J)*DUM 18 CONTINUE ENDIF 19 CONTINUE IF(A(N,N).EQ.0.)A(N,N)=TINY RETURN END SUBROUTINE RAHAOS(N,L,BETA,X,FL) C C THIS SUBROUTINE CALCULATES THE RADIAL HARMONIC OSCILLATOR WAVE C FUNCTIONS FOR ONE VALUE OF ANGULAR MOMENTUM AND FOR RADIAL QUANTUM C NUMBERS FROM 0 TO N-1. C C EXPLANATION OF PARAMETERS : C C N - NO. OF FUNCTIONS EVALUATED (QU.NO. 0 TO N-1) C L - ANGULAR MOMENTUM C BETA - OSCILLATOR WIDTH C X - ARGUMENT OF WAVE FUNCTIONS C FL - ARRAY STORING FUNCTION VALUES AFTER CALCULATION C C METHOD C C THE ZERO-NODE AND THE ONE-NODE WAVE FUNCTIONS ARE CALCULATED VIA THE C EXPLICIT FORMULA C G(0,L) = (4PI)**0.5*(PI*BETA**2)**(-0.75)*(2**L/(2L+1)!!)**0.5 C * EXP(-X**2/(2*BETA**2)) * X**(L+1)/BETA**L C G(1,L) = G(0,L) * 2/(2L+3) * (L+3/2-X**2/BETA**2) C C THE HIGHER FUNCTIONS ARE CALCULATED VIA A RECURSION RELATION DEDUCED C FROM THE RECURSION RELATION FOR LAGUERRE POLYNOMIALS C C IMPLICIT REAL*8(A-H,O-Z) DIMENSION FL(N),IL(72) DATA PI/3.141592653589793D0/ IF(BETA.LE.1.D-6*X) GOTO 10 B2=BETA*BETA XB=X/BETA XB2=XB*XB XB22=0.5D0*XB2 R01=DFLOAT(L)+1.5D0 R02=DFLOAT(2*L+3) IF(0.25D0*XB2.GT.R02+DFLOAT(4*N)) GOTO 10 DLXB=DLOG(XB) FACTOR=DSQRT(4.D0*PI)*(PI*B2)**(-0.75D0) TERM=DLOG(X)-XB22 IL(1)=0 IF(L.EQ.0) GOTO 11 DL2=DLOG(2.D0) LF=1 DO 12 LL=1,L LF=LF+2 12 TERM=TERM+0.5D0*(DL2-DLOG(DFLOAT(LF))) + DLXB IF(XB2.LT.R01) GOTO 11 13 IF(TERM.GT.-20.D0) GOTO 11 TERM=TERM+20.D0 IL(1)=IL(1)-20 GOTO 13 11 FL(1)=FACTOR*DEXP(TERM) C IF(N.EQ.1) GOTO 14 C F1=0.5D0*R02 FL(2)=FL(1)*(F1-XB2)*DSQRT(2.D0/R02) IL(2)=IL(1) C IF(N.EQ.2) GOTO 14 C DO 15 I=3,N IN=I-2 IN2=IN+IN IN1=I-1 FA1=DFLOAT(IN2)+F1 FA2=FA1+F1 FA3=FA2-2.D0 IL(I)=IL(I-1) RFL=DFLOAT(IN1)*FA2 FL(I) = (FA1-XB2)*DSQRT(2.D0/RFL) * FL(I-1) * - DSQRT(DFLOAT(IN)*FA3/RFL) * FL(I-2) 16 IF(FL(I).LE.1.D10) GOTO 15 FL(I)=FL(I)*1.D-10 FL(I-1)=FL(I-1)*1.D-10 IL(I)=IL(I)+10 IL(I-1)=IL(I-1)+10 GOTO 16 15 CONTINUE C 14 DO 17 I=2,N 17 FL(I)=FL(I)*10.D0**DFLOAT(IL(I)-IL(1))*DEXP(DFLOAT(IL(1))) FL(1)=FL(1)*DEXP(DFLOAT(IL(1))) RETURN 10 DO 18 I=1,N 18 FL(I)=0.D0 RETURN END c 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