!     Last change:  DOS  29 Jul 2000    4:52 pm
!                                       *** copyright 2000 ***
! *** filename numdif.f95               *** John F. Monahan **
!                                       **********************
      program numdif
!  example of numerical differentiation
      implicit none
      real algama
      real x,h,fx,fxph,fxmh,fp1,fp2,fpp1,fpp2
      integer k
!
  21  format(i4,4f15.8)
  22  format('   k',8x,'f1''',12x,'f2''',12x,'f1"',12x,'f2"')
!
      open( unit=6, file='numdif.out' )
!
!      numerical differentiation of the log of the gamma function
!       use x = 1/2 so derivative is 1.96351002
!        second derivative is 4.9348022005
!
      x = .5
      fx = algama(x)
!                    start with h = 1/4 -- huge
      h = .5 
      write(6,22) 
      do  k = 2,30
      h = h/2.
      fxph = algama(x+h)
      fxmh = algama(x-h) 
!                    one derivative as forward difference
      fp1 = (fxph - fx)/h
!                    other derivative as central difference
      fp2 = (fxph - fxmh) / (2.*h)
!                    second derivative -- simple formula
      fpp1 = (fxph - 2.*fx + fxmh)/(h*h)
!                    other formula difference of differences
      fpp2 = ( (fxph - fx) - (fx - fxmh) )/(h*h)
!                    write them out -- 2's should be better
      write(6,21) k,fp1,fp2,fpp1,fpp2
      end do  ! loop on k
      stop
      end program numdif
      REAL FUNCTION ALGAMA(X)
!  COMPUTES NATURAL LOG OF GAMMA FUNCTION FOR POSITIVE ARGUMENTS
!   X .LE. 0.  RETURNS -1.0
      IMPLICIT NONE
      REAL, INTENT(IN) :: X
      REAL, DIMENSION(5) :: P,Q         ! COEFS FOR RATIONAL FUNCTION
      REAL, DIMENSION(3) :: R           ! COEFS FOR STIRLING CORRECTION
      REAL XN,A,V
      REAL, PARAMETER :: LSQR2PI = .918938533
      DATA P/ -.510046056E1, -.368309363E2, -.664664448E2,  &
     &        -.195692802E3, -.185072903E1 /
      DATA Q/ -.114321842E2,  .368848969E2,  .16269403E2,   &
     &        -.195692802E3, 1.0 /
      DATA R/ -.277765545E-2, .833333330E-1, .77783067E-3 /      
!   USES METHODS 5230,5401 FROM 'COMPUTER APPROXIMATIONS'
!   ED. BY JF HART,  SIAM SERIES,  WILEY(1968)
!   27 MAR 74    JFM W/RWS,RM
!                         ** MODIFICATIONS   JUNE 1988, 1993 **
!    RECODED OCTOBER 1999 FOR FORTRAN 95
      IF( X .LE. 0. ) THEN
      ALGAMA = -1.0            ! NOT DESIGNED FOR NEGATIVE ARGUMENTS
      RETURN
                        ENDIF
!                          FOR LARGE X (>8) USE STIRLING ROUTE
      IF( X .GT. 8. ) THEN
!                           LARGE X -- USE STIRLING WITH CORRECTION
      XN = X * X
!                           CORRECTION FOR STIRLING
      V = R(2) + ( R(3)/XN + R(1) )/XN                                  
      ALGAMA = (X-0.5)*LOG(X) - X + V/X + LSQR2PI
                      ELSE
!                          FOR SMALL X, USE REPRODUCTION FORMULA
!                           TO GET NUMBER BETWEEN 2 AND 3
      V = 1.0                                                             
      XN = X                                                              
      DO WHILE( XN .LT. 2. )
      V = V/XN
      XN = XN + 1.
      END DO  ! WHILE( XN .LT. 2. )
      DO WHILE( XN .GT. 3. )
      XN = XN - 1.
      V = V * XN
      END DO  ! WHILE( XN .GT. 3. )
      A = XN-2.0
!                           RATIONAL FUNCTION APPROXIMATION
      V = V*(P(4)+A*(P(3)+A*(P(2)+A*(P(1)+A*P(5))))) /  &
     &   (Q(4)+A*(Q(3)+A*(Q(2)+A*(Q(1)+A*Q(5)))))
      ALGAMA = LOG(V)
                      ENDIF
      RETURN
      END FUNCTION ALGAMA
! *** end of filename numdif.f95        *********************