c clsq.f c c
```      module arrays
real, allocatable :: t(:), z(:)
integer np
c
c   t   -  array containing times (seconds) at which positions are measured
c   z   -  measured distance of the falling object from the nominal point of
c          release (fit may let you estimate any offset in the actual release
c          point
c   np  -  number of data points
c
end module
c
program leastsq
use arrays
implicit none
c
c   Program designed to perform a least squares fit of a quadratic equation
c   to some data.  In this specific example the data is the result of
c   measuring the location of a falling object at various times.  The result
c   of interest is gravitational acceleration obtained from the coefficient
c   off t**2.
c
c       John Mahaffy   2/6/95
c
real coef(0:2),g
c
c   Variables
c
c   t   -  array containing times (seconds) at which positions are measured
c   z   -  measured distance of the falling object from the nominal point of
c          release (fit may let you estimate any offset in the actual release
c          point
c   np  -  number of data points
c   coef - Coefficients in the second order polynomial approximating the data
c          z = coef(0) + coef(1)*t + coef(2)*t**2
c   g   -  approximation to the gravitational acceleration deduced from data
c
c
call getdata
call quadfit(t,z,np,coef(0))
g=2*coef(2)
write(6,2000) g
2000 format (' Predicted value of g = ', f6.3,' m/s**2')
write(6,2001) coef(0),coef(1)
2001 format(' Predicted initial offset = ',f7.4,' m'/
#         ' Predicted initial velocity = ',f7.4,' m/s')
stop
end
subroutine getdata
use arrays
implicit none
c
c   Get Experimental Data
c
c   Input Arguments: NONE
c   Output Arguments:
c     t     -    time (sec) of measurement
c     s     -    distance of fall measured (meters)
c     np    -    number of data points
c
integer i, iend
open(11,file='fall.data')
c
c    Count Data Pairs in fall.data
c
np=-1
do while (iend.eq.0)
read(11,*,iostat=iend)
np = np + 1
enddo
if (iend.gt.0.or. np.lt.1) then
print *, 'Empty or Bad File, check fall.data'
stop
endif
c
c    Rewind, Allocate Arrays, and Read data
c
rewind (11)
print *, np, ' data points read.'
allocate (t(1:np),z(1:np))
do i = 1,np
read (11,*) t(i),z(i)
enddo
c
end
c
subroutine quadfit (xd,yd,ndata,c)
implicit none
c
c       Take corresponding data points from the arrays xd and yd and fit
c       them with the following equation:
c
c       y = c(1) + c(2) * x + c(3) * x**2
c
real xd(ndata),yd(ndata),aa(3,3),c(3),fs,suma,sumb,sqrt,rsid
integer ipvt(3),ir,is,ij,ndata, info
c
c    Input Arguments:
c      xd   -   x values for data points
c      yd   -   y valuess for data point pairs
c      ndata -  number of data points
c
c    Output Arguements:
c      c    -   array containing the three coefficients in the 2nd order
c               polynomial that provide the "best" fit to the data from a
c               least squares method.  Note that it also temporarily holds
c               the values of the right hand sides of each Least Squares
c               equation.
c    Other Key Variables:
c      aa   -   matrix containing coefficients of the system of equations
c               generated by the least squares method
c      suma -   a variable that tallies the sum that is needed to generate
c               each element of aa.
c      sumb -   a variable that tallies the sum that is needed for the right
c               hand side of each Least Squares equation.
c      ir   -   an index that is used to keep track of the equation number
c               within the system of equations .
c      is   -   an index used to track the coefficient number within a given
c               Least Squares equation.
c
c
c    DO loop 55 generates terms for each of the 3 Least Squares Equations
c
do 55 ir=1,3
c
c       DO loop 45 generates the right hand side for a given equation
c
sumb=0.
do 45 ij=1,ndata
45        sumb=sumb+yd(ij)*xd(ij)**(ir-1)
c(ir)=sumb
c
c       DO loop 50 generates the coefficients a given equation
c
do 50 is=1,3
suma=0.
do 52 ij=1,ndata
52           suma=suma+xd(ij)**(is-1)*xd(ij)**(ir-1)
50        aa(ir,is)=suma
55     continue
c
c    Solve the Least Squares Equations
c
call sgefa(aa,3,3,ipvt,info)
call sgesl(aa,3,3,ipvt,c ,0)
c
c    Calculate a measure of mean error between the data and the curve
c
rsid=0.
do 65 ir=1,ndata
fs=0.
do 60 is=1,3
60             fs=fs+c(is)*xd(ir)**(is-1)
65          rsid=rsid+(yd(ir)-fs)**2
c
rsid=sqrt(rsid/float(ndata-1))
c
write(6,2001) rsid
2001       format(' Fit to 2nd order polynomial has a mean error of',
\$             1p,e12.5)
return
end
c```
c c