!Determination of spring constant - Linear Least-Squares Fit of (T**2, m).
!Because the plot uses T**2 instead of T, errors are recalculated for each A(i,j) 
!and b(i).

program SPRING
implicit none
                                                                                
     integer, parameter :: dbl=SELECTED_REAL_KIND(14), m=2, n=3
     real(KIND=dbl) :: A(n,m), b(n), c(m), Q(n,m), sx=0.01, sx_new(n), & 
                       & sy=0.0001, T(n), mass(n), K, pi=3.141592654
     real(KIND=dbl), dimension(m,m) :: R, R_INV, R_INV_T, E, ID
     integer :: i, j

T(1)=0.78      ; T(2)=0.80      ; T(3)=1.00
mass(1)=0.0698 ; mass(2)=0.0743 ; mass(3)=0.1194

do i = 1, n
   sx_new(i)=2*T(i)*sx
   A(i, 1) = 1/SQRT(sx_new(i)**2 + sy**2)
   do j = 2, m
   A(i,j) = T(i)**2/SQRT(sx_new(i)**2 + sy**2)
   ID(m,m)=1
   enddo
enddo
                                                                                
do i = 1, n
   b(i) = mass(i)/SQRT(sx_new(i)**2 + sy**2)
enddo

call MATOUT("A is:", A, n, m)
call MATOUT("b is:", b, n, 1)
                                                                                
CALL QR_GS(A, n, m, Q, R)
CALL QR_BACK(Q, R, n, m, b, c)

print *, "The linear fit is:"
print *, "m =", c(2), "*   T**2  ", c(1)

!Calculation of the covariance (error) matrix, E = R_INV*(R_INV)_T :
!R is already upper triangular => Q=identity matrix, R=R
CALL QR_INVERSE(ID, R, m, R_INV)
R_INV_T = TRANSPOSE(R_INV)
E = MATMUL(R_INV, R_INV_T)

call MATOUT("The covariance matrix is:", E, m, m)

K = c(2)*(2*pi)**2
print *, "The spring constant:"
print *, "k =", K, "with error:", SQRT(E(2,2))

end program SPRING