text:
# module Ecto

module English

function pluralize(word::String)
    string(lowercase(word), "s")
end

function singularize(word::String)
    endswith(word, "es") ? word[1:end-2] : endswith(word, "s") ? word[1:end-1] : word
end

end # module Ecto.English


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Comparison" begin
    @testset "Refraction" begin
        λ = 0.550
        r1 = OpticalRay([0.0, 0.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r2 = OpticalRay([2.0, 2.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r3 = OpticalRay([5.0, 5.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ)
        r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ)

        # Test a number of rays under standard environmental conditions
        a = TestData.planoplano()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [2.0, 2.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [5.0, 5.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, 0.7192321011619232, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.234903851062, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [0.7200535362928893, -6.141047252697188, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.2448952769065, rtol = COMP_TOLERANCE)

        a = TestData.concaveplano()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [3.633723967570758, 3.633723967570758, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.0772722321026, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [9.153654757938119, 9.153654757938119, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.5695599263722, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, -3.401152468994745, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.1609946836496, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [-3.457199906282556, -10.32836827735406, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.5393056652617, rtol = COMP_TOLERANCE)

        a = TestData.doubleconcave()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [5.235579681684886, 5.235579681684886, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.2624909340242, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [13.42351905137007, 13.42351905137007, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.8131310483928, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, -6.904467939352202, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.3286918571586, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [-7.089764102262856, -14.61837033417989, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.4286280772859, rtol = COMP_TOLERANCE)

        a = TestData.convexplano()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [0.8867891519368289, 0.8867891519368289, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9698941263224, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [2.212067233969831, 2.212067233969831, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.8895474037682, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, 3.489759589087602, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.4310036252394, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [3.491858246934857, -3.331802834115089, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.344370622227, rtol = COMP_TOLERANCE)

        a = TestData.doubleconvex()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [-0.06191521590711035, -0.06191521590711035, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.987421926598, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [-0.2491105067897657, -0.2491105067897657, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.0110871422915, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, 5.639876913179362, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.6758349889372, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [5.705452331562673, -0.7679713587894854, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.6175821043825, rtol = COMP_TOLERANCE)
    end # testset Refraction

    @testset "Temperature/Pressure" begin
        λ = 0.550
        r1 = OpticalRay([0.0, 0.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r2 = OpticalRay([2.0, 2.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r3 = OpticalRay([5.0, 5.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ)
        r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ)

        # Test a component whose material does have a ΔT component
        a = TestData.doubleconvex(temperature = 40 * u"°C")
        @test isapprox(point(intersection(trace(a, r1, test = true))), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(point(intersection(trace(a, r2, test = true))), [-0.06214282132053373, -0.06214282132053373, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(point(intersection(trace(a, r3, test = true))), [-0.2497005807710933, -0.2497005807710933, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(point(intersection(trace(a, r4, test = true))), [0.0, 5.640416741927821, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(point(intersection(trace(a, r5, test = true))), [5.706006107804734, -0.7673622008537766, -67.8], rtol = COMP_TOLERANCE)

        # note that this material doesn't have a ΔT component
        λ2 = 0.533
        r1 = OpticalRay([0.0, 0.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ2)
        r2 = OpticalRay([2.0, 2.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ2)
        r3 = OpticalRay([5.0, 5.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ2)
        r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ2)
        r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ2)
    end # testset Temperature/Pressure

    @testset "Reflection" begin
        λ = 0.550
        r1 = OpticalRay([0.0, 0.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r2 = OpticalRay([2.0, 2.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r3 = OpticalRay([5.0, 5.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ)
        r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ)

        a = TestData.planoconcaverefl()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, 47.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 79.1704477524159, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [-6.19723306038804, -6.19723306038804, 47.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 80.0979229417755, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [-17.7546255175372, -17.7546255175372, 47.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 86.2063094599075, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, 19.5349836031196, 47.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 83.5364646376395, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [21.2477706541112, 17.3463704045055, 47.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 87.4033932087711, rtol = COMP_TOLERANCE)
    end # testset Reflection

    @testset "Complex Lenses" begin
        λ = 0.550
        r1 = OpticalRay([0.0, 0.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r2 = OpticalRay([2.0, 2.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r3 = OpticalRay([5.0, 5.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
        r4 = OpticalRay([0.0, -5.0, 1.0], [0.0, 0.08715574274765818, -0.9961946980917454], 1.0, λ)
        r5 = OpticalRay([-5.0, -5.0, 1.0], [0.08715574274765818, -0.01738599476176408, -0.9960429728140486], 1.0, λ)

        a = TestData.conicsystemZ()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [-1.80279270185495, -1.80279270185495, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.1000975813922, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [-2.8229241607807, -2.8229241607807, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.282094499601, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, 9.01421545841289, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.2211078071893, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [8.13946939328266, 2.23006981338816, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.1025133542546, rtol = COMP_TOLERANCE)

        a = TestData.asphericsystem()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [0.0735282671574837, 0.0735282671574837, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.0161411455851, rtol = COMP_TOLERANCE)
        track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0)
        res = trace(a, r3, test = true, trackrays = track)
        @test (res === nothing) # TIR
        @test isapprox(point(track[end]), [-5.0, -5.0, 0.0], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(track[end]), 46.5556716286238, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, 11.9748998399885, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 76.0760286320348, rtol = COMP_TOLERANCE)
        track = Vector{OpticSim.LensTrace{Float64,3}}(undef, 0)
        res = trace(a, r5, test = true, trackrays = track)
        @test (res === nothing) # TIR
        @test isapprox(point(track[end]), [5.49905367197174, 5.66882664623822, 0.0], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(track[end]), 47.6825931025333, rtol = COMP_TOLERANCE)

        a = TestData.zernikesystem()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, -8.9787010034042, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.5977029819277, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [-1.58696749235066, -9.71313213852721, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.7603266537023, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [0.790348081563859, -0.762155123619682, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.1933359669848, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true, trackrays = track)
        @test isapprox(point(res), [0.0, 10.0037899692172, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 76.8041236301678, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [35.5152289731456, 28.3819941055557, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 94.3059947823954, rtol = COMP_TOLERANCE)

        a = TestData.conicsystemQ()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [0.0, 0.0, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 73.9852238762079, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [-1.80279270185495, -1.80279270185495, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.1000975813922, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [-2.8229241607807, -2.8229241607807, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.282094499601, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [0.0, 9.01421545841289, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.2211078071893, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [8.13946939328266, 2.23006981338816, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.1025133542546, rtol = COMP_TOLERANCE)

        # No NSC qtype so have less precise SC values and no angled rays
        a = TestData.qtypesystem()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [-4.3551074306, -0.318112811010, -67.8], rtol = 1e-12)
        # TODO these are the values that I get for comparison - I'm pretty certain that they are just wrong...
        # res = trace(a, r2, test = true)
        # @test isapprox(point(res), [-0.091773202667, 3.9362695851, -67.8], rtol = 1e-12)
        # res = trace(a, r3, test = true)
        # @test isapprox(point(res), [-7.0993098547, -2.6848726242, -67.8], rtol = 1e-12)

        # No NSC chebyshev so have less precise SC values and no angled rays
        a = TestData.chebyshevsystem()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [-1.07963031980, 0.53981515992, -67.8], rtol = 1e-12)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [0.00851939011, 1.25870229260, -67.8], rtol = 1e-12)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [-1.20441411240, -0.07554605390, -67.8], rtol = 1e-12)

        a = TestData.gridsagsystem()
        res = trace(a, r1, test = true)
        @test isapprox(point(res), [21.0407756733608, 21.724638830759, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 82.0777022031378, rtol = COMP_TOLERANCE)
        res = trace(a, r2, test = true)
        @test isapprox(point(res), [-0.489183765274452, 0.405160352533666, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 75.536902213118, rtol = COMP_TOLERANCE)
        res = trace(a, r3, test = true)
        @test isapprox(point(res), [-12.0770793886528, -8.7705340321259, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 77.5651917266687, rtol = COMP_TOLERANCE)
        res = trace(a, r4, test = true)
        @test isapprox(point(res), [-1.20535362019229, 6.07973526939659, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.7349148204077, rtol = COMP_TOLERANCE)
        res = trace(a, r5, test = true)
        @test isapprox(point(res), [6.5112407340601, -0.22055440245024, -67.8], rtol = COMP_TOLERANCE)
        @test isapprox(pathlength(res), 74.6955888563913, rtol = COMP_TOLERANCE)
    end #testset complex lenses

    # @testset "Power" begin
    #     λ = 0.550
    #     r1 = OpticalRay([0.0, 0.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
    #     r2 = OpticalRay([2.0, 2.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
    #     r3 = OpticalRay([5.0, 5.0, 1.0], [0.0, 0.0, -1.0], 1.0, λ)
    #     a = TestData.doubleconvex()
    #     res = trace(a, r1, test = true)
    #     @test isapprox(power(res), 0.915508396)
    #     res = trace(a, r2, test = true)
    #     @test isapprox(power(res), 0.915526077)
    #     res = trace(a, r3, test = true)
    #     @test isapprox(power(res), 0.915577586)
    # end
end # testset Comparison

ext: jl


text:
"""
	REDCap.Config(url::String, key::String; ssl::Bool = true)

Struct to hold api url and key/superkey.
`APIConfigObj = Config("http...","ABCD...")`

This will be passed to all function calls to both orient and authorize the api_pusher() function. The REDCap API cannot
be accessed without this config object properly set-up. Always test your config object before automating a project.

#### NOTE: SSL should always be on by default.
If for whatever reason, your project will not SSL verify AND you must use it, disable ssl verification with `ssl=false`
Leaving SSL verification disabled leaves you open for Man-in-the-Middle attacks and is generally just bad practice.

#### Parameters:
* `url` - The url of the REDCap instance.
* `key` - Either the standard or super API key.
* `ssl` - Flag to enable ssl verification
"""
struct Config
	url::String
	key::String
	ssl::Bool
	#basic validation - checks that the url starts and ends 'properly', and then checks the key length
	function Config(url::String, key::String; ssl::Bool = true)
		if isequal(url[end-4:end], "/api/") && (isequal(url[1:7], "http://") || isequal(url[1:8], "https://"))
			if (length(key)==32 || length(key)==64)
				new(url, key, ssl)
			else
				@error("Invalid Key: $key \nMust be 32 characters long for a standard key, or 64 characters long for a super key.")
			end
		else
			@error("Invalid URL: $url \nMust be in format of http(s)://<redcap-hosting-url>/api/")
		end
	end
end


text:
using HDF5
using Healpix
using Random
using Random: MersenneTwister

"""
    read_halo_catalog_hdf5(filename)

Utility function to read an HDF5 table with x, y, z, M_h as the four rows.
The hdf5 record is "halos". This is a format we use to distribute Websky halos.

## Example
```julia-repl
julia> halo_pos, halo_mass = read_halo_catalog_hdf5(
    "/global/cfs/cdirs/sobs/www/users/Radio_WebSky/websky_halos-light.hdf5")
```
"""
function read_halo_catalog_hdf5(filename)
    hdata = h5open(filename, "r") do file
        read(file, "halos")
    end
    pos = hdata[1:3, :]
    halo_mass = hdata[4, :]
    return pos, halo_mass
end

"""
    chunk(arr_len, chunksize::Integer)

Generates a list of tuples which describe starting and ending chunk indices.
Useful for parallelizing an array operation.
"""
function chunk(arr_len, chunksize::Integer)
    return [
        (i, min(i + chunksize - 1, arr_len)) for
        i in range(1, arr_len, step = chunksize)
    ]
end

function getrange(n)
    tid = Threads.threadid()
    nt = Threads.nthreads()
    d, r = divrem(n, nt)
    from = (tid - 1) * d + min(r, tid - 1) + 1
    to = from + d - 1 + (tid ≤ r ? 1 : 0)
    from:to
end

function threaded_rand!(arr::Array{T,1}; chunksize = 4096) where {T}

    num = size(arr, 1)
    Threads.@threads for (i1, i2) in chunk(num, chunksize)
        @views rand!(arr[i1:i2])
    end
end

"""
Generate an array where the value at index i corresponds to the index of the
first source of halo i. Takes an array where the value at index i corresponds
to the number of subhalos that halo i has.
"""
function generate_subhalo_offsets(num_subhalos)
    result = cumsum(num_subhalos)
    prepend!(result, 0)
    return result
end

"""
    get_basic_halo_properties(halo_pos::Array{T,2}, model::AbstractForegroundModel,
                              cosmo::Cosmology.FlatLCDM{T}, res::Resolution) where T

Compute distance, redshift, and healpix indices from halo positions and a cosmology.
"""
function get_basic_halo_properties(
    halo_pos::Array{T,2},
    model::AbstractForegroundModel,
    cosmo::Cosmology.FlatLCDM{T},
    res::Resolution,
) where {T}
    N_halos = size(halo_pos, 2)
    hp_ind = Array{Int64}(undef, N_halos)  # healpix index of halo
    redshift = Array{T}(undef, N_halos)
    dist = Array{T}(undef, N_halos)

    r2z = build_r2z_interpolator(model.min_redshift, model.max_redshift, cosmo)
    Threads.@threads for i = 1:N_halos
        dist[i] = sqrt(halo_pos[1, i]^2 + halo_pos[2, i]^2 + halo_pos[3, i]^2)
        redshift[i] = r2z(dist[i])
        hp_ind[i] = Healpix.vec2pixRing(
            res,
            halo_pos[1, i],
            halo_pos[2, i],
            halo_pos[3, i],
        )
    end

    return dist, redshift, hp_ind
end

"""
Compute angles of halos
"""
function get_angles(halo_pos::Array{T,2}) where {T}
    N_halos = size(halo_pos, 2)
    θ = Array{T}(undef, N_halos)
    ϕ = Array{T}(undef, N_halos)

    Threads.@threads for i = 1:N_halos
        θ[i], ϕ[i] =
            Healpix.vec2ang(halo_pos[1, i], halo_pos[2, i], halo_pos[3, i])
    end

    return θ, ϕ
end

"""
Utility function which prepends some zeros to an array. It makes a copy instead
of modifying the input.
"""
function ellpad(arr::Array{T,N}; nzeros = 1) where {T,N}
    result = arr[:]
    pushfirst!(result, zeros(T, nzeros)...)
    return result
end

function catalog2map!(m::Map{T,RingOrder}, flux, theta, phi) where {T}
    res = m.resolution
    pixel_array = m.pixels
    N_halo = length(flux)

    # try to prevent thread issues by sorting by theta
    perm = sortperm(theta, rev = true, alg = ThreadsX.MergeSort)
    Threads.@threads for i_perm = 1:N_halo
        i_halo = perm[i_perm]
        hp_ind = Healpix.ang2pixRing(res, theta[i_halo], phi[i_halo])
        pixel_array[hp_ind] += flux[i_halo]
    end

    # divide by healpix pixel size
    per_pixel_steradian = 1 / nside2pixarea(res.nside)
    pixel_array .*= per_pixel_steradian
end

import PhysicalConstants.CODATA2018: c_0, k_B, PlanckConstant

function BlackBody(
    ν::Unitful.Frequency,
    Tbb::Unitful.Temperature;
    unit_out::Unitful.FreeUnits = u"Jy/sr",
)
    hP = uconvert(u"erg/GHz", PlanckConstant)
    kB = uconvert(u"erg/K", k_B)
    c = uconvert(u"cm*GHz", c_0)
    x = hP * ν / (kB * Tbb)
    emission = 2 * hP * ν^3 / c^2 / (exp(x) - 1) / 1u"sr"
    emission = uconvert(unit_out, emission)
    return emission
end

function BlackBody(
    λ::Unitful.Length,
    Tbb::Unitful.Temperature;
    unit_out::Unitful.FreeUnits = u"Jy/sr",
)
    hP = uconvert(u"erg*s", PlanckConstant)
    kB = uconvert(u"erg/K", k_B)
    c = uconvert(u"μm/s", c_0)
    x = hP * c / (λ * kB * Tbb)
    ν = uconvert(u"GHz", c / λ)

    emission = (2hP * c^2 / λ^5 / (exp(x) - 1))c / ν^2 / 1u"sr"
    emission = uconvert(unit_out, emission)
    return emission

end

"""
Fixsen et al. (1998) modeled the Galactic and interplanetary foreground emission in the FIRAS data
to determine the spectrum of the cosmic infrared background.
More details can be found in:
"The Spectrum of the Extragalactic Far-Infrared Background from the COBE FIRAS Observations"
D.J. Fixsen, et al., 1998, ApJ 508, 123

We implemented multiple definition of the model depending whether we are estimating it in units
of frequency or wavelength.
"""

function fixsen_model(
    x::Unitful.Frequency;
    Tbb::Unitful.Temperature = 18.5u"K",
    unit_out::Unitful.FreeUnits = u"MJy/sr",
)
    B = BlackBody(x, Tbb)
    z0 = 100 / u"cm"
    x2 = uconvert(u"cm*GHz", c_0) / x
    z = uconvert(unit(z0), 1.0 / x2)
    model = (1.3E-5 * (z / z0)^0.64 * B)
    model = uconvert(unit_out, model)
    return model
end

function fixsen_model(
    x::Unitful.Length;
    Tbb::Unitful.Temperature = 18.5u"K",
    unit_out::Unitful.FreeUnits = u"MJy/sr",
)
    B = BlackBody(x, Tbb)
    z0 = 100 / u"cm"
    z = uconvert(unit(z0), 1.0 / x)
    model = (1.3E-5 * (z / z0)^0.64 * B)
    model = uconvert(unit_out, model)

    return model
end

function write_to_hdf5(filename, fluxes_cen, fluxes_sat)
    hdata = h5open(filename, "w")

    write(hdata, "fluxes_cen", fluxes_cen)

    write(hdata, "fluxes_sat", fluxes_sat)
    close(hdata)
end

function write_to_hdf5(filename, fluxes_cen, fluxes_sat, sources)
    hdata = h5open(filename, "w")

    write(hdata, "fluxes_cen", fluxes_cen)

    write(hdata, "fluxes_sat", fluxes_sat)
    write(hdata, "redshift_cen", sources.redshift_cen)

    write(hdata, "redshift_sat", sources.redshift_sat)
    write(hdata, "lum_cen", sources.lum_cen)

    write(hdata, "lum_sat", sources.lum_sat)
    close(hdata)
end

function write_to_hdf5(filename, fluxes_cen, fluxes_sat, sources, m)
    hdata = h5open(filename, "w")

    write(hdata, "fluxes_cen", fluxes_cen)

    write(hdata, "fluxes_sat", fluxes_sat)
    write(hdata, "redshift_cen", sources.redshift_cen)
    write(hdata, "redshift_sat", sources.redshift_sat)

    write(hdata, "healpix_index_cen", sources.hp_ind_cen)
    write(hdata, "healpix_index_sat", sources.hp_ind_sat)

    write(hdata, "lum_cen", sources.lum_cen)
    write(hdata, "lum_sat", sources.lum_sat)
    write(hdata, "healpix_map", m)

    close(hdata)
end

function read_fluxes_hdf5(filename)
    hdata = h5open(filename, "r")

    fcen = hdata["fluxes_cen"][:]
    fsat = hdata["fluxes_cen"][:]
    close(hdata)
    return fcen, fsat
end

function read_from_hdf5(filename, datalist)
    hdata = h5open(filename, "r")

    datadic = Dict(string => hdata[string][:] for string in datalist)
    close(hdata)
    return datadic
end


text:
"""
Slit cylinder under torsional loads.

Reference
Morris, A. J., NAFEMS: shell finite element evaluation tests
Reference P12
Section 3.3 In Extensional Bending
Slit Cylinder Under Applied Torque T = 4*pi*D*(1-nu)*C/R
x = axial coordinate
displacements: axial u = -C phi, circumferential v = C x, radial w = 0,
twisting moment M_x,phi =  D*(1-nu)*C/R^2
all other stress resultants (section forces and moments) are zero.

Note: our coordinates are different. The axial coordinate is y.
"""
module slit_cylinder_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

E = 3e6;
nu = 0.3;
R = 300.0;
L = 3*R;
thickness = 3.0;

# T = 4*pi*D*(1-nu)*C/R  -> C = T*R/(4*pi*D*(1-nu))
T = 1000.0
D = E*thickness^3/12/(1-nu^2)
C = T*R/(4*pi*D*(1-nu))
# The maximum displacement axially is:
vmax = C * pi
# Twisting moment magnitude
M_x_phi =  D*(1-nu)*C/R^2

cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
    r = vec(XYZ); r[2] = 0.0;
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 2] .= (0.0, 1.0, 0.0) #  this is along the axis
    cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
    return csmatout
end

function _execute(n = 2, visualize = true)
    formul = FEMMShellT3FFModule

    n = Int(round(n/2)*2)
    @info "Mesh: $n elements per side"

    # Mesh
    tolerance = min(R, L) / n  / 100

    tolerance = R/n/1000
    fens, fes = T3block(360.0, L, n, n);
    bfes = meshboundary(fes)
    al0 = selectelem(fens, bfes, box = [0.0 0.0 -Inf Inf], inflate = tolerance)
    al360 = selectelem(fens, bfes, box = [360.0 360.0 -Inf Inf], inflate = tolerance)
    ll0 = selectelem(fens, bfes, box = [-Inf Inf 0.0 0.0], inflate = tolerance)
    llL2 = selectelem(fens, bfes, box = [-Inf Inf L L], inflate = tolerance)
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a/180*pi), y-L/2, R*cos(a/180*pi))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)
    ocsys = CSys(3, 3, cylindrical!)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), ocsys, mater)
    femm.drilling_stiffness_scale = 0.1
    associategeometry! = formul.associategeometry!
    stiffness = formul.stiffness

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Fix one of the nodes to remove rigid body displacements
    @show     nl = selectnode(fens, nearestto = Float64[0.0, L/2, -R])
    for i in 1:6
        setebc!(dchi, nl, true, i)
    end

    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    loadbdry = subset(bfes, ll0)
    lfemm = FEMMBase(IntegDomain(loadbdry, GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0, 0, 0, T/(2*pi*R), 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 1);
    loadbdry = subset(bfes, llL2)
    lfemm = FEMMBase(IntegDomain(loadbdry, GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0, 0, 0, -T/(2*pi*R), 0]);
    F += distribloads(lfemm, geom0, dchi, fi, 1);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])

    @show vmax, maximum(dchi.values[connectednodes(subset(bfes, llL2)), 2])

    fld = fieldfromintegpoints(femm, geom0, dchi, :moment, 3, outputcsys = ocsys)
    @show -M_x_phi, minimum(fld.values[:]), maximum(fld.values[:])

        # Generate a graphical display of resultants
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
        push!(scalars, ("m$nc", fld.values))
        fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
        push!(scalars, ("em$nc", fld.values))
    end
    vtkwrite("slit_cylinder-$(n)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
        push!(scalars, ("n$nc", fld.values))
        fld = elemfieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
        push!(scalars, ("en$nc", fld.values))
    end
    vtkwrite("slit_cylinder-$(n)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:2
        fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
        push!(scalars, ("q$nc", fld.values))
        fld = elemfieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
        push!(scalars, ("eq$nc", fld.values))
    end
    vtkwrite("slit_cylinder-o-$(n)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    vtkwrite("slit_cylinder-$(n)-uur.vtu", fens, fes; vectors = [("u", dchi.values[:, 1:3]), ("ur", dchi.values[:, 4:6])])

    return true
end

function test_convergence()
    @info "Slit cylinder"
    for n in [8, 16, 32] # 3:64 #
        _execute(n, !false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
@testset "misc-html" begin
    fs()
    λ = "blah/blah.ext"
    set_curpath("cpB/blah.md")
    @test F.html_ahref(λ, 1)           == "1"
    @test F.html_ahref(λ, "bb")        == "bb"
    @test F.html_ahref_key("cc", "dd") == "dd"
end

@testset "misc-html 2" begin
   h = "
blah
"
   @test !F.is_html_escaped(h)
   @test F.html_code(h, "html") == """
<div class="foo">blah</div>
"""
   he = Markdown.htmlesc(h)
   @test F.is_html_escaped(he)
   @test F.html_code(h, "html") == F.html_code(h, "html")
end

@testset "html/hide" begin
    c = """
        a=5
        b=7
        """
    @test F.html_skip_hidden(c, "foo") == c
    c = """
        #hideall
        a=5
        b=7
        """
    @test F.html_skip_hidden(c, "julia") == ""
    c = """
        a=5 #hide
        b=7
        """
    @test F.html_skip_hidden(c, "julia") == "b=7"
end

@testset "html_code" begin
    c = """
        using Random
        Random.seed!(555) # hide
        a = randn()
        b = a + 5
        """
    @test F.html_code(c, "julia") ==
        """
$(F.htmlesc("""using Random
        a = randn()
        b = a + 5"""))
"""
end

@testset "html_content" begin
    h = F.html_content("div", "foo bar"; class="container", id="body")
    @test h == "
foo bar
"
end

@testset "fd2text" begin
    h = F.fd2text(
        """
        A **B** `C` [D](https://julialang.org)
        - E
        - F
        G < H
        """
    )
    @test isapproxstr(h, """
        A B C D E F G < H
        """)
end


text:
"""
### Nuclear attraction term
"""
function Aterm(i::Int64,r::Int64,u::Int64,l1::Int64,l2::Int64,ax::Float64,bx::Float64,cx::Float64,gamma::Float64)
    term1 = (-1)^(i+u)*binomial_prefactor(i,l1,l2,ax,bx)
    term2 = factorial(i)*cx^(i-2r-2u)
    term3 = (1/4/gamma)^(r+u)/factorial(r)/factorial(u)/factorial(i-2r-2u)
    return term1*term2*term3
end

function Aarray(l1::Int64,l2::Int64,a::Float64,b::Float64,c::Float64,g::Float64)
    Imax = l1+l2+1
    A = zeros(Float64,Imax)
    for i in 0:(Imax-1)
        for r in 0:div(i,2)
            for u in 0:div(i-2r,2)
                I = i-2r-u+1
                A[I] += Aterm(i,r,u,l1,l2,a,b,c,g)
            end
        end
    end
    return A
end

function nuclear_attraction(aexpn::Float64,ax::Float64,ay::Float64,az::Float64,
                            aI::Int64,aJ::Int64,aK::Int64,
                            bexpn::Float64,bx::Float64,by::Float64,bz::Float64,
                            bI::Int64,bJ::Int64,bK::Int64,
                            cx::Float64,cy::Float64,cz::Float64)
    px,py,pz = gaussian_product_center(aexpn,ax,ay,az,bexpn,bx,by,bz)
    gamma = aexpn+bexpn
    rab2 = dist2(ax-bx,ay-by,az-bz)
    rcp2 = dist2(cx-px,cy-py,cz-pz)
    Ax = Aarray(aI,bI,px-ax,px-bx,px-cx,gamma)
    Ay = Aarray(aJ,bJ,py-ay,py-by,py-cy,gamma)
    Az = Aarray(aK,bK,pz-az,pz-bz,pz-cz,gamma)
    total = 0.0
    for I in 0:(aI+bI)
        for J in 0:(aJ+bJ)
            for K in 0:(aK+bK)
                total += Ax[I+1]*Ay[J+1]*Az[K+1]*Fgamma(I+J+K,rcp2*gamma)
            end
        end
    end
    val=-2pi*exp(-aexpn*bexpn*rab2/gamma)*total/gamma
    #println(val)
    #println((Ax,Ay,Az,rcp2*gamma,Fgamma(0,rcp2*gamma)))
    return val
end

function nuclear_attraction(a::PGBF,b::PGBF,cx::Float64,cy::Float64,cz::Float64)
    return a.norm*b.norm*nuclear_attraction(a.expn,a.x,a.y,a.z,a.I,a.J,a.K,
                                            b.expn,b.x,b.y,b.z,b.I,b.J,b.K,cx,cy,cz)
end
nuclear_attraction(a::PGBF,b::PGBF,c::Atom) = c.atno*nuclear_attraction(a,b,c.x,c.y,c.z)
nuclear_attraction(a::PGBF,b::PGBF,m::Molecule) = sum([nuclear_attraction(a,b,c) for c in m.atomlist])

# Need a nested scope to squeeze this into the contract function
function nuclear_attraction(a::CGBF,b::CGBF,cx::Float64,cy::Float64,cz::Float64)
    na(a,b) = nuclear_attraction(a,b,cx,cy,cz)
    contract(na,a,b)
end
function nuclear_attraction(a::CGBF,b::CGBF,c::Atom)
    na(a,b) = nuclear_attraction(a,b,c)
    contract(na,a,b)
end
function nuclear_attraction(a::CGBF,b::CGBF,m::Molecule)
    na(a,b) = nuclear_attraction(a,b,m)
    contract(na,a,b)
end

function test_a_terms()
    @testset "a_terms" begin
    @test Aterm(0,0,0,0,0,0.,0.,0.,0.) == 1.0
    @test Aarray(0,0,0.,0.,0.,1.) == [1.0]
    @test Aarray(0,1,1.,1.,1.,1.) == [1.0, -1.0]
    @test Aarray(1,1,1.,1.,1.,1.) == [1.5, -2.5, 1.0]
    @test Aterm(0,0,0,0,0,0.,0.,0.,1.) == 1.0
    @test Aterm(0,0,0,0,1,1.,1.,1.,1.) == 1.0
    @test Aterm(1,0,0,0,1,1.,1.,1.,1.) == -1.0
    @test Aterm(0,0,0,1,1,1.,1.,1.,1.) == 1.0
    @test Aterm(1,0,0,1,1,1.,1.,1.,1.) == -2.0
    @test Aterm(2,0,0,1,1,1.,1.,1.,1.) == 1.0
    @test Aterm(2,0,1,1,1,1.,1.,1.,1.) == -0.5
    @test Aterm(2,1,0,1,1,1.,1.,1.,1.) == 0.5
    end #testset
end

function test_na()
    @testset "test_na" begin
    s = pgbf(1.0)
    c = cgbf(0.0,0.0,0.0)
    push!(c,1,1)
    @test amplitude(c,0,0,0) ≈ 0.71270547
    @test nuclear_attraction(s,s,0.,0.,0.) ≈ -1.59576912
    @test nuclear_attraction(c,c,0.,0.,0.) ≈ -1.59576912
    end #testset
end
#todo make into a test
function test_one()
    s1 = pgbf(1)
    s2 = pgbf(1,0,1,0)
    x=y=0.
    println("S: $(overlap(s1,s2))")
    println("T: $(kinetic(s1,s2))")
    for z in range(0,stop=1,length=5)
        println("V: $z $(nuclear_attraction(s1,s2,x,y,z))")
    end
end

function test_na2()
    @testset "test_na2" begin
    li,h = lih.atomlist
    bfs = build_basis(lih)
    s1,s2,x,y,z,h1s = bfs.bfs
    @test nuclear_attraction(s1,s1,lih) ≈ -8.307532656
    end #testset
end


text:
#############################################################################
# SCS.jl
# Wrapper around the SCS solver https://github.com/cvxgrp/scs
# See http://github.com/JuliaOpt/SCS.jl
#############################################################################
# test/options.jl
# Tests the ability to pass options
#############################################################################

using Base.Test
using MathProgBase.SolverInterface
using SCS

# The normal test
A = [1.0 1.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 1.0;
     0.0 0.0 1.0 1.0 1.0]
collb = [0.0, 0.0, 0.0, 0.0, 0.0]
obj   = [3.0, 4.0, 4.0, 9.0, 5.0]
rowub = [ 5.0,  3.0,  9.0]
for s in [SCSSolver(), SCSSolver(linearsolver=SCS.Indirect)]
    m = MathProgBase.ConicModel(s)
    MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
    MathProgBase.optimize!(m)
    @test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-3)
end

# With eps = 1e-8, solution should be far more accurate
for s in [SCSSolver(eps=1e-8), SCSSolver(eps=1e-8, linearsolver=SCS.Indirect)]
    m = MathProgBase.ConicModel(s)
    MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
    MathProgBase.optimize!(m)
    @test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-5)

    # With a warmstart from the eps = 1e-8 solution, solution should be extremely accurate even after 1 iteration
    push!(m.options, (:warm_start, true))
    push!(m.options, (:max_iters, 1))
    MathProgBase.optimize!(m)
    @test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-5)

    # Now let's do the same warmstart, but on a new instance of the same problem
    primal_sol = m.primal_sol
    dual_sol = m.dual_sol
    slack = m.slack

    s = SCSSolver(max_iters=1)
    m = MathProgBase.ConicModel(s)
    MathProgBase.loadproblem!(m, -obj, A, rowub, [(:NonNeg,1:3)],[(:NonNeg,1:5)])
    MathProgBase.setwarmstart!(m, primal_sol; dual_sol = dual_sol, slack = slack)
    MathProgBase.optimize!(m)
    @test isapprox(MathProgBase.getobjval(m), -99.0, atol=1e-5)
end


text:
using DataFrames
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = test_page_model()
    include("../src/components/TotalAbatementCosts.jl")

    add_comp!(m, TotalAbatementCosts)

    set_param!(m, :TotalAbatementCosts, :pop_population, readpagedata(m, "test/validationdata/$valdir/pop_population.csv"))
    set_param!(m, :TotalAbatementCosts, :tc_totalcosts_co2, readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_co2.csv"))
    set_param!(m, :TotalAbatementCosts, :tc_totalcosts_ch4, readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_ch4.csv"))
    set_param!(m, :TotalAbatementCosts, :tc_totalcosts_n2o, readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_n2o.csv"))
    set_param!(m, :TotalAbatementCosts, :tc_totalcosts_linear, readpagedata(m, "test/validationdata/$valdir/tc_totalcosts_linear.csv"))

    run(m)

    # Generated data
    abate_cost = m[:TotalAbatementCosts, :tct_totalcosts]
    abate_cost_per_cap = m[:TotalAbatementCosts, :tct_per_cap_totalcostspercap]

    # Recorded data
    cost_compare = readpagedata(m, "test/validationdata/$valdir/tct_totalcosts.csv")
    cost_cap_compare = readpagedata(m, "test/validationdata/$valdir/tct_per_cap_totalcostspercap.csv")

    @test abate_cost ≈ cost_compare rtol = 1e-4
    @test abate_cost_per_cap ≈ cost_cap_compare rtol = 1e-7
end


text:
module LIFNeuron
export LIFStates, LIFParams, evolve, update_dc
using ...Signals
using Printf

mutable struct LIFStates{T <: AbstractFloat}
    v::T
    t_refractory::T
    dc::T
    v_equilibrium::T
end

function gen_v_eqlbrm(dc::T, tau_refractory::T, v_steady::T) where {T <: AbstractFloat}
    dc * tau_refractory + v_steady
end

function LIFStates{T}(v::T,
                      v_steady::T) where {T <: AbstractFloat}
    LIFStates(v, zero(T), zero(T), v_steady)
end

## seems useless?
# function LIFStates{T}(v::T,
#                       t_refractory::T,
#                       dc::T,
#                       tau_refractory::T,
#                       v_steady::T) where {T <: AbstractFloat}
#     LIFStates(v, t_refractory, dc, gen_v_eqlbrm(dc, tau_refractory, v_steady))
# end

struct LIFParams{T <: AbstractFloat}
    v_steady::T
    v_reset::T
    tau_leak::T
    tau_refractory::T
    lazy_threshold::T
end

function update_dc(states::LIFStates{T},
                   params::LIFParams{T},
                   instruction::DCInstruction{T}) where {T <: AbstractFloat}
    states.dc = states.dc - instruction.previous + instruction.new
    states.v_equilibrium = gen_v_eqlbrm(states.dc, params.tau_refractory, params.v_steady)
end

# return (fired::Bool, triggered::Bool, states::LIFStates) ?
function evolve(dt::T,
                states::LIFStates,
                params::LIFParams,
                delta_v::T) where {T <: AbstractFloat}
    if states.t_refractory <= zero(T)
        new_v = states.v + dt * ((states.v_equilibrium - states.v) / params.tau_leak) + delta_v
        if new_v >= 30.0
            true, true, LIFStates(params.v_reset, params.tau_refractory, states.dc, states.v_equilibrium)
        else
            new_states = LIFStates(new_v, zero(T), states.dc, states.v_equilibrium)
            if abs(new_v - states.v_equilibrium) > params.lazy_threshold
                false, true, new_states
            else
                false, false, new_states
            end
        end
    else # in refractory period.
        false, true, LIFStates(states.v, states.t_refractory - dt, states.dc, states.v_equilibrium)
    end
end

end # module end


text:
export applytrap1, applytrap2, applytrap3, squaresub
"""
    applytrap1(square)

Returns the next square after trap 1, meaning going back to square 1,
has been applied when being at `square`.
"""
function applytrap1(square::Int64)
    return squaresub(square, 15)
end

"""
    applytrap2(square)

Returns the next square after trap 2, meaning going back 3 squares,
has been applied when being at `square`.
"""
function applytrap2(square::Int64)
    return squaresub(square, 3)
end

"""
    applytrap3(square)

Returns the next square after trap 3 (the magic trap),
meaning that one can go to every one of the 15 other squares
with equal probability, has been applied when being at `square`.
"""
function applytrap3(square::Int64)
    return rand(1:15)
end

"""
    squaresub(square, k)

Returns the next square after going back `k` squares when being at `square`.
"""
function squaresub(square::Int64, k::Int64)
    if k == 0
        return square
    elseif square <= k+1
        return 1
    end
    if square == 11
        return squaresub(3, k-1)
    else
        return squaresub(square-1, k-1)
    end
end

max_stars_repo_name: qlete/markov-decision

max_issues_repo_name: qlete/markov-decision

max_forks_repo_name: qlete/markov-decision


text:
module M2_2
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    1.3806522621004e+17*l1.^18 - 8.21485647116353e+17*l1.^17.*l2 - 1.60065545335102e+18*l1.^17.*l3 + 9.45494569700184e+15*l1.^17 + 2.95872947164765e+18*l1.^16.*l2.^2 + 7.83577706326508e+18*l1.^16.*l2.*l3 - 3.18029991626426e+16*l1.^16.*l2 + 9.37824823540379e+18*l1.^16.*l3.^2 - 9.28303759341999e+16*l1.^16.*l3 + 6.98362956528783e+15*l1.^16 - 7.18076310961731e+18*l1.^15.*l2.^3 - 2.54077238493916e+19*l1.^15.*l2.^2.*l3 + 5.84487552178296e+16*l1.^15.*l2.^2 - 3.65908602424017e+19*l1.^15.*l2.*l3.^2 + 2.06289724298222e+17*l1.^15.*l2.*l3 - 3.41077122580316e+16*l1.^15.*l2 - 3.68347640474837e+19*l1.^15.*l3.^3 + 4.43522907241177e+17*l1.^15.*l3.^2 - 7.91889848342627e+16*l1.^15.*l3 - 233793607153361.0*l1.^15 + 1.32045993751291e+19*l1.^14.*l2.^4 + 5.4915970224221e+19*l1.^14.*l2.^3.*l3 + 3.43816207163712e+15*l1.^14.*l2.^3 + 1.07152524912614e+20*l1.^14.*l2.^2.*l3.^2 - 3.09434586447333e+17*l1.^14.*l2.^2.*l3 + 1.06766736869272e+17*l1.^14.*l2.^2 + 1.10579520304018e+20*l1.^14.*l2.*l3.^3 - 4.22893934811355e+17*l1.^14.*l2.*l3.^2 + 3.02134632655923e+17*l1.^14.*l2.*l3 + 607863378598739.0*l1.^14.*l2 + 1.0820483759454e+20*l1.^14.*l3.^4 - 1.36838850451154e+18*l1.^14.*l3.^3 + 4.49463089102775e+17*l1.^14.*l3.^2 + 1.98390575212995e+15*l1.^14.*l3 - 132468346574870.0*l1.^14 - 1.8585176087238e+19*l1.^13.*l2.^5 - 9.34627626673762e+19*l1.^13.*l2.^4.*l3 - 2.28637777763863e+17*l1.^13.*l2.^4 - 1.98448951614853e+20*l1.^13.*l2.^3.*l3.^2 - 2.88805614017511e+17*l1.^13.*l2.^3.*l3 - 2.13391482686037e+17*l1.^13.*l2.^3 - 3.0095437439064e+20*l1.^13.*l2.^2.*l3.^3 + 2.88805614017511e+17*l1.^13.*l2.^2.*l3.^2 - 8.60794364302447e+17*l1.^13.*l2.^2.*l3 - 935174428613445.0*l1.^13.*l2.^2 - 2.37353371685459e+20*l1.^13.*l2.*l3.^4 - 4.81342690029185e+17*l1.^13.*l2.*l3.^3 - 1.26287635895607e+18*l1.^13.*l2.*l3.^2 - 2.61848840011765e+15*l1.^13.*l2.*l3 + 499081570986659.0*l1.^13.*l2 - 2.51507726600378e+20*l1.^13.*l3.^5 + 2.99635824543167e+18*l1.^13.*l3.^4 - 1.6947195767173e+18*l1.^13.*l3.^3 - 8.18277625036764e+15*l1.^13.*l3.^2 + 1.3789396743782e+15*l1.^13.*l3 + 3519658997501.61*l1.^13 + 2.0769731372755e+19*l1.^12.*l2.^6 + 1.21240761326154e+20*l1.^12.*l2.^5.*l3 + 6.37779064288669e+17*l1.^12.*l2.^5 + 3.09219480922851e+20*l1.^12.*l2.^4.*l3.^2 + 1.44402807008755e+18*l1.^12.*l2.^4.*l3 + 3.05276339280265e+17*l1.^12.*l2.^4 + 4.54712508554322e+20*l1.^12.*l2.^3.*l3.^3 + 3.17686175419262e+18*l1.^12.*l2.^3.*l3.^2 + 1.54986666970204e+18*l1.^12.*l2.^3.*l3 - 467587214306723.0*l1.^12.*l2.^3 + 6.39393825022315e+20*l1.^12.*l2.^2.*l3.^4 + 1.39589380108463e+18*l1.^12.*l2.^2.*l3.^3 + 3.27784475884175e+18*l1.^12.*l2.^2.*l3.^2 + 2.24441862867227e+15*l1.^12.*l2.^2.*l3 - 1.30279030323097e+15*l1.^12.*l2.^2 + 3.68436973916698e+20*l1.^12.*l2.*l3.^5 + 5.30680315757176e+18*l1.^12.*l2.*l3.^4 + 3.2871304746453e+18*l1.^12.*l2.*l3.^3 - 46758721430676.2*l1.^12.*l2.*l3.^2 - 3.93806995623569e+15*l1.^12.*l2.*l3 - 8338520980658.15*l1.^12.*l2 + 4.78232475724381e+20*l1.^12.*l3.^6 - 4.76529263128893e+18*l1.^12.*l3.^5 + 4.72167056515959e+18*l1.^12.*l3.^4 + 2.13687356938172e+16*l1.^12.*l3.^3 - 7.07202508447977e+15*l1.^12.*l3.^2 - 23385653741789.2*l1.^12.*l3 + 3998444404438.65*l1.^12 - 1.84438294242929e+19*l1.^11.*l2.^7 - 1.26547873036733e+20*l1.^11.*l2.^6.*l3 - 9.62685380058369e+17*l1.^11.*l2.^6 - 3.59353173407429e+20*l1.^11.*l2.^5.*l3.^2 - 3.46566736821013e+18*l1.^11.*l2.^5.*l3 - 2.76001641166066e+17*l1.^11.*l2.^5 - 6.5127311576983e+20*l1.^11.*l2.^4.*l3.^3 - 5.99271649086335e+18*l1.^11.*l2.^4.*l3.^2 - 2.26055195443658e+18*l1.^11.*l2.^4.*l3 + 3.41338666443908e+15*l1.^11.*l2.^4 - 7.35327068761165e+20*l1.^11.*l2.^3.*l3.^4 - 1.34775953208172e+19*l1.^11.*l2.^3.*l3.^3 - 4.7369484344801e+18*l1.^11.*l2.^3.*l3.^2 + 1.23443024576975e+16*l1.^11.*l2.^3.*l3 + 1.95942473343314e+15*l1.^11.*l2.^3 - 1.09999146791926e+21*l1.^11.*l2.^2.*l3.^5 - 4.98189684180206e+18*l1.^11.*l2.^2.*l3.^4 - 8.36258505812447e+18*l1.^11.*l2.^2.*l3.^3 + 5.33049424309664e+15*l1.^11.*l2.^2.*l3.^2 + 9.86024803137302e+15*l1.^11.*l2.^2.*l3 + 15755788935124.7*l1.^11.*l2.^2 - 3.89142937668125e+20*l1.^11.*l2.*l3.^6 - 1.70395312270331e+19*l1.^11.*l2.*l3.^5 - 5.62835257649039e+18*l1.^11.*l2.*l3.^4 + 3.10477910299664e+16*l1.^11.*l2.*l3.^3 + 1.38086438475336e+16*l1.^11.*l2.*l3.^2 + 12141642447757.2*l1.^11.*l2.*l3 - 13986463157622.5*l1.^11.*l2 - 7.59238457459379e+20*l1.^11.*l3.^7 + 5.24663532131811e+18*l1.^11.*l3.^6 - 1.02572123300015e+19*l1.^11.*l3.^5 - 3.81551166874286e+16*l1.^11.*l3.^4 + 2.38591285298109e+16*l1.^11.*l3.^3 + 78259263481361.3*l1.^11.*l3.^2 - 34308597689336.4*l1.^11.*l3 - 47059199054.2633*l1.^11 + 1.44307302506765e+19*l1.^10.*l2.^8 + 9.36708155517112e+19*l1.^10.*l2.^7.*l3 + 9.35180083485273e+17*l1.^10.*l2.^7 + 3.84531098732033e+20*l1.^10.*l2.^6.*l3.^2 + 4.90969543829768e+18*l1.^10.*l2.^6.*l3 + 1.04950401704453e+17*l1.^10.*l2.^6 + 5.71831001714585e+20*l1.^10.*l2.^5.*l3.^3 + 1.04692035081348e+19*l1.^10.*l2.^5.*l3.^2 + 2.35799829712165e+18*l1.^10.*l2.^5.*l3 - 7.24760182175421e+15*l1.^10.*l2.^5 + 1.10332795442878e+21*l1.^10.*l2.^4.*l3.^4 + 1.62693829229864e+19*l1.^10.*l2.^4.*l3.^3 + 6.44000255399745e+18*l1.^10.*l2.^4.*l3.^2 - 2.94112357798928e+16*l1.^10.*l2.^4.*l3 - 1.97254682230363e+15*l1.^10.*l2.^4 + 7.79244446676221e+20*l1.^10.*l2.^3.*l3.^5 + 3.31404442085094e+19*l1.^10.*l2.^3.*l3.^4 + 8.63770835934457e+18*l1.^10.*l2.^3.*l3.^3 - 6.49946227886344e+16*l1.^10.*l2.^3.*l3.^2 - 1.37266339176104e+16*l1.^10.*l2.^3.*l3 - 9472370859047.87*l1.^10.*l2.^3 + 1.63553252203771e+21*l1.^10.*l2.^2.*l3.^6 + 7.65334877146403e+18*l1.^10.*l2.^2.*l3.^5 + 1.64597557889265e+19*l1.^10.*l2.^2.*l3.^4 - 3.54431108444496e+16*l1.^10.*l2.^2.*l3.^3 - 3.36051930916874e+16*l1.^10.*l2.^2.*l3.^2 + 10866061334568.8*l1.^10.*l2.^2.*l3 + 30865253407485.6*l1.^10.*l2.^2 + 1.7726752237352e+20*l1.^10.*l2.*l3.^7 + 3.47048079511042e+19*l1.^10.*l2.*l3.^6 + 5.63600238197717e+18*l1.^10.*l2.*l3.^5 - 1.21853228048332e+17*l1.^10.*l2.*l3.^4 - 2.80574721550909e+16*l1.^10.*l2.*l3.^3 + 76841951133374.2*l1.^10.*l2.*l3.^2 + 92886084223492.2*l1.^10.*l2.*l3 + 114824445692.402*l1.^10.*l2 + 1.01940491607017e+21*l1.^10.*l3.^8 - 2.86742716774529e+18*l1.^10.*l3.^7 + 1.79155354191358e+19*l1.^10.*l3.^6 + 4.49818900163067e+16*l1.^10.*l3.^5 - 5.88628791294064e+16*l1.^10.*l3.^4 - 158786226719502.0*l1.^10.*l3.^3 + 142956995885481.0*l1.^10.*l3.^2 + 210825211763.1*l1.^10.*l3 - 35900414104.5541*l1.^10 - 1.22674308389357e+19*l1.^9.*l2.^9 - 3.64160135387586e+19*l1.^9.*l2.^8.*l3 - 3.77338287362164e+17*l1.^9.*l2.^8 - 3.43318996033344e+20*l1.^9.*l2.^7.*l3.^2 - 6.18869172894666e+18*l1.^9.*l2.^7.*l3 + 1.29671752315144e+17*l1.^9.*l2.^7 - 4.97860560973349e+20*l1.^9.*l2.^6.*l3.^3 - 7.72555017496841e+18*l1.^9.*l2.^6.*l3.^2 - 1.9082654881531e+18*l1.^9.*l2.^6.*l3 + 7.34111926461555e+15*l1.^9.*l2.^6 - 7.33167197716163e+20*l1.^9.*l2.^5.*l3.^4 - 3.02283209338328e+19*l1.^9.*l2.^5.*l3.^3 - 5.87441943771181e+18*l1.^9.*l2.^5.*l3.^2 + 5.72326750311429e+16*l1.^9.*l2.^5.*l3 + 588156984881232.0*l1.^9.*l2.^5 - 1.49851853102291e+21*l1.^9.*l2.^4.*l3.^5 - 1.87723649111382e+19*l1.^9.*l2.^4.*l3.^4 - 1.15149182988773e+19*l1.^9.*l2.^4.*l3.^3 + 7.2476018217542e+16*l1.^9.*l2.^4.*l3.^2 + 1.67972816306122e+16*l1.^9.*l2.^4.*l3 + 1889749797316.19*l1.^9.*l2.^4 - 3.70004423189357e+20*l1.^9.*l2.^3.*l3.^6 - 6.32484294698348e+19*l1.^9.*l2.^3.*l3.^5 - 9.72086906851566e+18*l1.^9.*l2.^3.*l3.^4 + 2.15090118581092e+17*l1.^9.*l2.^3.*l3.^3 + 3.43795560850134e+16*l1.^9.*l2.^3.*l3.^2 - 197715072544218.0*l1.^9.*l2.^3.*l3 - 40625021191171.1*l1.^9.*l2.^3 - 2.18053763151336e+21*l1.^9.*l2.^2.*l3.^7 - 2.76772046766781e+18*l1.^9.*l2.^2.*l3.^6 - 2.70093706591686e+19*l1.^9.*l2.^2.*l3.^5 + 5.82146081811869e+16*l1.^9.*l2.^2.*l3.^4 + 7.73949738647062e+16*l1.^9.*l2.^2.*l3.^3 - 109251160157348.0*l1.^9.*l2.^2.*l3.^2 - 189131685366031.0*l1.^9.*l2.^2.*l3 - 353132229703.191*l1.^9.*l2.^2 + 2.80882413602424e+20*l1.^9.*l2.*l3.^8 - 5.12011095708187e+19*l1.^9.*l2.*l3.^7 - 8.30609045744765e+16*l1.^9.*l2.*l3.^6 + 2.76437561098134e+17*l1.^9.*l2.*l3.^5 + 3.06729708299722e+16*l1.^9.*l2.*l3.^4 - 463933575241152.0*l1.^9.*l2.*l3.^3 - 274543278742612.0*l1.^9.*l2.*l3.^2 + 226637102645.332*l1.^9.*l2.*l3 + 100600003375.191*l1.^9.*l2 - 1.16583148796443e+21*l1.^9.*l3.^9 - 2.63277260635606e+18*l1.^9.*l3.^8 - 2.56054027372355e+19*l1.^9.*l3.^7 - 2.48288810796869e+16*l1.^9.*l3.^6 + 1.11856044022742e+17*l1.^9.*l3.^5 + 198659947442877.0*l1.^9.*l3.^4 - 386212772258539.0*l1.^9.*l3.^3 - 541745499512.68*l1.^9.*l3.^2 + 264176808169.723*l1.^9.*l3 + 865900574.905909*l1.^9 + 1.44307302506765e+19*l1.^8.*l2.^10 - 3.64160135387586e+19*l1.^8.*l2.^9.*l3 - 3.77338287362164e+17*l1.^8.*l2.^9 + 3.48476239140799e+20*l1.^8.*l2.^8.*l3.^2 + 6.21963518759139e+18*l1.^8.*l2.^8.*l3 - 2.3891532244108e+17*l1.^8.*l2.^8 + 1.10932740031378e+20*l1.^8.*l2.^7.*l3.^3 + 4.61057533806526e+18*l1.^8.*l2.^7.*l3.^2 + 7.20897447977008e+17*l1.^8.*l2.^7.*l3 - 3.50690410730042e+15*l1.^8.*l2.^7 + 9.64561382337457e+20*l1.^8.*l2.^6.*l3.^4 + 2.31044491214009e+19*l1.^8.*l2.^6.*l3.^3 + 5.83820900425927e+18*l1.^8.*l2.^6.*l3.^2 - 6.77533873530441e+16*l1.^8.*l2.^6.*l3 + 1.06539466904603e+15*l1.^8.*l2.^6 + 2.17677974975506e+20*l1.^8.*l2.^5.*l3.^5 + 4.15880084185215e+19*l1.^8.*l2.^5.*l3.^4 + 5.71304669746353e+18*l1.^8.*l2.^5.*l3.^3 - 1.05207123219013e+17*l1.^8.*l2.^5.*l3.^2 - 1.16045216369227e+16*l1.^8.*l2.^5.*l3 + 12188886192690.1*l1.^8.*l2.^5 + 2.08033842206562e+21*l1.^8.*l2.^4.*l3.^6 + 1.25630442097617e+19*l1.^8.*l2.^4.*l3.^5 + 1.84960572163466e+19*l1.^8.*l2.^4.*l3.^4 - 1.62720350578739e+17*l1.^8.*l2.^4.*l3.^3 - 4.54048183711123e+16*l1.^8.*l2.^4.*l3.^2 + 351493462300831.0*l1.^8.*l2.^4.*l3 + 38056227104539.1*l1.^8.*l2.^4 - 6.26471918973075e+20*l1.^8.*l2.^3.*l3.^7 + 9.09015670120115e+19*l1.^8.*l2.^3.*l3.^6 + 2.14803046294038e+18*l1.^8.*l2.^3.*l3.^5 - 4.1451606548291e+17*l1.^8.*l2.^3.*l3.^4 - 4.1243802343065e+16*l1.^8.*l2.^3.*l3.^3 + 609798637721503.0*l1.^8.*l2.^3.*l3.^2 + 216909016183870.0*l1.^8.*l2.^3.*l3 + 526310082222.881*l1.^8.*l2.^3 + 2.68014122311302e+21*l1.^8.*l2.^2.*l3.^8 - 9.46869834528839e+18*l1.^8.*l2.^2.*l3.^7 + 3.94002958440667e+19*l1.^8.*l2.^2.*l3.^6 - 3.57704218944643e+16*l1.^8.*l2.^2.*l3.^5 - 1.42623560281992e+17*l1.^8.*l2.^2.*l3.^4 + 349367493778849.0*l1.^8.*l2.^2.*l3.^3 + 524396343354119.0*l1.^8.*l2.^2.*l3.^2 - 399814955165.021*l1.^8.*l2.^2.*l3 - 200376207393.209*l1.^8.*l2.^2 - 8.3483893965458e+20*l1.^8.*l2.*l3.^9 + 5.62010567634968e+19*l1.^8.*l2.*l3.^8 - 1.15569684761763e+19*l1.^8.*l2.*l3.^7 - 4.26860367940607e+17*l1.^8.*l2.*l3.^6 + 3.41790744940337e+15*l1.^8.*l2.*l3.^5 + 1.18487312291732e+15*l1.^8.*l2.*l3.^4 + 470410598032541.0*l1.^8.*l2.*l3.^3 - 1629754181647.25*l1.^8.*l2.*l3.^2 - 514354254133.687*l1.^8.*l2.*l3 - 1455995781.50845*l1.^8.*l2 + 1.13901490129235e+21*l1.^8.*l3.^10 + 9.42400223835711e+18*l1.^8.*l3.^9 + 3.02279042622145e+19*l1.^8.*l3.^8 - 2.73538520369433e+16*l1.^8.*l3.^7 - 1.68455921243885e+17*l1.^8.*l3.^6 - 91983571384370.7*l1.^8.*l3.^5 + 752422684660545.0*l1.^8.*l3.^4 + 627957952180.093*l1.^8.*l3.^3 - 940983064062.558*l1.^8.*l3.^2 - 2809366309.69472*l1.^8.*l3 + 251120682.558624*l1.^8 - 1.84438294242929e+19*l1.^7.*l2.^11 + 9.36708155517112e+19*l1.^7.*l2.^10.*l3 + 9.35180083485273e+17*l1.^7.*l2.^10 - 3.43318996033344e+20*l1.^7.*l2.^9.*l3.^2 - 6.18869172894666e+18*l1.^7.*l2.^9.*l3 + 1.4707702187496e+17*l1.^7.*l2.^9 + 1.10932740031378e+20*l1.^7.*l2.^8.*l3.^3 + 4.61057533806526e+18*l1.^7.*l2.^8.*l3.^2 + 5.64250021938661e+17*l1.^7.*l2.^8.*l3 - 3.50690410730042e+15*l1.^7.*l2.^8 - 4.99291659203167e+20*l1.^7.*l2.^7.*l3.^4 - 3.33364194465927e+19*l1.^7.*l2.^7.*l3.^3 - 4.87553472732424e+18*l1.^7.*l2.^7.*l3.^2 + 7.95165457015318e+16*l1.^7.*l2.^7.*l3 - 1.9027669589023e+15*l1.^7.*l2.^7 - 8.36203331663008e+20*l1.^7.*l2.^6.*l3.^5 + 6.06491789436772e+18*l1.^7.*l2.^6.*l3.^4 - 4.15739646893326e+18*l1.^7.*l2.^6.*l3.^3 + 1.57109304007058e+16*l1.^7.*l2.^6.*l3.^2 + 4.77778629942405e+15*l1.^7.*l2.^6.*l3 - 2763759078575.12*l1.^7.*l2.^6 + 4.57398976730297e+20*l1.^7.*l2.^5.*l3.^6 - 7.9710349468833e+19*l1.^7.*l2.^5.*l3.^5 - 4.62515333664467e+18*l1.^7.*l2.^5.*l3.^4 + 3.38907212929513e+17*l1.^7.*l2.^5.*l3.^3 + 3.00886794260032e+16*l1.^7.*l2.^5.*l3.^2 - 598814467024602.0*l1.^7.*l2.^5.*l3 - 21386026905643.3*l1.^7.*l2.^5 - 2.67811314828077e+21*l1.^7.*l2.^4.*l3.^7 + 2.74365333316635e+19*l1.^7.*l2.^4.*l3.^6 - 2.49273905346508e+19*l1.^7.*l2.^4.*l3.^5 + 3.08607561442437e+16*l1.^7.*l2.^4.*l3.^4 + 7.01939267962585e+16*l1.^7.*l2.^4.*l3.^3 - 332005417516006.0*l1.^7.*l2.^4.*l3.^2 - 200870297280760.0*l1.^7.*l2.^4.*l3 - 502592245899.531*l1.^7.*l2.^4 + 1.87228553773066e+21*l1.^7.*l2.^3.*l3.^8 - 1.11891546459356e+20*l1.^7.*l2.^3.*l3.^7 + 1.45148282744523e+19*l1.^7.*l2.^3.*l3.^6 + 6.84547681745042e+17*l1.^7.*l2.^3.*l3.^5 + 9.38110101475783e+15*l1.^7.*l2.^3.*l3.^4 - 1.74754612506825e+15*l1.^7.*l2.^3.*l3.^3 - 459947187413500.0*l1.^7.*l2.^3.*l3.^2 + 1784484828137.67*l1.^7.*l2.^3.*l3 + 270246110575.56*l1.^7.*l2.^3 - 2.98968762329269e+21*l1.^7.*l2.^2.*l3.^9 + 2.26196682693e+19*l1.^7.*l2.^2.*l3.^8 - 5.13238297940285e+19*l1.^7.*l2.^2.*l3.^7 - 7.79935473463612e+16*l1.^7.*l2.^2.*l3.^6 + 2.22265258562414e+17*l1.^7.*l2.^2.*l3.^5 - 227832959938944.0*l1.^7.*l2.^2.*l3.^4 - 936349868356883.0*l1.^7.*l2.^2.*l3.^3 + 941936928270.134*l1.^7.*l2.^2.*l3.^2 + 827078392504.098*l1.^7.*l2.^2.*l3 + 5221059762.766*l1.^7.*l2.^2 + 1.2356537029059e+21*l1.^7.*l2.*l3.^10 - 4.52324602144568e+19*l1.^7.*l2.*l3.^9 + 2.47844398824155e+19*l1.^7.*l2.*l3.^8 + 4.75816749278521e+17*l1.^7.*l2.*l3.^7 - 8.08388415722778e+16*l1.^7.*l2.*l3.^6 - 1.97389090704181e+15*l1.^7.*l2.*l3.^5 - 464205723944267.0*l1.^7.*l2.*l3.^4 + 5064322765423.61*l1.^7.*l2.*l3.^3 + 1204511268611.23*l1.^7.*l2.*l3.^2 - 7748206625.82473*l1.^7.*l2.*l3 - 538955484.704588*l1.^7.*l2 - 9.49445811342512e+20*l1.^7.*l3.^11 - 1.44196517284457e+19*l1.^7.*l3.^10 - 2.95546141808913e+19*l1.^7.*l3.^9 + 9.24219528163974e+16*l1.^7.*l3.^8 + 2.03888245281702e+17*l1.^7.*l3.^7 - 190274182717285.0*l1.^7.*l3.^6 - 1.11073848947718e+15*l1.^7.*l3.^5 + 222495893128.556*l1.^7.*l3.^4 + 2119956909954.35*l1.^7.*l3.^3 + 5477622896.07145*l1.^7.*l3.^2 - 1570353858.04327*l1.^7.*l3 - 15051582.3418304*l1.^7 + 2.0769731372755e+19*l1.^6.*l2.^12 - 1.26547873036733e+20*l1.^6.*l2.^11.*l3 - 9.62685380058369e+17*l1.^6.*l2.^11 + 3.84531098732033e+20*l1.^6.*l2.^10.*l3.^2 + 4.90969543829768e+18*l1.^6.*l2.^10.*l3 + 7.41559996812465e+16*l1.^6.*l2.^10 - 4.97860560973349e+20*l1.^6.*l2.^9.*l3.^3 - 7.72555017496841e+18*l1.^6.*l2.^9.*l3.^2 - 1.72215835483975e+18*l1.^6.*l2.^9.*l3 + 7.34111926461555e+15*l1.^6.*l2.^9 + 9.64561382337457e+20*l1.^6.*l2.^8.*l3.^4 + 2.31044491214009e+19*l1.^6.*l2.^8.*l3.^3 + 5.54899289548341e+18*l1.^6.*l2.^8.*l3.^2 - 6.77533873530441e+16*l1.^6.*l2.^8.*l3 + 821113508760744.0*l1.^6.*l2.^8 - 8.36203331663008e+20*l1.^6.*l2.^7.*l3.^5 + 6.06491789436772e+18*l1.^6.*l2.^7.*l3.^4 - 3.38615351219764e+18*l1.^6.*l2.^7.*l3.^3 + 1.57109304007059e+16*l1.^6.*l2.^7.*l3.^2 + 6.73203558170636e+15*l1.^6.*l2.^7.*l3 - 2763759078575.17*l1.^6.*l2.^7 + 1.95114110721369e+21*l1.^6.*l2.^6.*l3.^6 - 4056.99047619048*l1.^6.*l2.^6.*l3.^5 + 1.29000863009655e+19*l1.^6.*l2.^6.*l3.^4 - 1.80675699608118e+17*l1.^6.*l2.^6.*l3.^3 - 3.82831315485985e+16*l1.^6.*l2.^6.*l3.^2 + 600373510607387.0*l1.^6.*l2.^6.*l3 + 9424132774611.62*l1.^6.*l2.^6 - 2.05086543889153e+21*l1.^6.*l2.^5.*l3.^7 + 8.28872112230256e+19*l1.^6.*l2.^5.*l3.^6 - 9.42908401067517e+18*l1.^6.*l2.^5.*l3.^5 - 3.18146340614294e+17*l1.^6.*l2.^5.*l3.^4 + 7.29776269004023e+15*l1.^6.*l2.^5.*l3.^3 + 422288214082793.0*l1.^6.*l2.^5.*l3.^2 + 91990323555287.8*l1.^6.*l2.^5.*l3 - 175060220481.859*l1.^6.*l2.^5 + 3.57748697773999e+21*l1.^6.*l2.^4.*l3.^8 - 5.71835115754671e+19*l1.^6.*l2.^4.*l3.^7 + 3.70399010213265e+19*l1.^6.*l2.^4.*l3.^6 + 9.03378498040588e+16*l1.^6.*l2.^4.*l3.^5 - 1.30552530563379e+17*l1.^6.*l2.^4.*l3.^4 + 845426815574378.0*l1.^6.*l2.^4.*l3.^3 + 466865946791125.0*l1.^6.*l2.^4.*l3.^2 - 1861285440994.22*l1.^6.*l2.^4.*l3 - 318363297536.372*l1.^6.*l2.^4 - 2.96839922900912e+21*l1.^6.*l2.^3.*l3.^9 + 1.0505304209887e+20*l1.^6.*l2.^3.*l3.^8 - 3.64323662415298e+19*l1.^6.*l2.^3.*l3.^7 - 7.8947425263547e+17*l1.^6.*l2.^3.*l3.^6 + 9.45423173706997e+16*l1.^6.*l2.^3.*l3.^5 + 2.42537575549314e+15*l1.^6.*l2.^3.*l3.^4 + 441922732142840.0*l1.^6.*l2.^3.*l3.^3 - 3959749245221.93*l1.^6.*l2.^3.*l3.^2 - 649313723013.393*l1.^6.*l2.^3.*l3 - 4656620869.494*l1.^6.*l2.^3 + 2.96793022843935e+21*l1.^6.*l2.^2.*l3.^10 - 2.30563148523979e+19*l1.^6.*l2.^2.*l3.^9 + 5.87121673362836e+19*l1.^6.*l2.^2.*l3.^8 + 1.70014711121924e+17*l1.^6.*l2.^2.*l3.^7 - 3.06678579039327e+17*l1.^6.*l2.^2.*l3.^6 - 88156828044805.5*l1.^6.*l2.^2.*l3.^5 + 1.30464784744811e+15*l1.^6.*l2.^2.*l3.^4 - 1513423841585.11*l1.^6.*l2.^2.*l3.^3 - 1906949380045.76*l1.^6.*l2.^2.*l3.^2 + 4156322759.54836*l1.^6.*l2.^2.*l3 + 1448456607.42789*l1.^6.*l2.^2 - 1.31168351457004e+21*l1.^6.*l2.*l3.^11 + 2.35376575424271e+19*l1.^6.*l2.*l3.^10 - 3.25959667226206e+19*l1.^6.*l2.*l3.^9 - 3.672429981165e+17*l1.^6.*l2.*l3.^8 + 1.70537790298301e+17*l1.^6.*l2.*l3.^7 + 2.18032245052615e+15*l1.^6.*l2.*l3.^6 + 116263321831980.0*l1.^6.*l2.*l3.^5 - 8824164533063.04*l1.^6.*l2.*l3.^4 - 1477965131898.66*l1.^6.*l2.*l3.^3 + 27567708673.6711*l1.^6.*l2.*l3.^2 + 905053987.919753*l1.^6.*l2.*l3 + 4391821.47651142*l1.^6.*l2 + 6.71734294136156e+20*l1.^6.*l3.^12 + 1.54511003499368e+19*l1.^6.*l3.^11 + 2.38557430901485e+19*l1.^6.*l3.^10 - 1.38125263106206e+17*l1.^6.*l3.^9 - 1.99267433879131e+17*l1.^6.*l3.^8 + 542405435574710.0*l1.^6.*l3.^7 + 1.27398958574552e+15*l1.^6.*l3.^6 - 2245288505277.01*l1.^6.*l3.^5 - 3342143730334.15*l1.^6.*l3.^4 - 1026252533.22183*l1.^6.*l3.^3 + 5163456156.24532*l1.^6.*l3.^2 + 35923394.1161248*l1.^6.*l3 - 4129655.42028746*l1.^6 - 1.8585176087238e+19*l1.^5.*l2.^13 + 1.21240761326154e+20*l1.^5.*l2.^12.*l3 + 6.3777906428867e+17*l1.^5.*l2.^12 - 3.59353173407429e+20*l1.^5.*l2.^11.*l3.^2 - 3.46566736821013e+18*l1.^5.*l2.^11.*l3 - 2.42508351877539e+17*l1.^5.*l2.^11 + 5.71831001714586e+20*l1.^5.*l2.^10.*l3.^3 + 1.04692035081348e+19*l1.^5.*l2.^10.*l3.^2 + 2.17433852708709e+18*l1.^5.*l2.^10.*l3 - 7.24760182175421e+15*l1.^5.*l2.^10 - 7.33167197716163e+20*l1.^5.*l2.^9.*l3.^4 - 3.02283209338328e+19*l1.^5.*l2.^9.*l3.^3 - 5.51444198747906e+18*l1.^5.*l2.^9.*l3.^2 + 5.72326750311429e+16*l1.^5.*l2.^9.*l3 + 974850745513828.0*l1.^5.*l2.^9 + 2.17677974975506e+20*l1.^5.*l2.^8.*l3.^5 + 4.15880084185215e+19*l1.^5.*l2.^8.*l3.^4 + 5.21154656431699e+18*l1.^5.*l2.^8.*l3.^3 - 1.05207123219013e+17*l1.^5.*l2.^8.*l3.^2 - 1.36190785209043e+16*l1.^5.*l2.^8.*l3 + 12188886192690.1*l1.^5.*l2.^8 + 4.57398976730297e+20*l1.^5.*l2.^7.*l3.^6 - 7.9710349468833e+19*l1.^5.*l2.^7.*l3.^5 - 4.779017505455e+18*l1.^5.*l2.^7.*l3.^4 + 3.38907212929513e+17*l1.^5.*l2.^7.*l3.^3 + 3.22841591150828e+16*l1.^5.*l2.^7.*l3.^2 - 598814467024602.0*l1.^5.*l2.^7.*l3 - 5852383258024.16*l1.^5.*l2.^7 - 2.05086543889153e+21*l1.^5.*l2.^6.*l3.^7 + 8.28872112230256e+19*l1.^5.*l2.^6.*l3.^6 - 9.2136741743407e+18*l1.^5.*l2.^6.*l3.^5 - 3.18146340614294e+17*l1.^5.*l2.^6.*l3.^4 + 2.17497674885439e+15*l1.^5.*l2.^6.*l3.^3 + 422288214082793.0*l1.^5.*l2.^6.*l3.^2 - 8364879810463.27*l1.^5.*l2.^6.*l3 + 1081985104655.62*l1.^5.*l2.^6 + 2.97232813729099e+21*l1.^5.*l2.^5.*l3.^8 - 1.24764025255565e+20*l1.^5.*l2.^5.*l3.^7 + 1.94806744730535e+19*l1.^5.*l2.^5.*l3.^6 + 7.54124659233882e+17*l1.^5.*l2.^5.*l3.^5 - 1.78747513276862e+16*l1.^5.*l2.^5.*l3.^4 - 2.19938530160657e+15*l1.^5.*l2.^5.*l3.^3 - 245342169425692.0*l1.^5.*l2.^5.*l3.^2 + 1660248542634.41*l1.^5.*l2.^5.*l3 + 292007735674.317*l1.^5.*l2.^5 - 3.97114301738214e+21*l1.^5.*l2.^4.*l3.^9 + 8.79413094683321e+19*l1.^5.*l2.^4.*l3.^8 - 4.61645075641977e+19*l1.^5.*l2.^4.*l3.^7 - 5.26316168423647e+17*l1.^5.*l2.^4.*l3.^6 + 1.75813841554379e+17*l1.^5.*l2.^4.*l3.^5 + 957040162978372.0*l1.^5.*l2.^4.*l3.^4 - 519399528669863.0*l1.^5.*l2.^4.*l3.^3 + 2442937141304.91*l1.^5.*l2.^4.*l3.^2 + 436757063392.546*l1.^5.*l2.^4.*l3 - 2007606518.11519*l1.^5.*l2.^4 + 3.31603561625243e+21*l1.^5.*l2.^3.*l3.^10 - 8.01916921588621e+19*l1.^5.*l2.^3.*l3.^9 + 5.10798540518479e+19*l1.^5.*l2.^3.*l3.^8 + 8.07990706322017e+17*l1.^5.*l2.^3.*l3.^7 - 2.15653720925927e+17*l1.^5.*l2.^3.*l3.^6 - 3.51195826707754e+15*l1.^5.*l2.^3.*l3.^5 - 133169193111702.0*l1.^5.*l2.^3.*l3.^4 + 6189225859616.72*l1.^5.*l2.^3.*l3.^3 + 1040047582392.16*l1.^5.*l2.^3.*l3.^2 - 10313837958.8793*l1.^5.*l2.^3.*l3 - 2478969400.63314*l1.^5.*l2.^3 - 2.51628679378898e+21*l1.^5.*l2.^2.*l3.^11 + 1.14800231571961e+19*l1.^5.*l2.^2.*l3.^10 - 5.63030397506682e+19*l1.^5.*l2.^2.*l3.^9 - 1.8221873741533e+17*l1.^5.*l2.^2.*l3.^8 + 3.55970349187962e+17*l1.^5.*l2.^2.*l3.^7 + 713971095298563.0*l1.^5.*l2.^2.*l3.^6 - 1.48381461562971e+15*l1.^5.*l2.^2.*l3.^5 + 1395964080745.67*l1.^5.*l2.^2.*l3.^4 + 2756575926723.46*l1.^5.*l2.^2.*l3.^3 - 7773862939.15528*l1.^5.*l2.^2.*l3.^2 - 1517767121.69319*l1.^5.*l2.^2.*l3 - 39590351.8537946*l1.^5.*l2.^2 + 1.07345796638641e+21*l1.^5.*l2.*l3.^12 - 3.17686175419262e+18*l1.^5.*l2.*l3.^11 + 3.09309076748349e+19*l1.^5.*l2.*l3.^10 + 1.69079536693311e+17*l1.^5.*l2.*l3.^9 - 2.18245683813449e+17*l1.^5.*l2.*l3.^8 - 1.6197045512798e+15*l1.^5.*l2.*l3.^7 + 388769868555016.0*l1.^5.*l2.*l3.^6 + 9385486659382.29*l1.^5.*l2.*l3.^5 + 752026916473.451*l1.^5.*l2.*l3.^4 - 51979690807.6851*l1.^5.*l2.*l3.^3 - 695106910.092239*l1.^5.*l2.*l3.^2 + 150686379.592149*l1.^5.*l2.*l3 + 5223539.06825786*l1.^5.*l2 - 3.99090686595394e+20*l1.^5.*l3.^13 - 1.26472791805168e+19*l1.^5.*l3.^12 - 1.57381731698543e+19*l1.^5.*l3.^11 + 1.41819202099229e+17*l1.^5.*l3.^10 + 1.56532485750635e+17*l1.^5.*l3.^9 - 758687299877560.0*l1.^5.*l3.^8 - 1.14325611305333e+15*l1.^5.*l3.^7 + 4487188748222.12*l1.^5.*l3.^6 + 3838333570461.99*l1.^5.*l3.^5 - 16375141983.2206*l1.^5.*l3.^4 - 10177219102.7171*l1.^5.*l3.^3 - 36904092.1157342*l1.^5.*l3.^2 + 20686659.92788*l1.^5.*l3 + 169839.971161965*l1.^5 + 1.32045993751291e+19*l1.^4.*l2.^14 - 9.34627626673762e+19*l1.^4.*l2.^13.*l3 - 2.28637777763863e+17*l1.^4.*l2.^13 + 3.09219480922851e+20*l1.^4.*l2.^12.*l3.^2 + 1.44402807008755e+18*l1.^4.*l2.^12.*l3 + 2.65509857803053e+17*l1.^4.*l2.^12 - 6.5127311576983e+20*l1.^4.*l2.^11.*l3.^3 - 5.99271649086335e+18*l1.^4.*l2.^11.*l3.^2 - 1.95082062315268e+18*l1.^4.*l2.^11.*l3 + 3.41338666443908e+15*l1.^4.*l2.^11 + 1.10332795442878e+21*l1.^4.*l2.^10.*l3.^4 + 1.62693829229864e+19*l1.^4.*l2.^10.*l3.^3 + 5.19562965702239e+18*l1.^4.*l2.^10.*l3.^2 - 2.94112357798929e+16*l1.^4.*l2.^10.*l3 - 2.19096518531071e+15*l1.^4.*l2.^10 - 1.49851853102291e+21*l1.^4.*l2.^9.*l3.^5 - 1.87723649111382e+19*l1.^4.*l2.^9.*l3.^4 - 7.96697105934834e+18*l1.^4.*l2.^9.*l3.^3 + 7.2476018217542e+16*l1.^4.*l2.^9.*l3.^2 + 1.70479964575201e+16*l1.^4.*l2.^9.*l3 + 1889749797316.23*l1.^4.*l2.^9 + 2.08033842206562e+21*l1.^4.*l2.^8.*l3.^6 + 1.25630442097617e+19*l1.^4.*l2.^8.*l3.^5 + 1.11400510938396e+19*l1.^4.*l2.^8.*l3.^4 - 1.62720350578739e+17*l1.^4.*l2.^8.*l3.^3 - 4.14966428822435e+16*l1.^4.*l2.^8.*l3.^2 + 351493462300831.0*l1.^4.*l2.^8.*l3 + 10582128194115.0*l1.^4.*l2.^8 - 2.67811314828077e+21*l1.^4.*l2.^7.*l3.^7 + 2.74365333316635e+19*l1.^4.*l2.^7.*l3.^6 - 1.30039165698294e+19*l1.^4.*l2.^7.*l3.^5 + 3.08607561442436e+16*l1.^4.*l2.^7.*l3.^4 + 5.61129926774757e+16*l1.^4.*l2.^7.*l3.^3 - 332005417516007.0*l1.^4.*l2.^7.*l3.^2 - 66493550767756.5*l1.^4.*l2.^7.*l3 - 1664766225743.62*l1.^4.*l2.^7 + 3.57748697773999e+21*l1.^4.*l2.^6.*l3.^8 - 5.71835115754671e+19*l1.^4.*l2.^6.*l3.^7 + 2.29496731987561e+19*l1.^4.*l2.^6.*l3.^6 + 9.03378498040586e+16*l1.^4.*l2.^6.*l3.^5 - 9.95074134290266e+16*l1.^4.*l2.^6.*l3.^4 + 845426815574379.0*l1.^4.*l2.^6.*l3.^3 + 254550693306996.0*l1.^4.*l2.^6.*l3.^2 + 273319828107.16*l1.^4.*l2.^6.*l3 - 204613400968.647*l1.^4.*l2.^6 - 3.97114301738214e+21*l1.^4.*l2.^5.*l3.^9 + 8.7941309468332e+19*l1.^4.*l2.^5.*l3.^8 - 3.4087169430566e+19*l1.^4.*l2.^5.*l3.^7 - 5.26316168423647e+17*l1.^4.*l2.^5.*l3.^6 + 1.38559700993156e+17*l1.^4.*l2.^5.*l3.^5 + 957040162978372.0*l1.^4.*l2.^5.*l3.^4 - 110897150401482.0*l1.^4.*l2.^5.*l3.^3 - 4814720773639.79*l1.^4.*l2.^5.*l3.^2 - 134137353025.932*l1.^4.*l2.^5.*l3 + 14733137930.0657*l1.^4.*l2.^5 + 3.93038933629301e+21*l1.^4.*l2.^4.*l3.^10 - 6.11305216337064e+19*l1.^4.*l2.^4.*l3.^9 + 5.0554747253888e+19*l1.^4.*l2.^4.*l3.^8 + 4.76938958592858e+17*l1.^4.*l2.^4.*l3.^7 - 2.62940149728079e+17*l1.^4.*l2.^4.*l3.^6 - 865977844620196.0*l1.^4.*l2.^4.*l3.^5 + 791221463811845.0*l1.^4.*l2.^4.*l3.^4 + 2010368983598.13*l1.^4.*l2.^4.*l3.^3 - 778955789558.545*l1.^4.*l2.^4.*l3.^2 - 6721954092.6029*l1.^4.*l2.^4.*l3 + 2970434790.20821*l1.^4.*l2.^4 - 2.91503601505726e+21*l1.^4.*l2.^3.*l3.^11 + 3.54508891206494e+19*l1.^4.*l2.^3.*l3.^10 - 5.12169062347554e+19*l1.^4.*l2.^3.*l3.^9 - 5.0008452570104e+17*l1.^4.*l2.^3.*l3.^8 + 3.06831262242348e+17*l1.^4.*l2.^3.*l3.^7 + 2.66844482317295e+15*l1.^4.*l2.^3.*l3.^6 - 512758669517507.0*l1.^4.*l2.^3.*l3.^5 - 9713018684799.96*l1.^4.*l2.^3.*l3.^4 - 603411428422.506*l1.^4.*l2.^3.*l3.^3 + 41127070268.8643*l1.^4.*l2.^3.*l3.^2 + 307729782.994391*l1.^4.*l2.^3.*l3 - 30806708.9007718*l1.^4.*l2.^3 + 1.78001966686822e+21*l1.^4.*l2.^2.*l3.^12 + 5.05409824530644e+18*l1.^4.*l2.^2.*l3.^11 + 4.35491669883574e+19*l1.^4.*l2.^2.*l3.^10 + 4.88628638950526e+16*l1.^4.*l2.^2.*l3.^9 - 3.37950078336592e+17*l1.^4.*l2.^2.*l3.^8 - 488264103881596.0*l1.^4.*l2.^2.*l3.^7 + 1.49360681427581e+15*l1.^4.*l2.^2.*l3.^6 + 643769843062.338*l1.^4.*l2.^2.*l3.^5 - 3013344303574.14*l1.^4.*l2.^2.*l3.^4 - 12109779892.0174*l1.^4.*l2.^2.*l3.^3 + 3379726077.54879*l1.^4.*l2.^2.*l3.^2 + 14070884.342221*l1.^4.*l2.^2.*l3 - 14413508.5309955*l1.^4.*l2.^2 - 6.90685619791023e+20*l1.^4.*l2.*l3.^13 - 7.73758374221914e+18*l1.^4.*l2.*l3.^12 - 2.19729146551759e+19*l1.^4.*l2.*l3.^11 + 7.62167159319957e+15*l1.^4.*l2.*l3.^10 + 1.96762314621126e+17*l1.^4.*l2.*l3.^9 + 559838377454955.0*l1.^4.*l2.*l3.^8 - 710376169343323.0*l1.^4.*l2.*l3.^7 - 6235532111486.12*l1.^4.*l2.*l3.^6 + 540687656055.598*l1.^4.*l2.*l3.^5 + 65532638324.5457*l1.^4.*l2.*l3.^4 - 1892699964.99567*l1.^4.*l2.*l3.^3 - 359127343.552596*l1.^4.*l2.*l3.^2 + 6886615.50979856*l1.^4.*l2.*l3 + 626594.06836453*l1.^4.*l2 + 1.95717816717947e+20*l1.^4.*l3.^14 + 8.0865571924903e+18*l1.^4.*l3.^13 + 8.32814778525283e+18*l1.^4.*l3.^12 - 1.08901062212036e+17*l1.^4.*l3.^11 - 9.74550763688724e+16*l1.^4.*l3.^10 + 731096952836744.0*l1.^4.*l3.^9 + 798923586058843.0*l1.^4.*l3.^8 - 5290206920884.06*l1.^4.*l3.^7 - 3249528772694.82*l1.^4.*l3.^6 + 35277430829.5*l1.^4.*l3.^5 + 13072999933.0996*l1.^4.*l3.^4 - 73019361.9274348*l1.^4.*l3.^3 - 55808537.3997485*l1.^4.*l3.^2 - 194573.947544775*l1.^4.*l3 + 91719.6490082086*l1.^4 - 7.18076310961731e+18*l1.^3.*l2.^15 + 5.4915970224221e+19*l1.^3.*l2.^14.*l3 + 3.43816207163708e+15*l1.^3.*l2.^14 - 1.98448951614853e+20*l1.^3.*l2.^13.*l3.^2 - 2.88805614017511e+17*l1.^3.*l2.^13.*l3 - 1.78691956675214e+17*l1.^3.*l2.^13 + 4.54712508554322e+20*l1.^3.*l2.^12.*l3.^3 + 3.17686175419262e+18*l1.^3.*l2.^12.*l3.^2 + 1.25783875747019e+18*l1.^3.*l2.^12.*l3 - 467587214306723.0*l1.^3.*l2.^12 - 7.35327068761165e+20*l1.^3.*l2.^11.*l3.^4 - 1.34775953208172e+19*l1.^3.*l2.^11.*l3.^3 - 3.60424362365678e+18*l1.^3.*l2.^11.*l3.^2 + 1.23443024576975e+16*l1.^3.*l2.^11.*l3 + 2.01318189266048e+15*l1.^3.*l2.^11 + 7.79244446676222e+20*l1.^3.*l2.^10.*l3.^5 + 3.31404442085094e+19*l1.^3.*l2.^10.*l3.^4 + 6.14362124895914e+18*l1.^3.*l2.^10.*l3.^3 - 6.49946227886345e+16*l1.^3.*l2.^10.*l3.^2 - 1.34442892170828e+16*l1.^3.*l2.^10.*l3 - 9472370859047.9*l1.^3.*l2.^10 - 3.70004423189357e+20*l1.^3.*l2.^9.*l3.^6 - 6.32484294698348e+19*l1.^3.*l2.^9.*l3.^5 - 7.03359853208103e+18*l1.^3.*l2.^9.*l3.^4 + 2.15090118581092e+17*l1.^3.*l2.^9.*l3.^3 + 3.24664029285597e+16*l1.^3.*l2.^9.*l3.^2 - 197715072544218.0*l1.^3.*l2.^9.*l3 - 11258550775256.1*l1.^3.*l2.^9 - 6.26471918973076e+20*l1.^3.*l2.^8.*l3.^7 + 9.09015670120115e+19*l1.^3.*l2.^8.*l3.^6 + 3.19574839536361e+18*l1.^3.*l2.^8.*l3.^5 - 4.1451606548291e+17*l1.^3.*l2.^8.*l3.^4 - 4.0715243525529e+16*l1.^3.*l2.^8.*l3.^3 + 609798637721503.0*l1.^3.*l2.^8.*l3.^2 + 68673815526401.7*l1.^3.*l2.^8.*l3 + 1451305698833.48*l1.^3.*l2.^8 + 1.87228553773066e+21*l1.^3.*l2.^7.*l3.^8 - 1.11891546459356e+20*l1.^3.*l2.^7.*l3.^7 + 5.16888838800707e+18*l1.^3.*l2.^7.*l3.^6 + 6.84547681745042e+17*l1.^3.*l2.^7.*l3.^5 + 2.24049174984686e+16*l1.^3.*l2.^7.*l3.^4 - 1.74754612506825e+15*l1.^3.*l2.^7.*l3.^3 - 146828201427198.0*l1.^3.*l2.^7.*l3.^2 - 1536012257131.15*l1.^3.*l2.^7.*l3 + 74326973207.8083*l1.^3.*l2.^7 - 2.96839922900912e+21*l1.^3.*l2.^6.*l3.^9 + 1.0505304209887e+20*l1.^3.*l2.^6.*l3.^8 - 1.903486759933e+19*l1.^3.*l2.^6.*l3.^7 - 7.89474252635471e+17*l1.^3.*l2.^6.*l3.^6 + 5.14728805860228e+16*l1.^3.*l2.^6.*l3.^5 + 2.42537575549314e+15*l1.^3.*l2.^6.*l3.^4 + 13900101575060.4*l1.^3.*l2.^6.*l3.^3 + 4815850194417.08*l1.^3.*l2.^6.*l3.^2 + 155515689302.095*l1.^3.*l2.^6.*l3 - 21397365317.6749*l1.^3.*l2.^6 + 3.31603561625242e+21*l1.^3.*l2.^5.*l3.^10 - 8.01916921588621e+19*l1.^3.*l2.^5.*l3.^9 + 3.19930610033822e+19*l1.^3.*l2.^5.*l3.^8 + 8.07990706322017e+17*l1.^3.*l2.^5.*l3.^7 - 1.47748190433768e+17*l1.^3.*l2.^5.*l3.^6 - 3.51195826707755e+15*l1.^3.*l2.^5.*l3.^5 + 84392364680531.3*l1.^3.*l2.^5.*l3.^4 - 3392780015016.14*l1.^3.*l2.^5.*l3.^3 - 109501517730.736*l1.^3.*l2.^5.*l3.^2 + 26631253237.1062*l1.^3.*l2.^5.*l3 - 1949941564.82257*l1.^3.*l2.^5 - 2.91503601505726e+21*l1.^3.*l2.^4.*l3.^11 + 3.54508891206494e+19*l1.^3.*l2.^4.*l3.^10 - 4.06131323189411e+19*l1.^3.*l2.^4.*l3.^9 - 5.0008452570104e+17*l1.^3.*l2.^4.*l3.^8 + 2.58327311890806e+17*l1.^3.*l2.^4.*l3.^7 + 2.66844482317294e+15*l1.^3.*l2.^4.*l3.^6 - 711504077159884.0*l1.^3.*l2.^4.*l3.^5 + 4399093927592.51*l1.^3.*l2.^4.*l3.^4 + 1088010996806.71*l1.^3.*l2.^4.*l3.^3 - 59894663470.1584*l1.^3.*l2.^4.*l3.^2 - 1262905500.83431*l1.^3.*l2.^4.*l3 + 130368875.382851*l1.^3.*l2.^4 + 1.94079018235802e+21*l1.^3.*l2.^3.*l3.^12 - 4.62088982428017e+18*l1.^3.*l2.^3.*l3.^11 + 3.68841371285331e+19*l1.^3.*l2.^3.*l3.^10 + 2.28182560581681e+17*l1.^3.*l2.^3.*l3.^9 - 2.84324943441848e+17*l1.^3.*l2.^3.*l3.^8 - 2.09478765032512e+15*l1.^3.*l2.^3.*l3.^7 + 818940990444009.0*l1.^3.*l2.^3.*l3.^6 + 2489243393174.34*l1.^3.*l2.^3.*l3.^5 - 432089647101.024*l1.^3.*l2.^3.*l3.^4 + 31505952769.9098*l1.^3.*l2.^3.*l3.^3 + 1665034805.63263*l1.^3.*l2.^3.*l3.^2 - 73594989.0141621*l1.^3.*l2.^3.*l3 + 16267568.9381165*l1.^3.*l2.^3 - 1.00196623479682e+21*l1.^3.*l2.^2.*l3.^13 - 1.24667756717559e+19*l1.^3.*l2.^2.*l3.^12 - 2.58562635547151e+19*l1.^3.*l2.^2.*l3.^11 + 6.45270355743277e+16*l1.^3.*l2.^2.*l3.^10 + 2.45040894741202e+17*l1.^3.*l2.^2.*l3.^9 + 88936349836200.5*l1.^3.*l2.^2.*l3.^8 - 1.2028142883742e+15*l1.^3.*l2.^2.*l3.^7 + 1732531472381.74*l1.^3.*l2.^2.*l3.^6 + 2626672686844.63*l1.^3.*l2.^2.*l3.^5 - 32442408206.4747*l1.^3.*l2.^2.*l3.^4 - 5842731347.42278*l1.^3.*l2.^2.*l3.^3 + 147872202.723705*l1.^3.*l2.^2.*l3.^2 + 5366584.85925438*l1.^3.*l2.^2.*l3 - 522711.367556726*l1.^3.*l2.^2 + 3.44602576540088e+20*l1.^3.*l2.*l3.^14 + 8.66416842052532e+18*l1.^3.*l2.*l3.^13 + 1.15896250377631e+19*l1.^3.*l2.*l3.^12 - 7.93027915464201e+16*l1.^3.*l2.*l3.^11 - 1.27582554581526e+17*l1.^3.*l2.*l3.^10 + 195825322746902.0*l1.^3.*l2.*l3.^9 + 656780134443356.0*l1.^3.*l2.*l3.^8 + 1066173213773.39*l1.^3.*l2.*l3.^7 - 1265663295056.12*l1.^3.*l2.*l3.^6 - 26682565863.7672*l1.^3.*l2.*l3.^5 + 5304556768.6566*l1.^3.*l2.*l3.^4 + 294294677.969726*l1.^3.*l2.*l3.^3 - 26318165.4690627*l1.^3.*l2.*l3.^2 - 2153504.87706336*l1.^3.*l2.*l3 - 44251.8187558517*l1.^3.*l2 - 7.69390145230899e+19*l1.^3.*l3.^15 - 3.9917061651706e+18*l1.^3.*l3.^14 - 3.42176510076508e+18*l1.^3.*l3.^13 + 6.28904803242542e+16*l1.^3.*l3.^12 + 4.67292386574698e+16*l1.^3.*l3.^11 - 507185223727228.0*l1.^3.*l3.^10 - 424904564012460.0*l1.^3.*l3.^9 + 4250010384988.64*l1.^3.*l3.^8 + 1997170086069.78*l1.^3.*l3.^7 - 42538167502.0443*l1.^3.*l3.^6 - 11152416544.2412*l1.^3.*l3.^5 + 275171066.977344*l1.^3.*l3.^4 + 81516851.7998442*l1.^3.*l3.^3 - 750263.950278582*l1.^3.*l3.^2 - 329513.643826071*l1.^3.*l3 - 1033.02998236332*l1.^3 + 2.95872947164765e+18*l1.^2.*l2.^16 - 2.54077238493915e+19*l1.^2.*l2.^15.*l3 + 5.84487552178296e+16*l1.^2.*l2.^15 + 1.07152524912614e+20*l1.^2.*l2.^14.*l3.^2 - 3.09434586447333e+17*l1.^2.*l2.^14.*l3 + 7.84631518274452e+16*l1.^2.*l2.^14 - 3.00954374390641e+20*l1.^2.*l2.^13.*l3.^3 + 2.88805614017511e+17*l1.^2.*l2.^13.*l3.^2 - 5.68642751749348e+17*l1.^2.*l2.^13.*l3 - 935174428613445.0*l1.^2.*l2.^13 + 6.39393825022315e+20*l1.^2.*l2.^12.*l3.^4 + 1.39589380108463e+18*l1.^2.*l2.^12.*l3.^3 + 1.81690114559439e+18*l1.^2.*l2.^12.*l3.^2 + 2.24441862867226e+15*l1.^2.*l2.^12.*l3 - 1.14021791655302e+15*l1.^2.*l2.^12 - 1.09999146791926e+21*l1.^2.*l2.^11.*l3.^5 - 4.98189684180206e+18*l1.^2.*l2.^11.*l3.^4 - 3.65151541595834e+18*l1.^2.*l2.^11.*l3.^3 + 5.33049424309663e+15*l1.^2.*l2.^11.*l3.^2 + 7.7481079135556e+15*l1.^2.*l2.^11.*l3 + 15755788935124.6*l1.^2.*l2.^11 + 1.63553252203771e+21*l1.^2.*l2.^10.*l3.^6 + 7.65334877146404e+18*l1.^2.*l2.^10.*l3.^5 + 5.37487960575866e+18*l1.^2.*l2.^10.*l3.^4 - 3.54431108444495e+16*l1.^2.*l2.^10.*l3.^3 - 2.2411939494483e+16*l1.^2.*l2.^10.*l3.^2 + 10866061334568.6*l1.^2.*l2.^10.*l3 + 8182072343480.51*l1.^2.*l2.^10 - 2.18053763151336e+21*l1.^2.*l2.^9.*l3.^7 - 2.76772046766781e+18*l1.^2.*l2.^9.*l3.^6 - 6.45198061469378e+18*l1.^2.*l2.^9.*l3.^5 + 5.82146081811869e+16*l1.^2.*l2.^9.*l3.^4 + 4.19972816971451e+16*l1.^2.*l2.^9.*l3.^3 - 109251160157348.0*l1.^2.*l2.^9.*l3.^2 - 53192720029585.8*l1.^2.*l2.^9.*l3 - 772147338082.351*l1.^2.*l2.^9 + 2.68014122311302e+21*l1.^2.*l2.^8.*l3.^8 - 9.46869834528839e+18*l1.^2.*l2.^8.*l3.^7 + 8.04035181114279e+18*l1.^2.*l2.^8.*l3.^6 - 3.57704218944643e+16*l1.^2.*l2.^8.*l3.^5 - 6.33751720181316e+16*l1.^2.*l2.^8.*l3.^4 + 349367493778850.0*l1.^2.*l2.^8.*l3.^3 + 141190437770387.0*l1.^2.*l2.^8.*l3.^2 + 1450176278056.17*l1.^2.*l2.^8.*l3 + 12243801149.4422*l1.^2.*l2.^8 - 2.98968762329269e+21*l1.^2.*l2.^7.*l3.^9 + 2.26196682693e+19*l1.^2.*l2.^7.*l3.^8 - 1.14784909479247e+19*l1.^2.*l2.^7.*l3.^7 - 7.79935473463614e+16*l1.^2.*l2.^7.*l3.^6 + 8.76535474500108e+16*l1.^2.*l2.^7.*l3.^5 - 227832959938946.0*l1.^2.*l2.^7.*l3.^4 - 238654763070797.0*l1.^2.*l2.^7.*l3.^3 - 4038808699633.08*l1.^2.*l2.^7.*l3.^2 - 220700664740.965*l1.^2.*l2.^7.*l3 + 15034599611.6996*l1.^2.*l2.^7 + 2.96793022843935e+21*l1.^2.*l2.^6.*l3.^10 - 2.30563148523979e+19*l1.^2.*l2.^6.*l3.^9 + 1.73234338551178e+19*l1.^2.*l2.^6.*l3.^8 + 1.70014711121925e+17*l1.^2.*l2.^6.*l3.^7 - 1.28096864349185e+17*l1.^2.*l2.^6.*l3.^6 - 88156828044806.3*l1.^2.*l2.^6.*l3.^5 + 417031661011432.0*l1.^2.*l2.^6.*l3.^4 + 5886541091299.67*l1.^2.*l2.^6.*l3.^3 + 451049356349.811*l1.^2.*l2.^6.*l3.^2 - 26438830887.1271*l1.^2.*l2.^6.*l3 + 202688924.603185*l1.^2.*l2.^6 - 2.51628679378898e+21*l1.^2.*l2.^5.*l3.^11 + 1.1480023157196e+19*l1.^2.*l2.^5.*l3.^10 - 2.23482864137358e+19*l1.^2.*l2.^5.*l3.^9 - 1.8221873741533e+17*l1.^2.*l2.^5.*l3.^8 + 1.73797937814058e+17*l1.^2.*l2.^5.*l3.^7 + 713971095298564.0*l1.^2.*l2.^5.*l3.^6 - 561830015505064.0*l1.^2.*l2.^5.*l3.^5 - 8067452612270.45*l1.^2.*l2.^5.*l3.^4 - 659783012961.096*l1.^2.*l2.^5.*l3.^3 + 59766381903.5056*l1.^2.*l2.^5.*l3.^2 + 3563877546.40519*l1.^2.*l2.^5.*l3 - 160472040.066512*l1.^2.*l2.^5 + 1.78001966686823e+21*l1.^2.*l2.^4.*l3.^12 + 5.05409824530644e+18*l1.^2.*l2.^4.*l3.^11 + 2.33906581451469e+19*l1.^2.*l2.^4.*l3.^10 + 4.88628638950524e+16*l1.^2.*l2.^4.*l3.^9 - 2.05903339846073e+17*l1.^2.*l2.^4.*l3.^8 - 488264103881599.0*l1.^2.*l2.^4.*l3.^7 + 830509417619546.0*l1.^2.*l2.^4.*l3.^6 + 5909129506845.74*l1.^2.*l2.^4.*l3.^5 - 102863064130.427*l1.^2.*l2.^4.*l3.^4 - 66950149636.0584*l1.^2.*l2.^4.*l3.^3 - 6296867732.43903*l1.^2.*l2.^4.*l3.^2 + 215540364.696749*l1.^2.*l2.^4.*l3 - 8946550.23595855*l1.^2.*l2.^4 - 1.00196623479682e+21*l1.^2.*l2.^3.*l3.^13 - 1.24667756717559e+19*l1.^2.*l2.^3.*l3.^12 - 1.85258967026386e+19*l1.^2.*l2.^3.*l3.^11 + 6.45270355743277e+16*l1.^2.*l2.^3.*l3.^10 + 1.86353455412083e+17*l1.^2.*l2.^3.*l3.^9 + 88936349836200.0*l1.^2.*l2.^3.*l3.^8 - 829955185651246.0*l1.^2.*l2.^3.*l3.^7 - 4623848662275.7*l1.^2.*l2.^3.*l3.^6 + 429134893275.437*l1.^2.*l2.^3.*l3.^5 + 65692990282.8616*l1.^2.*l2.^3.*l3.^4 + 5805805854.79946*l1.^2.*l2.^3.*l3.^3 - 416242342.268974*l1.^2.*l2.^3.*l3.^2 - 11884239.3961021*l1.^2.*l2.^3.*l3 + 169839.971161965*l1.^2.*l2.^3 + 4.34777633458962e+20*l1.^2.*l2.^2.*l3.^14 + 1.12634189466829e+19*l1.^2.*l2.^2.*l3.^13 + 1.1029766129526e+19*l1.^2.*l2.^2.*l3.^12 - 1.18393082662462e+17*l1.^2.*l2.^2.*l3.^11 - 1.29039725169923e+17*l1.^2.*l2.^2.*l3.^10 + 567869814093551.0*l1.^2.*l2.^2.*l3.^9 + 761904549695259.0*l1.^2.*l2.^2.*l3.^8 - 460803677139.349*l1.^2.*l2.^2.*l3.^7 - 1551037232579.74*l1.^2.*l2.^2.*l3.^6 - 20820098267.7377*l1.^2.*l2.^2.*l3.^5 + 524865366.225963*l1.^2.*l2.^2.*l3.^4 + 298302748.055086*l1.^2.*l2.^2.*l3.^3 + 14690298.4466565*l1.^2.*l2.^2.*l3.^2 + 1474144.9924155*l1.^2.*l2.^2.*l3 + 20146.893983931*l1.^2.*l2.^2 - 1.29196727131924e+20*l1.^2.*l2.*l3.^15 - 5.18818656610028e+18*l1.^2.*l2.*l3.^14 - 4.33571147161073e+18*l1.^2.*l2.*l3.^13 + 6.9623736210271e+16*l1.^2.*l2.*l3.^12 + 5.79921161251011e+16*l1.^2.*l2.*l3.^11 - 437075506246793.0*l1.^2.*l2.*l3.^10 - 387578749936886.0*l1.^2.*l2.*l3.^9 + 891112993291.535*l1.^2.*l2.*l3.^8 + 999092946496.164*l1.^2.*l2.*l3.^7 + 14187941271.7917*l1.^2.*l2.*l3.^6 - 3819838145.07941*l1.^2.*l2.*l3.^5 - 320368453.046296*l1.^2.*l2.*l3.^4 + 21812207.6780169*l1.^2.*l2.*l3.^3 + 2379408.52802636*l1.^2.*l2.*l3.^2 + 85157.4158991442*l1.^2.*l2.*l3 - 1095.6378600823*l1.^2.*l2 + 2.31139995715999e+19*l1.^2.*l3.^16 + 1.47153336666065e+18*l1.^2.*l3.^15 + 1.03017060414577e+18*l1.^2.*l3.^14 - 2.67927473797752e+16*l1.^2.*l3.^13 - 1.63731235407023e+16*l1.^2.*l3.^12 + 255352441362365.0*l1.^2.*l3.^11 + 164745436278446.0*l1.^2.*l3.^10 - 2225711878470.44*l1.^2.*l3.^9 - 851215774465.643*l1.^2.*l3.^8 + 27131551347.0518*l1.^2.*l3.^7 + 5905775784.97353*l1.^2.*l3.^6 - 276364960.194259*l1.^2.*l3.^5 - 63314813.0084641*l1.^2.*l3.^4 + 1965526.656554*l1.^2.*l3.^3 + 451065.679012345*l1.^2.*l3.^2 - 187.823633156966*l1.^2.*l3 - 318.735802469136*l1.^2 - 8.21485647116353e+17*l1.*l2.^17 + 7.83577706326508e+18*l1.*l2.^16.*l3 - 3.18029991626426e+16*l1.*l2.^16 - 3.65908602424017e+19*l1.*l2.^15.*l3.^2 + 2.06289724298222e+17*l1.*l2.^15.*l3 - 2.10662094153716e+16*l1.*l2.^15 + 1.10579520304018e+20*l1.*l2.^14.*l3.^3 - 4.22893934811355e+17*l1.*l2.^14.*l3.^2 + 1.63119260099674e+17*l1.*l2.^14.*l3 + 607863378598739.0*l1.*l2.^14 - 2.37353371685459e+20*l1.*l2.^13.*l3.^4 - 4.81342690029185e+17*l1.*l2.^13.*l3.^3 - 5.81920363615431e+17*l1.*l2.^13.*l3.^2 - 2.61848840011764e+15*l1.*l2.^13.*l3 + 381740453330136.0*l1.*l2.^13 + 3.68436973916698e+20*l1.*l2.^12.*l3.^5 + 5.30680315757176e+18*l1.*l2.^12.*l3.^4 + 1.3102835703341e+18*l1.*l2.^12.*l3.^3 - 46758721430675.3*l1.*l2.^12.*l3.^2 - 2.73778020005678e+15*l1.*l2.^12.*l3 - 8338520980658.15*l1.*l2.^12 - 3.89142937668124e+20*l1.*l2.^11.*l3.^6 - 1.70395312270331e+19*l1.*l2.^11.*l3.^5 - 2.05334668463985e+18*l1.*l2.^11.*l3.^4 + 3.10477910299664e+16*l1.*l2.^11.*l3.^3 + 8.71904542827763e+15*l1.*l2.^11.*l3.^2 + 12141642447757.2*l1.*l2.^11.*l3 - 3243737356422.3*l1.*l2.^11 + 1.7726752237352e+20*l1.*l2.^10.*l3.^7 + 3.47048079511042e+19*l1.*l2.^10.*l3.^6 + 2.20493989317309e+18*l1.*l2.^10.*l3.^5 - 1.21853228048332e+17*l1.*l2.^10.*l3.^4 - 1.68577803492885e+16*l1.*l2.^10.*l3.^3 + 76841951133374.2*l1.*l2.^10.*l3.^2 + 20873057082411.9*l1.*l2.^10.*l3 + 233413627309.146*l1.*l2.^10 + 2.80882413602424e+20*l1.*l2.^9.*l3.^8 - 5.12011095708187e+19*l1.*l2.^9.*l3.^7 - 1.21708677139263e+18*l1.*l2.^9.*l3.^6 + 2.76437561098134e+17*l1.*l2.^9.*l3.^5 + 2.03725873992467e+16*l1.*l2.^9.*l3.^4 - 463933575241152.0*l1.*l2.^9.*l3.^3 - 53474073944221.1*l1.*l2.^9.*l3.^2 - 690452568524.15*l1.*l2.^9.*l3 - 30619783730.6774*l1.*l2.^9 - 8.3483893965458e+20*l1.*l2.^8.*l3.^9 + 5.62010567634968e+19*l1.*l2.^8.*l3.^8 - 1.1389512094318e+18*l1.*l2.^8.*l3.^7 - 4.26860367940607e+17*l1.*l2.^8.*l3.^6 - 9.74075768083499e+15*l1.*l2.^8.*l3.^5 + 1.18487312291732e+15*l1.*l2.^8.*l3.^4 + 71528246253905.9*l1.*l2.^8.*l3.^3 + 2354842320675.32*l1.*l2.^8.*l3.^2 + 210596236737.168*l1.*l2.^8.*l3 - 6074132181.00664*l1.*l2.^8 + 1.2356537029059e+21*l1.*l2.^7.*l3.^10 - 4.52324602144568e+19*l1.*l2.^7.*l3.^9 + 4.40508790414506e+18*l1.*l2.^7.*l3.^8 + 4.75816749278521e+17*l1.*l2.^7.*l3.^7 - 1.8423384361159e+16*l1.*l2.^7.*l3.^6 - 1.97389090704181e+15*l1.*l2.^7.*l3.^5 - 32207296274029.3*l1.*l2.^7.*l3.^4 - 5086911180969.64*l1.*l2.^7.*l3.^3 - 555148909637.159*l1.*l2.^7.*l3.^2 + 19960611771.1644*l1.*l2.^7.*l3 + 576611267.117847*l1.*l2.^7 - 1.31168351457004e+21*l1.*l2.^6.*l3.^11 + 2.35376575424271e+19*l1.*l2.^6.*l3.^10 - 7.54786329926233e+18*l1.*l2.^6.*l3.^9 - 3.672429981165e+17*l1.*l2.^6.*l3.^8 + 5.70989860328557e+16*l1.*l2.^6.*l3.^7 + 2.18032245052615e+15*l1.*l2.^6.*l3.^6 - 111129988624206.0*l1.*l2.^6.*l3.^5 + 7944345747544.48*l1.*l2.^6.*l3.^4 + 854507400105.762*l1.*l2.^6.*l3.^3 - 63063218166.4807*l1.*l2.^6.*l3.^2 - 3740812974.67279*l1.*l2.^6.*l3 + 105126561.653776*l1.*l2.^6 + 1.07345796638641e+21*l1.*l2.^5.*l3.^12 - 3.17686175419262e+18*l1.*l2.^5.*l3.^11 + 9.11109495343344e+18*l1.*l2.^5.*l3.^10 + 1.69079536693311e+17*l1.*l2.^5.*l3.^9 - 8.76234777708895e+16*l1.*l2.^5.*l3.^8 - 1.6197045512798e+15*l1.*l2.^5.*l3.^7 + 289544472093328.0*l1.*l2.^5.*l3.^6 - 7975969529308.92*l1.*l2.^5.*l3.^5 - 765262674720.594*l1.*l2.^5.*l3.^4 + 111964151374.5*l1.*l2.^5.*l3.^3 + 6683406453.65381*l1.*l2.^5.*l3.^2 - 292546477.187814*l1.*l2.^5.*l3 - 1974797.57357581*l1.*l2.^5 - 6.90685619791023e+20*l1.*l2.^4.*l3.^13 - 7.73758374221914e+18*l1.*l2.^4.*l3.^12 - 8.38963453759149e+18*l1.*l2.^4.*l3.^11 + 7.62167159319958e+15*l1.*l2.^4.*l3.^10 + 9.48507027106993e+16*l1.*l2.^4.*l3.^9 + 559838377454954.0*l1.*l2.^4.*l3.^8 - 438479218352182.0*l1.*l2.^4.*l3.^7 + 5908000086068.41*l1.*l2.^4.*l3.^6 + 639787249499.758*l1.*l2.^4.*l3.^5 - 124965488154.754*l1.*l2.^4.*l3.^4 - 5455838920.56256*l1.*l2.^4.*l3.^3 + 668367006.255499*l1.*l2.^4.*l3.^2 + 25018868.7417061*l1.*l2.^4.*l3 - 412232.939713507*l1.*l2.^4 + 3.44602576540088e+20*l1.*l2.^3.*l3.^14 + 8.66416842052532e+18*l1.*l2.^3.*l3.^13 + 5.84013718988892e+18*l1.*l2.^3.*l3.^12 - 7.93027915464202e+16*l1.*l2.^3.*l3.^11 - 7.50804219515317e+16*l1.*l2.^3.*l3.^10 + 195825322746901.0*l1.*l2.^3.*l3.^9 + 424296386181778.0*l1.*l2.^3.*l3.^8 - 2918423288549.17*l1.*l2.^3.*l3.^7 - 481198051508.662*l1.*l2.^3.*l3.^6 + 79534571324.6909*l1.*l2.^3.*l3.^5 + 1297358293.37313*l1.*l2.^3.*l3.^4 - 753346619.873821*l1.*l2.^3.*l3.^3 - 36097829.8279387*l1.*l2.^3.*l3.^2 + 616700.477811406*l1.*l2.^3.*l3 - 41121.4248699022*l1.*l2.^3 - 1.29196727131924e+20*l1.*l2.^2.*l3.^15 - 5.18818656610028e+18*l1.*l2.^2.*l3.^14 - 3.00890658363976e+18*l1.*l2.^2.*l3.^13 + 6.9623736210271e+16*l1.*l2.^2.*l3.^12 + 4.36733526805572e+16*l1.*l2.^2.*l3.^11 - 437075506246793.0*l1.*l2.^2.*l3.^10 - 309293572638915.0*l1.*l2.^2.*l3.^9 + 1958415627842.22*l1.*l2.^2.*l3.^8 + 615612636123.816*l1.*l2.^2.*l3.^7 - 35457025022.8138*l1.*l2.^2.*l3.^6 - 314152567.304578*l1.*l2.^2.*l3.^5 + 505656416.40727*l1.*l2.^2.*l3.^4 + 22357854.1873104*l1.*l2.^2.*l3.^3 - 1775899.50428579*l1.*l2.^2.*l3.^2 + 97678.991442942*l1.*l2.^2.*l3 + 2191.27572016461*l1.*l2.^2 + 3.30670379789799e+19*l1.*l2.*l3.^16 + 1.84285487039745e+18*l1.*l2.*l3.^15 + 1.03462461681744e+18*l1.*l2.*l3.^14 - 3.0673721258521e+16*l1.*l2.*l3.^13 - 1.68024610730944e+16*l1.*l2.*l3.^12 + 255210710127566.0*l1.*l2.*l3.^11 + 138667196126174.0*l1.*l2.*l3.^10 - 1060902583479.31*l1.*l2.*l3.^9 - 396535388444.66*l1.*l2.*l3.^8 + 8569208652.40218*l1.*l2.*l3.^7 + 1233734807.08401*l1.*l2.*l3.^6 - 129835887.339585*l1.*l2.*l3.^5 - 24140699.192035*l1.*l2.*l3.^4 - 1233400.95562281*l1.*l2.*l3.^3 - 181792.942713436*l1.*l2.*l3.^2 + 313.039388594945*l1.*l2.*l3 + 178.184126984127*l1.*l2 - 4.8090190002e+18*l1.*l3.^17 - 3.68742882183072e+17*l1.*l3.^16 - 2.02922863756319e+17*l1.*l3.^15 + 7.66843031463025e+15*l1.*l3.^14 + 3.7679869572501e+15*l1.*l3.^13 - 84731656537169.6*l1.*l3.^12 - 42113760260152.3*l1.*l3.^11 + 736005873208.679*l1.*l3.^10 + 225146312236.476*l1.*l3.^9 - 10897519087.1492*l1.*l3.^8 - 1690761779.10642*l1.*l3.^7 + 181493088.492921*l1.*l3.^6 + 28192935.5365956*l1.*l3.^5 - 1470847.12889779*l1.*l3.^4 - 233086.809937509*l1.*l3.^3 + 3161.69782480894*l1.*l3.^2 + 473.667724867725*l1.*l3 + 2.15703703703704*l1 + 1.3806522621004e+17*l2.^18 - 1.60065545335102e+18*l2.^17.*l3 + 9.45494569700184e+15*l2.^17 + 9.37824823540379e+18*l2.^16.*l3.^2 - 9.28303759341999e+16*l2.^16.*l3 + 2.26430371962935e+15*l2.^16 - 3.68347640474837e+19*l2.^15.*l3.^3 + 4.43522907241177e+17*l2.^15.*l3.^2 - 1.67212741463872e+16*l2.^15.*l3 - 233793607153361.0*l2.^15 + 1.0820483759454e+20*l2.^14.*l3.^4 - 1.36838850451154e+18*l2.^14.*l3.^3 + 5.04629452218332e+16*l2.^14.*l3.^2 + 1.98390575212995e+15*l2.^14.*l3 - 58551401016888.8*l2.^14 - 2.51507726600378e+20*l2.^13.*l3.^5 + 2.99635824543167e+18*l2.^13.*l3.^4 - 5.97042370531174e+16*l2.^13.*l3.^3 - 8.18277625036765e+15*l2.^13.*l3.^2 + 461443554222984.0*l2.^13.*l3 + 3519658997501.61*l2.^13 + 4.78232475724381e+20*l2.^12.*l3.^6 - 4.76529263128893e+18*l2.^12.*l3.^5 - 9.7917562671199e+16*l2.^12.*l3.^4 + 2.13687356938172e+16*l2.^12.*l3.^3 - 1.7084451815603e+15*l2.^12.*l3.^2 - 23385653741789.2*l2.^12.*l3 + 770996718188.735*l2.^12 - 7.5923845745938e+20*l2.^11.*l3.^7 + 5.24663532131811e+18*l2.^11.*l3.^6 + 5.94797998422274e+17*l2.^11.*l3.^5 - 3.81551166874286e+16*l2.^11.*l3.^4 + 4.10134172455169e+15*l2.^11.*l3.^3 + 78259263481361.2*l2.^11.*l3.^2 - 6058419740833.74*l2.^11.*l3 - 7529471848.68216*l2.^11 + 1.01940491607017e+21*l2.^10.*l3.^8 - 2.86742716774529e+18*l2.^10.*l3.^7 - 1.40797310150712e+18*l2.^10.*l3.^6 + 4.49818900163067e+16*l2.^10.*l3.^5 - 7.32888836639421e+15*l2.^10.*l3.^4 - 158786226719502.0*l2.^10.*l3.^3 + 22691523589263.4*l2.^10.*l3.^2 - 57976933234.8516*l2.^10.*l3 + 13150132201.2872*l2.^10 - 1.16583148796442e+21*l2.^9.*l3.^9 - 2.63277260635606e+18*l2.^9.*l3.^8 + 2.16161313037133e+18*l2.^9.*l3.^7 - 2.48288810796869e+16*l2.^9.*l3.^6 + 1.08481391828629e+16*l2.^9.*l3.^5 + 198659947442876.0*l2.^9.*l3.^4 - 57541823054935.2*l2.^9.*l3.^3 + 122353917541.085*l2.^9.*l3.^2 - 102805761781.986*l2.^9.*l3 - 288633524.968636*l2.^9 + 1.13901490129234e+21*l2.^8.*l3.^10 + 9.4240022383571e+18*l2.^8.*l3.^9 - 2.31224074718622e+18*l2.^8.*l3.^8 - 2.73538520369434e+16*l2.^8.*l3.^7 - 1.47509531290267e+16*l2.^8.*l3.^6 - 91983571384371.1*l2.^8.*l3.^5 + 111843042352983.0*l2.^8.*l3.^4 + 153601225713.116*l2.^8.*l3.^3 + 369129713782.403*l2.^8.*l3.^2 + 3540571239.61527*l2.^8.*l3 - 465619164.455507*l2.^8 - 9.49445811342512e+20*l2.^7.*l3.^11 - 1.44196517284457e+19*l2.^7.*l3.^10 + 1.63433160282822e+18*l2.^7.*l3.^9 + 9.24219528163973e+16*l2.^7.*l3.^8 + 1.93089306917034e+16*l2.^7.*l3.^7 - 190274182717284.0*l2.^7.*l3.^6 - 174203045570290.0*l2.^7.*l3.^5 - 1271727795242.42*l2.^7.*l3.^4 - 808931637880.382*l2.^7.*l3.^3 - 4335916952.86219*l2.^7.*l3.^2 + 3182311153.15036*l2.^7.*l3 + 5095365.69362246*l2.^7 + 6.71734294136156e+20*l2.^6.*l3.^12 + 1.54511003499368e+19*l2.^6.*l3.^11 - 4.81329300790995e+17*l2.^6.*l3.^10 - 1.38125263106206e+17*l2.^6.*l3.^9 - 2.35806830797019e+16*l2.^6.*l3.^8 + 542405435574711.0*l2.^6.*l3.^7 + 225008773284054.0*l2.^6.*l3.^6 + 2925199813213.0*l2.^6.*l3.^5 + 1187785642626.88*l2.^6.*l3.^4 - 8530724182.40639*l2.^6.*l3.^3 - 9416563875.44215*l2.^6.*l3.^2 - 44664398.0256864*l2.^6.*l3 + 7958513.40098424*l2.^6 - 3.99090686595394e+20*l2.^5.*l3.^13 - 1.26472791805168e+19*l2.^5.*l3.^12 - 4.79787072850844e+17*l2.^5.*l3.^11 + 1.41819202099229e+17*l2.^5.*l3.^10 + 2.48944067685093e+16*l2.^5.*l3.^9 - 758687299877561.0*l2.^5.*l3.^8 - 238086475331268.0*l2.^5.*l3.^7 - 4051232328183.41*l2.^5.*l3.^6 - 1221772955260.3*l2.^5.*l3.^5 + 37887960710.8829*l2.^5.*l3.^4 + 16568087501.7102*l2.^5.*l3.^3 - 16757144.0802814*l2.^5.*l3.^2 - 45798268.5891143*l2.^5.*l3 - 3297.86351770806*l2.^5 + 1.95717816717947e+20*l2.^4.*l3.^14 + 8.0865571924903e+18*l2.^4.*l3.^13 + 8.38546901702678e+17*l2.^4.*l3.^12 - 1.08901062212036e+17*l2.^4.*l3.^11 - 2.14448756867667e+16*l2.^4.*l3.^10 + 731096952836743.0*l2.^4.*l3.^9 + 206637032323311.0*l2.^4.*l3.^8 + 3627699536695.04*l2.^4.*l3.^7 + 842312142678.177*l2.^4.*l3.^6 - 62280700609.8991*l2.^4.*l3.^5 - 18930596635.9307*l2.^4.*l3.^4 + 249331806.63981*l2.^4.*l3.^3 + 102491909.507981*l2.^4.*l3.^2 + 151701.72181457*l2.^4.*l3 - 25148.3894005705*l2.^4 - 7.69390145230899e+19*l2.^3.*l3.^15 - 3.9917061651706e+18*l2.^3.*l3.^14 - 6.75004225220867e+17*l2.^3.*l3.^13 + 6.28904803242542e+16*l2.^3.*l3.^12 + 1.43162445339773e+16*l2.^3.*l3.^11 - 507185223727228.0*l2.^3.*l3.^10 - 140560566597960.0*l2.^3.*l3.^9 - 2130087585992.17*l2.^3.*l3.^8 - 355230676660.777*l2.^3.*l3.^7 + 60215367386.7902*l2.^3.*l3.^6 + 14257984302.9574*l2.^3.*l3.^5 - 470266010.334411*l2.^3.*l3.^4 - 123982018.161775*l2.^3.*l3.^3 + 288563.057799454*l2.^3.*l3.^2 + 137958.509809045*l2.^3.*l3 - 563.470899470899*l2.^3 + 2.31139995715999e+19*l2.^2.*l3.^16 + 1.47153336666065e+18*l2.^2.*l3.^15 + 3.46807210245514e+17*l2.^2.*l3.^14 - 2.67927473797752e+16*l2.^2.*l3.^13 - 7.07664472278389e+15*l2.^2.*l3.^12 + 255352441362365.0*l2.^2.*l3.^11 + 73521687490755.2*l2.^2.*l3.^10 + 636240371213.633*l2.^2.*l3.^9 + 27877627610.3987*l2.^2.*l3.^8 - 34058755946.2991*l2.^2.*l3.^7 - 6622598642.47112*l2.^2.*l3.^6 + 408631273.011137*l2.^2.*l3.^5 + 86055918.9106354*l2.^2.*l3.^4 - 804678.698320764*l2.^2.*l3.^3 - 253272.945326278*l2.^2.*l3.^2 + 751.294532627866*l2.^2.*l3 + 103.867372134039*l2.^2 - 4.8090190002e+18*l2.*l3.^17 - 3.68742882183072e+17*l2.*l3.^16 - 1.11807744569618e+17*l2.*l3.^15 + 7.66843031463025e+15*l2.*l3.^14 + 2.33775944680111e+15*l2.*l3.^13 - 84731656537169.6*l2.*l3.^12 - 25978023432975.1*l2.*l3.^11 - 7152998256.24688*l2.*l3.^10 + 45613616106.7476*l2.*l3.^9 + 10461361760.5299*l2.*l3.^8 + 1543397519.55572*l2.*l3.^7 - 181151976.14523*l2.*l3.^6 - 27630899.6449716*l2.*l3.^5 + 1472496.06065665*l2.*l3.^4 + 221863.768153809*l2.*l3.^3 - 3412.12933568489*l2.*l3.^2 - 371.538624338624*l2.*l3 + 2.15703703703704*l2 + 5.53436344864651e+17*l3.^18 + 5.15724310745555e+16*l3.^17 + 1.915297717596e+16*l3.^16 - 1.24244602658643e+15*l3.^15 - 428385767786077.0*l3.^14 + 16251848256920.2*l3.^13 + 5627465260702.14*l3.^12 - 84330084705.2401*l3.^11 - 27041784416.5483*l3.^10 + 51312626.6610925*l3.^9 + 18242869.7627493*l3.^8 - 213195.217306373*l3.^7 - 86344.0630090864*l3.^6 + 6595.72703541596*l3.^5 + 2597.06752019505*l3.^4 - 125.215755437979*l3.^3 - 46.9559082892417*l3.^2 + 1.23259259259259*l3 + 0.457851851851852
end
end