PROGRAM SubEvilError

! ------------------------------------------------------------------------------
!
! Copyright 2006 Matthew Norton
!
! ------------------------------------------------------------------------------
!    This program is free software; you can redistribute it and/or modify
!    it under the terms of the GNU General Public License as published by
!    the Free Software Foundation; either version 2 of the License, or
!    (at your option) any later version.
!
!    This program is distributed in the hope that it will be useful,
!    but WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!    GNU General Public License for more details.
!
!    You should have received a copy of the GNU General Public License
!    along with this program; if not, write to the Free Software
!    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
!
! NOTE: This is in FORTRAN 90 and you will need a FORTRAN 90 complier to run 
! this code
! ------------------------------------------------------------------------------
!
! This program calculates the machine epislon, the smallest number that, when 
! added to 1, gives a number greater than 1. This is useful in finding out the 
! precision of the computer you run your FORTRAN program on. 
!
! ------------------------------------------------------------------------------
!
! PROGRAM HISTORY
!
! AUTHOR: Matthew S. Norton
! CREATION DATE: January 25, 2007
!
! ------------------------------------------------------------------------------

REAL*8 :: YEvil, YGood, X, Error
CHARACTER*20 :: FileName

PRINT *, "What do you want the output file to be (include extension)?"
READ *, FileName
OPEN (20, File = TRIM(FileName))
WRITE(20, '(1X, A, T21, A, T34, A, T51, A)') "X", "YGood", "YEvil", "Error"  
DO I = 1, 100
   X = 10 ** (1 + ((I - 1.) / 99.) * 5)
   PRINT *, X
   YEvil = SQRT(X ** 2 + 1.) - X
   YGood = 1. / (SQRT(X ** 2 + 1.) + X)
   Error = ABS(YEvil - YGood) / YGood
   WRITE(20, *) X, YGood, YEvil, Error
   PRINT *, X, YGood, YEvil, Error
END DO

!100 FORMAT (1X, F20.10, 1X, F15.13, 1X, F15.13, 1X, E20.15E4)

END PROGRAM