text:
@testset "4.2.10 (c+d x)^m (a+b cos)^n" begin
    (a, b, c, d, e, f, m, n, x, ) = @variables a b c d e f m n x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*cos(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(b*cos(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*cos(e+f*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [cos(a + b*x)*(c + d*x)^4, x, 5, -((24*d^3*(c + d*x)*cos(a + b*x))/b^4) + (4*d*(c + d*x)^3*cos(a + b*x))/b^2 + (24*d^4*sin(a + b*x))/b^5 - (12*d^2*(c + d*x)^2*sin(a + b*x))/b^3 + ((c + d*x)^4*sin(a + b*x))/b]
    @test_int [cos(a + b*x)*(c + d*x)^3, x, 4, -((6*d^3*cos(a + b*x))/b^4) + (3*d*(c + d*x)^2*cos(a + b*x))/b^2 - (6*d^2*(c + d*x)*sin(a + b*x))/b^3 + ((c + d*x)^3*sin(a + b*x))/b]
    @test_int [cos(a + b*x)*(c + d*x)^2, x, 3, (2*d*(c + d*x)*cos(a + b*x))/b^2 - (2*d^2*sin(a + b*x))/b^3 + ((c + d*x)^2*sin(a + b*x))/b]
    @test_int [cos(a + b*x)*(c + d*x)^1, x, 2, (d*cos(a + b*x))/b^2 + ((c + d*x)*sin(a + b*x))/b]
    @test_int [cos(a + b*x)/(c + d*x)^1, x, 3, (cos(a - (b*c)/d)*CosIntegral((b*c)/d + b*x))/d - (sin(a - (b*c)/d)*SinIntegral((b*c)/d + b*x))/d]
    @test_int [cos(a + b*x)/(c + d*x)^2, x, 4, -(cos(a + b*x)/(d*(c + d*x))) - (b*CosIntegral((b*c)/d + b*x)*sin(a - (b*c)/d))/d^2 - (b*cos(a - (b*c)/d)*SinIntegral((b*c)/d + b*x))/d^2]
    @test_int [cos(a + b*x)/(c + d*x)^3, x, 5, -(cos(a + b*x)/(2*d*(c + d*x)^2)) - (b^2*cos(a - (b*c)/d)*CosIntegral((b*c)/d + b*x))/(2*d^3) + (b*sin(a + b*x))/(2*d^2*(c + d*x)) + (b^2*sin(a - (b*c)/d)*SinIntegral((b*c)/d + b*x))/(2*d^3)]
    @test_int [cos(a + b*x)/(c + d*x)^4, x, 6, -(cos(a + b*x)/(3*d*(c + d*x)^3)) + (b^2*cos(a + b*x))/(6*d^3*(c + d*x)) + (b^3*CosIntegral((b*c)/d + b*x)*sin(a - (b*c)/d))/(6*d^4) + (b*sin(a + b*x))/(6*d^2*(c + d*x)^2) + (b^3*cos(a - (b*c)/d)*SinIntegral((b*c)/d + b*x))/(6*d^4)]

    @test_int [cos(a + b*x)^2*(c + d*x)^4, x, 6, (3*d^4*x)/(4*b^4) - (d*(c + d*x)^3)/(2*b^2) + (c + d*x)^5/(10*d) - (3*d^3*(c + d*x)*cos(a + b*x)^2)/(2*b^4) + (d*(c + d*x)^3*cos(a + b*x)^2)/b^2 + (3*d^4*cos(a + b*x)*sin(a + b*x))/(4*b^5) - (3*d^2*(c + d*x)^2*cos(a + b*x)*sin(a + b*x))/(2*b^3) + ((c + d*x)^4*cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [cos(a + b*x)^2*(c + d*x)^3, x, 4, -((3*c*d^2*x)/(4*b^2)) - (3*d^3*x^2)/(8*b^2) + (c + d*x)^4/(8*d) - (3*d^3*cos(a + b*x)^2)/(8*b^4) + (3*d*(c + d*x)^2*cos(a + b*x)^2)/(4*b^2) - (3*d^2*(c + d*x)*cos(a + b*x)*sin(a + b*x))/(4*b^3) + ((c + d*x)^3*cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [cos(a + b*x)^2*(c + d*x)^2, x, 4, -((d^2*x)/(4*b^2)) + (c + d*x)^3/(6*d) + (d*(c + d*x)*cos(a + b*x)^2)/(2*b^2) - (d^2*cos(a + b*x)*sin(a + b*x))/(4*b^3) + ((c + d*x)^2*cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [cos(a + b*x)^2*(c + d*x)^1, x, 2, (c*x)/2 + (d*x^2)/4 + (d*cos(a + b*x)^2)/(4*b^2) + ((c + d*x)*cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [cos(a + b*x)^2/(c + d*x)^1, x, 5, (cos(2*a - (2*b*c)/d)*CosIntegral((2*b*c)/d + 2*b*x))/(2*d) + log(c + d*x)/(2*d) - (sin(2*a - (2*b*c)/d)*SinIntegral((2*b*c)/d + 2*b*x))/(2*d)]
    @test_int [cos(a + b*x)^2/(c + d*x)^2, x, 5, -(cos(a + b*x)^2/(d*(c + d*x))) - (b*CosIntegral((2*b*c)/d + 2*b*x)*sin(2*a - (2*b*c)/d))/d^2 - (b*cos(2*a - (2*b*c)/d)*SinIntegral((2*b*c)/d + 2*b*x))/d^2]
    @test_int [cos(a + b*x)^2/(c + d*x)^3, x, 7, -(cos(a + b*x)^2/(2*d*(c + d*x)^2)) - (b^2*cos(2*a - (2*b*c)/d)*CosIntegral((2*b*c)/d + 2*b*x))/d^3 + (b*cos(a + b*x)*sin(a + b*x))/(d^2*(c + d*x)) + (b^2*sin(2*a - (2*b*c)/d)*SinIntegral((2*b*c)/d + 2*b*x))/d^3]

    @test_int [cos(a + b*x)^3*(c + d*x)^4, x, 12, -((160*d^3*(c + d*x)*cos(a + b*x))/(9*b^4)) + (8*d*(c + d*x)^3*cos(a + b*x))/(3*b^2) - (8*d^3*(c + d*x)*cos(a + b*x)^3)/(27*b^4) + (4*d*(c + d*x)^3*cos(a + b*x)^3)/(9*b^2) + (488*d^4*sin(a + b*x))/(27*b^5) - (80*d^2*(c + d*x)^2*sin(a + b*x))/(9*b^3) + (2*(c + d*x)^4*sin(a + b*x))/(3*b) - (4*d^2*(c + d*x)^2*cos(a + b*x)^2*sin(a + b*x))/(9*b^3) + ((c + d*x)^4*cos(a + b*x)^2*sin(a + b*x))/(3*b) - (8*d^4*sin(a + b*x)^3)/(81*b^5)]
    @test_int [cos(a + b*x)^3*(c + d*x)^3, x, 8, -((40*d^3*cos(a + b*x))/(9*b^4)) + (2*d*(c + d*x)^2*cos(a + b*x))/b^2 - (2*d^3*cos(a + b*x)^3)/(27*b^4) + (d*(c + d*x)^2*cos(a + b*x)^3)/(3*b^2) - (40*d^2*(c + d*x)*sin(a + b*x))/(9*b^3) + (2*(c + d*x)^3*sin(a + b*x))/(3*b) - (2*d^2*(c + d*x)*cos(a + b*x)^2*sin(a + b*x))/(9*b^3) + ((c + d*x)^3*cos(a + b*x)^2*sin(a + b*x))/(3*b)]
    @test_int [cos(a + b*x)^3*(c + d*x)^2, x, 6, (4*d*(c + d*x)*cos(a + b*x))/(3*b^2) + (2*d*(c + d*x)*cos(a + b*x)^3)/(9*b^2) - (14*d^2*sin(a + b*x))/(9*b^3) + (2*(c + d*x)^2*sin(a + b*x))/(3*b) + ((c + d*x)^2*cos(a + b*x)^2*sin(a + b*x))/(3*b) + (2*d^2*sin(a + b*x)^3)/(27*b^3)]
    @test_int [cos(a + b*x)^3*(c + d*x)^1, x, 3, (2*d*cos(a + b*x))/(3*b^2) + (d*cos(a + b*x)^3)/(9*b^2) + (2*(c + d*x)*sin(a + b*x))/(3*b) + ((c + d*x)*cos(a + b*x)^2*sin(a + b*x))/(3*b)]
    @test_int [cos(a + b*x)^3/(c + d*x)^1, x, 8, (3*cos(a - (b*c)/d)*CosIntegral((b*c)/d + b*x))/(4*d) + (cos(3*a - (3*b*c)/d)*CosIntegral((3*b*c)/d + 3*b*x))/(4*d) - (3*sin(a - (b*c)/d)*SinIntegral((b*c)/d + b*x))/(4*d) - (sin(3*a - (3*b*c)/d)*SinIntegral((3*b*c)/d + 3*b*x))/(4*d)]
    @test_int [cos(a + b*x)^3/(c + d*x)^2, x, 8, -(cos(a + b*x)^3/(d*(c + d*x))) - (3*b*CosIntegral((3*b*c)/d + 3*b*x)*sin(3*a - (3*b*c)/d))/(4*d^2) - (3*b*CosIntegral((b*c)/d + b*x)*sin(a - (b*c)/d))/(4*d^2) - (3*b*cos(a - (b*c)/d)*SinIntegral((b*c)/d + b*x))/(4*d^2) - (3*b*cos(3*a - (3*b*c)/d)*SinIntegral((3*b*c)/d + 3*b*x))/(4*d^2)]
    @test_int [cos(a + b*x)^3/(c + d*x)^3, x, 12, -(cos(a + b*x)^3/(2*d*(c + d*x)^2)) - (3*b^2*cos(a - (b*c)/d)*CosIntegral((b*c)/d + b*x))/(8*d^3) - (9*b^2*cos(3*a - (3*b*c)/d)*CosIntegral((3*b*c)/d + 3*b*x))/(8*d^3) + (3*b*cos(a + b*x)^2*sin(a + b*x))/(2*d^2*(c + d*x)) + (3*b^2*sin(a - (b*c)/d)*SinIntegral((b*c)/d + b*x))/(8*d^3) + (9*b^2*sin(3*a - (3*b*c)/d)*SinIntegral((3*b*c)/d + 3*b*x))/(8*d^3)]

    @test_int [x^3*cos(a + b*x)^4, x, 8, -((45*x^2)/(128*b^2)) + (3*x^4)/32 - (45*cos(a + b*x)^2)/(128*b^4) + (9*x^2*cos(a + b*x)^2)/(16*b^2) - (3*cos(a + b*x)^4)/(128*b^4) + (3*x^2*cos(a + b*x)^4)/(16*b^2) - (45*x*cos(a + b*x)*sin(a + b*x))/(64*b^3) + (3*x^3*cos(a + b*x)*sin(a + b*x))/(8*b) - (3*x*cos(a + b*x)^3*sin(a + b*x))/(32*b^3) + (x^3*cos(a + b*x)^3*sin(a + b*x))/(4*b)]
    @test_int [x^2*cos(a + b*x)^4, x, 8, -((15*x)/(64*b^2)) + x^3/8 + (3*x*cos(a + b*x)^2)/(8*b^2) + (x*cos(a + b*x)^4)/(8*b^2) - (15*cos(a + b*x)*sin(a + b*x))/(64*b^3) + (3*x^2*cos(a + b*x)*sin(a + b*x))/(8*b) - (cos(a + b*x)^3*sin(a + b*x))/(32*b^3) + (x^2*cos(a + b*x)^3*sin(a + b*x))/(4*b)]
    @test_int [x^1*cos(a + b*x)^4, x, 3, (3*x^2)/16 + (3*cos(a + b*x)^2)/(16*b^2) + cos(a + b*x)^4/(16*b^2) + (3*x*cos(a + b*x)*sin(a + b*x))/(8*b) + (x*cos(a + b*x)^3*sin(a + b*x))/(4*b)]
    @test_int [cos(a + b*x)^4/x^1, x, 8, (1/2)*cos(2*a)*CosIntegral(2*b*x) + (1/8)*cos(4*a)*CosIntegral(4*b*x) + (3*log(x))/8 - (1/2)*sin(2*a)*SinIntegral(2*b*x) - (1/8)*sin(4*a)*SinIntegral(4*b*x)]
    @test_int [cos(a + b*x)^4/x^2, x, 8, -(cos(a + b*x)^4/x) - b*CosIntegral(2*b*x)*sin(2*a) - (1/2)*b*CosIntegral(4*b*x)*sin(4*a) - b*cos(2*a)*SinIntegral(2*b*x) - (1/2)*b*cos(4*a)*SinIntegral(4*b*x)]
    @test_int [cos(a + b*x)^4/x^3, x, 14, -(cos(a + b*x)^4/(2*x^2)) - b^2*cos(2*a)*CosIntegral(2*b*x) - b^2*cos(4*a)*CosIntegral(4*b*x) + (2*b*cos(a + b*x)^3*sin(a + b*x))/x + b^2*sin(2*a)*SinIntegral(2*b*x) + b^2*sin(4*a)*SinIntegral(4*b*x)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [sec(a + b*x)*(c + d*x)^3, x, 9, -((2*im*(c + d*x)^3*atan(ℯ^(im*(a + b*x))))/b) + (3*im*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(im*(a + b*x))))/b^2 - (3*im*d*(c + d*x)^2*PolyLog(2, im*ℯ^(im*(a + b*x))))/b^2 - (6*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(im*(a + b*x))))/b^3 + (6*d^2*(c + d*x)*PolyLog(3, im*ℯ^(im*(a + b*x))))/b^3 - (6*im*d^3*PolyLog(4, (-im)*ℯ^(im*(a + b*x))))/b^4 + (6*im*d^3*PolyLog(4, im*ℯ^(im*(a + b*x))))/b^4]
    @test_int [sec(a + b*x)*(c + d*x)^2, x, 7, -((2*im*(c + d*x)^2*atan(ℯ^(im*(a + b*x))))/b) + (2*im*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(im*(a + b*x))))/b^2 - (2*im*d*(c + d*x)*PolyLog(2, im*ℯ^(im*(a + b*x))))/b^2 - (2*d^2*PolyLog(3, (-im)*ℯ^(im*(a + b*x))))/b^3 + (2*d^2*PolyLog(3, im*ℯ^(im*(a + b*x))))/b^3]
    @test_int [sec(a + b*x)*(c + d*x)^1, x, 5, -((2*im*(c + d*x)*atan(ℯ^(im*(a + b*x))))/b) + (im*d*PolyLog(2, (-im)*ℯ^(im*(a + b*x))))/b^2 - (im*d*PolyLog(2, im*ℯ^(im*(a + b*x))))/b^2]
    @test_int [sec(a + b*x)/(c + d*x)^1, x, 0, Unintegrable(sec(a + b*x)/(c + d*x), x)]

    @test_int [sec(a + b*x)^2*(c + d*x)^3, x, 6, -((im*(c + d*x)^3)/b) + (3*d*(c + d*x)^2*log(1 + ℯ^(2*im*(a + b*x))))/b^2 - (3*im*d^2*(c + d*x)*PolyLog(2, -ℯ^(2*im*(a + b*x))))/b^3 + (3*d^3*PolyLog(3, -ℯ^(2*im*(a + b*x))))/(2*b^4) + ((c + d*x)^3*tan(a + b*x))/b]
    @test_int [sec(a + b*x)^2*(c + d*x)^2, x, 5, -((im*(c + d*x)^2)/b) + (2*d*(c + d*x)*log(1 + ℯ^(2*im*(a + b*x))))/b^2 - (im*d^2*PolyLog(2, -ℯ^(2*im*(a + b*x))))/b^3 + ((c + d*x)^2*tan(a + b*x))/b]
    @test_int [sec(a + b*x)^2*(c + d*x)^1, x, 2, (d*log(cos(a + b*x)))/b^2 + ((c + d*x)*tan(a + b*x))/b]
    @test_int [sec(a + b*x)^2/(c + d*x)^1, x, 0, Unintegrable(sec(a + b*x)^2/(c + d*x), x)]

    @test_int [sec(a + b*x)^3*(c + d*x)^3, x, 15, -((6*im*d^2*(c + d*x)*atan(ℯ^(im*(a + b*x))))/b^3) - (im*(c + d*x)^3*atan(ℯ^(im*(a + b*x))))/b + (3*im*d^3*PolyLog(2, (-im)*ℯ^(im*(a + b*x))))/b^4 + (3*im*d*(c + d*x)^2*PolyLog(2, (-im)*ℯ^(im*(a + b*x))))/(2*b^2) - (3*im*d^3*PolyLog(2, im*ℯ^(im*(a + b*x))))/b^4 - (3*im*d*(c + d*x)^2*PolyLog(2, im*ℯ^(im*(a + b*x))))/(2*b^2) - (3*d^2*(c + d*x)*PolyLog(3, (-im)*ℯ^(im*(a + b*x))))/b^3 + (3*d^2*(c + d*x)*PolyLog(3, im*ℯ^(im*(a + b*x))))/b^3 - (3*im*d^3*PolyLog(4, (-im)*ℯ^(im*(a + b*x))))/b^4 + (3*im*d^3*PolyLog(4, im*ℯ^(im*(a + b*x))))/b^4 - (3*d*(c + d*x)^2*sec(a + b*x))/(2*b^2) + ((c + d*x)^3*sec(a + b*x)*tan(a + b*x))/(2*b)]
    @test_int [sec(a + b*x)^3*(c + d*x)^2, x, 9, -((im*(c + d*x)^2*atan(ℯ^(im*(a + b*x))))/b) + (d^2*atanh(sin(a + b*x)))/b^3 + (im*d*(c + d*x)*PolyLog(2, (-im)*ℯ^(im*(a + b*x))))/b^2 - (im*d*(c + d*x)*PolyLog(2, im*ℯ^(im*(a + b*x))))/b^2 - (d^2*PolyLog(3, (-im)*ℯ^(im*(a + b*x))))/b^3 + (d^2*PolyLog(3, im*ℯ^(im*(a + b*x))))/b^3 - (d*(c + d*x)*sec(a + b*x))/b^2 + ((c + d*x)^2*sec(a + b*x)*tan(a + b*x))/(2*b)]
    @test_int [sec(a + b*x)^3*(c + d*x)^1, x, 6, -((im*(c + d*x)*atan(ℯ^(im*(a + b*x))))/b) + (im*d*PolyLog(2, (-im)*ℯ^(im*(a + b*x))))/(2*b^2) - (im*d*PolyLog(2, im*ℯ^(im*(a + b*x))))/(2*b^2) - (d*sec(a + b*x))/(2*b^2) + ((c + d*x)*sec(a + b*x)*tan(a + b*x))/(2*b)]
    @test_int [sec(a + b*x)^2/(c + d*x)^1, x, 0, Unintegrable(sec(a + b*x)^2/(c + d*x), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^(m/2)*cos(e+f*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [cos(a + b*x)*(c + d*x)^(5/2), x, 8, (5*d*(c + d*x)^(3/2)*cos(a + b*x))/(2*b^2) + (15*d^(5/2)*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(4*b^(7/2)) + (15*d^(5/2)*sqrt(pi/2)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(4*b^(7/2)) - (15*d^2*sqrt(c + d*x)*sin(a + b*x))/(4*b^3) + ((c + d*x)^(5/2)*sin(a + b*x))/b]
    @test_int [cos(a + b*x)*(c + d*x)^(3/2), x, 7, (3*d*sqrt(c + d*x)*cos(a + b*x))/(2*b^2) - (3*d^(3/2)*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(2*b^(5/2)) + (3*d^(3/2)*sqrt(pi/2)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(2*b^(5/2)) + ((c + d*x)^(3/2)*sin(a + b*x))/b]
    @test_int [cos(a + b*x)*(c + d*x)^(1/2), x, 6, -((sqrt(d)*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/b^(3/2)) - (sqrt(d)*sqrt(pi/2)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/b^(3/2) + (sqrt(c + d*x)*sin(a + b*x))/b]
    @test_int [cos(a + b*x)/(c + d*x)^(1/2), x, 5, (sqrt(2*pi)*cos(a - (b*c)/d)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(sqrt(b)*sqrt(d)) - (sqrt(2*pi)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(sqrt(b)*sqrt(d))]
    @test_int [cos(a + b*x)/(c + d*x)^(3/2), x, 6, -((2*cos(a + b*x))/(d*sqrt(c + d*x))) - (2*sqrt(b)*sqrt(2*pi)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/d^(3/2) - (2*sqrt(b)*sqrt(2*pi)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/d^(3/2)]
    @test_int [cos(a + b*x)/(c + d*x)^(5/2), x, 7, -((2*cos(a + b*x))/(3*d*(c + d*x)^(3/2))) - (4*b^(3/2)*sqrt(2*pi)*cos(a - (b*c)/d)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(3*d^(5/2)) + (4*b^(3/2)*sqrt(2*pi)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(3*d^(5/2)) + (4*b*sin(a + b*x))/(3*d^2*sqrt(c + d*x))]
    @test_int [cos(a + b*x)/(c + d*x)^(7/2), x, 8, -((2*cos(a + b*x))/(5*d*(c + d*x)^(5/2))) + (8*b^2*cos(a + b*x))/(15*d^3*sqrt(c + d*x)) + (8*b^(5/2)*sqrt(2*pi)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(15*d^(7/2)) + (8*b^(5/2)*sqrt(2*pi)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(15*d^(7/2)) + (4*b*sin(a + b*x))/(15*d^2*(c + d*x)^(3/2))]

    @test_int [cos(a + b*x)^2*(c + d*x)^(5/2), x, 10, -((5*d*(c + d*x)^(3/2))/(16*b^2)) + (c + d*x)^(7/2)/(7*d) + (5*d*(c + d*x)^(3/2)*cos(a + b*x)^2)/(8*b^2) + (15*d^(5/2)*sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/(128*b^(7/2)) + (15*d^(5/2)*sqrt(pi)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/(128*b^(7/2)) + ((c + d*x)^(5/2)*cos(a + b*x)*sin(a + b*x))/(2*b) - (15*d^2*sqrt(c + d*x)*sin(2*a + 2*b*x))/(64*b^3)]
    @test_int [cos(a + b*x)^2*(c + d*x)^(3/2), x, 9, -((3*d*sqrt(c + d*x))/(16*b^2)) + (c + d*x)^(5/2)/(5*d) + (3*d*sqrt(c + d*x)*cos(a + b*x)^2)/(8*b^2) - (3*d^(3/2)*sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/(32*b^(5/2)) + (3*d^(3/2)*sqrt(pi)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/(32*b^(5/2)) + ((c + d*x)^(3/2)*cos(a + b*x)*sin(a + b*x))/(2*b)]
    @test_int [cos(a + b*x)^2*(c + d*x)^(1/2), x, 8, (c + d*x)^(3/2)/(3*d) - (sqrt(d)*sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/(8*b^(3/2)) - (sqrt(d)*sqrt(pi)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/(8*b^(3/2)) + (sqrt(c + d*x)*sin(2*a + 2*b*x))/(4*b)]
    @test_int [cos(a + b*x)^2/(c + d*x)^(1/2), x, 7, sqrt(c + d*x)/d + (sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/(2*sqrt(b)*sqrt(d)) - (sqrt(pi)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/(2*sqrt(b)*sqrt(d))]
    @test_int [cos(a + b*x)^2/(c + d*x)^(3/2), x, 7, -((2*cos(a + b*x)^2)/(d*sqrt(c + d*x))) - (2*sqrt(b)*sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/d^(3/2) - (2*sqrt(b)*sqrt(pi)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/d^(3/2)]
    @test_int [cos(a + b*x)^2/(c + d*x)^(5/2), x, 9, -((2*cos(a + b*x)^2)/(3*d*(c + d*x)^(3/2))) - (8*b^(3/2)*sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/(3*d^(5/2)) + (8*b^(3/2)*sqrt(pi)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/(3*d^(5/2)) + (8*b*cos(a + b*x)*sin(a + b*x))/(3*d^2*sqrt(c + d*x))]
    @test_int [cos(a + b*x)^2/(c + d*x)^(7/2), x, 9, -((16*b^2)/(15*d^3*sqrt(c + d*x))) - (2*cos(a + b*x)^2)/(5*d*(c + d*x)^(5/2)) + (32*b^2*cos(a + b*x)^2)/(15*d^3*sqrt(c + d*x)) + (32*b^(5/2)*sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/(15*d^(7/2)) + (32*b^(5/2)*sqrt(pi)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/(15*d^(7/2)) + (8*b*cos(a + b*x)*sin(a + b*x))/(15*d^2*(c + d*x)^(3/2))]
    @test_int [cos(a + b*x)^2/(c + d*x)^(9/2), x, 11, -((16*b^2)/(105*d^3*(c + d*x)^(3/2))) - (2*cos(a + b*x)^2)/(7*d*(c + d*x)^(7/2)) + (32*b^2*cos(a + b*x)^2)/(105*d^3*(c + d*x)^(3/2)) + (128*b^(7/2)*sqrt(pi)*cos(2*a - (2*b*c)/d)*FresnelC((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi))))/(105*d^(9/2)) - (128*b^(7/2)*sqrt(pi)*FresnelS((2*sqrt(b)*sqrt(c + d*x))/(sqrt(d)*sqrt(pi)))*sin(2*a - (2*b*c)/d))/(105*d^(9/2)) + (8*b*cos(a + b*x)*sin(a + b*x))/(35*d^2*(c + d*x)^(5/2)) - (128*b^3*cos(a + b*x)*sin(a + b*x))/(105*d^4*sqrt(c + d*x))]

    @test_int [cos(a + b*x)^3*(c + d*x)^(5/2), x, 23, (5*d*(c + d*x)^(3/2)*cos(a + b*x))/(3*b^2) + (5*d*(c + d*x)^(3/2)*cos(a + b*x)^3)/(18*b^2) + (45*d^(5/2)*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(16*b^(7/2)) + (5*d^(5/2)*sqrt(pi/6)*cos(3*a - (3*b*c)/d)*FresnelS((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d)))/(144*b^(7/2)) + (5*d^(5/2)*sqrt(pi/6)*FresnelC((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d))*sin(3*a - (3*b*c)/d))/(144*b^(7/2)) + (45*d^(5/2)*sqrt(pi/2)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(16*b^(7/2)) - (45*d^2*sqrt(c + d*x)*sin(a + b*x))/(16*b^3) + (2*(c + d*x)^(5/2)*sin(a + b*x))/(3*b) + ((c + d*x)^(5/2)*cos(a + b*x)^2*sin(a + b*x))/(3*b) - (5*d^2*sqrt(c + d*x)*sin(3*a + 3*b*x))/(144*b^3)]
    @test_int [cos(a + b*x)^3*(c + d*x)^(3/2), x, 20, (d*sqrt(c + d*x)*cos(a + b*x))/b^2 + (d*sqrt(c + d*x)*cos(a + b*x)^3)/(6*b^2) - (9*d^(3/2)*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(8*b^(5/2)) - (d^(3/2)*sqrt(pi/6)*cos(3*a - (3*b*c)/d)*FresnelC((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d)))/(24*b^(5/2)) + (d^(3/2)*sqrt(pi/6)*FresnelS((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d))*sin(3*a - (3*b*c)/d))/(24*b^(5/2)) + (9*d^(3/2)*sqrt(pi/2)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(8*b^(5/2)) + (2*(c + d*x)^(3/2)*sin(a + b*x))/(3*b) + ((c + d*x)^(3/2)*cos(a + b*x)^2*sin(a + b*x))/(3*b)]
    @test_int [cos(a + b*x)^3*(c + d*x)^(1/2), x, 14, -((3*sqrt(d)*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(4*b^(3/2))) - (sqrt(d)*sqrt(pi/6)*cos(3*a - (3*b*c)/d)*FresnelS((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d)))/(12*b^(3/2)) - (sqrt(d)*sqrt(pi/6)*FresnelC((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d))*sin(3*a - (3*b*c)/d))/(12*b^(3/2)) - (3*sqrt(d)*sqrt(pi/2)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(4*b^(3/2)) + (3*sqrt(c + d*x)*sin(a + b*x))/(4*b) + (sqrt(c + d*x)*sin(3*a + 3*b*x))/(12*b)]
    @test_int [cos(a + b*x)^3/(c + d*x)^(1/2), x, 12, (3*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(2*sqrt(b)*sqrt(d)) + (sqrt(pi/6)*cos(3*a - (3*b*c)/d)*FresnelC((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d)))/(2*sqrt(b)*sqrt(d)) - (sqrt(pi/6)*FresnelS((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d))*sin(3*a - (3*b*c)/d))/(2*sqrt(b)*sqrt(d)) - (3*sqrt(pi/2)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(2*sqrt(b)*sqrt(d))]
    @test_int [cos(a + b*x)^3/(c + d*x)^(3/2), x, 12, -((2*cos(a + b*x)^3)/(d*sqrt(c + d*x))) - (3*sqrt(b)*sqrt(pi/2)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/d^(3/2) - (sqrt(b)*sqrt((3*pi)/2)*cos(3*a - (3*b*c)/d)*FresnelS((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d)))/d^(3/2) - (sqrt(b)*sqrt((3*pi)/2)*FresnelC((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d))*sin(3*a - (3*b*c)/d))/d^(3/2) - (3*sqrt(b)*sqrt(pi/2)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/d^(3/2)]
    @test_int [cos(a + b*x)^3/(c + d*x)^(5/2), x, 18, -((2*cos(a + b*x)^3)/(3*d*(c + d*x)^(3/2))) - (b^(3/2)*sqrt(2*pi)*cos(a - (b*c)/d)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/d^(5/2) - (b^(3/2)*sqrt(6*pi)*cos(3*a - (3*b*c)/d)*FresnelC((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d)))/d^(5/2) + (b^(3/2)*sqrt(6*pi)*FresnelS((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d))*sin(3*a - (3*b*c)/d))/d^(5/2) + (b^(3/2)*sqrt(2*pi)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/d^(5/2) + (4*b*cos(a + b*x)^2*sin(a + b*x))/(d^2*sqrt(c + d*x))]
    @test_int [cos(a + b*x)^3/(c + d*x)^(7/2), x, 19, -((16*b^2*cos(a + b*x))/(5*d^3*sqrt(c + d*x))) - (2*cos(a + b*x)^3)/(5*d*(c + d*x)^(5/2)) + (24*b^2*cos(a + b*x)^3)/(5*d^3*sqrt(c + d*x)) + (2*b^(5/2)*sqrt(2*pi)*cos(a - (b*c)/d)*FresnelS((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d)))/(5*d^(7/2)) + (6*b^(5/2)*sqrt(6*pi)*cos(3*a - (3*b*c)/d)*FresnelS((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d)))/(5*d^(7/2)) + (6*b^(5/2)*sqrt(6*pi)*FresnelC((sqrt(b)*sqrt(6/pi)*sqrt(c + d*x))/sqrt(d))*sin(3*a - (3*b*c)/d))/(5*d^(7/2)) + (2*b^(5/2)*sqrt(2*pi)*FresnelC((sqrt(b)*sqrt(2/pi)*sqrt(c + d*x))/sqrt(d))*sin(a - (b*c)/d))/(5*d^(7/2)) + (4*b*cos(a + b*x)^2*sin(a + b*x))/(5*d^2*(c + d*x)^(3/2))]

    @test_int [x^(3/2)*cos(x), x, 4, (3/2)*sqrt(x)*cos(x) - (3/2)*sqrt(pi/2)*FresnelC(sqrt(2/pi)*sqrt(x)) + x^(3/2)*sin(x)]
    @test_int [x^(1/2)*cos(x), x, 3, (-sqrt(pi/2))*FresnelS(sqrt(2/pi)*sqrt(x)) + sqrt(x)*sin(x)]
    @test_int [cos(x)/x^(1/2), x, 2, sqrt(2*pi)*FresnelC(sqrt(2/pi)*sqrt(x))]
    @test_int [cos(x)/x^(3/2), x, 3, -((2*cos(x))/sqrt(x)) - 2*sqrt(2*pi)*FresnelS(sqrt(2/pi)*sqrt(x))]

    #= ::Subsubsection:: =#
    #=n<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^(m/3)*cos(e+f*x)^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [cos(a + b*x)*(c + d*x)^(4/3), x, 5, (4*d*(c + d*x)^(1/3)*cos(a + b*x))/(3*b^2) + (2*im*d^2*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(2/3)*Gamma(1/3, -((im*b*(c + d*x))/d)))/(9*b^3*(c + d*x)^(2/3)) - (2*im*d^2*((im*b*(c + d*x))/d)^(2/3)*Gamma(1/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(9*b^3*(c + d*x)^(2/3))) + ((c + d*x)^(4/3)*sin(a + b*x))/b]
    @test_int [cos(a + b*x)*(c + d*x)^(2/3), x, 4, (d*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(1/3)*Gamma(2/3, -((im*b*(c + d*x))/d)))/(3*b^2*(c + d*x)^(1/3)) + (d*((im*b*(c + d*x))/d)^(1/3)*Gamma(2/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(3*b^2*(c + d*x)^(1/3))) + ((c + d*x)^(2/3)*sin(a + b*x))/b]
    @test_int [cos(a + b*x)*(c + d*x)^(1/3), x, 4, (d*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(2/3)*Gamma(1/3, -((im*b*(c + d*x))/d)))/(6*b^2*(c + d*x)^(2/3)) + (d*((im*b*(c + d*x))/d)^(2/3)*Gamma(1/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(6*b^2*(c + d*x)^(2/3))) + ((c + d*x)^(1/3)*sin(a + b*x))/b]
    @test_int [cos(a + b*x)/(c + d*x)^(1/3), x, 3, -((im*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(1/3)*Gamma(2/3, -((im*b*(c + d*x))/d)))/(2*b*(c + d*x)^(1/3))) + (im*((im*b*(c + d*x))/d)^(1/3)*Gamma(2/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(2*b*(c + d*x)^(1/3)))]
    @test_int [cos(a + b*x)/(c + d*x)^(2/3), x, 3, -((im*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(2/3)*Gamma(1/3, -((im*b*(c + d*x))/d)))/(2*b*(c + d*x)^(2/3))) + (im*((im*b*(c + d*x))/d)^(2/3)*Gamma(1/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(2*b*(c + d*x)^(2/3)))]
    @test_int [cos(a + b*x)/(c + d*x)^(4/3), x, 4, -((3*cos(a + b*x))/(d*(c + d*x)^(1/3))) + (3*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(1/3)*Gamma(2/3, -((im*b*(c + d*x))/d)))/(2*d*(c + d*x)^(1/3)) + (3*((im*b*(c + d*x))/d)^(1/3)*Gamma(2/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(2*d*(c + d*x)^(1/3)))]
    @test_int [cos(a + b*x)/(c + d*x)^(5/3), x, 4, -((3*cos(a + b*x))/(2*d*(c + d*x)^(2/3))) + (3*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(2/3)*Gamma(1/3, -((im*b*(c + d*x))/d)))/(4*d*(c + d*x)^(2/3)) + (3*((im*b*(c + d*x))/d)^(2/3)*Gamma(1/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(4*d*(c + d*x)^(2/3)))]
    @test_int [cos(a + b*x)/(c + d*x)^(7/3), x, 5, -((3*cos(a + b*x))/(4*d*(c + d*x)^(4/3))) + (9*im*b*ℯ^(im*(a - (b*c)/d))*(-((im*b*(c + d*x))/d))^(1/3)*Gamma(2/3, -((im*b*(c + d*x))/d)))/(8*d^2*(c + d*x)^(1/3)) - (9*im*b*((im*b*(c + d*x))/d)^(1/3)*Gamma(2/3, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*(8*d^2*(c + d*x)^(1/3))) + (9*b*sin(a + b*x))/(4*d^2*(c + d*x)^(1/3))]

    #= ::Subsubsection:: =#
    #=n<0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*cos(e+f*x)^(n/2)=#

    @test_int [x^1*cos(a + b*x)^(1/2), x, 0, Unintegrable(x*sqrt(cos(a + b*x)), x)]
    @test_int [x^0*cos(a + b*x)^(1/2), x, 1, (2*Elliptic.E((1/2)*(a + b*x), 2))/b]
    @test_int [cos(a + b*x)^(1/2)/x^1, x, 0, Unintegrable(sqrt(cos(a + b*x))/x, x)]

    @test_int [x^1*cos(a + b*x)^(3/2), x, 1, (4*cos(a + b*x)^(3/2))/(9*b^2) + (1/3)*Unintegrable(x/sqrt(cos(a + b*x)), x) + (2*x*sqrt(cos(a + b*x))*sin(a + b*x))/(3*b)]
    @test_int [x^0*cos(a + b*x)^(3/2), x, 2, (2*Elliptic.F((1/2)*(a + b*x), 2))/(3*b) + (2*sqrt(cos(a + b*x))*sin(a + b*x))/(3*b)]
    @test_int [cos(a + b*x)^(3/2)/x^1, x, 0, Unintegrable(cos(a + b*x)^(3/2)/x, x)]

    @test_int [x^1*cos(a + b*x)^(3/2) - x/(3*sqrt(cos(a + b*x))), x, 2, (4*cos(a + b*x)^(3/2))/(9*b^2) + (2*x*sqrt(cos(a + b*x))*sin(a + b*x))/(3*b)]

    @test_int [cos(x)^(3/2)/x^3, x, 1, -(cos(x)^(3/2)/(2*x^2)) + (3/8)*Unintegrable(1/(x*sqrt(cos(x))), x) - (9/8)*Unintegrable(cos(x)^(3/2)/x, x) + (3*sqrt(cos(x))*sin(x))/(4*x)]

    @test_int [x^1/cos(a + b*x)^(1/2), x, 0, Unintegrable(x/sqrt(cos(a + b*x)), x)]
    @test_int [x^0/cos(a + b*x)^(1/2), x, 1, (2*Elliptic.F((1/2)*(a + b*x), 2))/b]
    @test_int [1/(x^1*cos(a + b*x)^(1/2)), x, 0, Unintegrable(1/(x*sqrt(cos(a + b*x))), x)]

    @test_int [x^1/cos(a + b*x)^(3/2), x, 1, (4*sqrt(cos(a + b*x)))/b^2 + (2*x*sin(a + b*x))/(b*sqrt(cos(a + b*x))) - Unintegrable(x*sqrt(cos(a + b*x)), x)]
    @test_int [x^0/cos(a + b*x)^(3/2), x, 2, -((2*Elliptic.E((1/2)*(a + b*x), 2))/b) + (2*sin(a + b*x))/(b*sqrt(cos(a + b*x)))]
    @test_int [1/(x^1*cos(a + b*x)^(3/2)), x, 0, Unintegrable(1/(x*cos(a + b*x)^(3/2)), x)]

    @test_int [x^1/cos(a + b*x)^(3/2) + x*sqrt(cos(a + b*x)), x, 2, (4*sqrt(cos(a + b*x)))/b^2 + (2*x*sin(a + b*x))/(b*sqrt(cos(a + b*x)))]

    @test_int [x/cos(x)^(3/2) + x*sqrt(cos(x)), x, 2, 4*sqrt(cos(x)) + (2*x*sin(x))/sqrt(cos(x))]
    @test_int [x/cos(x)^(5/2) - x/(3*sqrt(cos(x))), x, 2, -(4/(3*sqrt(cos(x)))) + (2*x*sin(x))/(3*cos(x)^(3/2))]
    @test_int [x/cos(x)^(7/2) + (3/5)*x*sqrt(cos(x)), x, 3, -(4/(15*cos(x)^(3/2))) + (12*sqrt(cos(x)))/5 + (2*x*sin(x))/(5*cos(x)^(5/2)) + (6*x*sin(x))/(5*sqrt(cos(x)))]
    @test_int [x^2/cos(x)^(3/2) + x^2*sqrt(cos(x)), x, 3, 8*x*sqrt(cos(x)) - 16*Elliptic.E(x/2, 2) + (2*x^2*sin(x))/sqrt(cos(x))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*x^m*sec(e+f*x)^(n/2)=#

    @test_int [x/sec(x)^(3/2) - (1/3)*x*sqrt(sec(x)), x, 4, 4/(9*sec(x)^(3/2)) + (2*x*sin(x))/(3*sqrt(sec(x)))]
    @test_int [x/sec(x)^(5/2) - (3/5)*x/sqrt(sec(x)), x, 4, 4/(25*sec(x)^(5/2)) + (2*x*sin(x))/(5*sec(x)^(3/2))]
    @test_int [x/sec(x)^(7/2) - (5/21)*x*sqrt(sec(x)), x, 5, 4/(49*sec(x)^(7/2)) + 20/(63*sec(x)^(3/2)) + (2*x*sin(x))/(7*sec(x)^(5/2)) + (10*x*sin(x))/(21*sqrt(sec(x)))]
    @test_int [x^2/sec(x)^(3/2) - (1/3)*x^2*sqrt(sec(x)), x, 7, (8*x)/(9*sec(x)^(3/2)) - (16/27)*sqrt(cos(x))*Elliptic.F(x/2, 2)*sqrt(sec(x)) - (16*sin(x))/(27*sqrt(sec(x))) + (2*x^2*sin(x))/(3*sqrt(sec(x)))]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*cos(e+f*x))^n*with*m*symbolic=#

    @test_int [(c + d*x)^m*(b*cos(e + f*x))^n, x, 0, Unintegrable((c + d*x)^m*(b*cos(e + f*x))^n, x)]

    @test_int [cos(a + b*x)^3*(c + d*x)^m, x, 8, -((3*im*ℯ^(im*(a - (b*c)/d))*(c + d*x)^m*Gamma(1 + m, -((im*b*(c + d*x))/d)))/((-((im*b*(c + d*x))/d))^m*(8*b))) + (3*im*(c + d*x)^m*Gamma(1 + m, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*((im*b*(c + d*x))/d)^m*(8*b)) - (im*3^(-1 - m)*ℯ^(3*im*(a - (b*c)/d))*(c + d*x)^m*Gamma(1 + m, -((3*im*b*(c + d*x))/d)))/((-((im*b*(c + d*x))/d))^m*(8*b)) + (im*3^(-1 - m)*(c + d*x)^m*Gamma(1 + m, (3*im*b*(c + d*x))/d))/(ℯ^(3*im*(a - (b*c)/d))*((im*b*(c + d*x))/d)^m*(8*b))]
    @test_int [cos(a + b*x)^2*(c + d*x)^m, x, 5, (c + d*x)^(1 + m)/(2*d*(1 + m)) - (im*2^(-3 - m)*ℯ^(2*im*(a - (b*c)/d))*(c + d*x)^m*Gamma(1 + m, -((2*im*b*(c + d*x))/d)))/((-((im*b*(c + d*x))/d))^m*b) + (im*2^(-3 - m)*(c + d*x)^m*Gamma(1 + m, (2*im*b*(c + d*x))/d))/(ℯ^(2*im*(a - (b*c)/d))*((im*b*(c + d*x))/d)^m*b)]
    @test_int [cos(a + b*x)*(c + d*x)^m, x, 3, -((im*ℯ^(im*(a - (b*c)/d))*(c + d*x)^m*Gamma(1 + m, -((im*b*(c + d*x))/d)))/((-((im*b*(c + d*x))/d))^m*(2*b))) + (im*(c + d*x)^m*Gamma(1 + m, (im*b*(c + d*x))/d))/(ℯ^(im*(a - (b*c)/d))*((im*b*(c + d*x))/d)^m*(2*b))]
    @test_int [sec(a + b*x)*(c + d*x)^m, x, 0, Unintegrable((c + d*x)^m*sec(a + b*x), x)]
    @test_int [sec(a + b*x)^2*(c + d*x)^m, x, 0, Unintegrable((c + d*x)^m*sec(a + b*x)^2, x)]

    @test_int [x^(m + 3)*cos(a + b*x), x, 3, -((ℯ^(im*a)*x^m*Gamma(4 + m, (-im)*b*x))/(((-im)*b*x)^m*(2*b^4))) - (x^m*Gamma(4 + m, im*b*x))/(ℯ^(im*a)*(im*b*x)^m*(2*b^4))]
    @test_int [x^(m + 2)*cos(a + b*x), x, 3, (im*ℯ^(im*a)*x^m*Gamma(3 + m, (-im)*b*x))/(((-im)*b*x)^m*(2*b^3)) - (im*x^m*Gamma(3 + m, im*b*x))/(ℯ^(im*a)*(im*b*x)^m*(2*b^3))]
    @test_int [x^(m + 1)*cos(a + b*x), x, 3, (ℯ^(im*a)*x^m*Gamma(2 + m, (-im)*b*x))/(((-im)*b*x)^m*(2*b^2)) + (x^m*Gamma(2 + m, im*b*x))/(ℯ^(im*a)*(im*b*x)^m*(2*b^2))]
    @test_int [x^(m + 0)*cos(a + b*x), x, 3, -((im*ℯ^(im*a)*x^m*Gamma(1 + m, (-im)*b*x))/(((-im)*b*x)^m*(2*b))) + (im*x^m*Gamma(1 + m, im*b*x))/(ℯ^(im*a)*(im*b*x)^m*(2*b))]
    @test_int [x^(m - 1)*cos(a + b*x), x, 3, ((-(1/2))*ℯ^(im*a)*x^m*Gamma(m, (-im)*b*x))/((-im)*b*x)^m - ((1/2)*x^m*Gamma(m, im*b*x))/(ℯ^(im*a)*(im*b*x)^m)]
    @test_int [x^(m - 2)*cos(a + b*x), x, 3, ((1/2)*im*b*ℯ^(im*a)*x^m*Gamma(-1 + m, (-im)*b*x))/((-im)*b*x)^m - ((1/2)*im*b*x^m*Gamma(-1 + m, im*b*x))/(ℯ^(im*a)*(im*b*x)^m)]
    @test_int [x^(m - 3)*cos(a + b*x), x, 3, ((1/2)*b^2*ℯ^(im*a)*x^m*Gamma(-2 + m, (-im)*b*x))/((-im)*b*x)^m + ((1/2)*b^2*x^m*Gamma(-2 + m, im*b*x))/(ℯ^(im*a)*(im*b*x)^m)]

    @test_int [x^(m + 3)*cos(a + b*x)^2, x, 5, x^(4 + m)/(2*(4 + m)) - (2^(-6 - m)*ℯ^(2*im*a)*x^m*Gamma(4 + m, -2*im*b*x))/(((-im)*b*x)^m*b^4) - (2^(-6 - m)*x^m*Gamma(4 + m, 2*im*b*x))/(ℯ^(2*im*a)*(im*b*x)^m*b^4)]
    @test_int [x^(m + 2)*cos(a + b*x)^2, x, 5, x^(3 + m)/(2*(3 + m)) + (im*2^(-5 - m)*ℯ^(2*im*a)*x^m*Gamma(3 + m, -2*im*b*x))/(((-im)*b*x)^m*b^3) - (im*2^(-5 - m)*x^m*Gamma(3 + m, 2*im*b*x))/(ℯ^(2*im*a)*(im*b*x)^m*b^3)]
    @test_int [x^(m + 1)*cos(a + b*x)^2, x, 5, x^(2 + m)/(2*(2 + m)) + (2^(-4 - m)*ℯ^(2*im*a)*x^m*Gamma(2 + m, -2*im*b*x))/(((-im)*b*x)^m*b^2) + (2^(-4 - m)*x^m*Gamma(2 + m, 2*im*b*x))/(ℯ^(2*im*a)*(im*b*x)^m*b^2)]
    @test_int [x^(m + 0)*cos(a + b*x)^2, x, 5, x^(1 + m)/(2*(1 + m)) - (im*2^(-3 - m)*ℯ^(2*im*a)*x^m*Gamma(1 + m, -2*im*b*x))/(((-im)*b*x)^m*b) + (im*2^(-3 - m)*x^m*Gamma(1 + m, 2*im*b*x))/(ℯ^(2*im*a)*(im*b*x)^m*b)]
    @test_int [x^(m - 1)*cos(a + b*x)^2, x, 5, x^m/(2*m) - (2^(-2 - m)*ℯ^(2*im*a)*x^m*Gamma(m, -2*im*b*x))/((-im)*b*x)^m - (2^(-2 - m)*x^m*Gamma(m, 2*im*b*x))/(ℯ^(2*im*a)*(im*b*x)^m)]
    @test_int [x^(m - 2)*cos(a + b*x)^2, x, 5, -(x^(-1 + m)/(2*(1 - m))) + (im*2^(-1 - m)*b*ℯ^(2*im*a)*x^m*Gamma(-1 + m, -2*im*b*x))/((-im)*b*x)^m - (im*2^(-1 - m)*b*x^m*Gamma(-1 + m, 2*im*b*x))/(ℯ^(2*im*a)*(im*b*x)^m)]
    @test_int [x^(m - 3)*cos(a + b*x)^2, x, 5, -(x^(-2 + m)/(2*(2 - m))) + (b^2*ℯ^(2*im*a)*x^m*Gamma(-2 + m, -2*im*b*x))/(2^m*((-im)*b*x)^m) + (b^2*x^m*Gamma(-2 + m, 2*im*b*x))/(2^m*ℯ^(2*im*a)*(im*b*x)^m)]

    #= ::Section:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(b*sec(e+f*x))^n=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*cos(e+f*x))^n*with*a^2-b^2=0=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*cos(e+f*x))^n=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [(a + a*cos(e + f*x))*(c + d*x)^3, x, 6, (a*(c + d*x)^4)/(4*d) - (6*a*d^3*cos(e + f*x))/f^4 + (3*a*d*(c + d*x)^2*cos(e + f*x))/f^2 - (6*a*d^2*(c + d*x)*sin(e + f*x))/f^3 + (a*(c + d*x)^3*sin(e + f*x))/f]
    @test_int [(a + a*cos(e + f*x))*(c + d*x)^2, x, 5, (a*(c + d*x)^3)/(3*d) + (2*a*d*(c + d*x)*cos(e + f*x))/f^2 - (2*a*d^2*sin(e + f*x))/f^3 + (a*(c + d*x)^2*sin(e + f*x))/f]
    @test_int [(a + a*cos(e + f*x))*(c + d*x)^1, x, 4, (a*(c + d*x)^2)/(2*d) + (a*d*cos(e + f*x))/f^2 + (a*(c + d*x)*sin(e + f*x))/f]
    @test_int [(a + a*cos(e + f*x))/(c + d*x)^1, x, 5, (a*cos(e - (c*f)/d)*CosIntegral((c*f)/d + f*x))/d + (a*log(c + d*x))/d - (a*sin(e - (c*f)/d)*SinIntegral((c*f)/d + f*x))/d]
    @test_int [(a + a*cos(e + f*x))/(c + d*x)^2, x, 6, -(a/(d*(c + d*x))) - (a*cos(e + f*x))/(d*(c + d*x)) - (a*f*CosIntegral((c*f)/d + f*x)*sin(e - (c*f)/d))/d^2 - (a*f*cos(e - (c*f)/d)*SinIntegral((c*f)/d + f*x))/d^2]

    @test_int [(a + a*cos(e + f*x))^2*(c + d*x)^3, x, 10, -((3*a^2*c*d^2*x)/(4*f^2)) - (3*a^2*d^3*x^2)/(8*f^2) + (3*a^2*(c + d*x)^4)/(8*d) - (12*a^2*d^3*cos(e + f*x))/f^4 + (6*a^2*d*(c + d*x)^2*cos(e + f*x))/f^2 - (3*a^2*d^3*cos(e + f*x)^2)/(8*f^4) + (3*a^2*d*(c + d*x)^2*cos(e + f*x)^2)/(4*f^2) - (12*a^2*d^2*(c + d*x)*sin(e + f*x))/f^3 + (2*a^2*(c + d*x)^3*sin(e + f*x))/f - (3*a^2*d^2*(c + d*x)*cos(e + f*x)*sin(e + f*x))/(4*f^3) + (a^2*(c + d*x)^3*cos(e + f*x)*sin(e + f*x))/(2*f)]
    @test_int [(a + a*cos(e + f*x))^2*(c + d*x)^2, x, 9, -((a^2*d^2*x)/(4*f^2)) + (a^2*(c + d*x)^3)/(2*d) + (4*a^2*d*(c + d*x)*cos(e + f*x))/f^2 + (a^2*d*(c + d*x)*cos(e + f*x)^2)/(2*f^2) - (4*a^2*d^2*sin(e + f*x))/f^3 + (2*a^2*(c + d*x)^2*sin(e + f*x))/f - (a^2*d^2*cos(e + f*x)*sin(e + f*x))/(4*f^3) + (a^2*(c + d*x)^2*cos(e + f*x)*sin(e + f*x))/(2*f)]
    @test_int [(a + a*cos(e + f*x))^2*(c + d*x)^1, x, 6, (1/2)*a^2*c*x + (1/4)*a^2*d*x^2 + (a^2*(c + d*x)^2)/(2*d) + (2*a^2*d*cos(e + f*x))/f^2 + (a^2*d*cos(e + f*x)^2)/(4*f^2) + (2*a^2*(c + d*x)*sin(e + f*x))/f + (a^2*(c + d*x)*cos(e + f*x)*sin(e + f*x))/(2*f)]
    @test_int [(a + a*cos(e + f*x))^2/(c + d*x)^1, x, 9, (2*a^2*cos(e - (c*f)/d)*CosIntegral((c*f)/d + f*x))/d + (a^2*cos(2*e - (2*c*f)/d)*CosIntegral((2*c*f)/d + 2*f*x))/(2*d) + (3*a^2*log(c + d*x))/(2*d) - (2*a^2*sin(e - (c*f)/d)*SinIntegral((c*f)/d + f*x))/d - (a^2*sin(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/(2*d)]
    @test_int [(a + a*cos(e + f*x))^2/(c + d*x)^2, x, 9, -((4*a^2*cos(e/2 + (f*x)/2)^4)/(d*(c + d*x))) - (a^2*f*CosIntegral((2*c*f)/d + 2*f*x)*sin(2*e - (2*c*f)/d))/d^2 - (2*a^2*f*CosIntegral((c*f)/d + f*x)*sin(e - (c*f)/d))/d^2 - (2*a^2*f*cos(e - (c*f)/d)*SinIntegral((c*f)/d + f*x))/d^2 - (a^2*f*cos(2*e - (2*c*f)/d)*SinIntegral((2*c*f)/d + 2*f*x))/d^2]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [1/(a + a*cos(e + f*x))*(c + d*x)^3, x, 7, -((im*(c + d*x)^3)/(a*f)) + (6*d*(c + d*x)^2*log(1 + ℯ^(im*(e + f*x))))/(a*f^2) - (12*im*d^2*(c + d*x)*PolyLog(2, -ℯ^(im*(e + f*x))))/(a*f^3) + (12*d^3*PolyLog(3, -ℯ^(im*(e + f*x))))/(a*f^4) + ((c + d*x)^3*tan(e/2 + (f*x)/2))/(a*f)]
    @test_int [1/(a + a*cos(e + f*x))*(c + d*x)^2, x, 6, -((im*(c + d*x)^2)/(a*f)) + (4*d*(c + d*x)*log(1 + ℯ^(im*(e + f*x))))/(a*f^2) - (4*im*d^2*PolyLog(2, -ℯ^(im*(e + f*x))))/(a*f^3) + ((c + d*x)^2*tan(e/2 + (f*x)/2))/(a*f)]
    @test_int [1/(a + a*cos(e + f*x))*(c + d*x)^1, x, 3, (2*d*log(cos(e/2 + (f*x)/2)))/(a*f^2) + ((c + d*x)*tan(e/2 + (f*x)/2))/(a*f)]
    @test_int [1/(a + a*cos(e + f*x))/(c + d*x)^1, x, 0, Unintegrable(1/((c + d*x)*(a + a*cos(e + f*x))), x)]
    @test_int [1/(a + a*cos(e + f*x))/(c + d*x)^2, x, 0, Unintegrable(1/((c + d*x)^2*(a + a*cos(e + f*x))), x)]

    @test_int [1/(a + a*cos(e + f*x))^2*(c + d*x)^3, x, 10, -((im*(c + d*x)^3)/(3*a^2*f)) + (2*d*(c + d*x)^2*log(1 + ℯ^(im*(e + f*x))))/(a^2*f^2) + (4*d^3*log(cos(e/2 + (f*x)/2)))/(a^2*f^4) - (4*im*d^2*(c + d*x)*PolyLog(2, -ℯ^(im*(e + f*x))))/(a^2*f^3) + (4*d^3*PolyLog(3, -ℯ^(im*(e + f*x))))/(a^2*f^4) - (d*(c + d*x)^2*sec(e/2 + (f*x)/2)^2)/(2*a^2*f^2) + (2*d^2*(c + d*x)*tan(e/2 + (f*x)/2))/(a^2*f^3) + ((c + d*x)^3*tan(e/2 + (f*x)/2))/(3*a^2*f) + ((c + d*x)^3*sec(e/2 + (f*x)/2)^2*tan(e/2 + (f*x)/2))/(6*a^2*f)]
    @test_int [1/(a + a*cos(e + f*x))^2*(c + d*x)^2, x, 9, -((im*(c + d*x)^2)/(3*a^2*f)) + (4*d*(c + d*x)*log(1 + ℯ^(im*(e + f*x))))/(3*a^2*f^2) - (4*im*d^2*PolyLog(2, -ℯ^(im*(e + f*x))))/(3*a^2*f^3) - (d*(c + d*x)*sec(e/2 + (f*x)/2)^2)/(3*a^2*f^2) + (2*d^2*tan(e/2 + (f*x)/2))/(3*a^2*f^3) + ((c + d*x)^2*tan(e/2 + (f*x)/2))/(3*a^2*f) + ((c + d*x)^2*sec(e/2 + (f*x)/2)^2*tan(e/2 + (f*x)/2))/(6*a^2*f)]
    @test_int [1/(a + a*cos(e + f*x))^2*(c + d*x)^1, x, 4, (2*d*log(cos(e/2 + (f*x)/2)))/(3*a^2*f^2) - (d*sec(e/2 + (f*x)/2)^2)/(6*a^2*f^2) + ((c + d*x)*tan(e/2 + (f*x)/2))/(3*a^2*f) + ((c + d*x)*sec(e/2 + (f*x)/2)^2*tan(e/2 + (f*x)/2))/(6*a^2*f)]
    @test_int [1/(a + a*cos(e + f*x))^2/(c + d*x)^1, x, 0, Unintegrable(1/((c + d*x)*(a + a*cos(e + f*x))^2), x)]
    @test_int [1/(a + a*cos(e + f*x))^2/(c + d*x)^2, x, 0, Unintegrable(1/((c + d*x)^2*(a + a*cos(e + f*x))^2), x)]

    @test_int [1/(a - a*cos(e + f*x))*(c + d*x)^3, x, 7, -((im*(c + d*x)^3)/(a*f)) - ((c + d*x)^3*cot(e/2 + (f*x)/2))/(a*f) + (6*d*(c + d*x)^2*log(1 - ℯ^(im*(e + f*x))))/(a*f^2) - (12*im*d^2*(c + d*x)*PolyLog(2, ℯ^(im*(e + f*x))))/(a*f^3) + (12*d^3*PolyLog(3, ℯ^(im*(e + f*x))))/(a*f^4)]
    @test_int [1/(a - a*cos(e + f*x))*(c + d*x)^2, x, 6, -((im*(c + d*x)^2)/(a*f)) - ((c + d*x)^2*cot(e/2 + (f*x)/2))/(a*f) + (4*d*(c + d*x)*log(1 - ℯ^(im*(e + f*x))))/(a*f^2) - (4*im*d^2*PolyLog(2, ℯ^(im*(e + f*x))))/(a*f^3)]
    @test_int [1/(a - a*cos(e + f*x))*(c + d*x)^1, x, 3, -(((c + d*x)*cot(e/2 + (f*x)/2))/(a*f)) + (2*d*log(sin(e/2 + (f*x)/2)))/(a*f^2)]
    @test_int [1/(a - a*cos(e + f*x))/(c + d*x)^1, x, 0, Unintegrable(1/((c + d*x)*(a - a*cos(e + f*x))), x)]
    @test_int [1/(a - a*cos(e + f*x))/(c + d*x)^2, x, 0, Unintegrable(1/((c + d*x)^2*(a - a*cos(e + f*x))), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*cos(e+f*x))^(n/2)=#

    #= ::Subsubsection::Closed:: =#
    #=n>0=#

    @test_int [x^3*sqrt(a + a*cos(c + d*x)), x, 5, -((96*sqrt(a + a*cos(c + d*x)))/d^4) + (12*x^2*sqrt(a + a*cos(c + d*x)))/d^2 - (48*x*sqrt(a + a*cos(c + d*x))*tan(c/2 + (d*x)/2))/d^3 + (2*x^3*sqrt(a + a*cos(c + d*x))*tan(c/2 + (d*x)/2))/d]
    @test_int [x^2*sqrt(a + a*cos(c + d*x)), x, 4, (8*x*sqrt(a + a*cos(c + d*x)))/d^2 - (16*sqrt(a + a*cos(c + d*x))*tan(c/2 + (d*x)/2))/d^3 + (2*x^2*sqrt(a + a*cos(c + d*x))*tan(c/2 + (d*x)/2))/d]
    @test_int [x*sqrt(a + a*cos(c + d*x)), x, 3, (4*sqrt(a + a*cos(c + d*x)))/d^2 + (2*x*sqrt(a + a*cos(c + d*x))*tan(c/2 + (d*x)/2))/d]
    @test_int [sqrt(a + a*cos(c + d*x)), x, 1, (2*a*sin(c + d*x))/(d*sqrt(a + a*cos(c + d*x)))]
    @test_int [sqrt(a + a*cos(c + d*x))/x, x, 4, cos(c/2)*sqrt(a + a*cos(c + d*x))*CosIntegral((d*x)/2)*sec(c/2 + (d*x)/2) - sqrt(a + a*cos(c + d*x))*sec(c/2 + (d*x)/2)*sin(c/2)*SinIntegral((d*x)/2)]
    @test_int [sqrt(a + a*cos(c + d*x))/x^2, x, 5, -(sqrt(a + a*cos(c + d*x))/x) - (1/2)*d*sqrt(a + a*cos(c + d*x))*CosIntegral((d*x)/2)*sec(c/2 + (d*x)/2)*sin(c/2) - (1/2)*d*cos(c/2)*sqrt(a + a*cos(c + d*x))*sec(c/2 + (d*x)/2)*SinIntegral((d*x)/2)]
    @test_int [sqrt(a + a*cos(c + d*x))/x^3, x, 6, -(sqrt(a + a*cos(c + d*x))/(2*x^2)) - (1/8)*d^2*cos(c/2)*sqrt(a + a*cos(c + d*x))*CosIntegral((d*x)/2)*sec(c/2 + (d*x)/2) + (1/8)*d^2*sqrt(a + a*cos(c + d*x))*sec(c/2 + (d*x)/2)*sin(c/2)*SinIntegral((d*x)/2) + (d*sqrt(a + a*cos(c + d*x))*tan(c/2 + (d*x)/2))/(4*x)]

    @test_int [x^3*sqrt(a + a*cos(x)), x, 5, -96*sqrt(a + a*cos(x)) + 12*x^2*sqrt(a + a*cos(x)) - 48*x*sqrt(a + a*cos(x))*tan(x/2) + 2*x^3*sqrt(a + a*cos(x))*tan(x/2)]
    @test_int [x^2*sqrt(a + a*cos(x)), x, 4, 8*x*sqrt(a + a*cos(x)) - 16*sqrt(a + a*cos(x))*tan(x/2) + 2*x^2*sqrt(a + a*cos(x))*tan(x/2)]
    @test_int [x*sqrt(a + a*cos(x)), x, 3, 4*sqrt(a + a*cos(x)) + 2*x*sqrt(a + a*cos(x))*tan(x/2)]
    @test_int [sqrt(a + a*cos(x)), x, 1, (2*a*sin(x))/sqrt(a + a*cos(x))]
    @test_int [sqrt(a + a*cos(x))/x, x, 2, sqrt(a + a*cos(x))*CosIntegral(x/2)*sec(x/2)]
    @test_int [sqrt(a + a*cos(x))/x^2, x, 3, -(sqrt(a + a*cos(x))/x) - (1/2)*sqrt(a + a*cos(x))*sec(x/2)*SinIntegral(x/2)]
    @test_int [sqrt(a + a*cos(x))/x^3, x, 4, -(sqrt(a + a*cos(x))/(2*x^2)) - (1/8)*sqrt(a + a*cos(x))*CosIntegral(x/2)*sec(x/2) + (sqrt(a + a*cos(x))*tan(x/2))/(4*x)]

    @test_int [x^3*sqrt(a - a*cos(x)), x, 5, -96*sqrt(a - a*cos(x)) + 12*x^2*sqrt(a - a*cos(x)) + 48*x*sqrt(a - a*cos(x))*cot(x/2) - 2*x^3*sqrt(a - a*cos(x))*cot(x/2)]
    @test_int [x^2*sqrt(a - a*cos(x)), x, 4, 8*x*sqrt(a - a*cos(x)) + 16*sqrt(a - a*cos(x))*cot(x/2) - 2*x^2*sqrt(a - a*cos(x))*cot(x/2)]
    @test_int [x*sqrt(a - a*cos(x)), x, 3, 4*sqrt(a - a*cos(x)) - 2*x*sqrt(a - a*cos(x))*cot(x/2)]
    @test_int [sqrt(a - a*cos(x)), x, 1, -((2*a*sin(x))/sqrt(a - a*cos(x)))]
    @test_int [sqrt(a - a*cos(x))/x, x, 2, sqrt(a - a*cos(x))*csc(x/2)*SinIntegral(x/2)]
    @test_int [sqrt(a - a*cos(x))/x^2, x, 3, -(sqrt(a - a*cos(x))/x) + (1/2)*sqrt(a - a*cos(x))*CosIntegral(x/2)*csc(x/2)]
    @test_int [sqrt(a - a*cos(x))/x^3, x, 4, -(sqrt(a - a*cos(x))/(2*x^2)) - (sqrt(a - a*cos(x))*cot(x/2))/(4*x) - (1/8)*sqrt(a - a*cos(x))*csc(x/2)*SinIntegral(x/2)]

    @test_int [x^3*(a + a*cos(x))^(3/2), x, 9, (-(1280/9))*a*sqrt(a + a*cos(x)) + 16*a*x^2*sqrt(a + a*cos(x)) - (64/27)*a*cos(x/2)^2*sqrt(a + a*cos(x)) + (8/3)*a*x^2*cos(x/2)^2*sqrt(a + a*cos(x)) - (32/9)*a*x*cos(x/2)*sqrt(a + a*cos(x))*sin(x/2) + (4/3)*a*x^3*cos(x/2)*sqrt(a + a*cos(x))*sin(x/2) - (640/9)*a*x*sqrt(a + a*cos(x))*tan(x/2) + (8/3)*a*x^3*sqrt(a + a*cos(x))*tan(x/2)]
    @test_int [x^2*(a + a*cos(x))^(3/2), x, 7, (32/3)*a*x*sqrt(a + a*cos(x)) + (16/9)*a*x*cos(x/2)^2*sqrt(a + a*cos(x)) + (4/3)*a*x^2*cos(x/2)*sqrt(a + a*cos(x))*sin(x/2) - (224/9)*a*sqrt(a + a*cos(x))*tan(x/2) + (8/3)*a*x^2*sqrt(a + a*cos(x))*tan(x/2) + (32/27)*a*sqrt(a + a*cos(x))*sin(x/2)^2*tan(x/2)]
    @test_int [x*(a + a*cos(x))^(3/2), x, 4, (16/3)*a*sqrt(a + a*cos(x)) + (8/9)*a*cos(x/2)^2*sqrt(a + a*cos(x)) + (4/3)*a*x*cos(x/2)*sqrt(a + a*cos(x))*sin(x/2) + (8/3)*a*x*sqrt(a + a*cos(x))*tan(x/2)]
    @test_int [(a + a*cos(x))^(3/2)/x, x, 5, (3/2)*a*sqrt(a + a*cos(x))*CosIntegral(x/2)*sec(x/2) + (1/2)*a*sqrt(a + a*cos(x))*CosIntegral((3*x)/2)*sec(x/2)]
    @test_int [(a + a*cos(x))^(3/2)/x^2, x, 5, -((2*a*cos(x/2)^2*sqrt(a + a*cos(x)))/x) - (3/4)*a*sqrt(a + a*cos(x))*sec(x/2)*SinIntegral(x/2) - (3/4)*a*sqrt(a + a*cos(x))*sec(x/2)*SinIntegral((3*x)/2)]
    @test_int [(a + a*cos(x))^(3/2)/x^3, x, 7, -((a*cos(x/2)^2*sqrt(a + a*cos(x)))/x^2) - (3/16)*a*sqrt(a + a*cos(x))*CosIntegral(x/2)*sec(x/2) - (9/16)*a*sqrt(a + a*cos(x))*CosIntegral((3*x)/2)*sec(x/2) + (3*a*cos(x/2)*sqrt(a + a*cos(x))*sin(x/2))/(2*x)]

    #= ::Subsubsection::Closed:: =#
    #=n<0=#

    @test_int [x^3/sqrt(a + a*cos(c + d*x)), x, 10, -((4*im*x^3*atan(ℯ^((1/2)*im*(c + d*x)))*cos(c/2 + (d*x)/2))/(d*sqrt(a + a*cos(c + d*x)))) + (12*im*x^2*cos(c/2 + (d*x)/2)*PolyLog(2, (-im)*ℯ^((1/2)*im*(c + d*x))))/(d^2*sqrt(a + a*cos(c + d*x))) - (12*im*x^2*cos(c/2 + (d*x)/2)*PolyLog(2, im*ℯ^((1/2)*im*(c + d*x))))/(d^2*sqrt(a + a*cos(c + d*x))) - (48*x*cos(c/2 + (d*x)/2)*PolyLog(3, (-im)*ℯ^((1/2)*im*(c + d*x))))/(d^3*sqrt(a + a*cos(c + d*x))) + (48*x*cos(c/2 + (d*x)/2)*PolyLog(3, im*ℯ^((1/2)*im*(c + d*x))))/(d^3*sqrt(a + a*cos(c + d*x))) - (96*im*cos(c/2 + (d*x)/2)*PolyLog(4, (-im)*ℯ^((1/2)*im*(c + d*x))))/(d^4*sqrt(a + a*cos(c + d*x))) + (96*im*cos(c/2 + (d*x)/2)*PolyLog(4, im*ℯ^((1/2)*im*(c + d*x))))/(d^4*sqrt(a + a*cos(c + d*x)))]
    @test_int [x^2/sqrt(a + a*cos(c + d*x)), x, 8, -((4*im*x^2*atan(ℯ^((1/2)*im*(c + d*x)))*cos(c/2 + (d*x)/2))/(d*sqrt(a + a*cos(c + d*x)))) + (8*im*x*cos(c/2 + (d*x)/2)*PolyLog(2, (-im)*ℯ^((1/2)*im*(c + d*x))))/(d^2*sqrt(a + a*cos(c + d*x))) - (8*im*x*cos(c/2 + (d*x)/2)*PolyLog(2, im*ℯ^((1/2)*im*(c + d*x))))/(d^2*sqrt(a + a*cos(c + d*x))) - (16*cos(c/2 + (d*x)/2)*PolyLog(3, (-im)*ℯ^((1/2)*im*(c + d*x))))/(d^3*sqrt(a + a*cos(c + d*x))) + (16*cos(c/2 + (d*x)/2)*PolyLog(3, im*ℯ^((1/2)*im*(c + d*x))))/(d^3*sqrt(a + a*cos(c + d*x)))]
    @test_int [x/sqrt(a + a*cos(c + d*x)), x, 6, -((4*im*x*atan(ℯ^((1/2)*im*(c + d*x)))*cos(c/2 + (d*x)/2))/(d*sqrt(a + a*cos(c + d*x)))) + (4*im*cos(c/2 + (d*x)/2)*PolyLog(2, (-im)*ℯ^((1/2)*im*(c + d*x))))/(d^2*sqrt(a + a*cos(c + d*x))) - (4*im*cos(c/2 + (d*x)/2)*PolyLog(2, im*ℯ^((1/2)*im*(c + d*x))))/(d^2*sqrt(a + a*cos(c + d*x)))]
    @test_int [1/sqrt(a + a*cos(c + d*x)), x, 2, (sqrt(2)*atanh((sqrt(a)*sin(c + d*x))/(sqrt(2)*sqrt(a + a*cos(c + d*x)))))/(sqrt(a)*d)]
    @test_int [1/(x*sqrt(a + a*cos(c + d*x))), x, 0, Unintegrable(1/(x*sqrt(a + a*cos(c + d*x))), x)]

    @test_int [x^3/sqrt(a - a*cos(x)), x, 10, -((4*x^3*atanh(ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x))) + (12*im*x^2*PolyLog(2, -ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) - (12*im*x^2*PolyLog(2, ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) - (48*x*PolyLog(3, -ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) + (48*x*PolyLog(3, ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) - (96*im*PolyLog(4, -ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) + (96*im*PolyLog(4, ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x))]
    @test_int [x^2/sqrt(a - a*cos(x)), x, 8, -((4*x^2*atanh(ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x))) + (8*im*x*PolyLog(2, -ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) - (8*im*x*PolyLog(2, ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) - (16*PolyLog(3, -ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) + (16*PolyLog(3, ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x))]
    @test_int [x/sqrt(a - a*cos(x)), x, 6, -((4*x*atanh(ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x))) + (4*im*PolyLog(2, -ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x)) - (4*im*PolyLog(2, ℯ^((im*x)/2))*sin(x/2))/sqrt(a - a*cos(x))]
    @test_int [1/sqrt(a - a*cos(x)), x, 2, -((sqrt(2)*atanh((sqrt(a)*sin(x))/(sqrt(2)*sqrt(a - a*cos(x)))))/sqrt(a))]
    @test_int [1/(x*sqrt(a - a*cos(x))), x, 0, Unintegrable(1/(x*sqrt(a - a*cos(x))), x)]

    @test_int [x^3/(a + a*cos(x))^(3/2), x, 16, -((3*x^2)/(a*sqrt(a + a*cos(x)))) - (24*im*x*atan(ℯ^((im*x)/2))*cos(x/2))/(a*sqrt(a + a*cos(x))) - (im*x^3*atan(ℯ^((im*x)/2))*cos(x/2))/(a*sqrt(a + a*cos(x))) + (24*im*cos(x/2)*PolyLog(2, (-im)*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) + (3*im*x^2*cos(x/2)*PolyLog(2, (-im)*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) - (24*im*cos(x/2)*PolyLog(2, im*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) - (3*im*x^2*cos(x/2)*PolyLog(2, im*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) - (12*x*cos(x/2)*PolyLog(3, (-im)*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) + (12*x*cos(x/2)*PolyLog(3, im*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) - (24*im*cos(x/2)*PolyLog(4, (-im)*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) + (24*im*cos(x/2)*PolyLog(4, im*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) + (x^3*tan(x/2))/(2*a*sqrt(a + a*cos(x)))]
    @test_int [x^2/(a + a*cos(x))^(3/2), x, 10, -((2*x)/(a*sqrt(a + a*cos(x)))) - (im*x^2*atan(ℯ^((im*x)/2))*cos(x/2))/(a*sqrt(a + a*cos(x))) + (4*atanh(sin(x/2))*cos(x/2))/(a*sqrt(a + a*cos(x))) + (2*im*x*cos(x/2)*PolyLog(2, (-im)*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) - (2*im*x*cos(x/2)*PolyLog(2, im*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) - (4*cos(x/2)*PolyLog(3, (-im)*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) + (4*cos(x/2)*PolyLog(3, im*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) + (x^2*tan(x/2))/(2*a*sqrt(a + a*cos(x)))]
    @test_int [x^1/(a + a*cos(x))^(3/2), x, 7, -(1/(a*sqrt(a + a*cos(x)))) - (im*x*atan(ℯ^((im*x)/2))*cos(x/2))/(a*sqrt(a + a*cos(x))) + (im*cos(x/2)*PolyLog(2, (-im)*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) - (im*cos(x/2)*PolyLog(2, im*ℯ^((im*x)/2)))/(a*sqrt(a + a*cos(x))) + (x*tan(x/2))/(2*a*sqrt(a + a*cos(x)))]
    @test_int [1/(x*(a + a*cos(x))^(3/2)), x, 0, Unintegrable(1/(x*(a + a*cos(x))^(3/2)), x)]

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*cos(e+f*x))^(n/3)=#

    #= Used*to*hang*Rubi =#
    @test_int [(a + a*cos(c + d*x))^(1/3)/x, x, 0, Unintegrable((a + a*cos(c + d*x))^(1/3)/x, x)]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+a*cos(e+f*x))^n*with*m*symbolic=#

    #= ::Section::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*cos(e+f*x))^n=#

    #= ::Subsection::Closed:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*cos(e+f*x))^n=#

    #= [x^3/(a + b*cos(c + d*x)), x, 12, -((im*x^3*log(1 + (b*ℯ^(im*c + im*d*x))/(a - sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d)) + (im*x^3*log(1 + (b*ℯ^(im*c + im*d*x))/(a + sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d) - (3*x^2*PolyLog(2, -((b*ℯ^(im*c + im*d*x))/(a - sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^2) + (3*x^2*PolyLog(2, -((b*ℯ^(im*c + im*d*x))/(a + sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^2) - (6*im*x*PolyLog(3, -((b*ℯ^(im*c + im*d*x))/(a - sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^3) + (6*im*x*PolyLog(3, -((b*ℯ^(im*c + im*d*x))/(a + sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^3) + (6*PolyLog(4, -((b*ℯ^(im*c + im*d*x))/(a - sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^4) - (6*PolyLog(4, -((b*ℯ^(im*c + im*d*x))/(a + sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^4)] =#
    @test_int [x^3/(a + b*cos(x)), x, 12, -((im*x^3*log(1 + (b*ℯ^(im*x))/(a - sqrt(a^2 - b^2))))/sqrt(a^2 - b^2)) + (im*x^3*log(1 + (b*ℯ^(im*x))/(a + sqrt(a^2 - b^2))))/sqrt(a^2 - b^2) - (3*x^2*PolyLog(2, -((b*ℯ^(im*x))/(a - sqrt(a^2 - b^2)))))/sqrt(a^2 - b^2) + (3*x^2*PolyLog(2, -((b*ℯ^(im*x))/(a + sqrt(a^2 - b^2)))))/sqrt(a^2 - b^2) - (6*im*x*PolyLog(3, -((b*ℯ^(im*x))/(a - sqrt(a^2 - b^2)))))/sqrt(a^2 - b^2) + (6*im*x*PolyLog(3, -((b*ℯ^(im*x))/(a + sqrt(a^2 - b^2)))))/sqrt(a^2 - b^2) + (6*PolyLog(4, -((b*ℯ^(im*x))/(a - sqrt(a^2 - b^2)))))/sqrt(a^2 - b^2) - (6*PolyLog(4, -((b*ℯ^(im*x))/(a + sqrt(a^2 - b^2)))))/sqrt(a^2 - b^2)]
    @test_int [x^2/(a + b*cos(c + d*x)), x, 10, -((im*x^2*log(1 + (b*ℯ^(im*(c + d*x)))/(a - sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d)) + (im*x^2*log(1 + (b*ℯ^(im*(c + d*x)))/(a + sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d) - (2*x*PolyLog(2, -((b*ℯ^(im*(c + d*x)))/(a - sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^2) + (2*x*PolyLog(2, -((b*ℯ^(im*(c + d*x)))/(a + sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^2) - (2*im*PolyLog(3, -((b*ℯ^(im*(c + d*x)))/(a - sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^3) + (2*im*PolyLog(3, -((b*ℯ^(im*(c + d*x)))/(a + sqrt(a^2 - b^2)))))/(sqrt(a^2 - b^2)*d^3)]
    @test_int [x^1/(a + b*cos(c + d*x)), x, 8, -((im*x*log(1 + (b*ℯ^(im*(c + d*x)))/(a - sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d)) + (im*x*log(1 + (b*ℯ^(im*(c + d*x)))/(a + sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d) - PolyLog(2, -((b*ℯ^(im*(c + d*x)))/(a - sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d^2) + PolyLog(2, -((b*ℯ^(im*(c + d*x)))/(a + sqrt(a^2 - b^2))))/(sqrt(a^2 - b^2)*d^2)]
    @test_int [1/(x^1*(a + b*cos(x))), x, 0, Unintegrable(1/(x*(a + b*cos(x))), x)]

    @test_int [(e + f*x)/(a + b*cos(c + d*x))^2, x, 11, -((im*a*(e + f*x)*log(1 + (b*ℯ^(im*(c + d*x)))/(a - sqrt(a^2 - b^2))))/((a^2 - b^2)^(3/2)*d)) + (im*a*(e + f*x)*log(1 + (b*ℯ^(im*(c + d*x)))/(a + sqrt(a^2 - b^2))))/((a^2 - b^2)^(3/2)*d) - (f*log(a + b*cos(c + d*x)))/((a^2 - b^2)*d^2) - (a*f*PolyLog(2, -((b*ℯ^(im*(c + d*x)))/(a - sqrt(a^2 - b^2)))))/((a^2 - b^2)^(3/2)*d^2) + (a*f*PolyLog(2, -((b*ℯ^(im*(c + d*x)))/(a + sqrt(a^2 - b^2)))))/((a^2 - b^2)^(3/2)*d^2) - (b*(e + f*x)*sin(c + d*x))/((a^2 - b^2)*d*(a + b*cos(c + d*x)))]

    #= ::Subsection:: =#
    #=Integrands*of*the*form*(c+d*x)^m*(a+b*cos(e+f*x))^n*with*m*symbolic=#
    end

max_stars_repo_path: test/4 Trig functions/4.2 Cosine/4.2.10 (c+d x)^m (a+b cos)^n.jl


text:
module AtomShell

using Sockets
using WebIO

abstract type Shell end

_shell = nothing

include("install.jl")
include("process.jl")
include("window.jl")
include("webio.jl")

end

ext: jl


text:
# TODO: rewrite this file based on `benchmarks.jl`, using atomic interaction
# functions, like move, add, etc.


text:
# Based on this guide: https://ase.tufts.edu/chemistry/lin/images/FortranMD_TeachersGuide.pdf

using Atomistic
using AtomsBase
using NBodySimulator
using Plots
using Unitful
using UnitfulAtomic

N = 864
element = :Ar
box_size = 3.47786u"nm"
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
initial_system = DynamicSystem(initial_bodies, box_size)
eq_simulator = NBSimulator(
    Δt = Δt,
    steps = 2000,
    t₀ = 0.0u"s",
    thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
)
potential = LJPotential(
    ϵ = 1.657e-21u"J",
    σ = 0.34u"nm",
    R = 0.765u"nm"
)
eq_result = @time simulate(initial_system, eq_simulator, potential)

temp = @time plot_temperature(eq_result, 10)
energy = @time plot_energy(eq_result, 10)

prod_simulator = NBSimulator(
    Δt = Δt,
    steps = 5000,
    t₀ = eq_simulator.steps * Δt
)
prod_result = @time simulate(get_system(eq_result), prod_simulator, potential)

display(@time plot_temperature!(temp, prod_result, 10))
display(@time plot_energy!(energy, prod_result, 10))

rdf = @time plot_rdf(prod_result, potential.σ, 0.05)
display(rdf)
savefig(rdf, "artifacts/argon_lj_rdf.svg")

;

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
using JuMP
using DisjunctiveProgramming

m = Model()
@variable(m, -5 ≤ x ≤ 10)
@disjunction(
    m,
    begin
        exp(x) ≤ 2
        -3 ≤ x
    end,
    begin
        3 ≤ exp(x)
        5 ≤ x
    end,
    reformulation=:hull,
    name=:z
)
print(m)

# Feasibility
# Subject to
#  XOR(disj_z) : z[1] + z[2] == 1.0             <- XOR constraint
#  x_aggregation : x - x_z1 - x_z2 == 0.0       <- aggregation of disaggregated variables
#  x_z1_lb : -5 z[1] - x_z1 <= 0.0              <- lower-bound constraint on disaggregated variable x_z1 (x in 1st disjunct)
#  x_z1_ub : -10 z[1] + x_z1 <= 0.0             <- upper-bound constraint on disaggregated variable x_z1 (x in 1st disjunct)
#  x_z2_lb : -5 z[2] - x_z2 <= 0.0              <- lower-bound constraint on disaggregated variable x_z2 (x in 2nd disjunct)
#  x_z2_ub : -10 z[2] + x_z2 <= 0.0             <- upper-bound constraint on disaggregated variable x_z2 (x in 2nd disjunct)
#  x >= -5.0                                    <- lower-bound on x
#  x_z1 >= -5.0                                 <- lower-bound on x_z1 (disaggregated x in 1st disjunct)
#  x_z2 >= -5.0                                 <- lower-bound on x_z2 (disaggregated x in 2nd disjunct)
#  x <= 10.0                                    <- upper-bound on x
#  x_z1 <= 10.0                                 <- upper-bound on x_z1 (disaggregated x in 1st disjunct)
#  x_z2 <= 10.0                                 <- upper-bound on x_z2 (disaggregated x in 2nd disjunct)
#  z[1] >= 0.0                                  <- lower bound on binary
#  z[2] >= 0.0                                  <- lower bound on binary
#  z[1] <= 1.0                                  <- upper bound on binary
#  z[2] <= 1.0                                  <- upper bound on binary
#  z[1] binary                                  <- indicator variable (1st disjunct) is binary
#  z[2] binary                                  <- indicator variable (2nd disjunct) is binary
#  Perspective Functions:
#  (-1.0e-6 + -1.9999989999999999 * z[1]) + (1.0e-6 + 0.999999 * z[1]) * exp(x_z1 / (1.0e-6 + 0.999999 * z[1])) <= 0
#  (1.0000000000000002e-6 + 2.999999 * z[2]) + (-1.0e-6 + -0.999999 * z[2]) * exp(x_z2 / (1.0e-6 + 0.999999 * z[2])) <= 0
#  -1.0 * x_z1 + -3.0 * z[1] <= 0
#  -1.0 * x_z2 + 5.0 * z[2] <= 0


text:
# pick a nice default x range given a distribution
function default_range(dist::Distribution, alpha = 0.0001)
    minval = isfinite(minimum(dist)) ? minimum(dist) : quantile(dist, alpha)
    maxval = isfinite(maximum(dist)) ? maximum(dist) : quantile(dist, 1-alpha)
    minval, maxval
end

function default_range(m::Distributions.UnivariateMixture, alpha = 0.0001)
    mapreduce(_minmax, 1:Distributions.ncomponents(m)) do k
        default_range(Distributions.component(m, k), alpha)
    end
end

_minmax((xmin, xmax), (ymin, ymax)) = (min(xmin, ymin), max(xmax, ymax))

yz_args(dist) = default_range(dist)
function yz_args(dist::DiscreteUnivariateDistribution)
    minval, maxval = extrema(dist)
    if isfinite(minval) && isfinite(maxval)  # bounded
        sup = support(dist)
        return sup isa AbstractVector ? (sup,) : ([sup...],)
    else  # unbounded
        return (UnitRange(promote(default_range(dist)...)...),)
    end
end

# this "user recipe" adds a default x vector based on the distribution's μ and σ
@recipe function f(dist::Distribution)
    if dist isa DiscreteUnivariateDistribution
        seriestype --> :sticks
    end
    (dist, yz_args(dist)...)
end

@recipe function f(m::Distributions.UnivariateMixture; components = true)
    if m isa DiscreteUnivariateDistribution
        seriestype --> :sticks
    end
    if components
        for k in 1:Distributions.ncomponents(m)
            c = Distributions.component(m, k)
            @series begin
                (c, yz_args(c)...)
            end
        end
    else
        (m, yz_args(m)...)
    end
end

@recipe function f(distvec::AbstractArray{<:Distribution}, yz...)
    for di in distvec
        @series begin
            seriesargs = isempty(yz) ? yz_args(di) : yz
            if di isa DiscreteUnivariateDistribution
                seriestype --> :sticks
            end
            (di, seriesargs...)
        end
    end
end

# this "type recipe" replaces any instance of a distribution with a function mapping xi to yi
@recipe function f(::Type{T}, dist::T; func = pdf) where T<:Distribution
    xi -> func(dist, xi)
end

#-----------------------------------------------------------------------------
# qqplots

@recipe function f(h::QQPair; qqline = :identity)
    if qqline in (:fit, :quantile, :identity, :R)
        xs = [extrema(h.qx)...]
        if qqline == :identity
            ys = xs
        elseif qqline == :fit
            itc, slp = hcat(fill!(similar(h.qx), 1), h.qx) \ h.qy
            ys = slp .* xs .+ itc
        else # if qqline == :quantile || qqline == :R
            quantx, quanty = quantile(h.qx, [0.25, 0.75]), quantile(h.qy, [0.25, 0.75])
            slp = diff(quanty) ./ diff(quantx)
            ys = quanty .+ slp .* (xs .- quantx)
        end

        @series begin
            primary := false
            seriestype := :path
            xs, ys
        end
    end

    seriestype --> :scatter
    legend --> false
    h.qx, h.qy
end

loc(D::Type{T}, x) where T<:Distribution = fit(D, x), x
loc(D, x) = D, x

@userplot QQPlot
recipetype(::Val{:qqplot}, args...) = QQPlot(args)
@recipe f(h::QQPlot) = qqbuild(loc(h.args[1], h.args[2])...)

@userplot QQNorm
recipetype(::Val{:qqnorm}, args...) = QQNorm(args)
@recipe f(h::QQNorm) = QQPlot((Normal, h.args[1]))


text:
function compute_marginals(pw::Vector, pogw::Matrix, pagow)

    card_o = size(pogw,1)
    card_a = size(pagow,1)

    #compute p(o)
    #p(o) = ∑_w p(o|w)p(w)
    po = pogw * pw

    #add some small value to prevent NaNs in the KL-terms
    po += eps()
    po /= sum(po)

    #compute p(a)
    #p(a) = ∑_w,o p(w)p(o|w)p(a|o,w)
    #compute p(a|w)
    pagw = marginalizeo(pogw, pagow)
    #p(a) = ∑_w p(a|w)p(w)
    pa = pagw * pw

    #add some small value to prevent NaNs in the KL-terms
    pa += eps()
    pa /= sum(pa)

    return po, pa, pagw
end

function marginalizeo(pogw::Matrix, pagow)
    card_a = size(pagow,1)
    card_w = size(pogw,2)

    #compute p(a|w)
    pagw = zeros(card_a,card_w)

    for j in 1:card_w
        #p(a|w) = ∑_o p(o|w)p(a|o,w)
        pagw[:,j] = pagow[:,:,j] * pogw[:,j]
    end

    return pagw
end

function compute_pago_iteration(pogw::Matrix, pagow, β2, β3,
                                U_pre::Matrix, pa::Vector, po::Vector, pw::Vector)

    card_a = size(U_pre,1)
    card_w = size(U_pre,2)
    card_o = size(pogw,1)

    pago = zeros(card_a,card_o)

    for k in 1:card_o
        #compute p(w|o=k)
        pwgo_k = vec(pogw[k,:]).*pw / po[k]

        #compute p(a|o=k)
        if(β3==0)
            #sequential case
            #p(a|o) ∝ p(a)exp( β2 ∑_w p(w|o)U(a,w) )
            pago[:,k] = boltzmanndist(pa, β2, U_pre*pwgo_k)
        else
            #general case
            #p(a|o) = ∑_w p(w|o)p(a|o,w)   with p(w|o) = p(o|w)p(w)/p(o)
            pago[:,k] = squeeze(pagow[:,k,:],2) * pwgo_k

            pago[:,k] += eps() #add some small value to prevent NaNs in the KL-terms
            pago[:,k] = pago[:,k] / sum(pago[:,k])
        end
    end

    return pago

end

function compute_pagow_iteration(pago::Matrix, β2, β3, U_pre::Matrix, pa::Vector)

    card_a = size(U_pre,1)
    card_w = size(U_pre,2)
    card_o = size(pago,2)

    pagow = zeros(card_a,card_o,card_w)

    for j in 1:card_w
        if(β3==0)
            #sequential case: p(a|o,w)=p(a|o) ∀w
            #save some computation by implementing this directly rather than evaluating the
            #expression for the general case
            pagow[:,:,j] = pago
        else
            #general case
            #p(a|o,w) ∝ p(a|o) exp( β3 U(a,w) - β3/β2 log(p(a|o)/p(a)) )
            for k in 1:card_o
                pagow_util_kj = U_pre[:,j] - (1/β2)*log_bits(pago[:,k]./pa)
                pagow[:,k,j] = boltzmanndist(vec(pago[:,k]), β3, pagow_util_kj)
            end
        end
    end

    return pagow

end

function compute_pogw_iteration(pago::Matrix, pagow, β1, β2, β3, U_pre::Matrix, pa::Vector, po::Vector)

    card_w = size(U_pre,2)
    card_o = size(po,1)

    pogw = zeros(card_o,card_w)

    #------- compute the conditional EU and the two KL terms -----------#
    EU_cond = zeros(card_o,card_w)
    DKL_a = zeros(card_o,card_w)
    DKL_ago = zeros(card_o,card_w)

    for j in 1:card_w
        for k in 1:card_o
            if(β3==0)
                #sequential case
                EU_cond[k,j] = (pago[:,k]' * U_pre[:,j])[1]  # E[U(a,w=j) = ∑_a p(a|o=k) U(a,w=j)
                                                                #use the (...)[1] syntax to get a scalar
                DKL_a[k,j] = kl_divergence_bits(vec(pago[:,k]),pa)

            else
                EU_cond[k,j] = (pagow[:,k,j]' * U_pre[:,j])[1]  # E[U(a,w=j) = ∑_a p(a|o=k,w=j) U(a,w=j)
                                                                #use the (...)[1] syntax to get a scalar
                DKL_a[k,j] = kl_divergence_bits(vec(pagow[:,k,j]),pa)
            end
            DKL_ago[k,j] = kl_divergence_bits(vec(pagow[:,k,j]),vec(pago[:,k]))
        end
    end

    if(β3==0)
        #sequential case
        #we have β3=0 which would lead to ∞*0=NaN in the computation for pogw_util
        #DKL_ago must be zero in the sequential case, therefore β3 does not matter

        #check that KL is really zero
        if(sum(DKL_ago) > 0)
            error("Sequential case: D_KL( p(a|o,w)||p(a|o) )=$(sum(DKL_ago)) is nonzero which violates sequential case assumption!")
        end

        #p(o|w) ∝ exp( β1 E[U] - 1/β2 DKL(p(a|o)||p(a)) )
        pogw_util = EU_cond - 1/β2 * DKL_a
    else
        #general case
        #p(o|w) ∝ p(o) exp( β1 (E[U] - 1/β2 DKL(p(a|o,w)||p(a))) - β1 DKL(p(a|o,w)||p(a|o)) (1/β3-1/β2) )
        pogw_util = EU_cond - 1/β2 * DKL_a - ( (1/β3)-(1/β2) ) * DKL_ago
    end

    #comptue p(o|w)
    for j in 1:card_w
        pogw[:,j] = boltzmanndist(po, β1, vec(pogw_util[:,j]))
    end

    return pogw
end

#This function performs Blahut-Arimoto iterations for the three-variable general case
#and initializes p(o|w) and p(a|o,w) either uniformly or randomly, depending on the value
#of the optional argument 'init_uniform'
function threevarBAiterations(cardinality_obs::Integer, β1, β2, β3, U_pre::Matrix, pw::Vector,
                              ε_conv::Real, maxiter::Integer; compute_performance::Bool=false,
                              performance_per_iteration::Bool=false, performance_as_dataframe::Bool=false,
                              init_pogw_uniformly::Bool=false, init_pogw_sparse::Bool=true,
                              init_pagow_uniformly::Bool=true)

    num_acts = size(U_pre,1)
    num_worldstates = size(U_pre,2)
    num_obs = cardinality_obs

    if init_pogw_uniformly && init_pogw_sparse
        warn("Cannot initialize p(o|w) uniformly and sparse at the same time - choosing uniform initialization.")
        init_pogw_sparse = false
    end

    #initialize p(o|w)
    if init_pogw_uniformly
        #init uniform
        p_ogw_init = ones(cardinality_obs, num_worldstates)
    elseif init_pogw_sparse
        #init with a sparse, diagonal pattern (maximizing H(O) - inspired by autoencoder pre-training)
        if(num_worldstates<=num_obs)
            p_ogw_init = eye(num_obs,num_worldstates)
        else
            for b in num_obs:num_obs:num_worldstates
                if b==num_obs
                    p_ogw_init = eye(num_obs, num_obs)
                else
                    p_ogw_init = [p_ogw_init eye(num_obs,b-size(p_ogw_init,2))]
                end
            end
            extra = num_worldstates - size(p_ogw_init,2)
            if extra>0
                p_ogw_init = [p_ogw_init eye(num_obs,extra)]
            end
        end
        #make sure that all elements are nonzero
        p_ogw_init += rand(size(p_ogw_init)) * 0.01  #this factor depends on the number of rows in p_ogw_init...
    else
        #init random
        p_ogw_init = rand(cardinality_obs, num_worldstates)
    end

    #initialize p(a|o,w)
    if init_pagow_uniformly
        #uniform initialization
        p_agow_init = ones(num_acts, cardinality_obs, num_worldstates)
    else
        #random initialization
        p_agow_init = rand(num_acts, cardinality_obs, num_worldstates)
    end

    #normalize
    for j in 1:num_worldstates
        p_ogw_init[:,j] /=  sum(p_ogw_init[:,j])
        for k in 1:cardinality_obs
            p_agow_init[:,:,j] /= sum(p_agow_init[:,:,j])
        end
    end

    #------- Blahut-Arimoto call --------#
    #Blahut-Arimoto iterations for the three-variable general case
    return threevarBAiterations(p_ogw_init, p_agow_init, β1, β2, β3, U_pre, pw, ε_conv, maxiter,
                                compute_performance=compute_performance,
                                performance_per_iteration=performance_per_iteration,
                                performance_as_dataframe=performance_as_dataframe)

end

#This function performs Blahut-Arimoto iterations for the three-variable general case
function threevarBAiterations(pogw_init::Matrix, pagow_init, β1, β2, β3,
                              U_pre::Matrix, pw::Vector, ε_conv::Real, maxiter::Integer;
                              compute_performance::Bool=false, performance_per_iteration::Bool=false,
                              performance_as_dataframe::Bool=false)

    card_a = size(U_pre,1)
    card_w = size(U_pre,2)
    card_o = size(pogw_init,1)

    pogw_new = pogw_init
    pagow_new = pagow_init

    #Initialize the marginals consistent with the conditionals
    po_new, pa_new, pagw = compute_marginals(pw, pogw_init, pagow_init)

    pago_new = compute_pago_iteration(pogw_init, pagow_init, β2, β3, U_pre, pa_new, po_new, pw)
    if(β3==0)
        #sequential case - make sure that D_KL( p(a|o,w)||p(a|o) ) = 0
        pagow_new = compute_pagow_iteration(pago_new, β2, β3, U_pre, pa_new)
    end

    #if performance measures don't need to be returned, don't compute them per iteration
    if compute_performance==false
        performance_per_iteration = false
    end

    #preallocate if necessary
    if performance_per_iteration
        I_ow_i = zeros(maxiter)  #I(O;W)
        I_ao_i = zeros(maxiter)  #I(A;O)
        I_awgo_i = zeros(maxiter)  #I(A;W|O)
        I_aw_i = zeros(maxiter) #I(A;W)
        Ha_i = zeros(maxiter)  #H(A)
        Ho_i = zeros(maxiter)  #H(O)
        Hago_i = zeros(maxiter)  #H(A|O)
        Hogw_i = zeros(maxiter)  #H(O|W)
        Hagow_i = zeros(maxiter)  #H(A|O,W)
        Hagw_i = zeros(maxiter) #H(A|W)
        EU_i = zeros(maxiter)
        ThreeVarRDobj_i = zeros(maxiter)
    end

    if(β3 == 0)
        sequential_case = true
    else
        sequential_case = false
    end

    #main iteration
    iter = 0 #initialize counter, so it persists beyond the loop
    for iter in 1:maxiter
        pa = pa_new
        po = po_new
        pago = pago_new
        pagow = pagow_new
        pogw = pogw_new

        #
        #figure out best or correct order of equations
        #make sure to use old distributions in curr. iteration and update all of them after every iteration
        #check the sequential case by reproducing the MEU solution
        #   perhaps it needs "layer-wise" initialization? If so, why?

        #compute p(o|w)
        pogw_new = compute_pogw_iteration(pago, pagow, β1, β2, β3, U_pre, pa, po)

        #compute p(a|o,w) and p(a|o)
        if(β3==0)
            #sequential case - compute p(a|o) first and then p(a|o,w)
            pago_new = compute_pago_iteration(pogw, pagow, β2, β3, U_pre, pa, po, pw)
            pagow_new = compute_pagow_iteration(pago_new, β2, β3, U_pre, pa)
        else
            #general case - compute p(a|o,w) first and then p(a|o)
            pagow_new = compute_pagow_iteration(pago, β2, β3, U_pre, pa)
            pago_new = compute_pago_iteration(pogw, pagow_new, β2, β3, U_pre, pa, po, pw)
        end

        #update the marginals p(o), p(a), p(a|w)
        po_new, pa_new, pagw = compute_marginals(pw, pogw_new, pagow_new)

        #compute entropic quantities (if requested with additional parameter)
        if performance_per_iteration
            I_ow_i[iter], I_ao_i[iter], I_awgo_i[iter], I_aw_i[iter],
            Ho_i[iter], Ha_i[iter], Hogw_i[iter], Hago_i[iter],
            Hagow_i[iter], Hagw_i[iter],
            EU_i[iter], ThreeVarRDobj_i[iter] = analyze_three_var_BAsolution(pw, po_new, pa_new, pogw_new, pago_new,
                                                                             pagow_new, pagw, U_pre, β1, β2, β3)
        end

        #check for convergence
        #if (norm(pa-pa_new) + norm(po-po_new)) < ε_conv
        if (norm(pagow[:]-pagow_new[:]) + norm(pogw[:]-pogw_new[:])) < ε_conv
            break
        end

    end

    #check if iteration limit has been reached (before convergence)
    if iter == maxiter
        warn("[Three variable BAiterations] maximum iteration reached - returning... (results might be inaccurate)")
    end

    #return results
    if compute_performance == false
        return po_new, pa_new, pogw_new, pago_new, pagow_new, pagw
    else
        if performance_per_iteration == false
            #compute performance measures for final solution
            I_ow, I_ao, I_awgo, I_aw, Ho, Ha, Hogw, Hago, Hagow, Hagw, EU, ThreeVarRDobj = analyze_three_var_BAsolution(pw, po_new,
                                                                                pa_new, pogw_new, pago_new, pagow_new, pagw,
                                                                                U_pre, β1, β2, β3)
        else
            #"cut" valid results from preallocated vector
            I_ow = I_ow_i[1:iter]
            I_ao = I_ao_i[1:iter]
            I_awgo = I_awgo_i[1:iter]
            I_aw = I_aw_i[1:iter]
            Ho = Ho_i[1:iter]
            Ha = Ha_i[1:iter]
            Hogw = Hogw_i[1:iter]
            Hago = Hago_i[1:iter]
            Hagow = Hagow_i[1:iter]
            Hagw = Hagw_i[1:iter]
            EU = EU_i[1:iter]
            ThreeVarRDobj = ThreeVarRDobj_i[1:iter]
        end

        #if needed, transform to data frame
        if performance_as_dataframe == false
            return po_new, pa_new, pogw_new, pago_new, pagow_new, pagw, I_ow, I_ao, I_awgo, I_aw, Ho, Ha, Hogw, Hago, Hagow, Hagw, EU, ThreeVarRDobj
        else
            performance_df = performancemeasures2DataFrame(I_ow, I_ao, I_awgo, I_aw, Ho, Ha, Hogw, Hago, Hagow, Hagw, EU, ThreeVarRDobj)
            return po_new, pa_new, pogw_new, pago_new, pagow_new, pagw, performance_df
        end
    end

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_json: true
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.11.3
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# # Gridded Domains Using `MeshArrays.jl`
#
# A `MeshArray` contains an array of subdomain arrays that (1) are connected at their edges and (2) collectively form a global grid. Grid specifications are contained in `gcmgrid` data structures. These merely define array sizes and how e.g. grid variables are represented in memory -- it is only when variables are read from file that larger memory allocations occur.

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Grid Configuration
#
# 1. import `MeshArrays` and plotting tools
# 2. choose e.g. a standard `MITgcm` grid

# + {"cell_style": "center", "slideshow": {"slide_type": "subslide"}}
using MeshArrays, Plots

pth=MeshArrays.GRID_LLC90
γ=GridSpec("LatLonCap",pth);
#Γ=GridLoad(γ)

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Read Example
#
# A `MeshArray` variable, on the chosen grid, can be accessed from file via `γ.read`(argument #1). Format conversion occurs inside the `read` function based on a propotype (argument #2). Further `read` / `write` calls convert back and forth between `MeshArray` and `Array` formats if needed.

# + {"slideshow": {"slide_type": "subslide"}}
D=γ.read(γ.path*"Depth.data",MeshArray(γ,Float64))
tmp1=write(D); tmp2=read(tmp1,D)
show(D)

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Subdomain Arrays
#
# The heatmap method is specialized for `MeshArray`s below. It operates on each `inner-array` sequentially, one after the other.

# + {"slideshow": {"slide_type": "subslide"}}
import Plots: heatmap

function heatmap(x::MeshArray; args...)
    n=x.grid.nFaces
    p=()
    for i=1:n; p=(p...,heatmap(x[i]; args...)); end
    plot(p...)
end

heatmap(D,title="Ocean Depth",clims=(0.,6000.))

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## `MeshArray` Behaves Like `Array`
#
# Here are a few examples that would be coded similarly in both cases

# +
size(D)
eltype(D)
view(D,:)

D .* 1.0
D .* D
1000*D
D*1000

D[findall(D .> 300.)] .= NaN
D[findall(D .< 1.)] .= NaN

D[1]=0.0 .+ D[1]
tmp=cos.(D);
# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## Try Another Grid
#
# The `cube-sphere` grid in `MeshArrays.GRID_CS32` has 6 subdomains, each of size `32x32`.
# -

pth=MeshArrays.GRID_CS32
γ=GridSpec("CubeSphere",pth)
D=γ.read(γ.path*"Depth.data",MeshArray(γ,Float32))
show(D)

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# ## The `exchange` Function
#
# This function adds neighboring data points (columns or rows) at face edges to slightly extend the computational domain. Doing this is often needed e.g. to compute partial derivatives between subdomains of the climate system.
# -

Dexch=exchange(D,4)
show(Dexch)

# Here is a visualization for subdomain `#6`:

# + {"slideshow": {"slide_type": "subslide"}}
P=heatmap(D.f[6],title="Ocean Depth (D, Face 6)",lims=(-4,36))
Pexch=heatmap(Dexch.f[6],title="...(Dexch, Face 6)",lims=(0,40))
plot(P,Pexch)

# + [markdown] {"slideshow": {"slide_type": "slide"}}
# # Global Diffusion Example
#
# Regular testing of `MeshArrays.jl` uses the `smooth()` function as an example. Starting from a random noise field, diffusive smoothing efficiency is predictable and can be set via a scale parameter [(see Weaver and Courtier, 2001)](https://doi.org/10.1002/qj.49712757518).
#
# This example also illustrates the generality of the approach chosen in `MeshArrays.jl`, where the same code (see `demo2`) is pretty much readily applicable to any `PeriodicDomain`, `PeriodicChannel`, `CubeSphere`, or `LatLonCap` grid. Here the chosen grid crudely represents the Earth as a unit-cube with `16*16` points on each of the 6 faces. Grid cell areas and distances are all artificially set to `1.0` for simplicity.

# + {"slideshow": {"slide_type": "subslide"}}
pth=dirname(pathof(MeshArrays))
include(joinpath(pth,"../examples/Demos.jl"))
γ,Γ=GridOfOnes("CubeSphere",6,16;option="full")
Δ=demo2(Γ);

# + [markdown] {"slideshow": {"slide_type": "subslide"}}
# The initial noise field is `D[1]` while the smoothed one is `D[2]`. After `smooth()` has been applied via `demo2()`, the noise field is visibly smoother and more muted.

# + {"cell_style": "split"}
heatmap(Δ[1],title="initial noise",clims=(-0.5,0.5))

# + {"cell_style": "split"}
heatmap(Δ[2],title="smoothed noise",clims=(-0.5,0.5))

# + [markdown] {"slideshow": {"slide_type": "subslide"}}
# The computational cost of `smooth()` predictably  increases with the decorrelation scale. For more about how this works, please refer to **Weaver and Courtier, 2001** _Correlation modelling on the sphere using a generalized diffusion equation_ https://doi.org/10.1002/qj.49712757518
# -

Rini=Δ[1]
DXCsm=Δ[3]
DYCsm=Δ[4]
@time Rend=smooth(Rini,DXCsm,DYCsm,Γ);
@time Rend=smooth(Rini,2DXCsm,2DYCsm,Γ);

max_stars_repo_name: gaelforget/Notebooks

max_issues_repo_name: gaelforget/Notebooks

max_forks_repo_name: gaelforget/Notebooks


text:
######### Stan program example  ###########

using StanSample

ProjDir = @__DIR__

bernoullimodel = "
data {
  int<lower=0> N;
  array[N] int<lower=0,upper=1> y;
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

observed_data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Default for tmpdir is to create a new tmpdir location
# To prevent recompilation of a Stan progam, choose a fixed location,
tmpdir= ProjDir * "/tmp"

sm = SampleModel("bernoulli2", bernoullimodel, tmpdir);

rc = stan_sample(sm; data=observed_data,
  save_warmup=true, num_warmups=1000,
  num_samples=1000, thin=1, delta=0.85
);

if success(rc)
  df = read_summary(sm, true)
  df |> display
end


text:
using TheAlgorithms
using Test

using LinearAlgebra

@testset "TheAlgorithms" begin

include("conversions.jl")
include("data_structures.jl")
include("knapsack.jl")
include("math.jl")
include("matrix.jl")
include("project-rosalind.jl")
include("scheduling.jl")
include("searches.jl")
include("sorts.jl")
include("statistics.jl")
include("strings.jl")

end

max_stars_repo_name: mayurdahibhate/Julia

max_issues_repo_name: mayurdahibhate/Julia

max_forks_repo_name: mayurdahibhate/Julia


text:
function alexnet()
  layers = Chain(Conv((11, 11), 3=>64, stride=(4, 4), relu, pad=(2, 2)),
                 MaxPool((3, 3), stride=(2, 2)),
                 Conv((5, 5), 64=>192, relu, pad=(2, 2)),
                 MaxPool((3, 3), stride=(2, 2)),
                 Conv((3, 3), 192=>384, relu, pad=(1, 1)),
                 Conv((3, 3), 384=>256, relu, pad=(1, 1)),
                 Conv((3, 3), 256=>256, relu, pad=(1, 1)),
                 MaxPool((3, 3), stride=(2, 2)),
                 AdaptiveMeanPool((6,6)),
                 flatten,
                 Dropout(0.5),
                 Dense(256 * 6 * 6, 4096, relu),
                 Dropout(0.5),
                 Dense(4096, 4096, relu),
                 Dense(4096, 1000))
  Flux.testmode!(layers, false)
  return layers
end


text:
# This file is a part of SolidStateDetectors.jl, licensed under the MIT License (MIT).
using SolidStateDetectors
using Test

T = Float32

@testset "Package SolidStateDetectors" begin

    @testset "Simulate example detector: Inverted Coax" begin
        sim = Simulation(SSD_examples[:InvertedCoax])
        simulate!(sim, max_refinements = 1, verbose = true)
        evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 40e-3 )]))
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end
    @testset "Simulate example detector: Coax" begin
        sim = Simulation(SSD_examples[:Coax])
        simulate!(sim, max_refinements = 0, verbose = true)
        evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(30), 12e-3 )]))
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-2 )
    end
    @testset "Simulate example detector: BEGe" begin
        sim = Simulation(SSD_examples[:BEGe])
        simulate!(sim, max_refinements = 1, verbose = true)
        evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 20e-3 )]))
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 1e-2 )
    end
    @testset "Simulate example detector: CGD" begin
        sim = Simulation(SSD_examples[:CGD])
        simulate!(sim, max_refinements = 1, verbose = true)
        evt = Event([CartesianPoint{T}(5e-3, 5e-3, 5e-3)])
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 1e-2 )
    end
    @testset "Simulate example detector: Spherical" begin
        sim = Simulation(SSD_examples[:Spherical])
        simulate!(sim, max_refinements = 2, verbose = true)
        evt = Event([CartesianPoint{T}(0,0,0)])
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 1e-2 )
    end
    @testset "Simulate example detector: SigGen PPC" begin
        sim = Simulation(SSD_examples[:SigGen])
        simulate!(sim, max_refinements = 1, verbose = true)
        evt = Event(CartesianPoint.([CylindricalPoint{T}(20e-3, deg2rad(10), 40e-3 )]))
        simulate!(evt, sim)
        signalsum = T(0)
        for i in 1:length(evt.signals)
            signalsum += abs(evt.signals[i][end])
        end
        @test isapprox( signalsum, T(2), atol = 5e-3 )
    end

end


text:
import Dance
import DanceWebpack
import HTTP

include("utils.jl")

# Dance by default does not automatically add Project.toml
# => test here with no values specified
Dance.start_project("demo")
cd("demo")
project_settings()
Dance.populate_settings(joinpath(abspath(@__DIR__), "demo"))
@test_logs (:info, "Webpack files added to `static` dir in project.\nPlease cd into `static` dir and add dependencies by installing NodeJS and running `npm install`") DanceWebpack.setup()

# Todo: fix for TravisCI Windows
if !Sys.iswindows()
    cd("static")
    run(`npm install`)
    static_files_append_content()

    #= Test Prod mode content is properly served =#
    run(`npm run build`)
    bundle_id = get_webpack_bundle_id()
    Dance.Router.delete_routes!()

    @testset "HTTP.listen" begin
        @async include(joinpath(abspath(pwd()), "dance.jl"))
        sleep(1)

        r = HTTP.request("GET", "http://127.0.0.1:8000/")
        @test r.status==200
        compare_http_header(r.headers, "content_type", "text/html; charset=UTF-8")
        @test extract_html_body_content(r.body)=="Hello World"

        r = HTTP.request("GET", "http://127.0.0.1:8000/static/main.$bundle_id.js")
        @test r.status==200
        compare_http_header(r.headers, "content_type", "text/javascript")

        r = HTTP.request("GET", "http://127.0.0.1:8000/static/vendors~main.$bundle_id.js")
        @test r.status==200
        compare_http_header(r.headers, "content_type", "text/javascript")

        r = HTTP.request("GET", "http://127.0.0.1:8000/static/main.$bundle_id.css")
        @test r.status==200
        compare_http_header(r.headers, "content_type", "text/css")
    end

    cd("static")

    #= Test over-writting `base.html` with Dev mode =#
    @testset "npm run develop" begin
        @async run(`npm run develop`)
        sleep(10)
        html_file = read("../html/base.html", String)
        @test occursin("<script src=\"http://localhost:3000/main.js\"></script>", html_file)
    end

    #= Test over-writting `base.html` with Prod mode =#
    run(`npm run build`)
    html_file = read("../html/base.html", String)
    @test occursin("<link rel=\"stylesheet\" href=\"/static/main.", html_file)
    @test occursin("<script src=\"/static/main.", html_file)
    @test occursin("<script src=\"/static/vendors~main.", html_file)

    delete_project()
end


text:
"""
Linear Algebra related graph functions.
Especially, methods for handling incidence matrixes [`IncMat`](@ref).
"""
module LinA
import SparseArrays: SparseMatrixCSC, sparse

const IncMat = SparseMatrixCSC{Float64,Int64}
const Edges = Vector{Tuple{Int, Int}}

"""
    is_incmat(::IncMat)

Check whether a matrix is a graph (oriented) edge-vertex incidence matrix.
"""
function is_incmat(d::IncMat)::Bool
    dt = IncMat(d')
    any(diff(diff(dt.colptr)) .!= 0) && return false
    any(dt.nzval[1:2:end] + dt.nzval[2:2:end] .!= 0) && return false
    return true
end

function edges_from_incmat(d::IncMat)
    m, n = size(d)
    head = Vector{Int}(undef, m)
    tail = Vector{Int}(undef, m)
    lam = Vector{Float64}(undef, m)
    k = Int(1)
    for i = 1:n
        for ci = d.colptr[i]:d.colptr[i+1]-1
            ei = d.rowval[ci]
            v = d.nzval[ci]
            lam[ei] = abs(v)
            if v > 0
                head[ei] = i
            else
                tail[ei] = i
            end
        end
    end
    return head, tail, lam
end

"""
    IncMat(edges, λ, n [, T = Int])

Create an incidence matrix from a list of edges, weights `λ` having `n` nodes.
The resulting matrix index will have type `T`.
"""
function IncMat(E::Edges,
                lam::Vector{Float64},
                n::Int,
                T = Int)::IncMat
    m = length(E)
    I = Vector{T}(undef, 2m)
    J = Vector{T}(undef, 2m)
    V = Vector{eltype(lam)}(undef, 2m)
    for (i, (u, v)) in enumerate(E)
        I[2i-1] = i
        J[2i-1] = u
        V[2i-1] = lam[i]
        I[2i-0] = i
        J[2i-0] = v
        V[2i-0] = -lam[i]
    end
    sparse(I, J, V, m, n)
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra

# ## Raibert Policy
include("policy/3D.jl")

# ## Simulation
s = get_simulation("hopper_3D", "flat_3D_lc", "flat")
model_sim = s.model
nq = model_sim.nq
nu = model_sim.nu
nc = model_sim.nc
nw = model_sim.nw

# ## Setup

using ContactImplicitMPCH = 92
h = 0.01
N_sample = 5
h_sim = h / N_sample
H_sim = 5000

# ## Policy
v0 = [0.0; 0.2]
Tstance = 0.13 # measure using hop-in-place gait
Tflight = 0.62 # measure using hop-in-place gait
p = raibert_policy(model_sim, v0=v0, Tstance=Tstance, Tflight=Tflight, h=h);

# ## Initial conditions
off0 = ContactImplicitMPC.SVector{nq}([0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5])
off1 = ContactImplicitMPC.SVector{nq}([v0[1] * h_sim, v0[2] * h_sim, 0.5, 0.0, 0.0, 0.0, 0.5])
q_ref = ContactImplicitMPC.SVector{nq}([0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0])
q0_sim = copy(q_ref) + off0
q1_sim = copy(q_ref) + off1

# ## Simulator
sim = simulator(model_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    ip_opts = InteriorPointOptions(
        γ_reg = 0.0,
        undercut = Inf,
        r_tol = 1.0e-8,
        κ_tol = 1.0e-8,),
    sim_opts = SimulatorOptions(warmstart = true)
    );

# ## Simulate
@time status = simulate!(sim);

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

# ## Visualize
visualize_robot!(vis, model_sim, sim.traj, sample=20);


text:
"""
    RangeMatrix{T<:Range}(rs::AbstractVector{T})

A RangeMatrix is a simple matrix representation of a vector of ranges, with
each range representing one column. Construct a RangeMatrix with a vector of
ranges; the ranges must all have the same length.

Note that it is only efficient when all the ranges are of the same type and in
a concretely typed Vector.
"""
immutable RangeMatrix{T,A} <: AbstractArray{T,2}
    rs::A # A <: AbstractVector{_<:Range{T}}
    dims::Tuple{Int,Int}
end
RangeMatrix(rs::Range...) = RangeMatrix(collect(rs)) # TODO: use tuple storage?
function RangeMatrix{T<:Range}(rs::AbstractVector{T})
    n = length(rs)
    n == 0 && return RangeMatrix{T}(rs, (0, 0))
    m = length(rs[1])
    for j=2:n
        m == length(rs[j]) || throw(ArgumentError("all ranges must have the same length; expected $m, got $(length(rs[j]))"))
    end
    RangeMatrix{eltype(T), typeof(rs)}(rs, (m, n))
end

Base.size(R::RangeMatrix) = R.dims
Base.linearindexing{R<:RangeMatrix}(::Type{R}) = Base.LinearSlow()

# Scalar indexing
Base.getindex(R::RangeMatrix, i::Int, j::Int) = (checkbounds(R, i, j); Base.unsafe_getindex(R, i, j))
Base.unsafe_getindex(R::RangeMatrix, i::Int, j::Int) = @inbounds return R.rs[j][i]

# For non-scalar indexing, only specialize with inner Ranges and Colons to
# return Ranges or RangeMatrixes. For everything else, we can use the fallbacks.
Base.getindex(R::RangeMatrix, I::Union{Range, Colon}, J) = (checkbounds(R, I, J); Base.unsafe_getindex(R, I, J))
Base.unsafe_getindex(R::RangeMatrix, I::Union{Range, Colon}, j::Real) = @inbounds return R.rs[j][I]
Base.unsafe_getindex(R::RangeMatrix, I::Union{Range, Colon}, ::Colon) = @inbounds return RangeMatrix([R.rs[j][I] for j=1:length(R.rs)])
Base.unsafe_getindex(R::RangeMatrix, I::Union{Range, Colon}, J)       = @inbounds return RangeMatrix([R.rs[j][I] for j in J])


text:
function get_new_best(results::Array{RemoteChannel}, best::Result)
    for channel in results
        try
            partial = take!(channel)
            if partial.cost_minimum < best.cost_minimum
                best = deepcopy(partial)
            end
        catch
        end
    end
    best
end

function get_new_best(results::Array{RemoteChannel})
    try
        take!(results[rand(1:length(results))])
    catch
    end
end


text:
struct NDLCDRefSpace{T} <: RefSpace{T,4} end

function valuetype(ref::NDLCDRefSpace{T}, charttype::Type) where {T}
    SVector{universedimension(charttype),T}
end

function (ϕ::NDLCDRefSpace)(ndlc)

    u,v,w = parametric(ndlc)
    j = jacobian(ndlc)

    tu = tangents(ndlc,1)
    tv = tangents(ndlc,2)
    tw = tangents(ndlc,3)

    # B = [tu tv tw]

    # return SVector((
    #     (value=(2*B*[(u-1),v,w]/j),divergence=(6/j)),
    #     (value=(2*B*[u,(v-1),w]/j),divergence=(6/j)),
    #     (value=(2*B*[u,v,(w-1)]/j),divergence=(6/j)),
    #     (value=(2*B*[u,v,w]/j)    ,divergence=(6/j))
    # ))

    return SVector((
        (value=(2*((u-1)*tu + v*tv + w*tw)/j), divergence=(6/j)),
        (value=(2*(u*tu + (v-1)*tv + w*tw)/j), divergence=(6/j)),
        (value=(2*(u*tu + v*tv + (w-1)*tw)/j),divergence=(6/j)),
        (value=(2*(u*tu + v*tv + w*tw)/j)    ,divergence=(6/j))
    ))
end

function ntrace(x::NDLCDRefSpace, el, q, fc)
    t = zeros(scalartype(x),1,4)
    t[q] = 1 / volume(fc)
    return t
end

"""
Does not give the correct result for an imput basis with non-vanishing
input basis on the boundary
"""
function ttrace(x::NDLCDRefSpace, el, q, fc)
    # el: the supporting chart of the input
    # q: the local index of the face of el on which we compute the trace
    # fc: the chart on which we compute the trace
    t = zeros(scalartype(x),3,4)
    fa = SArray{Tuple{2},Int,1,2}
    te = SArray{Tuple{3},Int,1,3}
    j = 0
    y = BEAST.NDRefSpace{scalartype(x)}()
    n = normal(fc)
    for fac in faces(el)
        j += 1
        j == q && continue
        # fac != fc || continue
        #vind locale index (i) in fc van de edge die in zowel fc als fac ligt
        A = setdiff(fc.vertices,setdiff(fc.vertices,fac.vertices))
        # @show length(A)
        @assert length(A) == 2
        edg = fa([findfirst(isequal(A[1]), fc.vertices),findfirst(isequal(A[2]), fc.vertices)])
        i = abs(CompScienceMeshes.relorientation(edg,te([1,2,3])))
        #i = findfirst(isequal(setdiff(fcv,facv)[1]), fcv)
        edg_chart = simplex(fc.vertices[edg[1]], fc.vertices[edg[2]])
        edg_centr = cartesian(center(edg_chart))
        p = neighborhood(fc, carttobary(fc, edg_centr))
        r = neighborhood(el, carttobary(el, edg_centr))
        yp = y(p)[i].value
        xr = x(r)[j].value
        tgt = edg_chart.vertices[1] - edg_chart.vertices[2]
        ypt = dot(yp, tgt)
        xrt = dot(n × xr, tgt)
        # t[i,j] = volume(el)/(volume(fac)*norm(A[1]-A[2]))
        t[i,j] = xrt / ypt
    end
    return t
end

divergence(ref::NDLCDRefSpace, sh, el) = Shape(sh.cellid, 1, sh.coeff/volume(el))

function restrict(ϕ::NDLCDRefSpace{T}, dom1, dom2) where {T}
    # dom2 is the smaller of the domains

    K = numfunctions(ϕ)
    D = dimension(dom1)

    @assert K == 4
    @assert D == 3
    @assert D == dimension(dom2)

    Q = zeros(T,K,K)
    for (i,face) in enumerate(faces(dom2))

        p = center(face)
        c = cartesian(p)
        A = volume(face)

        m = normal(p)
        u = carttobary(dom1,p)

        u = carttobary(dom1, c)
        x = neighborhood(dom1, u)

        y = ϕ(x)

        for j in 1:K
            Q[j,i] = dot(y[j].value, m) * A
        end
    end

    return Q
end


text:
rmlines = @λ begin
    e :: Expr           -> Expr(e.head, map(rmlines, filter(x -> !(x isa LineNumberNode), e.args))...)
    a                   -> a
end
expr = quote
    struct S{T}
        a :: Int
        b :: T
    end
end |> rmlines

@test @match expr begin
    quote
        struct $name{$tvar}
            $f1 :: $t1
            $f2 :: $t2
        end
    end =>
    quote
        struct $name{$tvar}
            $f1 :: $t1
            $f2 :: $t2
        end
    end |> rmlines == expr
end


text:
module DynNets

using Base: FlatteningRF, Bool
using Distributions
using StatsBase
using Optim
using LineSearches
using StatsFuns
using Statistics
using Roots
using MLBase
using GLM
using LinearAlgebra
using JLD2
using DataFrames
using ForwardDiff
using NLSolversBase
using RCall
using PyPlot
using Logging
using Distributed
using SharedArrays
using DataStructures

using ..AReg
using ..Utilities
using ..Scalings
using ..ErgmRcall
using ..StaticNets

#constants
const targetErrValDynNets = 0.01

include("./DynNets_models/DynNets_Abstract_SdErgm.jl")

include("./DynNets_models/DynNets_Abstract_SdErgm_utilities.jl")

include("./DynNets_models/DynNets_Abstract_SdErgm_SD__based_conf_Intervals.jl")

include("./DynNets_models/DynNets_SdErgmDirBin1.jl")

include("./DynNets_models/DynNets_SdErgmPml.jl")

include("./DynNets_models/DynNets_SdErgmDirBin0Rec0.jl")

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module MeshViz

using Meshes
using Tables

using CategoricalArrays: CategoricalValue
using CategoricalArrays: levelcode

using Colors: distinguishable_colors

import Makie

"""
    viz(object)

Visualize Meshes.jl `object` with various options:

* `size`          - size of points in point set
* `color`         - color of geometries or points
* `colormap`      - colormap used for colors
* `boundarycolor` - color of the boundary (e.g. segments)
* `facetcolor`    - color of the facets (e.g. edges)
* `showboundary`  - tells whether or not to show the boundary
* `showfacets`    - tells whether or not to show the facets
* `variable`      - informs which variable to visualize
* `decimation`    - decimation tolerance for polygons

The option `color` can be a single scalar or a vector
of scalars. For meshes, the length of the vector of
colors determines if the colors should be assigned to
vertices or elements.

## Examples

```
# vertex coloring (i.e. linear interpolation)
viz(mesh, color = 1:nvertices(mesh))

# element coloring (i.e. discrete colors)
viz(mesh, color = 1:nelements(mesh))
```
"""
@Makie.recipe(Viz, object) do scene
  Makie.Attributes(;
    size          = Makie.theme(scene, :markersize),
    color         = :slategray3,
    colormap      = :viridis,
    boundarycolor = :gray30,
    facetcolor    = :gray30,
    showboundary  = true,
    showfacets    = false,
    variable      = nothing,
    decimation    = 0.0,
  )
end

# domain
include("simplemesh.jl")
include("cartesiangrid.jl")
include("collection.jl")
include("partition.jl")
include("fallbacks.jl")

# data
include("data.jl")

end


text:
using Pkg
Pkg.install(["LinearAlgebra", "Plots", "PlotThemes", "ScikitLearn"])
using LinearAlgebra

m = [1 2 3 4; 4 3 2 1; 1 1 1 1; 0 0 0 0]
transpose(m)

max_stars_repo_path: src/website/test data/julia.jl

max_issues_repo_path: src/website/test data/julia.jl

max_forks_repo_path: src/website/test data/julia.jl