!****************************************************
!*      Program to demonstrate least squares        *
!*         polynomial fitting subroutine            *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!*   by F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981.  *
!*                                                  *
!*               F90 version by J-P Moreau, Paris   *
!*                                                  *
!*               Modified by Pat Fitzpatrick,       *
!*               Mississippi State University       *
!* ------------------------------------------------ *
!*                                                  *
!****************************************************
      PROGRAM DEMO_LS_POLY

      integer SIZE

      parameter(SIZE=25)

      integer i,l,m,n
      real*8  dd,e1,vv

      real*8 c(0:SIZE),x(SIZE),y(SIZE) 

      real*8 fit,variance,mean,explained_var,total_var

      x(1)=1.0
      y(1)=.018

      x(2)=1.5
      y(2)=.05

      x(3)=2.0
      y(3)=.09

      x(4)=2.5
      y(4)=.13

      x(5)=3.0
      y(5)=.145

      x(6)=4.0
      y(6)=.175

      x(7)=5.0
      y(7)=.2

      x(8)=6.0
      y(8)=.225

      x(9)=8.0
      y(9)=.27

      x(10)=10.0
      y(10)=.3

      x(11)=15.0
      y(11)=.35

      x(12)=20.0
      y(12)=.4

      x(13)=30.0
      y(13)=.5

      x(14)=40.0
      y(14)=.57

c      data x/1.0,1.5,2.0,2.5,3.0,4.0,5.0,6.0,8.0,10.0,&
c     &       15.0,20.0,30.0,40.0/
c      data y/.018,.05,.09,.13,.145,.175,.2,.225,.27,.3,&
c     &       .35,.4,.5,.57/

      explained_var=0
      total_var=0
      mean=0
      m=5
      e1=0
      n=14


      do i=1,n
         mean=mean+y(i)
      enddo 
      mean=mean/n
  
      call LS_POLY(m,e1,n,l,x,y,c,dd)

      print *,'Coefficients are:'
      print *,' '
      do i=0, l
         write(*,60)  i, c(i) 
      end do
60    format('  ',i2,'   ',f9.6)

      print*,' '

      do i=1,n
         fit=0
         do k=o,l
            fit=fit+c(k)*x(i)**k
         enddo
         print*,'Polynomial fit= ',fit,' y=',y(i)
         explained_var=explained_var+(fit-mean)**2
         total_var=total_var+(y(i)-mean)**2
      enddo

      variance=100.*explained_var/total_var

      print*,'The variance explained is ',variance,' percent'
      print*,'The standard deviation is ',dd

      stop 

      end

!*****************************************************************
!*         LEAST SQUARES POLYNOMIAL FITTING PROCEDURE            *
!* ------------------------------------------------------------- *
!* This program least squares fits a polynomial to input data.   *
!* forsythe orthogonal polynomials are used in the fitting.      *
!* The number of data points is n.                               *
!* The data is input to the subroutine in x[i], y[i] pairs.      *
!* The coefficients are returned in c[i],                        *
!* the smoothed data is returned in v[i],                        *
!* the order of the fit is specified by m.                       *
!* The standard deviation of the fit is returned in d.           *
!* There are two options available by use of the parameter e1:    *
!*  1. if e1 = 0, the fit is to order m,                          *
!*  2. if e1 > 0, the order of fit increases towards m, but will  *
!*     stop if the relative standard deviation does not decrease *
!*     by more than e1 between successive fits.                   *
!* The order of the fit then obtained is l.                      *
!*****************************************************************
      Subroutine LS_POLY(m,e1,n,l,x,y,c,dd)
      integer SIZE
      parameter(SIZE=25)
!Labels: 10,15,20,30,50
      real*8 x(SIZE),y(SIZE),v(SIZE),a(SIZE),b(SIZE)
      real*8 c(0:SIZE),d(SIZE),c2(SIZE),e(SIZE),f(SIZE)
      integer i,l,l2,m,n,n1
      real*8 a1,a2,b1,b2,c1,dd,d1,e1,f1,f2,v1,v2,w
      n1 = m + 1; l=0
      v1 = 1.d7
! Initialize the arrays
      do i=1, n1
         a(i) = 0.d0; b(i) = 0.d0; f(i) = 0.d0
      end do
      do i=1, n
         v(i) = 0.d0; d(i) = 0.d0
      end do
      d1 = dsqrt(dfloat(n)); w = d1;
      do i=1, n
         e(i) = 1.d0 / w
      end do
      f1 = d1; a1 = 0.d0
      do i=1, n
        a1 = a1 + x(i) * e(i) * e(i)
      end do
      c1 = 0.d0
      do i=1, n
         c1 = c1 + y(i) * e(i)
      end do
      b(1) = 1.d0 / f1; f(1) = b(1) * c1
      do i=1, n
         v(i) = v(i) + e(i) * c1
      end do
      m = 1
! Save latest results
10    do i=1, l
         c2(i) = c(i)
      end do
      l2 = l; v2 = v1; f2 = f1; a2 = a1; f1 = 0.d0
      do i=1, n
         b1 = e(i)
         e(i) = (x(i) - a2) * e(i) - f2 * d(i)
         d(i) = b1
         f1 = f1 + e(i) * e(i)
      end do
      f1 = dsqrt(f1)
      do i=1, n
         e(i) = e(i) / f1
      end do
      a1 = 0.d0
      do i=1, n  
         a1 = a1 + x(i) * e(i) * e(i)
      end do
      c1 = 0.d0
      do i=1, n  
         c1 = c1 + e(i) * y(i)
      end do	 
      m = m + 1; i = 0
15    l = m - i; b2 = b(l); d1 = 0.d0
      if (l > 1)  d1 = b(l - 1)
      d1 = d1 - a2 * b(l) - f2 * a(l)
      b(l) = d1 / f1; a(l) = b2; i = i + 1
      if (i.ne.m) goto 15
      do i=1, n
         v(i) = v(i) + e(i) * c1 
      end do
      do i=1, n1
         f(i) = f(i) + b(i) * c1
         c(i) = f(i)
      end do
      vv = 0.d0
      do i=1, n
         vv = vv + (v(i) - y(i)) * (v(i) - y(i))
      end do 
!Note the division is by the number of degrees of freedom
      vv = dsqrt(vv / dfloat(n - l - 1)); l = m
      if (e1.eq.0.d0) goto 20
!Test for minimal improvement
      if (dabs(v1 - vv) / vv < e1) goto 50
!if error is larger, quit
      if (e1 * vv > e1 * v1) goto 50
      v1 = vv
20    if (m.eq.n1) goto 30
      goto 10
!Shift the c[i] down, so c(0) is the constant term
30    do i=1, l  
         c(i - 1) = c(i)
      end do
      c(l) = 0.d0
! l is the order of the polynomial fitted
      l = l - 1; dd = vv
      return
! Aborted sequence, recover last values
50    l = l2; vv = v2
      do i=1, l  
         c(i) = c2(i)
      end do
      goto 30
      end