PROGRAM OPTMOD C******************************************************************** C Calculates Wavefunctions, S-matrix, and phase-shifts for a C Woods-Saxon+spin-orbit optical potential C C.A. Bertulani - version 02/Sept/94 C----------------------------------------------------------- C Inputs at OPT.IN C - 1st row: V0R,R0R,AR,FS0,RS0,RC0 c Potential parameters: V0R,R0R,AR - Real part c FS0,RS0 - Spin-orbit c RC0 - Coulomb radius c C - 2nd row: Z1,A1,Z2,A2,EI,EF,DE,L,J c System (Z1,A1,Z2,A2), energy range and step (EI,EF,DE), C angular momentum (L), and total angular momentum (J) c C - 3rd row: DR,NAX c Radial step (DR) and maximum number of steps (NAX) c------------------------------------------------------------ C Units are fm and MeV. Outputs in OPT.OUT and WAVE.OUT C------------------------------------------------------------ IMPLICIT REAL*8(A-H,O-Z) COMMON/SNBAR/ SNBAR,PSHIFT COMMON/DER/AM2,AUX1 COMMON/INTE1/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) REAL*8 PSH(0:500),DPSH(500) c OPEN(11,FILE='OPT.IN',STATUS='OLD') OPEN(12,FILE='OPT.OUT',STATUS='UNKNOWN') OPEN(13,FILE='lixo.OUT',STATUS='UNKNOWN') OPEN(16,FILE='WAVE.OUT',STATUS='UNKNOWN') c C input parameters c READ(11,*)V0R,R0R,AR,FS0,RS0,RC0 READ(11,*)Z1,A1,Z2,A2,EI,EF,DE,L,J READ(11,*)DR,NAX c c general parameters c PI=3.141519 AM0=A1*A2/(A1+A2) A13=A1**.333333+A2**.333333 RC=RC0*A13 RC2=RC*2. AM=931.5*A1*A2/(A1+A2) AM2=2.*AM H2=197.33**2 AUX4=R0R*A13 AUX6=1.44*Z1*Z2 C C******************** Energy mesh ********************** C IE=0. PSH(0)=0. DO 2000 ECM=EI,EF,DE IE=IE+1 PSH(IE)=0. AUX1=AM2*ECM/H2 XK=DSQRT(AUX1) AUX1=AUX1/ECM ETA=1.44*Z1*Z2*AM/XK/H2 C ALL=DBLE(L*(L+1)) AJ=DFLOAT(J)/2.D0 ALS=AJ*(AJ+1.D0)-ALL-0.75D0 CALL WAVE(PSI) PSH(IE)=PSHIFT IF(PSH(IE-1)*PSH(IE).LT.0.) PSH(IE)=PSH(IE)+PI DPSH(IE)=0. WRITE(6,40)ECM,L,J,PSH(IE) WRITE(12,40)ECM,L,J,PSH(IE) C WRITE(13,*)ECM,PSH(IE) 40 FORMAT(' /E, L, J, Ph. Shift/ =', F8.3,1X,I2,1X, & I2,'/2',1X,D12.4) 2000 CONTINUE C C******************** End energy mesh ********************** C C Find position and width of resonance c IEMAX=INT((EF-EI)/DE) PRINT*,'Want position and width of resonance? [1=yes, 0=no]' PRINT*,'If yes, use a small energy step' READ(*,*)IOPT IF(IOPT.EQ.1) THEN C derivative using five-point formula DO 3000 IE=3,IEMAX-2 DPSH(IE)=(PSH(IE-2)-8.*PSH(IE-1)+8.*PSH(IE+1)-PSH(IE+2))/DE/12. WRITE(13,*)EI+IE*DE,DPSH(IE) 3000 CONTINUE C resonance position (maximum of derivative of phase-shift) AA=0. BB=0. DO 4000 IE=2,IEMAX ER=EI+(IE-1)*DE AA=ABS(DPSH(IE-1)) BB=ABS(DPSH(IE)) IF(BB.LT.AA) GOTO 5000 4000 CONTINUE 5000 WIDTH=2./AA WRITE(6,60)ER,WIDTH WRITE(13,60)ER,WIDTH 60 FORMAT(' Resonance energy=',F8.3,' MeV',2X,'Width=',F8.3) ENDIF C printout of wavefunction PRINT*,'Want output of wavefunction? [1=yes, 0=no]' READ(*,*)IOPT IF(IOPT.EQ.1) THEN PRINT*,'If yes, enter value of energy' READ(*,*)ECM ENDIF AUX1=AM2*ECM/H2 XK=DSQRT(AUX1) AUX1=AUX1/ECM ETA=1.44*Z1*Z2*AM/XK/H2 CALL WAVE(PSI) DO I=1,NAX R=I*DR WRITE(16,800)R,REAL(PSI(I)),DIMAG(PSI(I)) 800 FORMAT(3(F8.3,1X)) ENDDO STOP END C C******************************************************************** c Real part of the potential + spin-orbit potential 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 VNR=0. CR=1.E15 TEMP=(R-AUX4)/AR 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 of an extended charge distribution of radius RC 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/INTE1/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 normalizes wavefunction ANORM=.5*(HM-SNUC*HP)*DCMPLX(CEXP(CMPLX(0.,-1.57)))/YY(1) 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