program hermite
c
c    abcd:hermite <-------------- hermite interpolation of Airy function
c                                                         february, 2004
c
c    written by: A Busy Code Developer
c    modified 1: 
c
c  PURPOSE:  This program performs a Hermite interpolation to the Bessel
c  function J1(x).
c
c  EXTERNALS:
c    HERINT
c
c-----------------------------------------------------------------------
c
       implicit      none
c
       integer       ndim    , idim    
       real*8        err     
       parameter     ( ndim=11, idim=101, err=1.0d-6 )
c
       integer       n       , i
       real*8        dx      , xint    , j1int
c   
       real*8        x       (ndim), j0     (ndim), j1      (ndim)
     1             , j2      (ndim), j1p    (ndim)
c
       real*8        herint
       external      herint   
c
       data          x
     1             / 0.0000000000d0, 1.0000000000d0, 2.0000000000d0
     2             , 3.0000000000d0, 4.0000000000d0, 5.0000000000d0
     3             , 6.0000000000d0, 7.0000000000d0, 8.0000000000d0
     4             , 9.0000000000d0, 1.0000000000d1 /
       data          j0
     1             / 1.0000000000d0, 0.7651976866d0, 0.2238907791d0
     2             ,-0.2600519549d0,-0.3971498099d0,-0.1775967713d0
     3             , 0.1506452573d0, 0.3000792705d0, 0.1716508071d0
     4             ,-0.0903336112d0,-0.2459357645d0 /
       data          j1
     1             / 0.0000000000d0, 0.4400505857d0, 0.5767248078d0
     2             , 0.3390589585d0,-0.0660433280d0,-0.3275791376d0
     3             ,-0.2766838581d0,-0.0046828235d0, 0.2346363469d0
     4             , 0.2453117866d0, 0.0434727462d0 /
       data          j2
     1             / 0.0000000000d0, 0.1149034849d0, 0.3528340286d0
     2             , 0.4860912606d0, 0.3641281459d0, 0.0465651163d0
     3             ,-0.2428732100d0,-0.3014172201d0,-0.1129917204d0
     4             , 0.1448473415d0, 0.2546303137d0 /
c
c-----------------------------------------------------------------------
c
       dx   = 10.0d0 / dble ( idim - 1 )
c
c      Open disc files "herraw.dat" and "herint.dat", and write headers.
c
       open  ( 10, file='herraw.dat', status='unknown' )
       write ( 10, 2010 )
       open  ( 11, file='herint.dat', status='unknown' )
       write ( 11, 2020 )
c
c      Compute d(J1)/dx.
c
       do 10 n=1,ndim
         j1p(n) = 0.5d0 * ( j0(n) - j2(n) )
10     continue
c
c      Tabulate the uninterpolated data.
c
       do 20 n=1,ndim
         write ( 10, 2030 ) x(n), j1(n)
20     continue
c
c      Tabulate interpolated data.
c
       xint = 0.0d0
       do 30 i=1,idim
         j1int  = herint ( xint, x, j1, j1p, ndim )
         write ( 11, 2030 ) xint, j1int
         xint = xint + dx
30     continue
c
c      Close data files.
c
       close ( 10 )
       close ( 11 )
c
       stop
c
c-----------------------------------------------------------------------
c----------------------- Write format statements -----------------------
c-----------------------------------------------------------------------
c
2010   format('#  Datafile containing the uninterpolated Bessel '
     1       ,'function J1(x).',/,'#',/
     2       ,'#        x                 J1(x)',/,'#')
2020   format('#  Datafile containing the interpolated Bessel '
     1       ,'function J1(x).',/,'#',/
     2       ,'#        x                 J1(x)',/,'#')
2030   format(1pg17.10,2x,g17.10)
c
       end
c
c=======================================================================
c
       function herint ( xint, x, fofx, dfdx, ndim )
c
c    abcd:hermite.herint <----------------------- performs interpolation
c                                                            march, 2002
c
c    written by: A Busy Code Developer
c    modified 1: 
c
c  PURPOSE:  This routine performs a Hermite interpolation on the input 
c  function at the specified location.  A zeroth order extrapolation is 
c  returned if "xint" falls outside the domain of "x".
c
c  EXTERNALS: [none]
c
c-----------------------------------------------------------------------
c
       implicit      none
       integer       ndim    , n       , n0      , n0p1
       real*8        herint  , xint    , x1      , x2      , f1      
     1             , f2      , fp1     , fp2     
c
       real*8        x       (ndim), fofx    (ndim), dfdx    (ndim)
c
c-----------------------------------------------------------------------
c
c      If "xint" is outside domain of "x", return zeroth order 
c  extrapolation.
c
       if ( xint .lt. x(1) ) then
         write(6, 2010) xint, x(1), x(ndim)
         herint = fofx(1)
         return
       endif
       if ( xint .gt. x(ndim) ) then
         write(6, 2010) xint, x(1), x(ndim)
         herint = fofx(ndim)
         return
       endif
c
c      Determine in which interval "xint" falls.
c
       do 10 n=1,ndim-1
         if ( xint .ge. x(n) ) n0 = n
10     continue
c
c      Perform the Hermite interpolation.
c
       n0p1   = n0 + 1
       x1     = x   (n0  )
       x2     = x   (n0p1)
       f1     = fofx(n0  )
       f2     = fofx(n0p1)
       fp1    = dfdx(n0  )
       fp2    = dfdx(n0p1)
       herint = ( 3.0d0*x1 - 2.0d0*xint - x2 ) * ( xint - x2 )**2 * f1
     1        / ( x1 - x2 )**3
     2        + ( 3.0d0*x2 - 2.0d0*xint - x1 ) * ( xint - x1 )**2 * f2
     3        / ( x2 - x1 )**3
     4        + ( xint - x1 ) * ( xint - x2 )**2 * fp1 / ( x1 - x2 )**2
     5        + ( xint - x2 ) * ( xint - x1 )**2 * fp2 / ( x2 - x1 )**2
c
       return
c
c-----------------------------------------------------------------------
c----------------------- Write format statements -----------------------
c-----------------------------------------------------------------------
c
2010   format('HERINT  : *** WARNING *** Requested value for "xint" '
     1       ,1pg12.5,' is outside',/
     2       ,'HERINT  : the data domain (',g12.5,',',g12.5,').  A '
     3       ,'zeroth order',/
     3       ,'HERINT  : extrapolation will be returned for "fint".')
c
       end