       subroutine bcui1d(x,f,xi,fi,fxi,fxxi)
	  
c***********************************************************************
c      Barycentric Cubic Interpolation 1D 
c      Written by Orlando Camargo Rodriguez
c      orodrig@ualg.pt
c      Universidade do Algarve
c      Physics Department
c      Signal Processing Laboratory
c      Faro, 23/09/2009 at 13:00
c****&******************************************************************
       
       integer*4 i
       real*8 x(4),f(4)
       real*8 a(3),px(3),sx(3),qx(3)
       real*8 xi,fi,fxi,fxxi
       real*8 x1,x2,x3,x4
       
c***********************************************************************
       
       x1 = x(1)
       x2 = x(2)
       x3 = x(3)
       x4 = x(4)

       px(1) = ( x2 - x1 )*( x2 - x3 )*( x2 - x4 )
       px(2) = ( x3 - x1 )*( x3 - x2 )*( x3 - x4 )
       px(3) = ( x4 - x1 )*( x4 - x2 )*( x4 - x3 )
       
       do i = 1,3
       
          a(i) = ( f(i+1) - f(1) )/px(i)
       
       end do

       px(1) = ( xi - x1 )*( xi - x3 )*( xi - x4 )
       px(2) = ( xi - x1 )*( xi - x2 )*( xi - x4 )
       px(3) = ( xi - x1 )*( xi - x2 )*( xi - x3 )
       
       sx(1) = (xi-x1)*(xi-x3)+(xi-x1)*(xi-x4)+(xi-x3)*(xi-x4)
       sx(2) = (xi-x1)*(xi-x2)+(xi-x1)*(xi-x4)+(xi-x2)*(xi-x4)
       sx(3) = (xi-x1)*(xi-x2)+(xi-x1)*(xi-x3)+(xi-x2)*(xi-x3)
       
       qx(1) = 2*( 3*xi - x1 - x3 - x4 )
       qx(2) = 2*( 3*xi - x1 - x2 - x4 )
       qx(3) = 2*( 3*xi - x1 - x2 - x3 )
       
       fi   = f(1)
       fxi  = 0.0
       fxxi = 0.0
       
       do i = 1,3
       
          fi   = fi   + a(i)*px(i)
          fxi  = fxi  + a(i)*sx(i)
          fxxi = fxxi + a(i)*qx(i)
       
       end do

c***********************************************************************
c      Back to main:
c***********************************************************************        
       
       return 
       
       end
