PROGRAM KAP13 *--------------------------------------------------------------------- * Program KAP13 solves the radial SCHR™DINGER equation for the * energy E and the angular momentum quantum number l in the * interval [0,b] by the FOX-GOODWIN method, with the à-à system * as example. * --Fixed quantities-- * H2M: hý/(8ãým) in units of 1MeVfmý * (h Planck's constant, m (reduced) mass, * in the à-à example H2M=10.375D0) * {DOUBLE PRECISION constant} * IMAX: Dimension of the arrays R,U,W (IMAX=2000) * {INTEGER constant} * BMAX: Upper limit for interval boundaries b (BMAX=30.D0) * {DOUBLE PRECISION constant} * --Input quantities-- * L: Angular momentum quantum number l * {INTEGER variable} * E: Energy in units of 1MEV * {DOUBLE PRECISION variable} * IB: Number of integration steps in the interval [0,b] * (IB ó IMAX) * {INTEGER variable} * B: Upper integration limit in units of 1fm * {DOUBLE PRECISION variable} * --Output quantities-- * U: Radial wave function u(r), unnormalised, * U(I) is the function value at the point R(I) * {DOUBLE PRECISION array of dimension (0:IMAX)} * --Other quantities-- * R: Distance r in units of 1fm * R(I)=I*H, H=B/IB, I=0,1,2,...,IB * {DOUBLE PRECISION array of dimension (0:IMAX)} * V: Potential V in units of 1MeV * {DOUBLE PRECISION FUNCTION} * W: Auxiliary array * {DOUBLE PRECISION array of dimension (0:IMAX)} * --Subprograms called-- * FXGOOD, GOPEN, GMOVE, GDRAW, GCLOSE, CURSOR {SUBROUTINES} * V {DOUBLE PRECISION FUNCTION} *--------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (IMAX=2000, H2M=10.375D0, BMAX=30.D0) DIMENSION R(0:IMAX), U(0:IMAX), W(0:IMAX) * >>>>>Input/output CHARACTER CONT*1 PARAMETER (ILINE=25, IOMAX=105) OPEN(6,FILE=' ') CALL CURSOR(1) WRITE(*,'(T29,A)') 'ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿' WRITE(*,'(T29,A)') '³ C H A P T E R 1 3 ³' WRITE(*,'(T29,A/)') 'ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ' WRITE(*,'(T5,A/8(T2/))') & 'SOLUTION OF THE RADIAL SCHR™DINGER EQUATION WITH THE FOX GOODWIN & METHOD' WRITE(*,'(T2,A)') &'Please input values for the following quantities:' IFLAG=1 1000 CONTINUE WRITE(*,'(T2/T2,A)') &'Angular momentum quantum number L (0 ó L ó 10): L = ' READ(*,*,IOSTAT=IOS) L IF((L.LT.0).OR.(L.GT.10).OR.(IOS.NE.0)) THEN WRITE(*,'(T2,A/)') 'Input error !' GOTO 1000 ENDIF 2000 CONTINUE WRITE(*,'(T2,A,F5.1,A)') & 'Upper integration limit B in units of 1fm (0.0 < B ó ', & BMAX,'): B = ' READ(*,*,IOSTAT=IOS) B IF((B.LE.0.D0).OR.(B.GT.BMAX).OR.(IOS.NE.0)) THEN WRITE(*,'(T2,A/)') 'Input error !' GOTO 2000 ENDIF 3000 CONTINUE WRITE(*,'(T2,A)') 'Energy E in units of 1MeV: E = ' READ(*,*,IOSTAT=IOS) E IF(IOS.NE.0) THEN WRITE(*,'(T2,A/)') 'Input error !' GOTO 3000 ENDIF 4000 CONTINUE WRITE(*,'(T2,A,I5,A)') &'Number IB of integration steps (2 ó IB ó ',IMAX,'): IB = ' READ(*,*,IOSTAT=IOS) IB IF((IB.LT.2).OR.(IB.GT.IMAX).OR.(IOS.NE.0)) THEN WRITE(*,'(T2,A/)') 'Input error !' GOTO 4000 ENDIF * <<<<< H=B/IB R(0)=0.D0 W(0)=0.D0 DO 10 I=1,IB R(I)=I*H W(I)=(E-V(R(I)))/H2M-L*(L+1)/(R(I)*R(I)) 10 CONTINUE CALL FOX(IB,B,W,U) * >>>>>Input/output UMIN=U(0) UMAX=U(0) DO 5000 I=1,IB UMIN=MIN(UMIN,U(I)) UMAX=MAX(UMAX,U(I)) 5000 CONTINUE CALL GOPEN(R(0),R(IB),-1.D0,UMIN,UMAX,0.D0, & 'Distance [fm]','Radial wave function',0,1) IF(IFLAG.EQ.1) THEN CALL GCOLOR('Please choose colour for the wave function',ICOLOR) IFLAG=0 ENDIF CALL GMOVE(1,R(0),U(0),ICOLOR,1) DO 6000 I=1,IB CALL GDRAW(1,R(I),U(I)) 6000 CONTINUE CALL GCLOSE 7000 CONTINUE WRITE(*,'(T2,A)') 'Numerical output required (Y/N) ? ' READ(*,'(A1)') CONT IF(INDEX('Yy ',CONT).NE.0) THEN IS=(IB-1)/IOMAX+1 JLINE=0 WRITE(6,'(T7,I10)')L WRITE(6,'(T7,F10.5)')B WRITE(6,'(T7,F10.5)')E DO 8000 I=0,IB,IS C IF(JLINE.EQ.0) WRITE(6,'(T7,A,T30,A/)') C & 'Distance [fm]','Radial wave function (unnormalized)' C WRITE(6,'(T7,F10.5,T30,E13.5)') R(I),U(I) WRITE(6,'(T7,F15.5)')U(I) JLINE=JLINE+1 IF((JLINE.EQ.ILINE-4).OR.(I.GT.IB-IS)) THEN WRITE(*,'(T2/T2,A)') 'Continue ? Please press ENTER key ...' READ(*,'(A1)') CONT JLINE=0 CALL CURSOR(1) ENDIF 8000 CONTINUE ELSEIF(INDEX('Nn',CONT).EQ.0) THEN GOTO 7000 ENDIF 9000 CONTINUE WRITE(*,'(T2/T2,A,I2,A,F7.3,A/T2,A)') &'Rerun program for L = ',L,' and B = ',B,' fm with', &'different values for E and IB (Y/N) ? ' READ(*,'(A1)') CONT IF(INDEX('Yy ',CONT).NE.0) THEN GOTO 3000 ELSEIF(INDEX('Nn',CONT).EQ.0) THEN GOTO 9000 ENDIF 10000 CONTINUE WRITE(*,'(T2/T2,A)') &'Rerun program with completely new values (Y/N) ? ' READ(*,'(A1)') CONT IF(INDEX('Yy ',CONT).NE.0) THEN GOTO 1000 ELSEIF(INDEX('Nn',CONT).EQ.0) THEN GOTO 10000 ENDIF WRITE(*,'(T2//T2,A///)') 'Program ends !' * <<<<< END DOUBLE PRECISION FUNCTION V(R) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (V0=122.694D0, BETA=0.22D0) PARAMETER (EXPMAX=50.D0) IF(BETA*R*R.LE.EXPMAX) THEN V=-V0*EXP(-BETA*R*R) ELSE V=0.D0 ENDIF END SUBROUTINE FOX(IB,B,W,U) *--------------------------------------------------------------------- * Subroutine FOX integrates the differential equation * u"(r)+w(r)u(r)=0 in the interval [0,b] by the FOX-GOODWIN * method. * --Input parameters-- * IB: Number of integration steps in the interval [0,b] * {INTEGER variable} * B: Upper integration limit b * {DOUBLE PRECISION variable} * W: Function w(r), with W(I) as the function value at the * point I*H, H=B/IB, I=0,1,2,...,IB * {DOUBLE PRECISION array of dimension (0:IB)} * --Output parameter-- * U: Solution function u(r), with U(I) as the function value * at the point I*H, H=B/IB, I=0,1,2,...,IB * {DOUBLE PRECISION array of dimension (0:IB)} *--------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION W(0:IB), U(0:IB) H=B/IB U(0)=0.D0 U(1)=1.D0 DO 10 I=1,IB-1 U(I+1)=(2.D0*U(I)-U(I-1)-H*H/12.D0*(10.D0*W(I)*U(I)+W(I-1) & *U(I-1)))/(1.D0+H*H/12.D0*W(I+1)) 10 CONTINUE END