program lagrange				!!   04-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....... Y(x) is tabulated at n points x_i, and L_{n,k)(x) are the n 
c....... Lagrange polynomials

c....... 	Y(x) = Sum Y(x_i)L_{ni}(x)	

c....... files used:
c....... unit=15 -- table.dat (input function)
c....... unit=20 -- Lpoly.dat (Legendre Polynomials)
	open(unit=15,file='table.dat',status='old')
	open(unit=20,file='Lpoly.dat',status='unknown')
	
	call readinput
	
	call printpolys
	
	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 printpolys
	implicit real*8 (a-h,o-z)

c....... Plot the Legendre Polynomials at nptsplot points

	parameter (mxpts=100)
	common/bkinput/x(mxpts),y(mxpts),npts	
	
	nptsplot=300
	rdelta = (x(npts)-x(1))/nptsplot
	r = x(1)
	do 100 i=1,nptsplot
	 write(20,*)  r,( rLpoly(npts,k,x,r), k=1,npts )
	 r = r + rdelta
100	continue

	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
	do 200 i=1,npts
	   result = result + y(i)*rLpoly(npts,i,x,r)
200	continue

	print*,' the interpolated value at r=:',r,' is =:',result
	
	go to 100
	
	return
	end
c	
c************************************************************************
c
	double precision function rLpoly(n,k,x,r)
	implicit real*8 (a-h,o-z)

c....... Lagrange interpolating polynomial L_nk (r)
c....... Polynomial that is zero at every X_n points, 
c....... and one at X_k
c....... Should use for interpolate a function F(x) = Sum f(x_i)L_ni(x)	
	dimension x(n)
	data one/1.0d0/
	
	rLpoly = one
	do 100 i=1,n	
	  if (i.eq.k) go to 100
	  rLpoly = rLpoly*(r-x(i))/(x(k) - x(i))
100	continue
	
	return
	end	
c	
c************************************************************************
c