module quadspline
implicit none

type qspline
real,allocatable,dimension(:)::x,y,b,c
end type qspline

contains

function qspline_alloc(x,y) result(s)
real::x(:),y(:)
real,allocatable::p(:),h(:)
type(qspline)::s
integer::n,i
n=size(x)
allocate(s%x(n), s%y(n), s%b(n-1), s%c(n-1))
s%x=x
s%y=y
allocate(p(n-1),h(n-1))
do i=1,n-1
	h(i)=x(i+1)-x(i)
	p(i)=(y(i+1)-y(i))/h(i)
end do
s%c(1)=0
do i=1,n-2
	s%c(i+1)=(p(i+1)-p(i)-s%c(i)*h(i))/h(i+1)
end do
s%c(n-1)=s%c(n-1)/2
do i=n-2,1,-1
	s%c(i)=(p(i+1)-p(i)-s%c(i+1)*h(i+1))/h(i)
end do
do i=1,n-1
	s%b(i)=p(i)-s%c(i)*h(i)
end do
deallocate(p)
deallocate(h)
end function qspline_alloc

function qspline_eval(s,z) result(f)
type(qspline)::s
real::z,f,h
integer::i,j,m
i=1
j=size(s%x)
do while (j-i>1)
	m=(i+j)/2
	if (z>=s%x(m)) then
		i=m
	else
		j=m
	end if
end do
h=z-s%x(i)
f=s%y(i)+s%b(i)*h+s%c(i)*h*h
end function qspline_eval
end module quadspline