       subroutine bpai2d(x,y,f,xi,yi,fi,fxi,fyi,fxxi,fyyi,fxyi)
	  
c***********************************************************************
c      Barycentric Parabolic 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, 12/04/2009 at 01:00
c****&******************************************************************
       
       integer*4 i,j
       real*8 f(3,3),a(3,3)
       real*8 x(3),y(3),px(3),py(3),sx(3),sy(3)
       real*8 fi,fxi,fyi,fxxi,fyyi,fxyi
       real*8 xi,x1,x2,x3
       real*8 yi,y1,y2,y3
       
c***********************************************************************
       
       x1 = x(1)
       x2 = x(2)
       x3 = x(3)
       
       y1 = y(1)
       y2 = y(2)
       y3 = y(3)

       px(1) = ( x1 - x2 )*( x1 - x3 )
       px(2) = ( x2 - x1 )*( x2 - x3 )
       px(3) = ( x3 - x1 )*( x3 - x2 )
       
       py(1) = ( y1 - y2 )*( y1 - y3 )
       py(2) = ( y2 - y1 )*( y2 - y3 )
       py(3) = ( y3 - y1 )*( y3 - y2 )       
       
       do i = 1,3
       
          do j = 1,3
       
          a(i,j) = f(i,j)/( px(j)*py(i) )
       
          end do
       
       end do

       px(1) = ( xi - x2 )*( xi - x3 )
       px(2) = ( xi - x1 )*( xi - x3 )
       px(3) = ( xi - x1 )*( xi - x2 )
       
       py(1) = ( yi - y2 )*( yi - y3 )
       py(2) = ( yi - y1 )*( yi - y3 )
       py(3) = ( yi - y1 )*( yi - y2 ) 
       
       sx(1) = 2.0*xi - x2 - x3
       sx(2) = 2.0*xi - x1 - x3
       sx(3) = 2.0*xi - x1 - x2
       
       sy(1) = 2.0*yi - y2 - y3
       sy(2) = 2.0*yi - y1 - y3
       sy(3) = 2.0*yi - y1 - y2     
      
       fi   = 0.0
       fxi  = 0.0
       fyi  = 0.0
       fxxi = 0.0
       fyyi = 0.0
       fxyi = 0.0
       
       do i = 1,3
       
          do j = 1,3
       
          fi   = fi   + a(i,j)*px(j)*py(i)
          fxi  = fxi  + a(i,j)*sx(j)*py(i)
          fyi  = fyi  + a(i,j)*px(j)*sy(i)
          fxxi = fxxi + a(i,j)* 2.0 *py(i)
          fyyi = fyyi + a(i,j)*px(j)* 2.0
          fxyi = fxyi + a(i,j)*sx(j)*sy(i) 
       
          end do
       
       end do

       return
       
       end

c***********************************************************************

