module interp
	interface
		real function pointerp(x, y, z, near)
			real, dimension(:), intent(in) :: x, y
			real, intent(in) :: z
			integer, intent(in) :: near
		end function pointerp

		real function qspline(x, y, z)
			real, dimension(:), intent(in) :: x, y
			real, intent(in) :: z
		end function qspline
	end interface
end module

! Polonomial interpolation of tabulated data (x,y) evaluated in z.
real function pointerp(x, y, z, near)
	implicit none
	real, dimension(:), intent(in) :: x, y
	real, intent(in) :: z
	integer, intent(in) :: near
	integer :: i, k, n, low, high, cen
	real :: prod
	logical :: dohigh
	n = size(x)

	! Check that z lines within the range of x:
	if (z > maxval(x) .or. z < minval(x)) then
		stop "z must be inside the range of x."
	endif

	! Find the interval in wich z lies:
	low = 0
	high = n-1
	do while (high-low > 1)
		cen = nint(real(high+low)/2.0)
		if (z > x(cen)) then
			low = cen
		else
			high = cen
		endif
	enddo

	! Add points on either side until reached near:
	dohigh = .true.
	do while (high-low < near)
		if (dohigh .and. high < n) then
			high = high+1
		elseif (low > 1) then
			low = low-1
		endif
		dohigh = .not. dohigh
	enddo
	
	pointerp = 0.0
	do i = 1,n
		prod = 1.0
		do k = low,high
			if (i/=k) then
				prod = prod * (z-x(k))/(x(i)-x(k))
			endif
		enddo
		pointerp = pointerp + y(i) * prod
	enddo
end function pointerp

! Quadratic spline interpolation of tabulated data (x,y) evaluated in z.
real function qspline(x, y, z)
	implicit none
	real, dimension(:), intent(in) :: x, y
	real, intent(in) :: z
	real, dimension(size(x)) :: b, c, b_up, c_up
	integer :: i, n, low, high, cen
	real :: h
	n = size(x)

	! Check that z lines within the range of x:
	if (z > maxval(x) .or. z < minval(x)) then
		stop "z must be inside the range of x."
	endif

	! Recursion UP:
	b(1) = (y(2)-y(1))/(x(2)-x(1))
	c(1) = 0.0
	do i = 2,n-1
		h = x(i+1)-x(i)
		b(i) = b(i-1) + 2*c(i-1)*(x(i)-x(i-1))
		c(i) = (y(i+1) - y(i) - b(i)*h) / h**2
	enddo
	b_up = b; b = 0.0;
	c_up = c; c = 0.0;

	! Recursion DOWN:
	b(n-1) = (y(n)-y(n-1))/(x(n)-x(n-1))
	c(n-1) = 0.0
	do i = n-2,1,-1
		h = x(i+1)-x(i)
		b(i) = 2*(y(i+1)-y(i))/h - b(i+1)
		c(i) = (b(i+1)-b(i))/(2*h)
	enddo

	! Take the average between the UP and DOWN recursions:
	b = (b_up+b)/2.
	c = (c_up+c)/2.

	! Find the interval in wich z lies:
	low = 1
	high = n
	do while (high-low > 1)
		cen = nint(real(high+low)/2.0)
		if (z > x(cen)) then
			low = cen
		else
			high = cen
		endif
	enddo

	! Calculate the polonomial:
	qspline = y(low) + b(low)*(z-x(low)) + c(low)*(z-x(low))**2
end function qspline