!___________________________________________________________________________________________________
! Aleksandar Donev
! PHY201 Solutions to Worksheet
! December, 2000, MSU
!___________________________________________________________________________________________________
PROGRAM Erf_Calculator
    USE Precision, ONLY : wp
    USE Erf_Series, ONLY: error, converged, ErfSeries
    USE Erf_Precise, ONLY: ErfPrecise

    ! We are in position to choose whether we want REAL or COMPLEX x:
    ! Comment one of these out:
    ! Real:
	! REAL(KIND=wp) :: x, erf1,erf2 ! x and erf(x) are REAL this time  
	! CHARACTER(*), PARAMETER :: format="(A,D20.15,A,D20.15)" ! We use this for dp REAL in g77	
	! Complex:
	COMPLEX(KIND=wp) :: x, erf1,erf2 ! x and erf(x) are COMPLEX this time  
	CHARACTER(*), PARAMETER :: format= &  ! We use this for dp COMPLEX in g77		
	            "(A,'(',D25.15,',',D25.15,')',A,'(',D25.15,',',D25.15,')')"
   INTEGER :: status ! I/O status
   
   MainLoop: DO
   
      WRITE(UNIT=*,FMT="(A)",ADVANCE="NO") "Enter the value of x and epsilon ('S' to stop): "
      READ(UNIT=*,FMT=*,IOSTAT=status) x,error      
      IF(status/=0) EXIT MainLoop   

      erf1=ErfSeries(x)
      IF(.NOT.converged) WRITE(UNIT=*,FMT=*) "Convergence was not achieved! Last known answer:"
      WRITE(UNIT=*,FMT=format) "Series calculation: Erf (",x,") = ",erf1     

      erf2=ErfPrecise(x)
      WRITE(UNIT=*,FMT=format)   "Precise result:     Erf (",x,") = ",erf2
      WRITE(UNIT=*,FMT=*)   "Error is: ",ABS(erf1-erf2)
      
   END DO MainLoop               

END PROGRAM Erf_Calculator
!