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

c....... by Dario Mitnik

c....... Interpolation of a function Y(x) with splines

c........uses free boundary conditions ( S''(0)=S''(n)=0 )

c....... files used:
c....... unit=15 -- curva.dat (input function)
c....... unit=20 -- neville.dat (Legendre interpolation Polynomials)
	open(unit=15,file='curva.dat',status='old')
	open(unit=20,file='splines.dat',status='unknown')
	
	call readinput
		
	call createsolvematrix
	
	call printpolynom
	
	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(0:mxpts),a(0:mxpts),npts
	
c....... here should be the data readed in a file, or a function generator
	npts = -1
	do 100 i=0,mxpts
	  read(15,*,end=2000) xx,yy
	  npts = npts + 1
	  x(i)=xx
	  a(i)=yy
100	continue

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

c....... create tridiagonal matrix 

c        rl1  rd1   ru1                                   hb1
c             rl2   rd2  ru2                              hb2
c                   rl3   rd3  ru3           =            hb3
c                         rl4   rd4  ru4                  hb3


c...... and solve by using Gauss Method

	parameter (mxpts=100)
	common/bkinput/x(0:mxpts),a(0:mxpts),npts
	common/bkspline/b(0:mxpts),c(0:mxpts),d(0:mxpts),h(0:mxpts)	
	common/bkmatrix/rupp(0:mxpts-1),rdiag(0:mxpts),HB(0:mxpts)

	data rzero,one,two,three/0.0d0,1.0d0,2.0d0,3.0d0/
		
c........set h_i = (x_i+1 - x_i)
	do 100 i=0,npts-1
		h(i) = x(i+1) - x(i)
100	continue

c....... set B matrix elements
	HB(0)=rzero
	HB(npts)=rzero
	do 200 i=1,npts-1
	    bb1 = three*(a(i+1) - a(i))/h(i)
	    bb2 = three*(a(i) - a(i-1))/h(i-1)
	    HB(i) = bb1 - bb2
200	continue

c........ set A matrix elements
c....... A at zero 
	rdiag(0) = one
	rupp(0) = rzero
	do 300 i=1,npts-1
		rdiag(i) = two*( h(i-1) + h(i) ) - h(i-1)*rupp(i-1)
		rupp(i) = h(i)/rdiag(i)
		HB(i) = ( HB(i) - h(i-1)*HB(i-1) ) / rdiag(i)
300	continue
c....... A at endpoints 
	rdiag(npts) = one
	c(npts) = rzero

c........ backsubstitution
	do 400 i=npts-1,0,-1
		c(i) = HB(i) - rupp(i)*c(i+1)
		bb1 = ( a(i+1) - a(i) )/h(i) 
		bb2 = h(i)*( c(i+1) + two*c(i) )/three 
		b(i) = bb1 - bb2
		d(i) = (c(i+1) - c(i) ) / (three* h(i) )
400	continue

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

c....... print polynom

c        S(x) = a_j + b_j*(x-x_j) + c_j*(x-x_j)^2 + d_j*(x-x_j)^3

	parameter (mxpts=100)
	common/bkinput/x(0:mxpts),a(0:mxpts),npts
	common/bkspline/b(0:mxpts),c(0:mxpts),d(0:mxpts),h(0:mxpts)	

	nptspernodes=10
	do 200 n=0,npts-1
		rsup = x(n+1)
		rinf = x(n)
		rdelta = (rsup-rinf)/nptspernodes
		r = rinf - rdelta
		NN= nptspernodes
		if (n.eq.npts-1) NN=NN+1
		do 100 i=1,NN
		  r = r + rdelta
		  rhh = r-x(n)
		  y = a(n) + b(n)*rhh + c(n)*rhh**2 + d(n)*rhh**3
		  write(20,*) r,' ',y
100		continue
200	continue

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