C********************************************************************************
C*        Spline interface routines                                             *
C*                                                                              *
C*	Date: 06/26/06                                                        *
C*	Washington, DC                                                        *
C*	Jan Vorberger                                                         *   
C********************************************************************************

	REAL*8 FUNCTION splint2(x,y,n,yp1,ypn,xx) 
	INTEGER n,NMAX
	PARAMETER (NMAX=500)
	REAL*8 yp1,ypn,x(n),y(n),y2(NMAX),xx,yy
	call spline(x,y,n,yp1,ypn,y2)
	call splint(x,y,y2,n,xx,yy)
	splint2 = yy
	RETURN
	END

	REAL*8 FUNCTION splintnat(x,y,n,xx) 
	INTEGER n,NMAX
	PARAMETER (NMAX=500)
	REAL*8 yp1,ypn,x(n),y(n),y2(NMAX),xx,yy
	yp1 = 1.d30
	ypn = 1.d30
	call spline(x,y,n,yp1,ypn,y2)
	call splint(x,y,y2,n,xx,yy)
	splintnat = yy
	RETURN
	END

C*****Spline calculation*********************************************************
C*	x(1:n)...argument array
C*	y(1:n)...function array
C*	n...number of field points
C*	yp1,ypn...first derivatives at beginning and end of array
C*	y2(1:n)...2nd derivates of interpolating function (output)

	SUBROUTINE spline(x,y,n,yp1,ypn,y2)
	INTEGER n,NMAX
	REAL*8 yp1,ypn,x(n),y(n),y2(n)
	PARAMETER (NMAX=500)
	INTEGER i,k
	REAL*8 p,qn,sig,un,u(NMAX)
	if (yp1.gt..99e30) then
	   y2(1)=0.
	   u(1)=0.
	else
	   y2(1)=-0.5
	   u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
	endif
	do 11 i=2,n-1
	   sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
	   p=sig*y2(i-1)+2.
	   y2(i)=(sig-1.)/p
	   u(i)=(6.*((y(i+1)-y(i))/(x(i+
     *  1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
     *  u(i-1))/p
 11	continue
	if (ypn.gt..99e30) then
	   qn=0.
	   un=0.
	else
	   qn=0.5
	   un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
	endif
	y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
	do 12 k=n-1,1,-1
	   y2(k)=y2(k)*y2(k+1)+u(k)
 12	continue
	return
	END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.

C*****Spline Interpolation*******************************************************
C*	xa(1:n)...argument array
C*	ya(1:n)...function array
C*	y2(1:n)...2nd derivates of interpolating function
C*	n...number of field points
C*	x...argument for desired function value
C*	y...(output) value at point x (by spline interpolation)

	SUBROUTINE splint(xa,ya,y2a,n,x,y)
	INTEGER n
	REAL*8 x,y,xa(n),y2a(n),ya(n),eps
	INTEGER k,khi,klo
	REAL*8 a,b,h

	eps = 0.1d0
	if (x.lt.xa(1)-eps*(xa(2)-xa(1))) then
	   write(*,*) " x-value below spline x-range ",xa(1),x,xa(n),n
	endif
	if (x.gt.xa(n)+eps*(xa(n)-xa(n-1))) then
	   write(*,*) " x-value above spline x-range ",xa(1),x,xa(n),n
	endif

	klo=1
	khi=n
 1	if (khi-klo.gt.1) then
	   k=(khi+klo)/2
	   if(xa(k).gt.x)then
	      khi=k
	   else
	      klo=k
	   endif
	   goto 1
	endif
	h=xa(khi)-xa(klo)
	if (h.eq.0.) then
	   write(*,*) 'bad xa input in splint'
	   stop
	endif
	a=(xa(khi)-x)/h
	b=(x-xa(klo))/h
	y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
     *  2)/6.
	return
	END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.
