program neville				!!   05-Apr-2006
	implicit real*8 (a-h,o-z)

c....... by Dario Mitnik

c....... Interpolation of a function Y(x) by using Lagrange polynomials
c........The polinomials are generated by Nevill''e method

c........           (x-x_{i-j+1})A_{i,j-1} - (x-x_i)A_{i-1,j-1}
c........A_{i,j} = ------------------------------------------
c........                       x_i - x_{i-j+1}


c....... files used:
c....... unit=15 -- table.dat (input function)
c....... unit=20 -- neville.dat (Legendre interpolation Polynomials)
	open(unit=15,file='table.dat',status='old')
	open(unit=20,file='neville.dat',status='unknown')
	
	call readinput
		
	call interpol
	
	stop
	end	
c	
c************************************************************************
c
	subroutine readinput
	implicit real*8 (a-h,o-z)
	
c....... input parameters, arrays, dimensions, etc.

	parameter (mxpts=100)
	common/bkinput/x(mxpts),y(mxpts),npts
	
	npts = 5
c....... here should be the data readed in a file, or a function generator
	do 100 i=1,npts
	  read(15,*,end=2000) x(i),y(i)
100	continue

2000	return
	end
c	
c************************************************************************
c
	subroutine interpol
	implicit real*8 (a-h,o-z)

	parameter (mxpts=100)
	common/bkinput/x(mxpts),y(mxpts),npts	
	data rzero/0.0d0/
	
100	print*,'give the r point for interpolation (le.x(1) for quit):'
	read*,r
	if (r.lt.x(1)) return

	result = rzero
	result = rNeville(npts,x,y,r)

	print777, r,result
777	format('the interpolated value at r=: ',1pg8.3,' is=: ',1pg14.7)		
	go to 100
	
	return
	end
c	
c************************************************************************
c
	double precision function rNeville(n,x,y,r)
	implicit real*8 (a-h,o-z)

c.......Neville''s method for interpolation

c........           (x-x_{i-j+1})A_{i,j-1} - (x-x_i)A_{i-1,j-1}
c........A_{i,j} = ------------------------------------------
c........                       x_i - x_{i-j+1}


	dimension A(n,n),x(n),y(n)
	data rzero,one/0.0d0,1.0d0/
	
c....... First column	
	do 100 irow=1,n	
		A(irow,1) = y(irow)
100	continue
	write(20,555) A(1,1)

c....... Neville algorithm
	 do 300 irow=2,n
	   do 200 jcol=2,irow
	     rnum1 = (r-x(irow-jcol+1))*A(irow,jcol-1)
	     rnum2 = (r-x(irow))*A(irow-1,jcol-1)
	     rden = x(irow)- x(irow-jcol+1)
	     A(irow,jcol) = (rnum1-rnum2)/rden	
200	   continue
	   write(20,555) (A(irow,j),j=1,irow)
300	continue
555	format(5(1pg14.7,2x))	


	rNeville = A(n,n)
	
	return
	end	
c	
c************************************************************************
c