program diff_sine4
C
C  a 4th program to compare numerical differentiation
C  (somehwat) better documentation
C  compare asymmetric 2-pt formula against 
C  symmetric 3-pt formula
C  CWJ --SDSU - August 2006
C  

      implicit none

      real x    ! argument (in radians)
      real dx   ! delta x for numerical derivative 
      real dfdx  ! derivative of sine

      integer i  ! counter
      
100   continue
      print*,' enter x '
      read*,x

 
      dx = 1.0

      do i = 1,6
        dx = dx / sqrt(10.)

C........... 2 -point formula      
         dfdx = (sin(x+dx)-sin(x))/dx
         print*,dfdx,cos(x)
         write(8,*)dx,abs(dfdx-cos(x))
C............3 - point formula
         dfdx = (sin(x+dx) - sin(x-dx))/(2*dx)
         write(10,*)dx,abs(dfdx-cos(x))

         write(9,*)x,cos(x)
      enddo
      print*,' Approx written to fort.8 '
      print*,' Exact written to fort.9 '
      print*,' Improved approx to fort.10 '

      end