!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
!
! The modules Precision and Erf_Series is defined are the files
! precision.f90 and erf_series.f90, which should be compiled previously.
! One should then link in the produced module (.vo) files (with -I, not needed if in same directory)
! and the object files (if there were any procedures in the module).
! Alternatively you can use the precompiled (double precision) ones in /classes/phy201
! f90-vast -I/classes/phy201 /classes/phy201/erf_series.o erf_file.f90 -o erf_file.x
!
PROGRAM Erf_File
      USE Precision
      USE Erf_Series, ONLY: error, converged, ErfSeries
!
      ! We are in position to choose whether we want REAL or COMPLEX x:
      ! Comment one of these out:
      REAL (KIND=wp) :: x, erf ! x and erf(x)
      REAL (KIND=wp) :: x_min, x_max ! Range of plot
      INTEGER :: i, status, n_points ! I/O status
!
      OPEN (FILE="Erf.dat", UNIT=10, STATUS="UNKNOWN", ACTION="READ", IOSTAT=STATUS)
      IF (status /= 0) THEN
         WRITE (*,*) "Error opening file:", status
         STOP
      END IF
!
      UserInput: DO
         WRITE (UNIT=*, FMT="(A)", ADVANCE="NO") "Enter: epsilon, n_points, x_min and x_max: "
         READ (UNIT=*, FMT=*, IOSTAT=STATUS) error, n_points, x_min, x_max
         IF (status == 0) EXIT UserInput ! Do not exit until the user enters acceptable numbers
      END DO UserInput
!
      FileWrite: DO i = 1, n_points
         x = REAL (i-1, wp) / REAL (n_points-1, wp) * (x_max-x_min) + x_min ! We will plot the interval [x_min,x_max]
         WRITE (UNIT=10, FMT=*, IOSTAT=STATUS) x, ErfSeries (x)
         IF (status /= 0) EXIT FileWrite
      END DO FileWrite
!
      IF ( .NOT. converged) WRITE (UNIT=*, FMT=*) "Convergence was not achieved for some points!"
!
      CLOSE (UNIT=10, IOSTAT=STATUS)
      IF (status /= 0) THEN
         WRITE (*,*) "Error closing file:", status
         STOP
      END IF
!
END PROGRAM Erf_File
!