quadfit.f

      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
            rsid=sqrt(rsid/float(ndata-1))
            write(6,2001) rsid
 2001       format(' Fit to 2nd order polynomial has a mean error of',
     $             1p,e12.5)
      return
      end