!****************************************************
!*      Program to demonstrate least squares        *
!*         polynomial fitting subroutine            *
!* ------------------------------------------------ *
!* Reference: BASIC Scientific Subroutines, Vol. II *
!* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].*
!*                                                  *
!*                F90 version by J-P Moreau, Paris  *
!*                        (www.jpmoreau.fr)         *
!* ------------------------------------------------ *
!* SAMPLE RUN:                                      *
!*                                                  *
!* LEAST SQUARES POLYNOMIAL FITTING                 *
!*                                                  *
!* What is the order of the fit      : 4            *
!* What is the error reduction factor: 0            *
!* How many data pooints are there   : 11           *
!*                                                  *
!* Input the data points as prompted:               *
!*                                                  *
!*  1   X, Y = 1,1                                  *
!*  2   X, Y = 2,2                                  *
!*  3   X, Y = 3,3                                  *
!*  4   X, Y = 4,4                                  *
!*  5   X, Y = 5,5                                  *
!*  6   X, Y = 6,6                                  *
!*  7   X, Y = 7,7                                  *
!*  8   X, Y = 8,8                                  *
!*  9   X, Y = 9,9                                  *
!*  10   X, Y = 0,0                                 *
!*  11   X, Y = 1,1                                 *
!*                                                  *
!* Coefficients are:                                *
!*  0   -0.000000                                   *
!*  1   1.000000                                    *
!*  2   -0.000000                                   *
!*  3   0.000000                                    *
!*  4   -0.000000                                   *
!*                                                  *
!* Standard deviation = 0.000000                    *
!*                                                  *
!****************************************************
PROGRAM DEMO_LS_POLY

parameter(SIZE=25)

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

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

  print *,' '
  print *,'LEAST SQUARES POLYNOMIAL FITTING'
  
  write(*,"(' What is the order of the fit      : ')",advance='no')
  read *, m
  write(*,"(' What is the error reduction factor: ')",advance='no')
  read *, e1
  write(*,"(' How many data pooints are there   : ')",advance='no')
  read *, n
  
  print *,' '
  print *,'Input the data points as prompted:'
  
  do i=1, n
    write(*,50,advance='no') i
    read *, x(i), y(i)
  end do
  print *,' '

  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
  
  write(*,70)  dd
  stop 

50 format('  ',i2,'   X  Y = ')
60 format('  ',i2,'   ',f9.6)
70 format(' Standard deviation: ',f9.6//) 

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 e:    *
!*  1. if e = 0, the fit is to order m,                          *
!*  2. if e > 0, the order of fit increases towards m, but will  *
!*     stop if the relative standard deviation does not decrease *
!*     by more than e between successive fits.                   *
!* The order of the fit then obtained is l.                      *
!*****************************************************************
Subroutine LS_POLY(m,e1,n,l,x,y,c,dd)
  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

! End of file lsqrpoly.f90