function bess1d(x,f1,f2,f3,f4)
c
c=======================================================================
c                                                                    ===
c  This function performs a 1-d  cubic Bessel interpolation at       ===
c  the point X given the value of function at four consecutive       ===
c  points (F1,F2,F3,F4).   The value of the function at X lies       ===
c  between F2 and F3.                                                ===
c                                                                    ===
c     On Input:                                                      ===
c                                                                    ===
c     F1,F2,F3,F4  field to interpolate from (real)                  ===
c     X            position at which the field is interpolated (grid ===
c                  units; real)                                      ===
c                                                                    ===
c     On Output:                                                     ===
c                                                                    ===
c     BESS1D       interpolated value (real)                         ===
c                                                                    ===
c     Calls:          none                                           ===
c                                                                    ===
c=======================================================================
c
c-----------------------------------------------------------------------
c  Define global data.
c-----------------------------------------------------------------------
c
#include <cdefs.h>
#include <pconst.h>
c
c-----------------------------------------------------------------------
c  Define local data.
c-----------------------------------------------------------------------
c
      FLOAT bess1d
      FLOAT
     *      a0,a1,a2,a3,f1,f2,f3,f4,fp1,x
c
c-----------------------------------------------------------------------
c  Begin executable code.
c-----------------------------------------------------------------------
c
      fp1=p5*(f4-f2)
      a0=f2
      a1=p5*(f3-f1)
      a2=c3*f3-fp1-c3*a0-c2*a1
      a3=f3-a0-a1-a2
c
c-----------------------------------------------------------------------
c  Interpolate F(X) with a cubic polynomial.
c-----------------------------------------------------------------------
c
      bess1d=a0+x*(a1+x*(a2+a3*x))
c
      return
      end