C
C  program deriv_prec
C
C  program to illustrate different dependence of the precision of 
C  numerical derivatives on step size (dx)
C
C  compares left-asymmetric, right-asymmetric, and symmetric 2-point formulas
C
C  By:  CWJ @ SDSU , September 2009
C
C  SUBROUTINES CALLED:
C
C  FUNCTIONS CALLED
C    df(x) = exact derivative (analytic)
 
      implicit none

      real x     !  variable
      real exactderiv    
      real approxderiv
      real dx
      real df
      real maxdx, mindx, dlogdx

      integer n,nsteps
C............................................

      print*,' Enter value of x '
      read*,x
      print*,' Enter max, min values of dx and number of points wanted '
      read*,maxdx,mindx,nsteps
C.............. ERROR TRAP....................
      if( maxdx .lt. mindx .or. mindx .le. 0.0)then
         print*,' wrong order ',maxdx,mindx
         stop
      endif
C................END ERROR TRAP..............

      dlogdx = (log(maxdx)-log(mindx))/float(nsteps-1)
      exactderiv = df(x)

      do n = 1,nsteps
         dx = mindx* exp( (n-1)* dlogdx)    ! logarithmic spacing
         call derivative_2pt_asym_left(x,dx,approxderiv)
         print*,dx,approxderiv
         write(17,*)dx,abs(approxderiv-exactderiv)  
         write(27,*)log(dx),log(abs(approxderiv-exactderiv))  
  

         call derivative_2pt_sym(x,dx,approxderiv)
         print*,dx,approxderiv
         write(18,*)dx,abs(approxderiv-exactderiv)      
         write(28,*)log(dx),log(abs(approxderiv-exactderiv))      

      enddo   ! loop over n
            

      end

C=====================================================
C
C  subroutine to compute numerical derivative 
C  via asymmetric 2-point function 
C  to the left
C
C  FUNCTIONS CALLED:
C    f(x) : the function whose derivative we need
C

       
      subroutine derivative_2pt_asym_left(x,dx,dfdx)

      implicit none

      real x   ! value where derivative is taken
      real dx  ! step size
     
      real dfdx  ! output derivative

      real f   ! function; must be defined outside

C------------ ERROR TRAP ------ check dx > 0
      if( dx .le. 0.0)then
        print*,' dx = ',dx
        stop
      endif
C------------END ERROR TRAP ---------------
      dfdx =( f(x+dx) - f(x) ) / dx

      return
      end   subroutine derivative_2pt_asym_left
C----------------------------------------------------
C
C  subroutine to compute numerical derivative 
C  via symmetric 2-point function 
C
C  FUNCTIONS CALLED:
C    f(x) : the function whose derivative we need
C

       
      subroutine derivative_2pt_sym(x,dx,dfdx)

      implicit none

      real x   ! value where derivative is taken
      real dx  ! step size
     
      real dfdx  ! output derivative

      real f   ! function; must be defined outside

C------------ ERROR TRAP ------ check dx > 0
      if( dx .le. 0.0)then
        print*,' dx = ',dx
        stop
      endif
C------------END ERROR TRAP ---------------
      dfdx =( f(x+dx) - f(x-dx) ) / (2.*dx)

      return
      end   subroutine derivative_2pt_sym

C function f
C
C 
      real function f(x)
      real x

      f = sin(x)   ! note we can change this! 
      return
      end
C----------------------------------------
      real function df(x)
      real x
      df = cos(x)
      return
      end