PROGRAM H1
C...
C...  H1 - ERF SOLUTION OF FOURIER'S SECOND LAW; MATLAB PLOTTING
C...       OF THE SOLUTION
C...
C...  DOUBLE PRECISION CODING IS USED
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C...
C...  SPATIAL GRID
      PARAMETER(NX=51)
      DIMENSION X(NX)
C...
C...  TIME GRID
      PARAMETER(NT=5)
      DIMENSION T(NT)
C...
C...  EXTERNAL FUNCTION TO CALCULATE INTEGRAND
C...
C...     FE   INTEGRAND OF ERF(X)
C...
      EXTERNAL FE
C...
C...  OPEN OUTPUT FILE
      NO=7
      OPEN(NO,FILE='output')
C...
C...  OPEN A FILE FOR MATLAB PLOTTING
      OPEN(1,FILE='h1.out')
C...
C...  HEAT CONDUCTION PROBLEM
C...
C...  DEFINE SPATIAL GRID
      XL=1.0D0
      DO 1 I=1,NX
         X(I)=XL*DFLOAT(I-1)/DFLOAT(NX-1)
1     CONTINUE
C...
C...  DEFINE TIME GRID
      TF=500.0D0
      DO 2 I=1,NT
         T(I)=TF*DFLOAT(I)/DFLOAT(NT)
2     CONTINUE
C...
C...  THERMAL DIFFUSIVITY, INITIAL, SURFACE TEMPERATURES
      ALPHA=0.0001D0
      T0= 25.0D0
      TS=250.0D0
C...
C...  PI
C...
C...  STEP THROUGH TIMES
      DO 10 I=1,NT
      PI=4.0D0*DATAN(1.0D0)
C...
C...  LABEL FOR TABULATED OUTPUT
      WRITE(NO,6)
6     FORMAT(4X,'I',4X,'J',6X,'T(I)',6X,'X(J)',
     +       5X,'ERFAB',5X,'ERFHA',8X,'T1')
C...
C...  STEP THROUGH SPACE
      DO 11 J=1,NX
C...
C...  LIMITS
      A=0.0D0
      B=X(J)/DSQRT(4.0D0*ALPHA*T(I))
C...
C...  EVALUATION OF ERF(X)
      NP=100
      ERFAB=SIMP(FE,A,B,NP)
C...
C...  INCLUDE FACTOR 2/SQRT(PI)
      ERFAB=2.0D0/DSQRT(PI)*ERFAB
C...
C...  TEMPERATURE
      T1=TS+(T0-TS)*ERFAB
C...
C...  ERF(X) FROM HASTINGS APPROXIMATION
      ERFHA=ERF(B)
C...
C...  MONITOR CALCULATION
      WRITE(*,*)I,J
C...
C...  WRITE NUMERICAL RESULTS
      WRITE(NO,4)I,J,T(I),X(J),ERFAB,ERFHA,T1
4     FORMAT(2I5,F10.1,F10.3,2F10.5,F10.2)
C...
C...  WRITE THE OUTPUT FOR MATLAB PLOTTING
      WRITE(1,5)X(J),T1
5     FORMAT(2F10.5)    
C...
C...  NEXT SPACE POINT
11    CONTINUE
C...
C...  NEXT TIME
10    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