      PROGRAM H3SS1
C...
C...  H3 - HEAT CONDUCTION IN A CIRCULAR FIN, STEADY STATE
C...       SOLUTION, VERSION 1
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C...
C...  ADDITIONAL VARIABLES ARE DECLARED DOUBLE PRECISION (TO 
C...  ACCOMMODATE THE MODIFIED BESSEL FUNCTIONS)
      DOUBLE PRECISION I00, I01, I11,
     +                 K00, K01, K11 
C...
C...  COMMON AREA WITH PROBLEM PARAMETERS
      PARAMETER(NR=6)
      COMMON/P/  ALPHA,   BETA,     UA,     U0,  U(NR),
     +              DR,     R0,     R1,  R(NR),      B,  
     +           RHOCP
C...
C...  OPEN AN OUTPUT FILE
      NO=6
      OPEN(NO,FILE='OUTPUT', STATUS='UNKNOWN')
C...
C...  CALL SUBROUTINE BFTEST TO CHECK THE PERFORMANCE OF FUNCTIONS
C...  BESSI0, BESSI1, BESSK0 AND BESSK1
      CALL BFTEST(NO)
C...
C...  PROBLEM PARAMETERS
C...
C...     FIN INSIDE RADIUS
         R0=0.5D0
C...
C...     FIN OUTSIDE RADIUS
         R1=2.0D0
C...
C...     FIN THICKNESS
         B=0.2D0
C...
C...     HEAT TRANSFER COEFFICIENT
         H=0.01D0
C...
C...     RHO*CP
         RHOCP=2.0D0
C...
C...     THERMAL DIFFUSIVITY
         ALPHA=0.01D0
C...
C...     AMBIENT TEMPERATURE
         UA=25.0D0
C...
C...     INSIDE TEMPERATURE
         U0=150.0D0
C...
C...  RECIPROCAL TIME FOR HEAT TRANSFER
      BETA=(2.0D0*H)/(B*RHOCP)
C...
C...  INCREMENT FOR RADIAL GRID
      DR=(R1-R0)/DFLOAT(NR-1)
C...
C...  RADIAL GRID
      DO 1 I=1,NR
         R(I)=R0+DFLOAT(I-1)*DR
1     CONTINUE
C...
C...  STEADY STATE SOLUTION WITH DIRICHLET BOUNDARY CONDITIONS
C...
C...     CONSTANTS C1, C2 
         ARG0=DSQRT(BETA/ALPHA)*R0
         ARG1=DSQRT(BETA/ALPHA)*R1
         I00=BESSI0(ARG0)
         I01=BESSI0(ARG1)
         K00=BESSK0(ARG0)
         K01=BESSK0(ARG1)
         DEN=I00*K01-I01*K00
         C1= K01/DEN
         C2=-I01/DEN
C...
C...     EVALUATE THE SOLUTION AT A SERIES OF RADIAL POSITIONS
         DO 2 I=1,NR
         ARG=DSQRT(BETA/ALPHA)*R(I)
         U(I)=C1*BESSI0(ARG)+C2*BESSK0(ARG)
C...
C...     CONVERSION FROM DIMENSIONLESS TO DIMENSIONAL TEMPERATURE
         U(I)=UA+(U0-UA)*U(I)
2        CONTINUE
C...
C...  PRINT THE SOLUTION
      WRITE(NO,3)(R(I),I=1,NR),(U(I),I=1,NR)
3     FORMAT(/,' Dirichlet boundary conditions',/,
     +         '      r ',6F9.2,/,'   u(r) ',6F9.3,/)
C...
C...  COMPUTE AND PRINT THE TOTAL HEAT TRANSFER TO THE FIN
      DTDR=(C1*BESSI1(ARG0)-C2*BESSK1(ARG0))*DSQRT(BETA/ALPHA)*(U0-UA)
      PI=4.0D0*DATAN(1.0D0)
      TOTAL=-2.0D0*PI*R0*B*ALPHA*RHOCP*DTDR
      WRITE(NO,4)TOTAL
4     FORMAT(' Total heat transfer rate = ',F7.3,//)
C...
C...  STEADY STATE SOLUTION WITH DIRICHLET-NEUMANN BOUNDARY 
C...  CONDITIONS
C...
C...     CONSTANTS C1, C2 
         ARG0=DSQRT(BETA/ALPHA)*R0
         ARG1=DSQRT(BETA/ALPHA)*R1
         I00=BESSI0(ARG0)
         I11=BESSI1(ARG1)
         K00=BESSK0(ARG0)
         K11=BESSK1(ARG1)
         DEN=I00*K11+I11*K00
         C1=K11/DEN
         C2=I11/DEN
C...
C...     EVALUATE THE SOLUTION AT A SERIES OF RADIAL POSITIONS
         DO 12 I=1,NR
         ARG=DSQRT(BETA/ALPHA)*R(I)
         U(I)=C1*BESSI0(ARG)+C2*BESSK0(ARG)
C...
C...     CONVERSION FROM DIMENSIONLESS TO DIMENSIONAL TEMPERATURE
         U(I)=UA+(U0-UA)*U(I)
12       CONTINUE
C...
C...  PRINT THE SOLUTION
      WRITE(NO,13)(R(I),I=1,NR),(U(I),I=1,NR)
13    FORMAT(/,' Dirichlet-Neumann boundary conditions',/,
     +         '      r ',6F9.2,/,'   u(r) ',6F9.3,/)
C...
C...  COMPUTE AND PRINT THE TOTAL HEAT TRANSFER TO THE FIN
      DTDR=(C1*BESSI1(ARG0)-C2*BESSK1(ARG0))*DSQRT(BETA/ALPHA)*(U0-UA)
      PI=4.0D0*DATAN(1.0D0)
      TOTAL=-2.0D0*PI*R0*B*ALPHA*RHOCP*DTDR
      WRITE(NO,14)TOTAL
14    FORMAT(' Total heat transfer rate = ',F7.3,//)
      STOP
      END
      SUBROUTINE BFTEST(NO)
C...
C...  PROGRAM BFTEST CALLS FUNCTIONS BESSI0, BESSI1, BESSK0 AND 
C...  BESSK1 TO CHECK THE CALCULATION OF I0(X), I1(X), K0(X) AND 
C...  K1(X).
C...
C...  THE FOLLOWING TABULATED VALUES CAN BE USED TO CHECK THE
C...  PERFORMANCE OF BESSI0 TO BESSK1 [SPIEGEL (1991), PP. 246-
C...  247]
C...
C...     I0(1) = 1.266,  I0(5) = 27.24
C...
C...     I1(1) = 0.5652, I1(5) = 24.34
C...
C...     K0(1) = 0.4210, K0(5) = 0.003691
C...
C...     K1(1) = 0.6019, K1(5) = 0.004045
C...
C...  SPIEGEL, M. R., (1991), MATHEMATICAL HANDBOOK OF FORMULAS AND 
C...  TABLES, MCGRAW-HILL, NEW YORK
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  STEP THROUGH A SERIES OF VALUES OF X
      DO 1 I=1,10
         X=0.5D0*DFLOAT(I)
         YI0=BESSI0(X)
         YI1=BESSI1(X)
         YK0=BESSK0(X)
         YK1=BESSK1(X)
C...
C...  WRITE THE COMPUTED BESSEL FUNCTIONS
      WRITE(NO,2)X, YI0, YI1, YK0, YK1
2     FORMAT(' X = ',F5.2,'   YI0 = ',F8.4,'   YI1 = ',F8.4, 
     +                    '   YK0 = ',F8.4,'   YK1 = ',F8.4,/) 
C...
C...  NEXT X
1     CONTINUE
C...
C...  END OF CALCULATION
      END
      DOUBLE PRECISION FUNCTION BESSI0(X)
C...
C...  FUNCTION BESSI0 COMPUTES THE BESSEL FUNCTION I0(X).  BESSI0 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            P1,                  P2,                  P3,
     +                P4,                  P5,                  P6,
     +                P7          
     +            /1.0D0,         3.5156229D0,         3.0899424D0,
     +       1.2067492D0,         0.2659732D0,         0.360768D-1,
     +        0.45813D-2/
      DATA            Q1,                  Q2,                  Q3,
     +                Q4,                  Q5,                  Q6,
     +                Q7,                  Q8,                  Q9
     +     /0.39894228D0,        0.1328592D-1,         0.225319D-2,
     +      -0.157565D-2,         0.916281D-2,       -0.2057706D-1,
     +      0.2635537D-1,       -0.1647633D-1,         0.392377D-2/
      IF(DABS(X).LT.3.75D0)THEN
         Y=(X/3.75D0)**2
         BESSI0=P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7)))))
      ELSE
         AX=DABS(X)
         Y=3.75D0/AX
         BESSI0=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4
     +         +Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION BESSI1(X)
C...
C...  FUNCTION BESSI1 COMPUTES THE BESSEL FUNCTION I1(X).  BESSI1 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            P1,                  P2,                  P3,
     +                P4,                  P5,                  P6,
     +                P7
     +            /0.5D0,        0.87890594D0,        0.51498869D0,
     +      0.15084934D0,        0.2658733D-1,         0.301532D-2,
     +        0.32411D-3/
      DATA            Q1,                  Q2,                  Q3,
     +                Q4,                  Q5,                  Q6,
     +                Q7,                  Q8,                  Q9
     +     /0.39894228D0,       -0.3988024D-1,        -0.362018D-2,
     +       0.163801D-2,       -0.1031555D-1,        0.2282967D-1,
     +     -0.2895312D-1,        0.1787654D-1,        -0.420059D-2/
      IF(DABS(X).LT.3.75D0)THEN
         Y=(X/3.75D0)**2
         BESSI1=X*(P1+Y*(P2+Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
         AX=DABS(X)
         Y=3.75D0/AX
         BESSI1=(DEXP(AX)/DSQRT(AX))*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+
     +         Y*(Q5+Y*(Q6+Y*(Q7+Y*(Q8+Y*Q9))))))))
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION BESSK0(X)
C...
C...  FUNCTION BESSK0 COMPUTES THE BESSEL FUNCTION K0(X).  BESSK0 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            P1,                  P2,                  P3,
     +                P4,                  P5,                  P6,
     +                P7
     +    /-0.57721566D0,        0.42278420D0,        0.23069756D0,
     +      0.3488590D-1,         0.262698D-2,          0.10750D-3,
     +           0.74D-5/
      DATA            Q1,                  Q2,                  Q3,
     +                Q4,                  Q5,                  Q6,
     +                Q7
     +     /1.25331414D0,       -0.7832358D-1,        0.2189568D-1,
     +     -0.1062446D-1,         0.587872D-2,        -0.251540D-2,
     +        0.53208D-3/
      IF(X.LE.2.0D0)THEN
         Y=X*X/4.0D0
         BESSK0=(-DLOG(X/2.0D0)*BESSI0(X))+(P1+Y*(P2+Y*(P3+
     +         Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
         Y=(2.0D0/X)
         BESSK0=(DEXP(-X)/DSQRT(X))*(Q1+Y*(Q2+Y*(Q3+
     +         Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION BESSK1(X)
C...
C...  FUNCTION BESSK1 COMPUTES THE BESSEL FUNCTION K1(X).  BESSK1 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            P1,                  P2,                  P3,
     +                P4,                  P5,                  P6,
     +                P7
     +            /1.0D0,        0.15443144D0,       -0.67278579D0,
     +     -0.18156897D0,       -0.1919402D-1,        -0.110404D-2,
     +        -0.4686D-4/
      DATA            Q1,                  Q2,                  Q3,
     +                Q4,                  Q5,                  Q6,
     +                Q7
     +     /1.25331414D0,        0.23498619D0,       -0.3655620D-1,
     +      0.1504268D-1,        -0.780353D-2,         0.325614D-2,
     +       -0.68245D-3/
      IF(X.LE.2.0D0)THEN
         Y=X*X/4.0D0
         BESSK1=(DLOG(X/2.0D0)*BESSI1(X))+(1.0/X)*(P1+Y*(P2+
     +         Y*(P3+Y*(P4+Y*(P5+Y*(P6+Y*P7))))))
      ELSE
         Y=2.0D0/X
         BESSK1=(DEXP(-X)/DSQRT(X))*(Q1+Y*(Q2+Y*(Q3+
     +         Y*(Q4+Y*(Q5+Y*(Q6+Y*Q7))))))
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION DFLOAT(I)
C...
C...  FUNCTION DFLOAT CONVERTS A SINGLE PRECISION INTEGER INTO A
C...  DOUBLE PRECISION FLOATING POINT.  THIS FUNCTION IS PROVIDED
C...  IN CASE THE USER'S FORTRAN COMPILER DOES NOT INCLUDE DFLOAT.
C...
      DFLOAT=DBLE(FLOAT(I))
      RETURN
      END
