       subroutine cbcu1d(x,z,xi,zi,zxi,zxxi)
	  
c***********************************************************************
c      Complex 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, 08/10/2010 at 18:30
c****&******************************************************************
       
       integer*4 i

       real*8 x(4)
       real*8 px(3),sx(3),qx(3)
       real*8 xi,x1,x2,x3,x4

       complex*16 z(4)
       complex*16 a(3)
       complex*16 zi,zxi,zxxi 

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) = ( z(i+1) - z(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 )
       
       zi   = z(1)
       zxi  = (0.0,0.0)
       zxxi = (0.0,0.0)
       
       do i = 1,3
       
          zi   = zi   + a(i)*px(i)
          zxi  = zxi  + a(i)*sx(i)
          zxxi = zxxi + a(i)*qx(i)
       
       end do
 
c***********************************************************************
c      Back to main:
c***********************************************************************        
       
       return 
       
       end
