Professor Diomar Cesar Lobao

Universidade Federal Fluminense-Volta Redonda, RJ, Brasil

Diomar Cesar


Dept. Ciências Exatas - Exact Science Dept.

Search

lsq.f

c<html>
c<head><title>lsq.f</title></head>
c<body>
c<pre>
      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<a name=1><font color=FF0000>
            rsid=sqrt(rsid/float(ndata-1))
c</font>
            write(6,2001) rsid
 2001       format(' Fit to 2nd order polynomial has a mean error of',
     $             1p,e12.5)
      return
      end
c</pre>
c</body>
c</html>
Skip to content