      PROGRAM M1SER
C...
C...  PROGRAM M1SER COMPUTES THE SOLUTION TO THE BLASIUS EQUATION
C...  FROM THE SERIES SOLUTION [1]
C...
C...         INF           (ALPHA**(N+1) CN)
C...     F = SUM (-1/2)**N ----------------- ETA**(3*N + 2)   (1)
C...         N=0               (3*N + 2)!
C...
C...  WHERE ALPHA =  0.3320, C0 = 1, C1 = 1, C2 = 11, C3 = 375,
C...  C4 = 27,897, C5 = 3,817,137
C...
C...  [1]  HUGHES, W. F., AND J. A. BRIGHTON, THEORY AND PROBLEMS
C...  OF FLUID DYNAMICS, 2ND ED., MCGRAW-HILL, INC., NEW YORK,
C...  1991
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  UPPER LIMIT OF THE SERIES (RATHER THAN INFINITY)
      PARAMETER (NU=5)
C...
C...  COMMON AREA FOR THE SERIES COEFFICIENTS
      COMMON/COEFF/ ALPHA, C(0:NU)
C...
C...  OPEN OUTPUT FILE
      NO=1
      OPEN(NO,FILE='OUTPUT',STATUS='UNKNOWN')
C...
C...  SET THE COEFFICIENTS IN THE SERIES
      CALL COEFF
C...
C...  WRITE A HEADING FOR THE COMPUTED OUTPUT
      WRITE(NO,5)
5     FORMAT(//,5X,' Series Solution',//,4X,
     +             'eta',9X,'f(eta)',/,
     +         15X,'(series)')
C...
C...  STEP THROUGH A SERIES OF VALUES OF ETA
      DO 1 I=1,11
         ETA=DFLOAT(I-1)
C...
C...     EVALUATE THE SERIES
         SUM=0.0D0
         DO 2 N=0,5
         SUM=SUM+TERM(N,ETA)
2        CONTINUE
C...
C...     PRINT THE SERIES SOLUTION
         WRITE(NO,3)ETA,SUM
3        FORMAT(F7.2,F15.4)
C...
C...     NEXT VALUE OF ETA
1     CONTINUE
C...
C...  PERFORM A HAND CALCULATION FOR ETA = 1 AS A CHECK OF THE
C...  PRECEDING PROGRAMMING
      CALL HAND(F)
      WRITE(NO,4)F
4     FORMAT(//,'  Series for eta = 1 by',/,
     +          ' hand calculation = ',F6.4)
C...
C...  CALCULATIONS ARE COMPLETE
      STOP
      END
      SUBROUTINE COEFF
C...
C...  SUBROUTINE COEFF SETS THE COEFFICIENTS IN THE SERIES
C...  SOLUTION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  UPPER LIMIT OF THE SERIES
      PARAMETER (NU=5)
C...
C...  COMMON AREA FOR THE SERIES COEFFICIENTS
      COMMON/COEFF/ ALPHA, C(0:NU)
C...
C...  SET THE COEFFICIENTS
      ALPHA=0.3320D0
      C(0)=1.0D0
      C(1)=1.0D0
      C(2)=11.0D0
      C(3)=375.0D0
      C(4)=27897.0D0
      C(5)=3817137.0D0
      RETURN
      END
      DOUBLE PRECISION FUNCTION TERM(N,ETA)
C...
C...  FUNCTION TERM COMPUTES A TERM IN THE SERIES SOLUTION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  UPPER LIMIT OF THE SERIES
      PARAMETER (NU=5)
C...
C...  COMMON AREA FOR THE SERIES COEFFICIENTS
      COMMON/COEFF/ ALPHA, C(0:NU)
C...
C...  TERM IN THE SERIES
      TERM=((-0.5D0)**N)*(ALPHA**(N+1))*C(N)*
     +     (1.0D0/FACT(3*N+2))*(ETA**(3*N+2))
      RETURN
      END
      DOUBLE PRECISION FUNCTION FACT(N)
C...
C...  FUNCTION COMPUTES THE FACTORIAL IN THE SERIES SOLUTION
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      FACT=1.0D0
      DO 2 I=1,N
         FACT=FACT*DFLOAT(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE HAND(F)
C...
C...  SUBROUTINE HAND DOES A 'HAND' CALCULATION OF THE SERIES
C...  FOR ETA = 1 TO CHECK THE PROGRAMMING OF THE SERIES
C...
C...         INF           (ALPHA**(N+1) CN)
C...     F = SUM (-1/2)**N ----------------- ETA**(3*N + 2)   (1)
C...         N=0               (3*N + 2)!
C...
C...  WHERE ALPHA =  0.3320, C0 = 1, C1 = 1, C2 = 11, C3 = 375,
C...  C4 = 27,897, C5 = 3,817,137
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  COMPUTE SIX TERMS IN THE SERIES
C...
C...  N=0
      TERM0=0.3320D0/2.0D0
C...
C...  N=1
      FACT5=1.0D0*2.0D0*3.0D0*4.0D0*5.0D0
      TERM1=-0.5D0*(0.3320**2)/FACT5
C...
C...  N=2
      FACT8=FACT5*6.0D0*7.0D0*8.0D0
      TERM2=(-0.5D0**2)*(0.3320**3)*11.0D0/FACT8
C...
C...  N=3
      FACT11=FACT8*9.0D0*10.0D0*11.0D0
      TERM3=(-0.5D0**3)*(0.3320**4)*375.0D0/FACT11
C...
C...  N=4
      FACT14=FACT11*12.0D0*13.0D0*14.0D0
      TERM4=(-0.5D0**4)*(0.3320**5)*27897.0D0/FACT14
C...
C...  N=5
      FACT17=FACT14*15.0D0*16.0D0*17.0D0
      TERM5=(-0.5D0**5)*(0.3320**6)*3817137.0D0/FACT17

C...  SUM SIX TERMS
      F=TERM0+TERM1+TERM2+TERM3+TERM4+TERM5
C...
C...  HAND CALCULATION IS COMPLETE
      RETURN
      END
      DOUBLE PRECISION FUNCTION DFLOAT(I)

C...  FUNCTION DFLOAT CONVERTS A SINGLE PRECISION INTEGER INTO
C...  A DOUBLE PRECISION FLOATING POINT.  THIS FUNCTION IS
C..   PROVIDED IN CASE THE USER'S FORTRAN COMPILER DOES NOT
C...  INCLUDE DFLOAT.
C...
      DFLOAT=DBLE(FLOAT(I))
      RETURN
      END
