      PROGRAM H1
C...
C...  H1 - ERF SOLUTION OF FOURIER'S SECOND LAW
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  EXTERNAL FUNCTIONS TO CALCULATION INTEGRANDS FOR VARIOUS CASES
C...
C...     F3    THIRD ORDER POLYNOMIAL   (CUBIC FUNCTION)
C...
C...     F4   FOURTH ORDER POLYNOMIAL (QUARTIC FUNCTION)
C...
C...     FE   INTEGRAND OF ERF(X)
C...
      EXTERNAL F3, F4, FE
C...
C...  OPEN OUTPUT FILE
      NO=7
      OPEN(NO,FILE='output')
C...
C...  INTEGRATION LIMITS
      A=0.0D0
      B=1.0D0
C...
C...  STEP THROUGH A SERIES OF RUNS FOR POLYNOMIALS WITH DIFFERING
C...  NUMBERS OF PANELS
      DO 2 NORUN=1,4
      IF(NORUN.EQ.1)NP=2
      IF(NORUN.EQ.2)NP=4
      IF(NORUN.EQ.3)NP=10
      IF(NORUN.EQ.4)NP=100
C...
C...  INTEGRATION OF THIRD ORDER POLYNOMIAL
      CI3=SIMP(F3,A,B,NP)
C...
C...  INTEGRATION OF FOURTH ORDER POLYNOMIAL
      CI4=SIMP(F4,A,B,NP)
C...
C...  PRINT NUMERICAL INTEGRALS
      WRITE (NO,1)NP,
     +           A,B,CI3,
     +           A,B,CI4
1     FORMAT(' NP = ',I3,//,
     +       ' A = ',F6.2,'   B = ',F6.2,'   I3 = ',F10.7,/,
     +       ' A = ',F6.2,'   B = ',F6.2,'   I4 = ',F10.7,/)
2     CONTINUE
C...
C...  HEAT CONDUCTION PROBLEM
C...
C...  THREE RUNS FOR T = 30, 60, 120 S (SECONDS)

      DO 3 NORUN=1,3
      IF(NORUN.EQ.1)T= 50.0D0
      IF(NORUN.EQ.2)T=100.0D0
      IF(NORUN.EQ.3)T=200.0D0
C...
C...  THERMAL DIFFUSIVITY, INITIAL, SURFACE TEMPERATURES
      ALPHA=0.0001D0
      T0= 25.0D0
      TS=250.0D0
C...
C...  PI
      PI=4.0D0*DATAN(1.0D0)
C...
C...  X COORDINATES
      X1=0.25D0
      X2=1.00D0
C...
C...  LIMITS
      A=0.0D0
      B1=X1/DSQRT(4.0D0*ALPHA*T)
      B2=X2/DSQRT(4.0D0*ALPHA*T)
C...
C...  EVALUATION OF ERF(X)
      NP=100
      ERFB1=SIMP(FE,A,B1,NP)
      ERFB2=SIMP(FE,A,B2,NP)
C...
C...  INCLUDE FACTOR 2/SQRT(PI)
      ERFB1=2.0D0/DSQRT(PI)*ERFB1
      ERFB2=2.0D0/DSQRT(PI)*ERFB2
C...
C...  TEMPERATURES AT X1, X2
      T1=TS+(T0-TS)*ERFB1
      T2=TS+(T0-TS)*ERFB2
C...
C...  PRINT NUMERICAL RESULTS
      WRITE(NO,4)NP,T,A,B1,ERFB1,T1,A,B2,ERFB2,T2
4     FORMAT(//,' NP = ',I3,' T = ',F6.1,//,
     + '   A = ',F6.2,'   B1 = ',F6.2,
     + '   ERF(B1) = ',F10.7,'   T1 = ',F10.7,/,
     + '   A = ',F6.2,'   B2 = ',F6.2,
     + '   ERF(B2) = ',F10.7,'   T2 = ',F10.7,/)
C...
C...  ERF(X) FROM HASTINGS APPROXIMATION
      ERFB1=ERF(B1)
      ERFB2=ERF(B2)
      WRITE(NO,5)B1,ERFB1,B2,ERFB2
5     FORMAT(//,' ERF FROM HASTINGS APPROXIMATION',/,
     +       '   B1 = ',F6.2,'   ERF(B1) = ',F10.7,/,
     +       '   B2 = ',F6.2,'   ERF(B2) = ',F10.7,/)
C...
C...  NEXT TIME
3     CONTINUE
C...
C...  CALCULATIONS COMPLETE
      STOP
      END
      DOUBLE PRECISION FUNCTION SIMP(F,A,B,NP)
C...
C...  FUNCTION SIMP COMPUTES A SIMPSON'S RULE QUADRATURE
C...
C...  ARGUMENT LIST
C...
C...     F      EXTERNAL FUNCTION TO DEFINE THE INTEGRAND (INPUT)
C...
C...     A      LOWER LIMIT OF THE INTEGRAL (INPUT)
C...
C...     B      UPPER LIMIT OF THE INTEGRAL (INPUT)
C...
C...     NP     NUMBER OF PANELS (MUST BE EVEN) (INPUT)
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  EXTERNAL FUNCTION FOR INTEGRAND
      EXTERNAL F
C...
C...  CHECK FOR AN EVEN NUMBER OF PANELS
      IF((NP/2*2).NE.NP)THEN
         WRITE(*,2)NP
2        FORMAT(//,' NP = ',I3,' (NOT EVEN)')
         SIMP=0.0D0
         RETURN
      END IF
C...
C...  INTEGRATION INTERVAL
      DX=(B-A)/DFLOAT(NP)
C...
C...  INITIALIZE SUM
      SIMP=(F(A)+4.0D0*F(A+DX)+F(B))
C...
C...  CONTINUE SUM
      IF(NP.GT.2)THEN
         DO 1 I=2,NP-1,2
            SIMP=SIMP
     +          +2.0D0*F(A+DFLOAT(I  )*DX)
     +          +4.0D0*F(A+DFLOAT(I+1)*DX)

1        CONTINUE
      END IF
C...
C...  INTEGRAL
      SIMP=(DX/3.0D0)*SIMP
      RETURN
      END
      DOUBLE PRECISION FUNCTION F3(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      F3=X**3
      RETURN
      END
      DOUBLE PRECISION FUNCTION F4(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      F4=X**4
      RETURN
      END
      DOUBLE PRECISION FUNCTION FE(X)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      FE=DEXP(-X**2)
      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
      DOUBLE PRECISION FUNCTION ERF(X)
C...
C...  FUNCTION ERF COMPUTES THE ERROR FUNCTION ACCORDING TO THE
C...  HASTINGS APPROXIMATION:
C...
C...  GREENBERG, M. D. 1978. FOUNDATIONS OF APPLIED MATHEMATICS.
C...  ENGLEWOOD CLIFFS, NJ: PRENTICE-HALL.
C...
C...  HASTINGS, C., JR. 1955. APPROXIMATIONS FOR DIGITAL COMPUTERS.
C...  PRINCETON, NJ: PRINCETON UNIVERSITY PRESS.
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  THE FOLLOWING CODING IS A STRAIGHTFORWARD IMPLEMENTATION OF
C...  THE HASTINGS APPROXIMATION CITED ABOVE
      T=  1.0D0/(1.0D0+0.3275911D0*X)
      A1= 0.254829592D0
      A2=-0.284496736D0
      A3= 1.421413741D0
      A4=-1.453152027D0
      A5= 1.061405429D0
      ERF=1.0D0-(A1*T+A2*T**2+A3*T**3+A4*T**4+A5*T**5)*DEXP(-X**2)
      RETURN
      END
