
const modlobP = [
    -1.0,
    -0.81649658092772603273242802490196,
    -.44721359549995793928183473374626,
    0,
    0.44721359549995793928183473374626,
    0.81649658092772603273242802490196,
    1.0,
]

const modlobWei = [
    0.052380952380952380952380952380952,
    0.29387755102040816326530612244898,
    0.42517006802721088435374149659864,
    0.45714285714285714285714285714286,
]

const modlobW = [
    0.052380952380952380952380952380952,
    0.29387755102040816326530612244898,
    0.42517006802721088435374149659864,
    0.45714285714285714285714285714286,
    0.42517006802721088435374149659864,
    0.29387755102040816326530612244898,
    0.052380952380952380952380952380952,
]

const modlobNulw = [
    0.1148787002745345 -0.2954023721345174 0.4102810724090519 -0.4595148010981382
    0.2168775557284140 -0.4553480127399748 0.3463949695733805 0.0
    0.2957100706945793 -0.4080798975585194 -0.1182840282778318 0.4613077102835437
    0.3560582128767142 -0.1825453130808602 -0.4628897488019905 0.0
]
modlobNw = zeros(Float64, (4, 7))

for j = 1:4
    for k = 1:4
        modlobNw[j, k] = modlobNulw[j, k]
    end
    for k = 1:3
        modlobNw[j, 8-k] = modlobNulw[j, k]
    end
end
for k = 5:7
    modlobNw[2, k] = -modlobNw[2, k]
end
for k = 5:7
    modlobNw[4, k] = -modlobNw[4, k]
end

#Espelid 'modlob'
function modlob(f, a::T, b::T; tol = sqrt(eps(T))) where {T}
    #  %MODLOB Numerically evaluates an integral using adaptively
    #  % a Lobatto 7 point rule developed by Gander and Gautschi.
    #  %
    #  % Q=MODLOB('F',A,B) approximates the integral of F(X) from
    # % A to B to a relative precision TOL= 10^-3. 'F' is a
    #  % string containing the name of the function. The
    # % function F must return a vector of output values if
    #  % given a vector of input values.
    # % Q=MODLOB('F',A,B,TOL) integrates to a relative
    #  % error of TOL.
    #  %
    #  % Q=MODLOB('F',A,B,TOL,TRACE) displays the left
    #  % end point of the current interval, the interval
    #  % length, and the partial integral.
    #  %
    # % Q=MODLOB('F',A,B,TOL,TRACE,P1,P2,...) allows
    #  % coefficients P1, ... to be passed directly to the
    #  % function F: G=F(X,P1,P2,...). To use default values
    #  % for TOL or TRACE, one may pass the empty matrix ([]).
    #  %
    #  % Terje O. Espelid 31/01/02
    #  % Reference: Gander and Gautschi: Bit 2000
    #  % Berntsen and Espelid, TOMS 1991
    #  % Espelid, in Espelid and Genz, 1992.
    #  % Espelid, Tech report May 2002.

    if tol <= 0
        tol = eps(T)
    end
    # % 7 points Lobatto rule in the interval [-1,1]; we give only 4 points due to symmetry.
    # % 4 Orthogonal Nullrules of degree 5,4,3 and 2 are given. The odd numbered are symmetric.
    # % All even numbered are antisymmetric. We given first the five elements.
    # % Then we use the symmetry
    # % Thus: p is the points, w is the weights and
    # % nw is the nullrules associated with the 7 points Lobatto rule.
    h = (b - a) / 2
    center = (b + a) / 2
    x = Vector{T}(undef, 7)
    y = Vector{T}(undef, 7)
    is = zero(T)
    isabs = zero(T)
    for j = 1:length(x)
        x[j] = center + h * modlobP[j]
        y[j] = f(x[j])
        is += h * modlobW[j] * y[j]
        isabs += h * modlobW[j] * abs(y[j])
    end
    nEval = 32 #not more than 32 levels
    fa = y[1]
    fb = y[7]
    c = x[4]
    fc = y[4]
    if isnan(fa)
        throw(DomainError(string("function is NaN at ", a)))
    elseif isnan(fb)
        throw(DomainError(string("function is NaN at ", b)))
    elseif isnan(fc)
        throw(DomainError(string("function is NaN at ", c)))
    end
    # % Apply the seven point rule to compute an estimate of the integral  to be used in the
    # % error estimator.
    # % Define the noise level
    noise = 50 * eps(T) * isabs
    # % Check that the estimate is not too small.
    absis = abs(is)
    if absis <= noise
        is = b - a
        absis = is
    end
    # % If tol is chosen too small, redefine tol too be on the noise level.
    tol = max(tol, noise / absis)
    return lobattor(f, a, c, fa, fc, is, tol, nEval) +
           lobattor(f, c, b, fc, fb, is, tol, nEval)
end

function lobattor(f, a::T, b::T, fa::T, fb::T, is::T, tol::T, nEval::Int)::T where {T}
    # %LOBATTOR Recursive function used by Lobatto.
    # %
    # % Q = LOBATTOR('F',A,B,FA,FB,IS,TOL) tries to
    # % approximate the integral of F(X) from A to B to
    # % an appropriate relative error. The argument 'F' is
    # % a string containing the name of f. The remaining
    # % arguments are generated by MODLOB or by recursion.
    # %
    #
    h = (b - a) / 2
    center = (b + a) / 2
    x = (
        center + h * modlobP[1],
        center + h * modlobP[2],
        center + h * modlobP[3],
        center + h * modlobP[4],
        center + h * modlobP[5],
        center + h * modlobP[6],
        center + h * modlobP[7],
    )
    z = (x[2], x[3], x[4], x[5], x[6])
    y = (fa, f(z[1]), f(z[2]), f(z[3]), f(z[4]), f(z[5]), fb)
    c = x[4]
    fc = y[4]
    # i = h* (w' * y)
    # iabs = h* (w' * abs.(y))
    # E = abs(h * (modelobNw' * y))
    # Emin = min(E[2:4])
    i = zero(T)
    iabs = zero(T)

    E = zeros(T, 4)
    for j = 1:length(y)
    	i += h * modlobW[j] * y[j]
    	iabs += h * modlobW[j] * abs(y[j])
    	for k = 1:length(E)
    		E[k] += h * modlobNw[k,j] * y[j]
    	end
    end
    @. E = abs(E)
    noise = 50 * iabs * eps(T)
    Emin = minimum(E[2:4])
    rmax = 2 * one(T)
    # Compute the local error estimate
    if Emin != 0
        # r=E(1:3)./E(2:4);
        # if (sum(isinf(r))>0)
        # rmax=2;
        # else
        # rmax=max(r);
        r1, r2, r3 = E[1] / E[2], E[2] / E[3], E[3] / E[4]
        if isinf(r1) || isinf(r2) || isinf(r3) || isnan(r1) || isnan(r2) || isnan(r3)
            #rmax = 2
        else
            rmax = max(r1, r2, r3)
        end
    end
    err = zero(T)
    C = 10
    if rmax > 1
        err = C * maximum(E)
    elseif one(T) / 2 < rmax
        err = C * rmax * E[2]
    else
        err = C * (2* rmax)^4 * rmax * E[2]
    end
    # % If the best null-rules are on the noise level, then put the error to zero.
    if (E[1] < noise) && (E[2] < noise)
        err = zero(T)
    end
    # % Check if the error is small enough or the interval is too small.
    if (err <= max(tol * abs(is), noise)) || (x[2] <= a) || (b <= x[6])
        if ((x[2] <= a) || (b <= x[6]))
            #("Interval too small: required tolerance may not be met.")
        end
        return i
    end
    if nEval == 0
        return i
    end
    if isnan(fa) || isnan(fb) || isnan(fc)
        return T(NaN)
    end
    return lobattor(f, a, c, fa, fc, is, tol, nEval-1) +
           lobattor(f, c, b, fc, fb, is, tol, nEval-1)
end
