def qwe(rtol, atol, maxint, inp, intervals, lambd=None, off=None,
        factAng=None):
    r"""Quadrature-With-Extrapolation.

    This is the kernel of the QWE method, used for the Hankel (``hqwe``) and
    the Fourier (``fqwe``) Transforms. See ``hqwe`` for an extensive
    description.

    This function is based on ``qwe.m`` from the source code distributed with
    [Key12]_.

    """
    def getweights(i, inpint):
        r"""Return weights for this interval."""
        return (np.atleast_2d(inpint)[:,  i+1] - np.atleast_2d(inpint)[:, i])/2

    # 1. Calculate the first interval for all offsets
    if hasattr(inp, '__call__'):  # Hankel and not spline
        EM0 = inp(0, lambd, off, factAng)
    else:                         # Fourier or Hankel with spline
        EM0 = inp[:, 0]
    EM0 *= getweights(0, intervals)

    # 2. Pre-allocate arrays and initialize
    EM = np.zeros(EM0.size, dtype=EM0.dtype)                # EM array
    om = np.ones(EM0.size, dtype=bool)                      # Convergence array
    S = np.zeros((EM0.size, maxint), dtype=EM0.dtype)  # Working arr. 4 recurs.
    relErr = np.zeros((EM0.size, maxint))                   # Relative error
    extrap = np.zeros((EM0.size, maxint), dtype=EM0.dtype)  # extrap. result
    kcount = 1  # Initialize kernel count (only important for Hankel)

    # 3. The extrapolation transformation loop
    for i in range(1, maxint):
        # 3.a Calculate the field for this interval
        if hasattr(inp, '__call__'):  # Hankel and not spline
            EMi = inp(i, lambd[om, :], off[om], factAng[om])
            kcount += 1  # Update count
        else:                         # Fourier or Hankel with spline
            EMi = inp[om, i]
        EMi *= getweights(i, intervals[om, :])

        # 3.b Compute Shanks transformation
        # Using the epsilon algorithm: structured after [Weni89]_, p26.
        S[:, i][om] = S[:, i-1][om] + EMi  # working array for transformation

        # Recursive loop
        aux2 = np.zeros(om.sum(), dtype=EM0.dtype)
        for k in range(i, 0, -1):
            aux1, aux2 = aux2, S[om, k-1]
            ddff = S[om, k] - aux2
            S[om, k-1] = np.where(np.abs(ddff) < np.finfo(np.double).tiny,
                                  np.finfo(np.double).max, aux1 + 1/ddff)

        # The extrapolated result plus the first interval term
        extrap[om, i-1] = S[om, np.mod(i, 2)] + EM0[om]

        # 3.c Analyze for convergence
        if i > 1:
            # Calculate relative and absolute error
            rErr = (extrap[om, i-1] - extrap[om, i-2])/extrap[om, i-1]
            relErr[om, i-1] = np.abs(rErr)
            abserr = atol/np.abs(extrap[om, i-1])

            # Update booleans
            om[om] *= relErr[om, i-1] >= rtol + abserr

            # Store in EM
            EM[om] = extrap[om, i-1]

        if (~om).all():
            break

    # 4. Cleaning up

    # Warning if maxint is potentially too small
    conv = i+1 != maxint

    # Catch the ones that did not converge
    EM[om] = extrap[om, i-1]

    # Set np.finfo(np.double).max to 0
    EM.real[EM.real == np.finfo(np.double).max] = 0

    return EM, kcount, conv
def qwe ( rtol , atol , maxint , inp , intervals , lambd = None , off = None , factAng = None ) : def getweights ( i , inpint ) : return ( np . atleast_2d ( inpint ) [ : , i + 1 ] - np . atleast_2d ( inpint ) [ : , i ] ) / 2 if hasattr ( inp , '__call__' ) : EM0 = inp ( 0 , lambd , off , factAng ) else : EM0 = inp [ : , 0 ] EM0 *= getweights ( 0 , intervals ) EM = np . zeros ( EM0 . size , dtype = EM0 . dtype ) om = np . ones ( EM0 . size , dtype = bool ) S = np . zeros ( ( EM0 . size , maxint ) , dtype = EM0 . dtype ) relErr = np . zeros ( ( EM0 . size , maxint ) ) extrap = np . zeros ( ( EM0 . size , maxint ) , dtype = EM0 . dtype ) kcount = 1 for i in range ( 1 , maxint ) : if hasattr ( inp , '__call__' ) : EMi = inp ( i , lambd [ om , : ] , off [ om ] , factAng [ om ] ) kcount += 1 else : EMi = inp [ om , i ] EMi *= getweights ( i , intervals [ om , : ] ) S [ : , i ] [ om ] = S [ : , i - 1 ] [ om ] + EMi aux2 = np . zeros ( om . sum ( ) , dtype = EM0 . dtype ) for k in range ( i , 0 , - 1 ) : aux1 , aux2 = aux2 , S [ om , k - 1 ] ddff = S [ om , k ] - aux2 S [ om , k - 1 ] = np . where ( np . abs ( ddff ) < np . finfo ( np . double ) . tiny , np . finfo ( np . double ) . max , aux1 + 1 / ddff ) extrap [ om , i - 1 ] = S [ om , np . mod ( i , 2 ) ] + EM0 [ om ] if i > 1 : rErr = ( extrap [ om , i - 1 ] - extrap [ om , i - 2 ] ) / extrap [ om , i - 1 ] relErr [ om , i - 1 ] = np . abs ( rErr ) abserr = atol / np . abs ( extrap [ om , i - 1 ] ) om [ om ] *= relErr [ om , i - 1 ] >= rtol + abserr EM [ om ] = extrap [ om , i - 1 ] if ( ~ om ) . all ( ) : break conv = i + 1 != maxint EM [ om ] = extrap [ om , i - 1 ] EM . real [ EM . real == np . finfo ( np . double ) . max ] = 0 return EM , kcount , conv
r"""Quadrature-With-Extrapolation.

    This is the kernel of the QWE method, used for the Hankel (``hqwe``) and
    the Fourier (``fqwe``) Transforms. See ``hqwe`` for an extensive
    description.

    This function is based on ``qwe.m`` from the source code distributed with
    [Key12]_.