      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
