subroutine qspline(x,y,n,a,b)
*
* quadratic spline y = y_i + a_i (x-x_i) + b_i (x-x_i)}^2
*
	dimension x(n),y(n),a(n-1),b(n-1)

	a(1)=(y(2)-y(1))/(x(2)-x(1))
	do i=1,n-1
		dx=x(i+1)-x(i)
		b(i)=( y(i+1)-y(i)-a(i)*dx )/dx**2
		a(i+1)=a(i)+2*b(i)*dx
	end do		
	pn=( a(n-1)+2*b(n-1)*(x(n)-x(n-1)) )
	a(n)=( pn+(y(n)-y(n-1))/(x(n)-x(n-1)) )/2
	do i=n-1,1,-1
		dx=x(i+1)-x(i)
		b(i)=(a(i+1)*dx-y(i+1)+y(i))/dx**2
		a(i)=a(i+1)-2*b(i)*dx
	end do
	return
	end