C---- COULOMB DRIVER (FULL TEST VERSION) FILE = COUL90_T.FOR C---- A.R.BARNETT MANCHESTER EASTER 1995 (REVISED) C---- COMPUTES COULOMBS, SPHERICAL BESSELS, AND CYLINDRICAL BESSELS C---- PREPARED FOR 'ATOMIC PHYSICS' BOOK 1995, EDITOR K. BARTSCHAT C---- THE DRIVER TESTS COUL90, SBESJY AND RICBES SUBROUTINES PROGRAM COUL90_TEST IMPLICIT NONE INTEGER M PARAMETER (M=1200) DOUBLE PRECISION FC(0:M), GC(0:M), FCP(0:M), GCP(0:M) DOUBLE PRECISION XJ(0:M), Y(0:M), XJP(0:M), YP(0:M) DOUBLE PRECISION PSI(0:M), CHI(0:M), PSIP(0:M), CHIP(0:M) DOUBLE PRECISION CF1, P, Q, F, GAMMA, WRONSK, PACCQ, ERR DOUBLE PRECISION X(21), ETA, XM, ZERO INTEGER LAMBDA(10), LRANGE, NL, NX, KFN INTEGER LMIN, IFAIL, K, N, L, NF, NFP, NPQ, IEXP, MINL CHARACTER TEXT*72, CHFN*4, INP*7 LOGICAL SPHRIC, FIRST COMMON /STEDE/ ERR, NF COMMON /STEED/ PACCQ, NFP, NPQ, IEXP, MINL COMMON /DESET/ CF1, P, Q, F, GAMMA, WRONSK ZERO = 0.D0 OPEN (1, FILE = 'COULTEST.IN') OPEN (2, FILE = 'COULTEST.OUT') WRITE (2, 1000) READ (1, 1001) TEXT 10 CONTINUE NL=0 NX=0 READ (1, 1002, ERR=41, END=50) CHFN, INP, NL, NX, ETA, XM IF ( NL.LE.0 ) STOP READ (1, *, ERR=42) (LAMBDA(K), K = 1, NL) ! number of L-values READ (1, *, ERR=43) (X(K), K = 1, ABS(NX)) ! number of X-values IF ( NX.LT.0 ) THEN READ (1, 1001) TEXT WRITE (2, 1003) TEXT ! annotate output END IF DO 40 K = 1, ABS(NX) ! main loop over x-values LRANGE = LAMBDA(NL) - LAMBDA(1) ! additional L-values FIRST = ( K.EQ.1 ) IF ( CHFN.EQ.'COUL' ) KFN = 0 IF ( CHFN.EQ.'SBES' ) KFN = 1 IF ( CHFN.EQ.'BESS' ) KFN = 2 IF ( CHFN.EQ.'BESI' .AND. FIRST ) THEN KFN = 2 ! Special case XM = 1.D0 / XM ! eg for 1/3 Bessels END IF C KFN = 0 (COULOMBS), = 1 (SPHERICAL BESSELS), = 2 (CYLINDER BESSELS) CALL COUL90( X(K), ETA, XM, LRANGE, FC, GC, FCP, GCP, KFN, IFAIL) IF ( IFAIL.NE.0 ) WRITE (2,1200) ' COUL90 ERROR! IFAIL=',IFAIL SPHRIC = ( KFN.NE.2 .AND. ETA.EQ.ZERO ) IF ( SPHRIC ) THEN CALL SBESJY(X(K), LRANGE, XJ, Y, XJP, YP, IFAIL) IF ( IFAIL.NE.0 ) WRITE (2,1200) ' SBESJY ERROR! IFAIL=',IFAIL CALL RICBES(X(K), LRANGE, PSI, CHI, PSIP, CHIP, IFAIL) IF ( IFAIL.NE.0 ) WRITE (2,1200) ' RICBES ERROR! IFAIL=',IFAIL END IF IF ( KFN.EQ.0 ) WRITE (2, 1004) ETA, X(K), XM, KFN,' (Coulomb)' IF ( KFN.EQ.1 ) WRITE (2, 1005) X(K), XM, KFN,' (SphBess)' IF ( KFN.EQ.2 ) WRITE (2, 1006) X(K), XM, KFN,' (CylBess)' IF ( IFAIL.NE.0 ) THEN WRITE (2, 1007) IFAIL, NL, KFN, X(K), ETA, XM, LAMBDA(NL) GOTO 10 ! abort & try next parameters END IF TEXT = + ' The 3 sets of results are COUL90(KFN), SBESJY & (1/X) RICBES' IF ( SPHRIC ) WRITE (2, 1003) TEXT LMIN = LAMBDA(1) - IDINT(XM) WRITE (6, 1201) X(K), ETA, XM, LRANGE, LMIN C compile these lines for diagnostics : C WRITE (2, 1008) NFP, NPQ, PACCQ C IF ( KFN.GT.0 ) WRITE (2, 1009) NF , KFN, ERR C WRITE (2, 1202) CF1, P, Q, F, GAMMA, WRONSK DO 30 N = 1, NL L = LMIN + LAMBDA(N) WRITE (2, 1010) LAMBDA(N), FC(L), GC(L), FCP(L), GCP(L) IF ( SPHRIC ) THEN WRITE (2, 1010) LAMBDA(N), XJ(L), Y(L), XJP(L), YP(L) WRITE (2, 1010) LAMBDA(N), PSI(L), CHI(L), PSIP(L), CHIP(L) END IF 30 CONTINUE IF ( IEXP.GT.1 ) WRITE (2, 1011) IEXP 40 CONTINUE GOTO 10 ! more data 41 WRITE(6,1002) CHFN, INP, NL, NX, ETA, XM 42 WRITE(6,1012) NL, NX, (REAL(LAMBDA(K)), K = 1, NL) 43 WRITE(6,1012) NL, NX, (X(K), K = 1, ABS(NX)) 50 CLOSE (1) CLOSE (2) C---- 1000 FORMAT (10X, + ' TEST OF THE CONTINUED-FRACTION COULOMB & BESSEL', + ' PROGRAM - COUL90'//12X, + ' WHEN LAMBDA IS REAL (IN GENERAL) = L + XM '// + ' F IS REGULAR AT THE ORIGIN ( X = 0 ) WHILE',/ + ' G IS IRREGULAR ( => INFINITY AT X = 0 )'//5X, + 'L', 2X, ' F(ETA,X,LAMBDA)', 2X, ' G(ETA,X,LAMBDA)', 1X, + ' FCP = D/DX (F)', 2X, ' GCP = D/DX (G)'/) 1001 FORMAT (A72) 1002 FORMAT (1X, A4, A7, 2I3, 2F10.2 ) 1003 FORMAT (/4X, A72) 1004 FORMAT (/1X, 'ETA =', F9.3, 4X, ' X = ', F10.3, 4X, 'XLMIN = ', + F6.2, 4X, 'KFN = ', I2, A10) 1005 FORMAT (/5X, 'sph Bessels', 3X, ' X = ', F10.3, 4X, 'XLMIN = ', + F6.2, 4X, 'KFN = ', I2, A10) 1006 FORMAT (/5X, 'CYL Bessels', 3X, ' X = ', F10.3, 4X, 'XLMIN = ', + F6.2, 4X, 'KFN = ', I2, A10) 1007 FORMAT (1X, 'IFAIL =', I4, ' NL,KFN =', 2I4, 3F12.4, I12) 1008 FORMAT (1X,': NFP,NPQ,PACCQ = ', 2I5, 1PD9.0, '---') 1009 FORMAT (1X,': NF ,KFN,ERR = ', 2I5, 1PD9.0, '---') 1010 FORMAT (1X, I5, 1P4D17.7 ) 1011 FORMAT (12X, ' **** IEXP = ', I6, ' F,FP *10**(-IEXP)', + ' G,GP *10**(+IEXP)') 1012 FORMAT (1X, 2I3, F10.3, 9F6.2) 1200 FORMAT (1X, A25, I5) 1201 FORMAT (3F10.3, 3I10, F10.3) 1202 FORMAT (6(1PD13.5)) END $INCLUDE: 'COUL90.FOR' $INCLUDE: 'RICBES.FOR' $INCLUDE: 'SBESJY.FOR'