       subroutine bcui2d(x,y,f,xi,yi,fi,fxi,fyi,fxxi,fyyi,fxyi)
	  
c***********************************************************************
c      Barycentric Cubic Interpolation 2D 
c      Written by Orlando Camargo Rodriguez
c      orodrig@ualg.pt
c      Universidade do Algarve
c      Physics Department
c      Signal Processing Laboratory
c      Faro, 16/06/2011 at 18:00
c***********************************************************************

       integer*4 i,j
       
       real*8 f(4,4),a(4,4)
       real*8 x(4),y(4)
       real*8 px(4),py(4)
       real*8 sx(4),sy(4)
       real*8 qx(4),qy(4)
       real*8 xi,x1,x2,x3,x4
       real*8 yi,y1,y2,y3,y4
       real*8 fi,fxi,fyi
       real*8 fxxi,fyyi,fxyi

       x1 = x(1)
       x2 = x(2)
       x3 = x(3)
       x4 = x(4)
       
       y1 = y(1)
       y2 = y(2)
       y3 = y(3)
       y4 = y(4)
       
       px(1) = (x1-x2)*(x1-x3)*(x1-x4)
       px(2) = (x2-x1)*(x2-x3)*(x2-x4)
       px(3) = (x3-x1)*(x3-x2)*(x3-x4)
       px(4) = (x4-x1)*(x4-x2)*(x4-x3)
       
       py(1) = (y1-y2)*(y1-y3)*(y1-y4)
       py(2) = (y2-y1)*(y2-y3)*(y2-y4)
       py(3) = (y3-y1)*(y3-y2)*(y3-y4)
       py(4) = (y4-y1)*(y4-y2)*(y4-y3)
      
c      Yes, a(1,1) is zero, but using the do loops 
c      shortens the code. 

       do j = 1,4
       
          do i = 1,4
       
             a(j,i) = (f(j,i) - f(1,1))/( py(j)*px(i) )

          end do

       end do
       
       px(1) = (xi-x2)*(xi-x3)*(xi-x4)
       px(2) = (xi-x1)*(xi-x3)*(xi-x4)
       px(3) = (xi-x1)*(xi-x2)*(xi-x4)
       px(4) = (xi-x1)*(xi-x2)*(xi-x3)
       
       py(1) = (yi-y2)*(yi-y3)*(yi-y4)
       py(2) = (yi-y1)*(yi-y3)*(yi-y4)
       py(3) = (yi-y1)*(yi-y2)*(yi-y4)
       py(4) = (yi-y1)*(yi-y2)*(yi-y3)

       sx(1) = (xi-x2)*(xi-x3) + (xi-x2)*(xi-x4) + (xi-x3)*(xi-x4)
       sx(2) = (xi-x1)*(xi-x3) + (xi-x1)*(xi-x4) + (xi-x3)*(xi-x4)
       sx(3) = (xi-x1)*(xi-x2) + (xi-x1)*(xi-x4) + (xi-x2)*(xi-x4)
       sx(4) = (xi-x1)*(xi-x2) + (xi-x1)*(xi-x3) + (xi-x2)*(xi-x3)
       
       sy(1) = (yi-y2)*(yi-y3) + (yi-y2)*(yi-y4) + (yi-y3)*(yi-y4)
       sy(2) = (yi-y1)*(yi-y3) + (yi-y1)*(yi-y4) + (yi-y3)*(yi-y4)
       sy(3) = (yi-y1)*(yi-y2) + (yi-y1)*(yi-y4) + (yi-y2)*(yi-y4)
       sy(4) = (yi-y1)*(yi-y2) + (yi-y1)*(yi-y3) + (yi-y2)*(yi-y3)

       qx(1) = 3*xi - x2 - x3 - x4
       qx(2) = 3*xi - x1 - x3 - x4
       qx(3) = 3*xi - x1 - x2 - x4
       qx(4) = 3*xi - x1 - x2 - x3
       
       qy(1) = 3*yi - y2 - y3 - y4
       qy(2) = 3*yi - y1 - y3 - y4
       qy(3) = 3*yi - y1 - y2 - y4
       qy(4) = 3*yi - y1 - y2 - y3
       
       fi   = f(1,1)
       fxi  = 0.0
       fyi  = 0.0
       fxxi = 0.0
       fyyi = 0.0
       fxyi = 0.0

c      No patience for third order derivatives...

       do j = 1,4
       
          do i = 1,4
       
              fi  =  fi + a(j,i)*px(i)*py(j)
	     fxi =  fxi + a(j,i)*sx(i)*py(j)
	    fxxi = fxxi + a(j,i)*qx(i)*py(j)
             fyi =  fyi + a(j,i)*px(i)*sy(j)
	    fyyi = fyyi + a(j,i)*px(i)*qy(j)
	    fxyi = fxyi + a(j,i)*sx(i)*sy(j)
	    
          end do

       end do

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