PROGRAM M2CL C... C... PROGRAM M2CL EVALUATES THE INFINITE SERIES C... C... INF 1 C... PHI(0) = 1 - 8 SUM ------------------*EXP(-E(N)**2*T) C... N=1 (E(N)**3)*J1(E(N)) C... C... WHICH IS EQ. (4.1-40), BIRD ET AL (1960), WITH R = 0 (FOR THE C... CENTERLINE DIMENSIONLESS VELOCITY). HERE C... C... E(N) NTH ROOT OF J0(X) C... C... J1(X) BESSEL FUNCTION OF THE FIRST KIND C... C... BIRD, R. B., W. E. STEWART AND E. N. LIGHTFOOT (1960), TRANSPORT C... PHENOMENA, JOHN WILEY AND SONS, NEW YORK, P 129 C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... INPUT/OUTPUT UNIT NUMBERS COMMON/IO/ NI, NO C... C... NUMBERS OF TERMS IN THE SERIES, NUMBER OF VALUES OF TAU PARAMETER(NTERMS=20, NTAU=9) C... C... DIMENSION ARRAYS DIMENSION TAU(NTAU), E(NTERMS), BESSJ1(NTERMS) C... C... DIMENSIONLESS TIMES DATA TAU/ 0.0D0, 0.05D0, 0.1D0, 0.15D0, + 0.2D0, 0.3D0, 0.4D0, 0.5D0, + 2.0D0/ C... C... OUTPUT FILE NO=1 OPEN(NO,FILE='OUTPUT',STATUS='UNKNOWN') C... C... COMPUTE E(N), J1(E(N)) CALL ROOTJ0(NTERMS,E,BESSJ1) C... C... PRINT HEADING FOR SERIES SOLUTION WRITE(NO,1) 1 FORMAT(//,' EQ. (4.1-40), BSL',/) C... C... EVALUATE THE SERIES FOR NTAU VALUES OF TAU DO 3 NT=1,NTAU C... C... INITIALIZE SERIES SERIES=1.0D0 C... C... SUM SERIES DO 4 N=1,NTERMS SERIES=SERIES-8.0D0/((E(N)**3)*BESSJ1(N)) + *DEXP(-(E(N)**2)*TAU(NT)) 4 CONTINUE C... C... WRITE SERIES WRITE(NO,2)NT,TAU(NT),SERIES 2 FORMAT(' NT = ',I5,' TAU = ',F6.2,' PHI(0) = ',F10.5) C... C... NEXT VALUE OF TAU 3 CONTINUE C... C... CALCULATIONS ARE COMPLETE END SUBROUTINE ROOTJ0(NROOTS,ROOTS,BFJ1) C... C... SUBROUTINE ROOTJ0 COMPUTES THE ROOTS OF J0(X), AND EVALUATES C... J1(X) AT THESE ROOTS C... C... ARGUMENTS C... C... NROOTS NUMBER OF SUCESSIVE ROOTS OF JO(X) TO BE COMPUTED C... STARTING WITH THE FIRST ROOT (INPUT) C... C... ROOTS ARRAY WITH THE NROOTS COMPUTED ROOTS (OUTPUT) C... C... BFJ1 ARRAY WITH NROOTS VALUES OF J1(X) EVALUATED AT C... THE NROOTS OF J0(X) (OUTPUT) C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... INPUT/OUTPUT UNIT NUMBERS COMMON/IO/ NI, NO C... C... ARRAYS FOR ROOTS AND J1 DIMENSION ROOTS(NROOTS), BFJ1(NROOTS) C... C... STEPS THROUGH NROOTS ROOTS DO 1 N=1,NROOTS C... C... COMPUTE ROOT OF J0(X) ROOTS(N)=ZERO(N) C... C... EVALUATE J1(X) AT THE ROOT BFJ1(N)=BESSJ1(ROOTS(N)) C... C... WRITE RESULTS TO OUTPUT FILE WITH LABELS WRITE(NO,2)N,ROOTS(N),BFJ1(N) 2 FORMAT(' N = ',I4, + ' Z FOR J0(Z) = 0',F15.10, + ' J1(Z)',F14.10,/) C... C... NEXT ROOT 1 CONTINUE RETURN END DOUBLE PRECISION FUNCTION BESSJ1(X) C... C... FUNCTION BESSJ1 COMPUTES THE BESSEL FUNCTION J1(X). BESSJ1 IS C... TAKEN FROM PRESS, W. H., B. P. FLANNERY, S. A. TEUKOLSKY AND C... W. T. VETTERLING (1985), NUMERICAL RECIPIES, CAMBRIDGE UNIVERSITY C... PRESS, CAMBRIDGE, CHAPTER 6, WITH CONVERSION TO DOUBLE PRECISION C... AND EDITING BY W. E. SCHIESSER. C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA R1, R2, R3, + R4, R5, R6 + /72362614232.D0, -7895059235.D0, 242396853.1D0, + -2972611.439D0, 15704.48260D0, -30.16036606D0/ DATA S1, S2, S3, + S4, S5, S6 + /144725228442.D0, 2300535178.D0, 18583304.74D0, + 99447.43394D0, 376.9991397D0, 1.0D0/ DATA P1, P2, P3, + P4, P5/ + 1.D0, 0.183105D-2, -0.3516396496D-4, + 0.2457520174D-5, -0.240337019D-6/ DATA Q1, Q2, Q3, + Q4, Q5/ + 0.04687499995D0, -0.2002690873D-3, 0.8449199096D-5, + -0.88228987D-6, 0.105787412D-6/ IF(DABS(X).LT.8.0D0)THEN Y=X**2 BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6))))) + /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) ELSE AX=DABS(X) Z=8.0D0/AX Y=Z**2 XX=AX-2.356194491D0 BESSJ1=DSQRT(0.636619772D0/AX)*(DCOS(XX)*(P1+Y*(P2+Y*(P3+Y* + (P4+Y*P5))))-Z*DSIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5))))) + *DSIGN(1.0D0,X) ENDIF RETURN END DOUBLE PRECISION FUNCTION ZERO(N) C... C... FUNCTION ZERO COMPUTES THE ZEROS OF J0(Z) FROM A FORMULA GIVEN BY C... WATSON, G. N., A TREATISE ON THE THEORY OF BESSEL FUNCTIONS, CAM- C... BRIDGE UNIVERSITY PRESS, SECOND EDITION, 1966, PAGE 506. C... C... THE FORMULA FROM WATSON IS FOR "LARGE ZEROS". FOR THE FIRST FIVE C... "SMALL ZEROS", THE VALUES ARE SET DIRECTLY (RATHER THAN BY THE C... APPROXIMATION) C... C... DOUBLE PRECISION CODING IS USED IMPLICIT DOUBLE PRECISION (A-H,O-Z) C... C... THE FIRST FIVE ROOTS OF J0(X) ARE SET EXPLICITLY IF(N.EQ.1)ZERO=2.4048255577D0 IF(N.EQ.2)ZERO=5.5200781103D0 IF(N.EQ.3)ZERO=8.6537279129D0 IF(N.EQ.4)ZERO=11.7915344391D0 IF(N.EQ.5)ZERO=14.9309177086D0 IF(N.LE.5)RETURN C... C... FOR N GT 5, THE WATSON FORMULA IS USED PI=4.0D0*DATAN(1.0D0) B=(DFLOAT(N)-0.25D0)*PI ZERO=B+1.0D0/((2.0D0**3)*B)-31.0D0/(3.0D0*(2.0D0**7)*(B**3)) 1 +3779.0D0/(15.0D0*(2.0D0**15)*(B**5)) 2 -6277237.0D0/(105.0D0*(2.0D0**15)*(B**7)) RETURN END DOUBLE PRECISION FUNCTION DFLOAT(I) C... C... FUNCTION DFLOAT CONVERTS A SINGLE PRECISION INTEGER INTO A DOUBLE C... PRECISION FLOATING POINT NUMBER C... IMPLICIT DOUBLE PRECISION (A-H,O-Z) DFLOAT=DBLE(FLOAT(I)) RETURN END