program main
	use interp
	implicit none
	integer, parameter :: N = 14, Ni = 1000, near = 5
	real, dimension(N) :: x, y
	real :: xmin, xmax, dx, z
	integer :: i
	
	! Data to interpolate:
	!x = (/ -5., -4., -3., -2., -1., 1., 2., 3., 4., 5. /)
	!y = (/ 0.05, 0.1, 0.15, 0.3, 0.8, 0.8, 0.3, 0.15, 0.1, 0.05 /)
	x = (/ -5., -4., -3., -2., -1., 1., 2., 3., 4., 5., 6., 7., 8., 9. /)
	y = (/ 0.05, 0.1, 0.15, 0.3, 0.8, 0.8, 0.3, 0.15, 0.1, 0.05, 0.4, 0.5, 0.8, 1.0 /)
	
	! Write points to file:
	open(51, file='points.dat', status='replace', action='write')
	do i = 1,N
		write(51, '(2f10.5)'), x(i), y(i)
	enddo
	close(51)
	
	! Variables used for interpolation:
	xmin = minval(x)+0.01
	xmax = maxval(x)-0.01
	dx = (xmax-xmin)/(Ni-1)
	
	! Do polonomial interpolation on dataset and write the results to file:
	print *, "Doing polonomial interpolation with the ", near
	print *, 'nearest points...'
	open(52, file='pointerp.dat', status='replace', action='write')
	do i = 1,Ni
		z = xmin+(i-1)*dx
		write(52, '(2f10.5)'), z, pointerp(x, y, z, near)
	enddo
	close(52)
	print *, "Done. Results can be fount in pointerp.dat."
	
	! Dp quadratic-spline interpolation on dataset and write results to file:
	print *, "Doing quadratic-spline interpolation..."
	open(52, file='qspline.dat', status='replace', action='write')
	do i = 1,Ni
		z = xmin+(i-1)*dx
		write(52, '(2f10.5)'), z, qspline(x, y, z)
	enddo
	close(52)
	print *, "Done. Results can be found in qspline.dat."
	
	print *, ""
	print *, "As it can be seen from the plot, the quadratic spline"
	print *, "will in this case begin to oscillate around the points."
	print *, "In this case the polonomial interpolation gives the "
	print *, "better result."
end program main