!-------------------------------------------------------------------------------
! Copyright (c) 2021, Whitman T. Dailey
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
! 1. Redistributions of source code must retain the above copyright notice, this
!    list of conditions and the following disclaimer.
! 2. Redistributions in binary form must reproduce the above copyright notice,
!    this list of conditions and the following disclaimer in the documentation
!    and/or other materials provided with the distribution.
! 3. Neither the name of the copyright holder nor the names of its
!    contributors may be used to endorse or promote products derived from
!    this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!-------------------------------------------------------------------------------
Module Quadratures

    Use Kinds, Only: dp
    Implicit None
    Private
    Public :: Romberg_Quad
    Public :: Composite_Quad
    Public :: Composite_Trapezoid
    Public :: Composite_Simpson
    Public :: Adaptive_Quad
    Public :: Adaptive_Trapezoid
    Public :: Adaptive_Simpson
    Public :: MidPoint_Quad
    Public :: Trapezoid_Quad
    Public :: Simpson_Quad
    Public :: GaussLegendre5
    Public :: GaussLegendre16
    Public :: GaussLegendre96
    Public :: GaussLegendreN
    Public :: Progressive_GaussLegendre

    Real(dp), Parameter :: one_third = 1._dp / 3._dp
    Real(dp), Parameter :: one_sixth = 1._dp / 6._dp
    Integer, Parameter :: bisect_limit = 53

    Interface Romberg_Quad
        Module Procedure Romb_Quad
        Module Procedure Romb_Quad_ranges
    End Interface

Contains

Function Romb_Quad_ranges(f,ab,aTol,rTol,p) Result(q)
    !Integrates f(x) on the intervals between entries in ab(:) by extrapolation on successive composite trapezoid
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: ab(:)  !list of limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerances for convergence
    Real(dp), Intent(Out), Optional :: p  !estimated precision achived in the final estimate
    Real(dp) :: h  !the total width of the integration window from ab(1) to ab(n)
    Real(dp) :: hi  !the width of the i-th integration window
    Integer :: n  !length of the list of integration limits, resulting in n-1 integrations to be performed
    Integer :: i  !loop counter
    Real(dp) :: p_i  !estimated precision of the integral from the i-th interval

    q = 0._dp
    n = Size(ab)
    h = ab(n) - ab(1)
    Do i = 1,n-1
        hi = ab(i+1) - ab(i)
        If (Present(p)) Then
            q = q + Romb_Quad(f,ab(i),ab(i+1),(hi/h)*aTol,rTol,p=p_i)
            p = Min(p,p_i)
        Else
            q = q + Romb_Quad(f,ab(i),ab(i+1),(hi/h)*aTol,rTol)
        End If
    End Do
End Function Romb_Quad_ranges

Function Romb_Quad(f,a,b,aTol,rTol,allow_GL_fallback,p,n_ord,n_ext) Result(q)
    !Integrates f(x) on (a,b) by extrapolation on successive composite trapezoid
    Use Kinds, Only: dp
    Use Utilities, Only: Converged
    Use Utilities, Only: Prec
    Implicit None
    Real(dp) :: q    !the result of the integration
    Interface
        Function f(x)    !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerances for convergence
    Logical, Intent(In), Optional :: allow_GL_fallback  !flag signaling whether fallback to progressive G-L should be used in the 
                                                        !event of a failure to converge during extrapolation stages
    Real(dp), Intent(Out), Optional :: p  !estimated precision achived in the final estimate
    Integer, Intent(Out), Optional :: n_ord  !total number of function evaluations (ordinates)
    Integer, Intent(Out), Optional :: n_ext  !total number of extrapolation stages (rows in the table)
    Integer, Parameter :: Tmax = 21  !maximum number of extrapolations in the table
    Real(dp) :: T(0:Tmax)  !Extrapolation table previous row
    Real(dp) :: Tk0,Tk  !Extrapolation table current row values
    Integer :: i,j,k  !counters: i for table row, j for quadrature ordinates, k for table column
    Integer :: n      !number of intervals
    Real(dp) :: h0,h  !spacing between quadrature ordinates
    Real(dp) :: fk    !multiplier for extrapolation steps
    Real(dp) :: s     !sum of function values at quadrature ordinates
    Real(dp) :: Prec0 !precision of the previous extrapolation (used to monitor convergence)
    Integer :: ifin   !final number of extrapolations performed
    Logical :: GL_fallback
    !    If the preprocessor directive ROMB_TABLES is set at compile-time, 
    ! then the tables of computed values are dumped to file in the working 
    ! directory as the routine runs.
#   if ROMB_TABLES
        Integer :: unit
#   endif

    If (a .EQ. b) Then !zero-width interval for integration
        q = 0._dp
        If (Present(p)) p = Prec(0._dp,q)
        If (Present(n_ord)) n_ord = 0
        If (Present(n_ext)) n_ext = 0
        Return
    End If
    !Initial trapezoid estimate
    n = 1
    s = 0.5_dp * (f(a) + f(b))
    h0 = b - a
    T(0) = h0 * s
    Prec0 = -1._dp
#   if ROMB_TABLES
        Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
        Do k = 0,Tmax
            Write(unit,'(I24)',ADVANCE='NO') k
        End Do
        Write(unit,'(ES24.15)') T(0)
        Close(unit)
#   endif
    Do i = 1,Tmax !up to Tmax rows in the table
        !Trapezoid estimate for the 0-th column of the i-th row of table
        n = n * 2
        h = h0 / Real(n,dp)
        Do j = 1,n-1,2  !Odd values of j are NEW points at which to evaluate f
            s = s + f(a + Real(j,dp)*h)
        End Do
        Tk0 = h * s
        !Fill i-th row, columns k = 1:i, with extrapolated estimates
        fk = 1._dp
        Do k = 1,i  !up to i columns this row
            fk = fk * 4._dp
            Tk = (fk * Tk0 - T(k-1)) / (fk - 1._dp)
            If (k .LT. i) Then
                T(k-1) = Tk0  !store Tk0 for next i
                Tk0 = Tk  !store Tk for next k
            End If !otherwise, skip storage steps if working final column
        End Do
        !Check for convergence
        If (Converged(T(i-1),Tk,rTol=Real(2**k-1)*rTol,aTol=aTol)) Then
            q = Tk
#           if ROMB_TABLES
                T(i-1) = Tk0
                T(i) = Tk
                Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
                Do k = 0,i
                    Write(unit,'(ES24.15)',ADVANCE='NO') T(k)
                End Do
                Write(unit,'(A)') '*'
                Write(unit,*)
                Close(unit)
#           endif
            If (Present(p)) p = Prec(T(i-1),Tk)
            If (Present(n_ord)) n_ord = n
            If (Present(n_ext)) n_ext = i
            Return  !Normal exit
        Else  !Check for failures and prep for the next time though the loop
            !check for exit conditions other than convergence (failures)
            If (i .EQ. Tmax) Then  !extrapolation limit has been reached without convergence
                ifin = Tmax
                Exit
            Else If (i .GT. Tmax/2) Then
                If (Prec(T(i-1),Tk) .LT. Prec0) Then  !precision was LOST on this extrapolation, convergence will not occur
                    ifin = i
                    Exit
                End If
            End If
            !store Tk0 and Tk for next i, and update precision monitor
            Prec0 = Prec(T(i-1),Tk)
            T(i-1) = Tk0
            T(i) = Tk
#           if ROMB_TABLES
                Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
                Do k = 0,i
                    Write(unit,'(ES24.15)',ADVANCE='NO') T(k)
                End Do
                Write(unit,*)
                Close(unit)
#           endif
        End If
    End Do
    !If we get this far, we did not converge
    If (Present(allow_GL_fallback)) Then
        GL_fallback = allow_GL_fallback
    Else
        GL_fallback = .TRUE.
    End If
    If (GL_fallback) Then
        If (Present(p)) Then
            q = Progressive_GaussLegendre(f,a,b,aTol=aTol,rTol=rTol,n_start=5,n_stride=3,n_done=n,p=p)
        Else
            q = Progressive_GaussLegendre(f,a,b,aTol=aTol,rTol=rTol,n_start=5,n_stride=3,n_done=n)
        End If
        If (Present(n_ord)) n_ord = n
        If (Present(n_ext)) n_ext = -1
    Else
        Write(*,*)
        Write(*,'(A,I0,A)')          'WARNING:  Quadratures: Romberg_Quad:  Failed to converge in ',ifin,' extrapolations.'
        If (ifin .LT. Tmax) Then
            Write(*,'(A10,A)')              '','Extrapolation was terminated for loss of precision.'
        End If
        Write(*,'(A10,A,ES23.15E2,A,F0.5,A)') '','Final estimated value: ',Tk,' (~',Prec(T(ifin-1),Tk),' good digits)'
        Write(*,'(A10,A,2(ES10.3E2,A))')      '','Final Extrapolation Error: ',Abs(Tk-T(ifin-1)),   ' (abs), ', & 
                                       &                                       Abs(Tk-T(ifin-1))/Tk,' (rel)'
        Write(*,'(A10,A,2(ES10.3E2,A))')      '','Convergence Criteria:      ',atol,                ' (abs), ', & 
                                       &                                       rtol,                ' (rel)'
        q = Tk
#       if ROMB_TABLES
            T(i-1) = Tk0
            T(i) = Tk
            Open(NEWUNIT=unit,FILE='Romb_Tables.tst',ACTION='WRITE',STATUS='UNKNOWN',POSITION='APPEND')
            Do k = 0,i
                Write(unit,'(ES24.15)',ADVANCE='NO') T(k)
            End Do
            Write(unit,'(A)') '*'
            Write(unit,*)
            Close(unit)
#       endif
        If (Present(p)) p = Prec(T(ifin-1),Tk)
        If (Present(n_ord)) n_ord = n
        If (Present(n_ext)) n_ext = ifin
    End If
    ! ERROR STOP
    !Currently, no STOP is implemented for a failure to converge.  The result of the quadrature is set to the final value from the 
    !Romberg table, and additional required values are set, and execution returns to the calling routine after the above warnings
    !are sent to the screen.
End Function Romb_Quad

Function Richardson_Extrapolation(fh,f2h,p) Result(r)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: r
    Real(dp), Intent(In) :: fh,f2h
    Integer, Intent(In) :: p

    r = ( fh*Real(2**p,dp) - f2h ) / Real(2**p-1,dp)
End Function Richardson_Extrapolation

Function Adaptive_Quad(f,a0,b0,Quad,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by adaptively applying composite quadrature
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
        Function Quad(f,a,b,p)  !the quadrature rule to use
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: Quad  !the result of the quadrature
            Interface
                Function f(x)    !the function to be integrated
                    Use Kinds,Only: dp
                    Implicit None
                    Real(dp) :: f
                    Real(dp), Intent(In) :: x
                End Function f
            End Interface
            Real(dp),Intent(In) :: a,b  !limits of integration
            Integer, Intent(Out), Optional :: p
        End Function Quad
    End Interface
    Real(dp), Intent(In) :: a0,b0  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: n,i
    Real(dp), Allocatable :: s(:,:)
    integer, Allocatable :: sb(:)
    Real(dp) :: a,b,ab
    Real(dp) :: q0,q1,q2
    Real(dp) :: h
    Real(dp) :: rTol_per_h,aTol_per_h
    Real(dp) :: r,r0,dr
    Integer :: p
    Integer :: bisects

    n = 16
    i = 0
    Allocate(s(1:3,1:n))
    s = 0._dp
    Allocate(sb(1:n))
    sb = 0
    q0 = Quad(f,a0,b0,p=p)
    Call Put_s(a0,b0,q0,0)
    rTol_per_h = rTol * Real(2**p-1,dp)
    aTol_per_h = aTol / (b0 - a0)
    q = 0._dp
    Do
        If (Get_s(a,b,q0,bisects)) Then  !got an interval to do
            ab = 0.5_dp * (a + b)
            q1 = Quad(f,a,ab)
            q2 = Quad(f,ab,b)
            If (p .GT. 1) Then
                r = Richardson_Extrapolation(q1+q2,q0,p)
                r0 = q1 + q2
            Else
                r = q1 + q2
                r0 = q0
            End If
            h = b - a
            dr = Abs( r - r0 )
            If (      dr .LT. h * rTol_per_h    & !relative convergence
               & .OR. dr .LT. h * aTol_per_h    & !absolute convergence
               & .OR. bisects .GT. bisect_limit & !bisection limit exceeded
               &      ) Then
                q = q + r  !add to result of the quadrature
            Else  !add the subintervals to the list to be refined again
                Call Put_s(a,ab,q1,bisects+1)
                Call Put_s(ab,b,q2,bisects+1)
            End If
        Else !no intervals left to quadrature
            Exit
        End If
    End Do
Contains
    Subroutine Put_s(c1,c2,c3,bs)
        Use Kinds, Only: dp
        Implicit None
        Real(dp), Intent(In) :: c1,c2,c3
        Integer, Intent(In) :: bs
        Real(dp), Allocatable :: swaps(:,:)
        Integer, Allocatable :: swapsb(:)

        i = i + 1
        If (i .GT. n) Then !need more space
            Allocate(swaps(1:3,1:n))
            swaps = s
            Deallocate(s)
            Allocate(s(1:3,1:2*n))
            s = 0._dp
            s(:,1:n) = swaps
            Deallocate(swaps)
            Allocate(swapsb(1:n))
            swapsb = sb
            Deallocate(sb)
            Allocate(sb(1:2*n))
            sb = 0
            sb(1:n) = swapsb
            Deallocate(swapsb)
            n = n * 2
        End If
        s(:,i) = (/ c1 , c2 , c3 /)
        sb(i) = bs
    End Subroutine Put_s
    Function Get_s(c1,c2,c3,bs) Result(bingo)
        Use Kinds, Only: dp
        Implicit None
        Logical :: bingo
        Real(dp), Intent(Out) :: c1,c2,c3
        Integer, Intent(Out) :: bs

        If (i .LT. 1) Then !there was no s to get
            bingo = .FALSE.
        Else
            bingo = .TRUE.
            c1 = s(1,i)
            c2 = s(2,i)
            c3 = s(3,i)
            bs = sb(i)
            i = i - 1
        End If
    End Function Get_s
End Function Adaptive_Quad

Function Adaptive_Trapezoid(f,a0,b0,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by adaptively applying composite trapezoid
    Use Kinds,Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a0,b0  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: n,i
    Real(dp), Allocatable :: s(:,:)
    Integer, Allocatable :: sb(:)
    Real(dp) :: a,b,ab
    Real(dp) :: fa,fb,fab
    Real(dp) :: q0,q1,q2,r,dr
    Real(dp) :: h
    Real(dp) :: rTol_per_h,aTol_per_h
    Integer :: bisects

    n = 16
    i = 0
    Allocate(s(1:5,1:n))
    s = 0._dp
    Allocate(sb(1:n))
    sb = 0
    i = 0
    !Initial Trapezoid estimate
    fa = f(a0)
    fb = f(b0)
    q0 = (b0 - a0) * (fa + fb) * 0.5_dp
    Call Put_s(a0,b0,fa,fb,q0,0)
    rTol_per_h = rTol * 3._dp
    aTol_per_h = aTol / (b0 - a0)
    q = 0._dp
    Do
        If (Get_s(a,b,fa,fb,q0,bisects)) Then  !got an interval to do
            ab = 0.5_dp * (a + b)
            fab = f(ab)
            h = b - a
            q1 = 0.25_dp * h * (fa + fab)
            q2 = 0.25_dp * h * (fab + fb)
            r = Richardson_Extrapolation(q1+q2,q0,2)
            dr = Abs( q1 + q2 - r )
            If (      dr .LT. h * rTol_per_h    & !relative convergence
               & .OR. dr .LT. h * aTol_per_h    & !absolute convergence
               & .OR. bisects .GT. bisect_limit & !bisection limit exceeded
               &      ) Then
                q = q + r
            Else  !add the subintervals to the list to be quadratured
                Call Put_s(a,ab,fa,fab,q1,bisects+1)
                Call Put_s(ab,b,fab,fb,q2,bisects+1)
            End If
        Else !no intervals left to quadrature
            Exit
        End If
    End Do
Contains
      Subroutine Put_s(c1,c2,c3,c4,c5,bs)
        Use Kinds, Only: dp
        Implicit None
        Real(dp), Intent(In) :: c1,c2,c3,c4,c5
        Integer, Intent(In) :: bs
        Real(dp), Allocatable :: swaps(:,:)
        Integer, Allocatable :: swapsb(:)

        i = i + 1
        If (i .GT. n) Then !need more space
            Allocate(swaps(1:5,1:n))
            swaps = s
            Deallocate(s)
            Allocate(s(1:5,1:2*n))
            s = 0._dp
            s(:,1:n) = swaps
            Deallocate(swaps)
            Allocate(swapsb(1:n))
            swapsb = sb
            Deallocate(sb)
            Allocate(sb(1:2*n))
            sb = 0
            sb(1:n) = swapsb
            Deallocate(swapsb)
            n = n * 2
        End If
        s(:,i) = (/ c1 , c2 , c3 , c4, c5 /)
        sb(i) = bs
    End Subroutine Put_s
    Function Get_s(c1,c2,c3,c4,c5,bs) Result(bingo)
        Use Kinds, Only: dp
        Implicit None
        Logical :: bingo
        Real(dp), Intent(Out) :: c1,c2,c3,c4,c5
        Integer, Intent(Out) :: bs

        If (i .LT. 1) Then !there was no s to get
            bingo = .FALSE.
        Else
            bingo = .TRUE.
            c1 = s(1,i)
            c2 = s(2,i)
            c3 = s(3,i)
            c4 = s(4,i)
            c5 = s(5,i)
            bs = sb(i)
            i = i - 1
        End If
    End Function Get_s
End Function Adaptive_Trapezoid

Function Adaptive_Simpson(f,a0,b0,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by adaptively applying composite Simpson
    Use Kinds,Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds,Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a0,b0  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: n,i
    Real(dp), Allocatable :: s(:,:)
    Integer, Allocatable :: sb(:)
    Real(dp) :: x1,x2,x3,x4,x5
    Real(dp) :: y1,y2,y3,y4,y5
    Real(dp) :: q0,q1,q2,r,dr
    Real(dp) :: h
    Real(dp) :: rTol_per_h,aTol_per_h
    Integer :: bisects

    n = 16
    i = 0
    Allocate(s(1:7,1:n))
    s = 0._dp
    Allocate(sb(1:n))
    sb = 0
    i = 0
    !Initial Simpson estimate
    x3 = 0.5_dp * (a0 + b0)
    y1 = f(a0)
    y3 = f(x3)
    y5 = f(b0)
    q0 = (b0 - a0) * (y1 + 4._dp*y3 + y5) * one_sixth
    Call Put_s(a0,x3,b0,y1,y3,y5,q0,0)
    rTol_per_h = rTol * 15._dp
    aTol_per_h = aTol / (b0 - a0)
    q = 0._dp
    Do
        If (Get_s(x1,x3,x5,y1,y3,y5,q0,bisects)) Then  !got an interval to do
            x2 = 0.5_dp * (x1 + x3)
            x4 = 0.5_dp * (x3 + x5)
            y2 = f(x2)
            y4 = f(x4)
            q1 = (x3-x1) * (y1 + 4._dp*y2 + y3) * one_sixth
            q2 = (x5-x3) * (y3 + 4._dp*y4 + y5) * one_sixth
            r = Richardson_Extrapolation(q1+q2,q0,4)
            h = x5 - x1
            dr = Abs( q1 + q2 - r )
            If (      dr .LT. h * rTol_per_h    & !relative convergence
               & .OR. dr .LT. h * aTol_per_h    & !absolute convergence
               & .OR. bisects .GT. bisect_limit & !bisection limit exceeded
               &      ) Then
                q = q + r
            Else  !add the subintervals to the list to be quadratured
                Call Put_s(x1,x2,x3,y1,y2,y3,q1,bisects+1)
                Call Put_s(x3,x4,x5,y3,y4,y5,q2,bisects+1)
            End If
        Else !no intervals left to quadrature
            Exit
        End If
    End Do
Contains
      Subroutine Put_s(c1,c2,c3,c4,c5,c6,c7,bs)
        Use Kinds, Only: dp
        Implicit None
        Real(dp), Intent(In) :: c1,c2,c3,c4,c5,c6,c7
        Integer, Intent(In) :: bs
        Real(dp), Allocatable :: swaps(:,:)
        Integer, Allocatable :: swapsb(:)

        i = i + 1
        If (i .GT. n) Then !need more space
            Allocate(swaps(1:7,1:n))
            swaps = s
            Deallocate(s)
            Allocate(s(1:7,1:2*n))
            s = 0._dp
            s(:,1:n) = swaps
            Deallocate(swaps)
            Allocate(swapsb(1:n))
            swapsb = sb
            Deallocate(sb)
            Allocate(sb(1:2*n))
            sb = 0
            sb(1:n) = swapsb
            Deallocate(swapsb)
            n = n * 2
        End If
        s(:,i) = (/ c1 , c2 , c3 , c4, c5 , c6 , c7 /)
        sb(i) = bs
    End Subroutine Put_s
    Function Get_s(c1,c2,c3,c4,c5,c6,c7,bs) Result(bingo)
        Use Kinds, Only: dp
        Implicit None
        Logical :: bingo
        Real(dp), Intent(Out) :: c1,c2,c3,c4,c5,c6,c7
        Integer, Intent(Out) :: bs

        If (i .LT. 1) Then !there was no s to get
            bingo = .FALSE.
        Else
            bingo = .TRUE.
            c1 = s(1,i)
            c2 = s(2,i)
            c3 = s(3,i)
            c4 = s(4,i)
            c5 = s(5,i)
            c6 = s(6,i)
            c7 = s(7,i)
            bs = sb(i)
            i = i - 1
        End If
    End Function Get_s
End Function Adaptive_Simpson

Function Composite_Quad(f,a,b,Quad,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by composite quadrature
    Use Kinds, Only: dp
    Use Utilities, Only: Converged
    Implicit None
    Real(dp) :: q    !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
        Function Quad(f,a,b,p)  !the quadrature rule to use
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: Quad  !the result of the quadrature
            Interface
                Function f(x)    !the function to be integrated
                    Use Kinds ,Only: dp
                    Implicit None
                    Real(dp) :: f
                    Real(dp), Intent(In) :: x
                End Function f
            End Interface
            Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
            Integer, Intent(Out), Optional :: p
        End Function Quad
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: i,j,n
    Real(dp) :: h,h0
    Real(dp) :: q0,r,dr
    Real(dp) :: eps
    Integer :: p

    n = 1
    q0 = Quad(f,a,b,p=p)
    h0 = b - a
    eps = Real(2**p-1,dp) * h0 * rTol
    Do i = 1,bisect_limit+1  !max intervals is bounded by splitting interval bisect_limit+1 times
        q = 0._dp
        n = n * 2
        h = h0 / Real(n,dp)
        Do j = 1,n
            q = q + Quad(f , a+h*Real((j-1),dp) , a+h*Real(j,dp))
        End Do
        If (p .GT. 1) Then !apply a stage of Richardson extrapolation for a low-cost refinement of the estimate
            r = Richardson_Extrapolation(q,q0,p)
            dr = Abs( r - q )
            If (      dr .LT. eps  & !relative convergence
               & .OR. dr .LT. aTol & !absolute convergence
               &      ) Then
                q = r
                Return
            End If
        Else !Richardson extrapolation does not provide a benefit for methods with p<2 (i.e. Midpoint, etc.)
            If (Converged(q0,q,rTol=rTol,aTol=aTol)) Return
        End If
        q0 = q
    End Do
    !if we get this far, we failed to converge
    Write(*,'(A,I0,A)') 'ERROR:  Quadratures: Composite_Quad:  Failed to converge on ',n,' intervals.'
    ERROR STOP
End Function Composite_Quad

Function MidPoint_Quad(f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the quadrature
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p

    q = (b - a) * f( 0.5_dp * (a + b) )
    If (Present(p)) p = 2
End Function MidPoint_Quad

Function Trapezoid_Quad(f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the quadrature
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p

    q = 0.5_dp * (b - a) * ( f(a) + f(b) )
    If (Present(p)) p = 2
End Function Trapezoid_Quad

Function Simpson_Quad(f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp):: q  !the result of the quadrature
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p

    q = one_sixth * (b - a) * ( f(a) + 4._dp * f( 0.5_dp * (a + b) ) + f(b) )
    If (Present(p)) p = 4
End Function Simpson_Quad

Function Composite_Trapezoid(f,a,b,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by composite trapezoid
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: i,j,n
    Real(dp) :: h,h0
    Real(dp) :: q0,s
    Real(dp) :: r,dr
    Real(dp) :: eps

    n = 1
    s = 0.5_dp * ( f(a) + f(b) )
    h0 = b - a
    eps = 3._dp * h0 * rTol
    q0 = h0 * s
    Do i = 1,bisect_limit+1  !max intervals is bounded by splitting interval bisect_limit+1 times
        n = n * 2
        h = h0 / Real(n,dp)
        Do j = 1,n-1,2  !only odd values of j, these are the NEW points at which to evaluate f
            s = s + f( a + h*Real(j,dp) )
        End Do
        q = h * s
        r = Richardson_Extrapolation(q,q0,2)
        dr = Abs( r - q )
        If (      dr .LT. eps  & !relative convergence
           & .OR. dr .LT. aTol & !absolute convergence
           &      ) Then
             q = r
            Return
        End If
        q0 = q
    End Do
    !if we get this far, we failed to converge
    Write(*,'(A,I0,A)') 'ERROR:  Quadratures: Progressive_Trapezoid:  Failed to converge on ',n,' intervals.'
    ERROR STOP
End Function Composite_Trapezoid

Function Composite_Simpson(f,a,b,aTol,rTol) Result(q)
    !Integrates f(x) on (a,b) by composite Simpson
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !the result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer :: i,j,n
    Real(dp) :: h,h0
    Real(dp) :: q0,s1,s2,s3
    Real(dp) :: r,dr,eps

    n = 2
    s1 = f(a) + f(b)
    s2 = 0._dp
    s3 = f( 0.5_dp * (a + b) )
    h0 = b - a
    eps = 15._dp * h0 * rTol
    q0 = one_sixth * h0 * (s1 + 4._dp*s3)
    Do i = 1,bisect_limit+1  !max intervals is bounded by splitting interval bisect_limit+1 times
        n = n * 2
        h = h0 / Real(n,dp)
        s2 = s2 + s3
        s3 = 0._dp
        Do j = 1,n-1,2  !only odd values of j, these are the NEW points at which to evaluate f
            s3 = s3 + f( a + h*Real(j,dp) )
        End Do
        q = one_third * h * (s1 + 2._dp*s2 + 4._dp*s3)
        r = Richardson_Extrapolation(q,q0,4)
        dr = Abs( r - q )
        If (      dr .LT. eps  & !relative convergence
           & .OR. dr .LT. aTol & !absolute convergence
           &      ) Then
             q = r
            Return
        End If
        q0 = q
    End Do
    !if we get this far, we failed to converge
    Write(*,'(A)') 'ERROR:  Quadratures: Progressive_Simpson:  Failed to converge on ',n,' intervals.'
    ERROR STOP
End Function Composite_Simpson

Function GaussLegendre(f,a,b,n,wi,xi) Result(q)
    Use Kinds, Only: dp
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(In) :: n
    Real(dp), Intent(In) :: wi(1:n),xi(1:n)
    Real(dp) :: c1,c2  !changes interval from (a,b) to (-1,1)
    Integer :: i
    Real(dp) :: fi(1:n)

    c1 = 0.5_dp * (b - a)
    c2 = 0.5_dp * (b + a)
    Do i = 1,n
        fi(i) = f(c1 * xi(i) + c2)
    End Do
    q = c1 * Dot_Product(wi,fi)
End Function GaussLegendre

Function GaussLegendre5(f,a,b) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a, Only: n => n005
    Use GaussLegendre_w_a, Only: wi => w005
    Use GaussLegendre_w_a, Only: xi => a005
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule

    q = GaussLegendre(f,a,b,n,wi,xi)
End Function GaussLegendre5

Function GaussLegendre16(f,a,b) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a, Only: n => n016
    Use GaussLegendre_w_a, Only: wi => w016
    Use GaussLegendre_w_a, Only: xi => a016
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule

    q = GaussLegendre(f,a,b,n,wi,xi)
End Function GaussLegendre16

Function GaussLegendre96(f,a,b) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a, Only: n => n096
    Use GaussLegendre_w_a, Only: wi => w096
    Use GaussLegendre_w_a, Only: xi => a096
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule

    q = GaussLegendre(f,a,b,n,wi,xi)
End Function GaussLegendre96

Function GaussLegendreN(n,f,a,b,p) Result(q)
    Use Kinds, Only: dp
    Use GaussLegendre_w_a
    Implicit None
    Real(dp) :: q  !result of the quadrature
    Integer, Intent(In) :: n  !number of quadrature points
    Interface
        Function f(x)  !the function to integrate
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function f
    End Interface
    Real(dp), Intent(In) :: a,b  !interval on which to apply the rule
    Integer, Intent(Out), Optional :: p
    Real(dp), Allocatable :: wi(:),ai(:)  !weights and abscissa

    Allocate(wi(1:n))
    Allocate(ai(1:n))
    Select Case (N)
        Case (1)
            wi = w001
            ai = a001
        Case (2)
            wi = w002
            ai = a002
        Case (3)
            wi = w003
            ai = a003
        Case (4)
            wi = w004
            ai = a004
        Case (5)
            wi = w005
            ai = a005
        Case (6)
            wi = w006
            ai = a006
        Case (7)
            wi = w007
            ai = a007
        Case (8)
            wi = w008
            ai = a008
        Case (9)
            wi = w009
            ai = a009
        Case (10)
            wi = w010
            ai = a010
        Case (11)
            wi = w011
            ai = a011
        Case (12)
            wi = w012
            ai = a012
        Case (13)
            wi = w013
            ai = a013
        Case (14)
            wi = w014
            ai = a014
        Case (15)
            wi = w015
            ai = a015
        Case (16)
            wi = w016
            ai = a016
        Case (17)
            wi = w017
            ai = a017
        Case (18)
            wi = w018
            ai = a018
        Case (19)
            wi = w019
            ai = a019
        Case (20)
            wi = w020
            ai = a020
        Case (21)
            wi = w021
            ai = a021
        Case (22)
            wi = w022
            ai = a022
        Case (23)
            wi = w023
            ai = a023
        Case (24)
            wi = w024
            ai = a024
        Case (25)
            wi = w025
            ai = a025
        Case (26)
            wi = w026
            ai = a026
        Case (27)
            wi = w027
            ai = a027
        Case (28)
            wi = w028
            ai = a028
        Case (29)
            wi = w029
            ai = a029
        Case (30)
            wi = w030
            ai = a030
        Case (31)
            wi = w031
            ai = a031
        Case (32)
            wi = w032
            ai = a032
        Case (33)
            wi = w033
            ai = a033
        Case (34)
            wi = w034
            ai = a034
        Case (35)
            wi = w035
            ai = a035
        Case (36)
            wi = w036
            ai = a036
        Case (37)
            wi = w037
            ai = a037
        Case (38)
            wi = w038
            ai = a038
        Case (39)
            wi = w039
            ai = a039
        Case (40)
            wi = w040
            ai = a040
        Case (41)
            wi = w041
            ai = a041
        Case (42)
            wi = w042
            ai = a042
        Case (43)
            wi = w043
            ai = a043
        Case (44)
            wi = w044
            ai = a044
        Case (45)
            wi = w045
            ai = a045
        Case (46)
            wi = w046
            ai = a046
        Case (47)
            wi = w047
            ai = a047
        Case (48)
            wi = w048
            ai = a048
        Case (49)
            wi = w049
            ai = a049
        Case (50)
            wi = w050
            ai = a050
        Case (51)
            wi = w051
            ai = a051
        Case (52)
            wi = w052
            ai = a052
        Case (53)
            wi = w053
            ai = a053
        Case (54)
            wi = w054
            ai = a054
        Case (55)
            wi = w055
            ai = a055
        Case (56)
            wi = w056
            ai = a056
        Case (57)
            wi = w057
            ai = a057
        Case (58)
            wi = w058
            ai = a058
        Case (59)
            wi = w059
            ai = a059
        Case (60)
            wi = w060
            ai = a060
        Case (61)
            wi = w061
            ai = a061
        Case (62)
            wi = w062
            ai = a062
        Case (63)
            wi = w063
            ai = a063
        Case (64)
            wi = w064
            ai = a064
        Case (65)
            wi = w065
            ai = a065
        Case (66)
            wi = w066
            ai = a066
        Case (67)
            wi = w067
            ai = a067
        Case (68)
            wi = w068
            ai = a068
        Case (69)
            wi = w069
            ai = a069
        Case (70)
            wi = w070
            ai = a070
        Case (71)
            wi = w071
            ai = a071
        Case (72)
            wi = w072
            ai = a072
        Case (73)
            wi = w073
            ai = a073
        Case (74)
            wi = w074
            ai = a074
        Case (75)
            wi = w075
            ai = a075
        Case (76)
            wi = w076
            ai = a076
        Case (77)
            wi = w077
            ai = a077
        Case (78)
            wi = w078
            ai = a078
        Case (79)
            wi = w079
            ai = a079
        Case (80)
            wi = w080
            ai = a080
        Case (81)
            wi = w081
            ai = a081
        Case (82)
            wi = w082
            ai = a082
        Case (83)
            wi = w083
            ai = a083
        Case (84)
            wi = w084
            ai = a084
        Case (85)
            wi = w085
            ai = a085
        Case (86)
            wi = w086
            ai = a086
        Case (87)
            wi = w087
            ai = a087
        Case (88)
            wi = w088
            ai = a088
        Case (89)
            wi = w089
            ai = a089
        Case (90)
            wi = w090
            ai = a090
        Case (91)
            wi = w091
            ai = a091
        Case (92)
            wi = w092
            ai = a092
        Case (93)
            wi = w093
            ai = a093
        Case (94)
            wi = w094
            ai = a094
        Case (95)
            wi = w095
            ai = a095
        Case (96)
            wi = w096
            ai = a096
        Case (97)
            wi = w097
            ai = a097
        Case (98)
            wi = w098
            ai = a098
        Case (99)
            wi = w099
            ai = a099
        Case (100)
            wi = w100
            ai = a100
        Case (101)
            wi = w101
            ai = a101
        Case (102)
            wi = w102
            ai = a102
        Case (103)
            wi = w103
            ai = a103
        Case (104)
            wi = w104
            ai = a104
        Case (105)
            wi = w105
            ai = a105
        Case (106)
            wi = w106
            ai = a106
        Case (107)
            wi = w107
            ai = a107
        Case (108)
            wi = w108
            ai = a108
        Case (109)
            wi = w109
            ai = a109
        Case (110)
            wi = w110
            ai = a110
        Case (111)
            wi = w111
            ai = a111
        Case (112)
            wi = w112
            ai = a112
        Case (113)
            wi = w113
            ai = a113
        Case (114)
            wi = w114
            ai = a114
        Case (115)
            wi = w115
            ai = a115
        Case (116)
            wi = w116
            ai = a116
        Case (117)
            wi = w117
            ai = a117
        Case (118)
            wi = w118
            ai = a118
        Case (119)
            wi = w119
            ai = a119
        Case (120)
            wi = w120
            ai = a120
        Case (121)
            wi = w121
            ai = a121
        Case (122)
            wi = w122
            ai = a122
        Case (123)
            wi = w123
            ai = a123
        Case (124)
            wi = w124
            ai = a124
        Case (125)
            wi = w125
            ai = a125
        Case (126)
            wi = w126
            ai = a126
        Case (127)
            wi = w127
            ai = a127
        Case (128)
            wi = w128
            ai = a128
        Case (129)
            wi = w129
            ai = a129
        Case (130)
            wi = w130
            ai = a130
        Case (131)
            wi = w131
            ai = a131
        Case (132)
            wi = w132
            ai = a132
        Case (133)
            wi = w133
            ai = a133
        Case (134)
            wi = w134
            ai = a134
        Case (135)
            wi = w135
            ai = a135
        Case (136)
            wi = w136
            ai = a136
        Case (137)
            wi = w137
            ai = a137
        Case (138)
            wi = w138
            ai = a138
        Case (139)
            wi = w139
            ai = a139
        Case (140)
            wi = w140
            ai = a140
        Case (141)
            wi = w141
            ai = a141
        Case (142)
            wi = w142
            ai = a142
        Case (143)
            wi = w143
            ai = a143
        Case (144)
            wi = w144
            ai = a144
        Case (145)
            wi = w145
            ai = a145
        Case (146)
            wi = w146
            ai = a146
        Case (147)
            wi = w147
            ai = a147
        Case (148)
            wi = w148
            ai = a148
        Case (149)
            wi = w149
            ai = a149
        Case (150)
            wi = w150
            ai = a150
        Case Default
            Write(*,'(A,I0)') 'ERROR:  Quadratures: GaussLegendreN:  Unsupported N, n = ',N
            ERROR STOP
    End Select
    q = GaussLegendre(f,a,b,n,wi,ai)
    If (Present(p)) p = 2*n - 1
End Function GaussLegendreN

Function Progressive_GaussLegendre(f,a,b,aTol,rTol,n_start,n_stride,n_done,p) Result(q)
    !Integrates f(x) on (a,b) by successive quadrature via Gauss-Legendre n-point rules
    Use Kinds, Only: dp
    Use Utilities, Only: Converged
    Use Utilities, Only: Prec
    Implicit None
    Real(dp) :: q  !result of the integration
    Interface
        Function f(x)  !the function to be integrated
            Use Kinds, Only: dp
            Implicit None
            Real(dp) :: f
            Real(dp), Intent(In) :: x
        End Function
    End Interface
    Real(dp), Intent(In) :: a,b  !limits of integration
    Real(dp), Intent(In) :: aTol,rTol  !absolute and relative tolerance to be used as criteria for convergence
    Integer, Intent(In), Optional :: n_start,n_stride
    Integer, Intent(Out), Optional :: n_done
    Real(dp), Intent(Out), Optional :: p
    Real(dp) :: q0
    Integer :: n,dn
    Integer :: i

    If (Present(n_start)) Then
        n = n_start
    Else
        n = 1
    End If
    If (Present(n_stride)) Then
        dn = n_stride
    Else
        dn = 1
    End If
    q0 = GaussLegendreN(n,f,a,b)
    Do i = n+dn,150,dn !Gauss-Legendre weights and abscissa are available up to 150 points
        q = GaussLegendreN(i,f,a,b)
        If (Converged(q,q0,rTol=rTol,aTol=aTol)) Then
            If (Present(n_done)) n_done = i
            If (Present(p)) Then
                If (i .LT. 150) Then !progress one more stage of GL points to estimate precision
                    q0 = q
                    q = GaussLegendreN(i+1,f,a,b)
                    p = Prec(q0,q)
                Else
                    p = Prec(q0,q)
                End If
            End If
            RETURN !NORMAL EXIT
        End If
        q0 = q
    End Do
    !If we get this far, we failed to converge on 150-point Gauss
    q0 = GaussLegendreN(149,f,a,b) !regress one stage of GL points to estimate precision
    Write(*,*)
    Write(*,'(A,I0,A)') 'WARNING:  Quadratures: Progressive_GaussLegendre:  Failed to converge with ',150,' Gauss-points.'
    Write(*,'(A10,A,ES23.15E2,A,F0.5,A)') '','Final estimated value: ',q,' (~',Prec(q,q0),' good digits)'
    Write(*,'(A10,A,2(ES10.3E2,A))')      '','Final estimated error:     ',Abs(q-q0),  ' (abs), ', & 
                                   &                                       Abs(q-q0)/q,' (rel)'
    Write(*,'(A10,A,2(ES10.3E2,A))')      '','Convergence Criteria:      ',atol,       ' (abs), ', & 
                                   &                                       rtol,       ' (rel)'
    If (Present(n_done)) n_done = -1
    If (Present(p)) p = Prec(q0,q)
    ! ERROR STOP
    !Currently, no STOP is implemented for a failure to converge.  The result of the quadrature is set to the final value from the 
    !quadrature, and additional required values are set, and execution returns to the calling routine after the above warnings
    !are sent to the screen.
End Function Progressive_GaussLegendre

End Module Quadratures
