C..............
C
      program extrap_power
C
C   Many algorithms have errors that go as (dx)^k,
C   where k is some known integer
C 
C  This program does a linear least squares fit, 
C  and extrapolates for dx -> 0
C
C  CWJ @ SDSU Sept 20101
C
C  MAIN PROGRAM
C  subroutines called:
C       read_from_keyboard(maxdim,npts,k,array_dx,array_ans)
C       extrapolate_pwrk(maxdim,npts,k,array_dx,array_ans,anslim)
C

      implicit none
      integer maxdim        ! maximum dimension of arrays
      parameter (maxdim = 100)
      real array_dx(maxdim)
      real array_ans(maxdim)
      real anslim
      integer npts  ! # of points
      integer k      ! assume error goes like (dx)^k

      call read_from_keyboard(maxdim,npts,k,array_dx,array_ans)

      call extrapolate_pwrk(maxdim,npts,k,array_dx,array_ans,anslim)

      print*,' Final answer is...',anslim,' ! '
      end
C----------------------------------------
C
C reads in data from keyboard
C
C  INPUT:
C      maxdim = max dim of arrays; set in main program
C  OUTPUT
C      ndim   = # of points
C      k      = assumed power of error
C      array_dx(i), i = 1,ndim   = array of dx
C      array_ans(i),i = 1,ndim   = array of ans(dx)
C  

      subroutine read_from_keyboard(maxdim,ndim,k,array_dx,array_ans)

      implicit none
      integer maxdim   ! = max dimension of arrays
      integer ndim     ! # of points actually used
      integer k        ! assume error goes like (dx)^k
      real array_dx(maxdim)    ! array of dxs
      real array_ans(maxdim)   ! array of ans(dx)

      integer i   ! dummy integer for looping

!............... ERROR TRAP...........
      if(maxdim < 1)then
           print*,' problem with max dim ',maxdim
           stop
      endif
!...............END ERROR TRAP.............
    
      print*,' Welcome to extrapolating with (dx)^k ! '
      print*,' Enter your choice of integer k '
      read*,k
      print*,' Enter  # of data points '
      read*,ndim

      if(ndim > maxdim .or. ndim < 1)then   ! ERROR TRAP
          print*,' problem ',ndim
          stop
      endif

      do i = 1,ndim
         print*,' Enter dx, ans '
         read*,array_dx(i), array_ans(i)
      end do  ! loop over i

      return
      end   ! subroutine read_from_keyboard

C----------------------------------------------------------
C
C  computes extrapolated answer
C
C  INPUT:
C      maxdim = max dim of arrays; set in main program
C      ndim   = # of points
C      k      = assumed power of error
C      array_dx(i), i = 1,ndim   = array of dx
C      array_ans(i),i = 1,ndim   = array of ans(dx)
C  OUTPUT
C      anslim  = limit answer
C

      subroutine extrapolate_pwrk(maxdim,ndim,k,arr_dx,arr_ans,anslim)

      implicit none
      integer maxdim   ! = max dimension of arrays
      integer ndim     ! # of points actually used
      integer k        ! assume error goes like (dx)^k
      real arr_dx(maxdim)    ! array of dxs
      real arr_ans(maxdim)   ! array of ans(dx)
      real anslim
      
      real sum0, sum1, sum2, sum3
      integer i   ! dummy integer for looping

      sum0 = 0.0
      sum1 = 0.0
      sum2 = 0.0
      sum3 = 0.0

      do i = 1,ndim
        sum0 = sum0+ arr_ans(i)
        sum1 = sum1+ arr_ans(i)*arr_dx(i)**k
        sum2 = sum2+ arr_dx(i)**k
        sum3 = sum3+ arr_dx(i)**(2*k)
      end do
      sum0 = sum0/float(ndim)
      sum1 = sum1/float(ndim)
      sum2 = sum2/float(ndim)
      sum3 = sum3/float(ndim)

      anslim = (sum3*sum0-sum1*sum2)/(sum3-sum2*sum2)

      return
      end