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

c....... by Dario Mitnik

c....... Interpolation of a function Y(x) by using Divided-Difference 

c....P_n(x) = f[x_1] + Sum_{k=1}^{n+1} f[x_1,x_2,...,x_k](x - x1)..(x - x_{k-1})


c.....              f(i,j-1) - f(i-1,j-1)
c....... f(i,j) = ------------------------
c.......              x(i) - x(i-j+1)


c....... files used:
c....... unit=15 -- table.dat (input function)
c....... unit=20 -- divided.dat (Legendre interpolation Polynomials)
	open(unit=15,file='table.dat',status='old')
	open(unit=20,file='divided.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),A(mxpts,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),A(mxpts,mxpts),npts	
	data rzero,one/0.0d0,1.0d0/
	
100	print*,'give the r point for interpolation (le.x(1) for quit):'
	read*,r
	if (r.lt.x(1)) return

	call rdividiff

	result = A(1,1)
	prod = one
	do 200 i=2,npts
	  prod = prod*(r - x(i-1) )
	  result = result + A(i,i)*prod
200	continue

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

c.......Divided-Difference method for interpolation

c.....              A(i,j-1) - A(i-1,j-1)
c....... A(i,j) = ------------------------
c.......              x(i) - x(i-j+1)

	parameter (mxpts=100)
	common/bkinput/x(mxpts),y(mxpts),A(mxpts,mxpts),npts	

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

c....... dividiff algorithm
	 do 300 irow=2,npts
	   do 200 jcol=2,irow
	     rnum1 = A(irow,jcol-1)
	     rnum2 = 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))	


	return
	end	
c	
c************************************************************************
c