program main
	use matrix
	use leastsqr
	implicit none
	real, external :: func, func2
	integer, parameter :: kmax = 2
	real, dimension(0:10) :: x, y, sigma
	real, dimension(kmax) :: sol
	real, dimension(kmax,kmax) :: cov
	real, dimension(3) :: sol2
	real, dimension(3,3) :: cov2
	integer :: i, N
	real :: xmax, xmin, dx, xi
	
	! Make data:
	open(51, file='points.dat', status='replace', action='write')
	do i = 0,10
		! Calculate the datapoints:
		x(i) = i+0.5
		y(i) = i+cos(2.0*i)
		sigma(i) = 1+sin(2.0*i)**2
		
		! Write the points to a file:
		write(51, '(3f10.5)') x(i), y(i), sigma(i)
	enddo
	close(51)
	
	! Parameters used for saving the data for the graphical plots:
	N = 100 ! The number of points in fit-plot.
	xmin = minval(x); xmax = maxval(x); dx = (xmax-xmin)/N;
	
	print *, "*********************************"
	print *, "Linear fit:"
	
	! Do the fit:
	call leastsqrfit(x, y, sigma, kmax, func, sol, cov)
	
	! Write the solution and error-matrix:
	call print_matrix('solution =', sol)
	call print_matrix('cov =', cov)

	! Write the solution to file for plot:
	open(52, file='fit1.dat', status='replace', action='write')
	do i = 1,N
		xi = dx*(i-1)+xmin
		write(52, '(2f10.5)') xi, sol(2)*xi+sol(1)
	enddo
	close(52)
	
	
	print *, "*********************************"
	print *, "Fit with parabola:"
	
	! Do fit with parabola-function instead:
	call leastsqrfit(x, y, sigma, 3, func2, sol2, cov2)
	
	! Print solution and error matrix:
	call print_matrix('solution =', sol2)
	call print_matrix('cov =', cov2)
	
	! Write the solution to file for plot:
	open(52, file='fit2.dat', status='replace', action='write')
	do i = 1,N
		xi = dx*(i-1)+xmin
		write(52, '(2f10.5)') xi, sol2(2)*xi+sol2(1)+sol2(3)*xi**2
	enddo
	close(52)
end program main

! Function to make linear fit:
real function func(k, x)
	implicit none
	integer, intent(in) :: k
	real, intent(in) :: x
	
	if (k == 1) then
		func = 1.0
	elseif (k == 2) then
		func = x
	endif
end function func

! Function to make parabola-fit:
real function func2(k, x)
	implicit none
	integer, intent(in) :: k
	real, intent(in) :: x
	
	if (k == 1) then
		func2 = 1.0
	elseif (k == 2) then
		func2 = x
	elseif (k == 3) then
		func2 = x**2
	endif
end function func2