! least square fit 

module realkinds
integer, parameter :: dpt = kind(1.0d0)
end module realkinds

program chpt5
use realkinds
implicit none

character(len=25) :: filename
integer, parameter :: lun = 18
integer :: ierror
integer :: i 
integer :: n = 0
real (kind=dpt), allocatable, dimension (:) :: x, y
real (kind=dpt) :: xum = 0
real (kind=dpt) :: xum2 = 0
real (kind=dpt) :: xyum = 0
real (kind=dpt) :: yum = 0
real (kind=dpt) :: yum2 = 0
real (kind=dpt) :: xavg, yavg, b, m , r, c

write(*,100) 
100 format (1x, 'if you enter the name of a compatable data file, this ' ,/, &
            1x, 'program will calculate a least-square fit of the given ' ,/, &
            1x, '(x, y) pairs')
read (*,120) filename
120 format (A)
write (*,*)

open (unit = lun, file = filename, status = 'old', action = 'read', iostat = ierror)

do while (ierror == 0)
  read(lun, *, iostat = ierror) c
  n = n + 1
end do

n = n - 1

allocate (x(n))
allocate (y(n))

rewind (lun)

do i = 1, n, 1
  read (lun, *, iostat = ierror) x(i), y(i)
end do

close (unit = lun)

if (ierror == 0) then
  do i = 1, n, 1
    xum = xum + x(i)
    yum = yum + y(i)
    xum2 = xum2 + x(i)**2_dpt
    yum2 = yum2 + y(i)**2_dpt
    xyum = xyum + x(i)*y(i)
  end do
 
 xavg = xum / (1.0_dpt*n)
 yavg = yum / (1.0_dpt*n)
 m = (xyum - xum*yavg)/(xum2 - xum*xavg)
 b = (yavg - m*xavg)
 r = ((n*xyum - xum*yum) / (sqrt((n*xum2 - xum**2_dpt)*(n*yum2 - yum**2_dpt))))

 write(*,*) "equation of lin for least...f it"
 write(*,110) 'y=', m, 'x+', b
 110 format (A2, 1F3.1, A2, 1F3.1)

 write(*,*) 'the correlation coefficient is =', r
   sub1: if (abs(r) < 0.3) then
             write(*,*) 'warning the correlation coefficient for the fit&
                         of (x,y) points is very low and shows a week   &
                         linear relationship between x and y'
            end if sub1
else if (ierror /= 0) then
  write(*,*) 'file could not open properly'
end if 

deallocate (x,y)

stop
end program chpt5