text:
"""
    $(TYPEDEF)
An abstract type for all VLBI interfermetry data types. See [EHTVisibilityDatum](@ref) for an example.
"""
abstract type AbstractInterferometryDatum{T} end

abstract type AbstractVisibilityDatum{T} <: AbstractInterferometryDatum{T} end
abstract type AbstractLinearPolDatum{S<:AbstractVisibilityDatum, T} <: AbstractInterferometryDatum{T} end
abstract type AbstractCrossPolDatum{S,T} <: AbstractInterferometryDatum{T} end

abstract type ClosureProducts{T} <: AbstractInterferometryDatum{T} end

abstract type Observation{T} end

using DelimitedFiles
using AstroTime: modified_julian

export uvpositions, stations, getdata, arrayconfig,
       getuv, baselines, rescaleuv!

"""
    $(TYPEDEF)
This defined the abstract type for an array configuration. Namely, baseline
times, SEFD's, bandwidth, observation frequencies, etc.
"""
abstract type ArrayConfiguration end

"""
    $(TYPEDEF)
Stores all the non-visibility data products for an EHT array. This is useful when evaluating
model visibilities.

# $(FIELDS)
"""
struct EHTArrayConfiguration{S,F,T<:AbstractArray} <: ArrayConfiguration
    """
    List of stations in the array
    """
    stations::S
    """
    Observing frequency (Hz)
    """
    frequency::F
    """
    Observing bandwith (Hz)
    """
    bandwidth::F
    """
    A struct array of `ArrayBaselineDatum` holding time, u, v, baselines.
    """
    uvsamples::T
end

"""
    $(SIGNATURES)
Get the u, v positions of the array.
"""
function getuv(ac::ArrayConfiguration)
    return ac.uvsamples.u, ac.uvsamples.v
end

"""
    $(SIGNATURES)
Get the u, v, time, freq of the array as a tuple.
"""
function uvtimefreq(ac::EHTArrayConfiguration)
    u,v = getuv(ac)
    t = ac.uvsamples.time
    ν = ac.frequency
    return u, v, t, fill(ν, length(u))
end

"""
    $(TYPEDEF)
A single datum of an observing array.
"""
struct ArrayBaselineDatum{T}
    time::T
    u::T
    v::T
    baseline::Tuple{Symbol, Symbol}
end

const ArrayQuadrangleDatum = NTuple{4, ArrayBaselineDatum{T}} where {T}
const ArrayTriangleDatum = NTuple{3, ArrayBaselineDatum{T}} where {T}

"""
    uvpositions(datum)
Get the uvp positions of an inferometric datum.
"""
uvpositions(D::AbstractVisibilityDatum) = D.u, D.v

"""
    $(SIGNATURES)
The main data product type in `Comrade` this stores the `data` which can be a StructArray
of any `AbstractInterferometryDatum` type.

# ($FIELDS)
"""
Base.@kwdef struct EHTObservation{F,T<:AbstractInterferometryDatum{F},S<:StructArray{T}, N} <: Observation{F}
    """
    StructArray of data productts
    """
    data::S
    """
    modified julia date of the observation
    """
    mjd::N
    """
    RA of the observation in J2000 (deg)
    """
    ra::F
    """
    DEC of the observation in J2000 (deg)
    """
    dec::F
    """
    bandwidth of the observation (Hz)
    """
    bandwidth::F
    """
    frequency of the observation (Hz)
    """
    frequency::F
    """
    Common source name
    """
    source::Symbol
    """
    Time zone used.
    """
    timetype::Symbol = :UTC
end

Base.getindex(data::EHTObservation, i::Int) = data.data[i]
Base.length(data::EHTObservation) = length(data.data)

"""
    $(SIGNATURES)
Get all the stations in a observation. The result is a vector of symbols.
"""
function stations(d::EHTObservation{T,A}) where {T,A<:AbstractInterferometryDatum}
    bl = getdata(d, :baseline)
    s1 = first.(bl)
    s2 = last.(bl)
    return unique([s1..., s2...])
end

function Base.show(io::IO, d::EHTObservation{F,D}) where {F,D}
    println(io, "EHTObservation{$F,$D, ...}")
    println(io, "  source: ", d.source)
    println(io, "  mjd: ", d.mjd)
    println(io, "  frequency: ", d.frequency)
    println(io, "  bandwidth: ", d.bandwidth)
    println(io, "  stations: ", stations(d))
    println(io, "  nsamples: ", length(d))
end

"""
    getdata(obs::EHTObservation, s::Symbol)
Pass-through function that gets the array of `s` from the EHTObservation. For example
say you want the times of all measurement then

```julia
getdata(obs, :time)
```
"""
getdata(obs::Observation, s::Symbol) = getproperty(obs.data, s)
Base.getindex(obs::Observation, i) = getindex(obs.data, i)

"""
    $(SIGNATURES)
A struct holding the information for a single measured visibility.

# $(FIELDS)

"""
Base.@kwdef struct EHTVisibilityDatum{T<:Number} <: AbstractVisibilityDatum{T}
    """
    real component of the visibility (Jy)
    """
    visr::T
    """
    imaginary component of the visibility (Jy)
    """
    visi::T
    """
    error of the visibility (Jy)
    """
    error::T
    """
    x-direction baseline length in λ
    """
    u::T
    """
    y-direction baseline length in λ
    """
    v::T
    """
    Time of the observation in hours
    """
    time::T
    """
    station baseline codes
    """
    baseline::NTuple{2,Symbol}
end

"""
    $(SIGNATURES)
Return the complex visibility of the visibility datum
"""
@inline function visibility(D::EHTVisibilityDatum{T}) where {T}
        return Complex{T}(D.visr, D.visi)
end

"""
    $(SIGNATURES)
Get the amplitude of a visibility datum
"""
function amplitude(D::EHTVisibilityDatum)
    amp = hypot(D.visr,D.visi)
    return EHTVisibilityAmplitudeDatum(amp, D.error,
                                       D.u, D.v,
                                       D.time,
                                       D.frequency,
                                       D.bandwidth,
                                       D.baseline
                                    )
end

"""
    $(SIGNATURES)
A struct holding the information for a single measured visibility amplitude.

# $(FIELDS)

"""
Base.@kwdef struct EHTVisibilityAmplitudeDatum{T<:Number} <: AbstractVisibilityDatum{T}
    """
    amplitude (Jy)
    """
    amp::T
    """
    error of the visibility amplitude (Jy)
    """
    error::T
    """
    x-direction baseline length in λ
    """
    u::T
    """
    y-direction baseline length in λ
    """
    v::T
    """
    Time of the observation in hours
    """
    time::T
    """
    station baseline codes
    """
    baseline::NTuple{2,Symbol}
end

function checktriangle(D1::EHTVisibilityDatum,
                       D2::EHTVisibilityDatum,
                       D3::EHTVisibilityDatum)
    b1 = D1.baseline
    b2 = D2.baseline
    b3 = D3.baseline
    l = length(unique([b1..., b2..., b3...]))
    @assert l == 3 "For a valid closure phase you need 3 unique stations not $l"
    @assert (D1.time == D2.time == D3.time) "For a valid closure phase the times need to match"

end

@inline function amplitude(D::EHTVisibilityAmplitudeDatum{T}) where {T}
    return D.amp
end

@inline function bispectrum(D1::EHTVisibilityDatum, D2::EHTVisibilityDatum, D3::EHTVisibilityDatum)
    visibility(D1)*visibility(D2)*visibility(D3)
end

"""
    rescaleuv!(data)
rescale the u-v lengths according to scale. This can be useful when you want the spatial scales
to be in 1/μas instead of 1/rad
"""
function rescaleuv!(data::EHTObservation{T,D}, scale) where {T,D<:Union{EHTVisibilityAmplitudeDatum,EHTVisibilityDatum}}
    data.data.u .= data.data.u.*scale
    data.data.v .= data.data.v.*scale
    return data
end

"""
    $(TYPEDEF)

A Datum for a single closure phase.

# $(FIELDS)

"""
Base.@kwdef struct EHTClosurePhaseDatum{T<:Number} <: ClosureProducts{T}
    """
    closure phase (rad)
    """
    phase::T
    """
    error of the closure phase assuming the high-snr limit
    """
    error::T
    """
    u (λ) of first station
    """
    u1::T
    """
    v (λ) of first station
    """
    v1::T
    """
    u (λ) of second station
    """
    u2::T
    """
    v (λ) of second station
    """
    v2::T
    """
    u (λ) of third station
    """
    u3::T
    """
    v (λ) of third station
    """
    v3::T
    """
    Measured time of closure phase in hours
    """
    time::T
    """
    station baselines used
    """
    triangle::NTuple{3,Symbol}
end

function rescaleuv!(data::EHTObservation{T,D}, scale) where {T,D<:EHTClosurePhaseDatum}
    data.data.u1 .= data.data.u1.*scale
    data.data.v1 .= data.data.v1.*scale
    data.data.u2 .= data.data.u2.*scale
    data.data.v2 .= data.data.v2.*scale
    data.data.u3 .= data.data.u3.*scale
    data.data.v3 .= data.data.v3.*scale
    return data
end

function stations(d::EHTObservation{T,A}) where {T,A<:EHTClosurePhaseDatum}
    bl = getdata(d, :triangle)
    return unique(vcat(collect.(bl)...))
end

"""
    $(SIGNATURES)
Computes the closure phase of the three visibility datums.

# Notes
We currently use the high SNR Gaussian error approximation for the closure phase.
In the future we may use the moment matching from Monte Carlo sampling.
"""
function closure_phase(D1::EHTVisibilityDatum,
                      D2::EHTVisibilityDatum,
                      D3::EHTVisibilityDatum)

    checktriangle(D1,D2,D3)

    amp1 = amplitude(D1).amp
    amp2 = amplitude(D2).amp
    amp3 = amplitude(D3).amp
    u1,v1 = uvpositions(D1)
    u2,v2 = uvpositions(D2)
    u3,v3 = uvpositions(D3)

    bis = bispectrum(D1, D2, D3)
    s12 = unique([D1.baseline..., D2.baseline...])
    s123 = unique([s12..., D3.baseline...])
    #Use the Gaussian approximation TODO hook this into Measurements.jl?
    err = sqrt((D1.error/amp1)^2 + (D2.error/amp2)^2 + (D3.error/amp3)^2)
    return EHTClosurePhaseDatum(angle(bis), err,
                                u1, v1, u2, v2, u3, v3,
                                time, s123)
end

function baselines(CP::EHTClosurePhaseDatum)
    tri = CP.triangle
    return ((tri[1],tri[2]), (tri[2], tri[3]), (tri[3], tri[1]))
end

uvpositions(datum::EHTClosurePhaseDatum) = (datum.u1, datum.v1, datum.u2, datum.v2, datum.u3, datum.v3)

"""
    $(TYPEDEF)

A Datum for a single log closure amplitude.

# $(FIELDS)

"""
Base.@kwdef struct EHTLogClosureAmplitudeDatum{T<:Number} <: ClosureProducts{T}
    """
    log-closure amplitude
    """
    amp::T
    """
    log-closure amplitude error in the high-snr limit
    """
    error::T
    """
    u (λ) of first station
    """
    u1::T
    """
    v (λ) of first station
    """
    v1::T
    """
    u (λ) of second station
    """
    u2::T
    """
    v (λ) of second station
    """
    v2::T
    """
    u (λ) of third station
    """
    u3::T
    """
    v (λ) of third station
    """
    v3::T
    """
    u (λ) of fourth station
    """
    u4::T
    """
    v (λ) of fourth station
    """
    v4::T
    """
    Observation time of the quadrangle
    """
    time::T
    """
    station codes for the quadrangle
    """
    quadrangle::NTuple{4,Symbol}
end

function rescaleuv!(data::EHTObservation{T,D}, scale) where {T,D<:EHTLogClosureAmplitudeDatum}
    data.data.u1 .= data.data.u1.*scale
    data.data.v1 .= data.data.v1.*scale
    data.data.u2 .= data.data.u2.*scale
    data.data.v2 .= data.data.v2.*scale
    data.data.u3 .= data.data.u3.*scale
    data.data.v3 .= data.data.v3.*scale
    data.data.u4 .= data.data.u4.*scale
    data.data.v4 .= data.data.v4.*scale
    return data
end

function baselines(CP::EHTLogClosureAmplitudeDatum)
    quad = CP.quadrangle
    return ((quad[1],quad[2]), (quad[3], quad[4]), (quad[1], quad[3]), (quad[2], quad[4]))
end

function stations(d::EHTObservation{T,A}) where {T,A<:EHTLogClosureAmplitudeDatum}
    bl = getdata(d, :quadrangle)
    return unique(vcat(collect.(bl)...))
end

uvpositions(datum::EHTLogClosureAmplitudeDatum) = (datum.u1, datum.v1, datum.u2, datum.v2, datum.u3, datum.v3, datum.u4, datum.v4)

# @doc raw"""
#     $(SIGNATURES)
# Compute the log-closure amplitude of 4 visibility datums.
# This uses the same ordering as eht-imaging i.e.

# ```math
# C_A = \log \frac{|V_{1}| |V_{2}|}{|V_{3}| |V_{4}|}
# ```

# # Notes

# Currently we use the high SNR Gaussian approximation to the thermal error.
# In the future we may use the moment matching approach using Monte Carlo sampling.
# """
# function logclosure_amplitude(D1::EHTVisibilityDatum,
#                               D2::EHTVisibilityDatum,
#                               D3::EHTVisibilityDatum,
#                               D4::EHTVisibilityDatum)
#     check_quadrangle(D1,D2,D3,D4)

#     amp1 = amplitude(D1)
#     amp2 = amplitude(D2)
#     amp3 = amplitude(D3)
#     amp4 = amplitude(D4)

#     u1,v1 = uvpositions(D1)
#     u2,v2 = uvpositions(D2)
#     u3,v3 = uvpositions(D3)
#     u4,v4 = uvpositions(D4)

#     #Construct the quadrangle
#     s12 = unique([D1.baseline..., D2.baseline...])
#     s123 = unique([s12..., D3.baseline...])
#     s1234 = unique([s123..., D4.baseline...])

#     # Now do the error propogation
#     err = hypot(D1.err/map1, D2.err/amp2, D3.err/amp3, D4.err/amp4)

#     #construct the lcamp
#     lcamp = log(amp1*amp2/(amp3*amp4))

#     return EHTLogClosureAmplitudeDatum(lcamp, err,
#                                        u1, v1,
#                                        u2, v2,
#                                        u3, v3,
#                                        u4, v4,
#                                        time, s1234)
# end

"""
    $(SIGNATURES)
Extract the array configuration from a visibility EHT observation.
"""
function arrayconfig(vis::EHTObservation{F,A}) where {F,A<:Union{EHTVisibilityDatum, EHTVisibilityAmplitudeDatum}}
    u = getdata(vis, :u)
    v = getdata(vis, :v)
    st = stations(vis)
    times = getdata(vis, :time)
    bandwidth = vis.bandwidth
    frequency = vis.frequency
    error = getdata(vis, :error)
    baseline = getdata(vis, :baseline)
    uvsamples = StructArray{ArrayBaselineDatum}(time=times,
                                        u=u,
                                        v=v,
                                        baseline=baseline,
                                        error_real=error,
                                        error_imag=error
                                    )
    return EHTArrayConfiguration(st, frequency, bandwidth, uvsamples)
end

include(joinpath(@__DIR__, "io.jl"))
include(joinpath(@__DIR__, "ehtim.jl"))


text:
using COESA
using Test

M0 = COESA.M0

# altitude conversions
@test isapprox(round(COESA.geopotential_altitude(0)), 0)
@test isapprox(round(COESA.geopotential_altitude(-5000)), -5004)
@test isapprox(round(COESA.geopotential_altitude(85_500)), 84_365)
@test isapprox(round(COESA.geopotential_altitude(1_000_000)), 864_071)

# out of bounds errors
@test_throws ErrorException atmosphere(prevfloat(-5_000.0))
@test_throws ErrorException atmosphere(nextfloat(1_000_000.0))

# b = 0, Z < 0
Z = -5_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 320.676)
@test isapprox(round(pressure(atmos), sigdigits=5), 1.7776e5)
@test isapprox(round(density(atmos), sigdigits=5), 1.9311)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 358.99)

# b = 0, Z = 0
Z = 0

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 288.150)
@test isapprox(round(pressure(atmos), sigdigits=6), 101325)
@test isapprox(round(density(atmos), sigdigits=5), 1.2250)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 340.29)

# b = 0, Z > 0
Z = 5_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 255.676)
@test isapprox(round(pressure(atmos), sigdigits=5), 5.4048e4)
@test isapprox(round(density(atmos), sigdigits=5), 7.3643e-1)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 320.55)

# b = 1
Z = 15_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 216.650)
@test isapprox(round(pressure(atmos), sigdigits=4), round(1.2111e4, sigdigits=4))
@test isapprox(round(density(atmos), sigdigits=5), 1.9476e-1)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 295.07)

# b = 2
Z = 25_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 221.552)
@test isapprox(round(pressure(atmos), sigdigits=5), 2.5492e3)
@test isapprox(round(density(atmos), sigdigits=5), 4.0084e-2)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 298.39)

# b = 3
Z = 40_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 250.350)
@test isapprox(round(pressure(atmos), sigdigits=5), 2.8714e2)
@test isapprox(round(density(atmos), sigdigits=5), 3.9957e-3)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 317.19)

# b = 4
Z = 50_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 270.650)
@test isapprox(round(pressure(atmos), sigdigits=5), 7.9779e1)
@test isapprox(round(density(atmos), sigdigits=5), 1.0269e-3)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 329.80)

# b = 5
Z = 60_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 247.021)
@test isapprox(round(pressure(atmos), sigdigits=4), round(2.1958e1, sigdigits=4))
@test isapprox(round(density(atmos), sigdigits=5), 3.0968e-4)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 315.07)

# b = 6
Z = 75_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=5), 28.964)
@test isapprox(round(temperature(atmos), sigdigits=6), 208.399)
@test isapprox(round(pressure(atmos), sigdigits=5), 2.3881)
@test isapprox(round(density(atmos), sigdigits=5), 3.9921e-5)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 289.40)

# b = 7
# The tables show values of M that have a discontinuity at Z = 86km, but our
# routine applies the suggested blending
Z = 85_000

atmos = atmosphere(Z)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=6), round(M0 * 0.999694, sigdigits=6))
@test isapprox(round(temperature(atmos), sigdigits=6), round(188.893 / M0 * mean_molecular_weight(atmos), sigdigits=6))
@test isapprox(round(pressure(atmos), sigdigits=5), 4.4568e-1)
@test isapprox(round(density(atmos), sigdigits=4), round(8.2196e-6, sigdigits=4))
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 275.52)

# b = 8, on a breakpoint
Z = 86_000

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=5), 186.87)
@test isapprox(round(pressure(atmos), sigdigits=5), 3.7338e-1)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=4), 28.95)
@test isapprox(round(density(atmos), sigdigits=3), round(6.958e-6, sigdigits=3))
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

# b = 8, in between breakpoints
Z = 86_500

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=5), 186.87)
@test isapprox(round(pressure(atmos), sigdigits=5), 3.4163e-1)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=4), 28.95)
@test isapprox(round(density(atmos), sigdigits=4), 6.366e-6)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

# b = 8, in between breakpoints
Z = 100_000

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=5), 195.08)
@test isapprox(round(pressure(atmos), sigdigits=2), round(3.2011e-2, sigdigits=2))
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=4), 28.40)
@test isapprox(round(density(atmos), sigdigits=2), round(5.604e-7, sigdigits=2))
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

# b = 9
Z = 115_000

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=5), 300.00)
@test isapprox(round(pressure(atmos), sigdigits=5), 4.0096e-3)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=4), 26.68)
@test isapprox(round(density(atmos), sigdigits=4), 4.289e-8)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

# b = 10
Z = 200_000

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=5), 854.56)
@test isapprox(round(pressure(atmos), sigdigits=5), 8.4736e-5)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=4), 21.30)
@test isapprox(round(density(atmos), sigdigits=3), round(2.541e-10, sigdigits=3))
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

# b = 11
Z = 750_000

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=5), 999.99)
@test isapprox(round(pressure(atmos), sigdigits=5), 2.2599e-8)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=3), 6.58)
@test isapprox(round(density(atmos), sigdigits=3), round(1.788e-14, sigdigits=3))
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

# b = 11, in between breakpoints
Z = 985_000

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=6), 1000.00)
@test isapprox(round(pressure(atmos), sigdigits=3), round(7.9185e-9, sigdigits=3))
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=3), 3.99)
@test isapprox(round(density(atmos), sigdigits=3), round(3.797e-15, sigdigits=3))
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

# b = 12
Z = 1_000_000

atmos = atmosphere(Z)
@test isapprox(round(temperature(atmos), sigdigits=6), 1000.00)
@test isapprox(round(pressure(atmos), sigdigits=5), 7.5138e-9)
@test isapprox(round(mean_molecular_weight(atmos), sigdigits=3), 3.94)
@test isapprox(round(density(atmos), sigdigits=4), 3.561e-15)
@test isapprox(round(speed_of_sound(atmos), sigdigits=5), 274.10)

ext: jl


text:
include("LibSCIP.jl")
using .LibSCIP

const SCIP_ = LibSCIP.SCIP
const TRUE = LibSCIP.TRUE
const FALSE = LibSCIP.FALSE

# SCIP_CALL: macro to check return codes, inspired by @assert
macro SC(ex)
    return :(@assert $(esc(ex)) == SCIP_OKAY)
end

const SCIP_STATISTICS_FUNCS = [
    :SCIPprintStatusStatistics,
    :SCIPprintTimingStatistics,
    :SCIPprintOrigProblemStatistics,
    :SCIPprintTransProblemStatistics,
    :SCIPprintPresolverStatistics,
    :SCIPprintConstraintStatistics,
    :SCIPprintConstraintTimingStatistics,
    :SCIPprintPropagatorStatistics,
    :SCIPprintConflictStatistics,
    :SCIPprintSeparatorStatistics,
    :SCIPprintPricerStatistics,
    :SCIPprintBranchruleStatistics,
    :SCIPprintHeuristicStatistics,
    :SCIPprintCompressionStatistics,
    :SCIPprintLPStatistics,
    :SCIPprintNLPStatistics,
    :SCIPprintRelaxatorStatistics,
    :SCIPprintTreeStatistics,
    :SCIPprintRootStatistics,
    :SCIPprintSolutionStatistics,
    :SCIPprintConcsolverStatistics,
    :SCIPprintBendersStatistics,
    :SCIPprintStatistics,
    :SCIPprintReoptStatistics,
    :SCIPprintBranchingStatistics
]


text:
int_rules_1_1_3_8 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.3.8*P(x)*(c*x)^m*(a+b*x^n)^p =#
    @apply_utils Antiderivative(((~x) ^ ~(m') * (~e + ~(f') * (~x) ^ ~(q') + ~(g') * (~x) ^ ~(r') + ~(h') * (~x) ^ ~(n'))) / (~a + ~(c') * (~x) ^ ~(n')) ^ (3 / 2), ~x) => -(((2 * ~a * ~g + 4 * ~a * ~h * (~x) ^ (~n / 4)) - 2 * ~c * ~f * (~x) ^ (~n / 2))) / (~a * ~c * ~n * sqrt(~a + ~c * (~x) ^ ~n)) <-- FreeQ([~a, ~c, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ(~q, ~n / 4) && (EqQ(~r, (3 * ~n) / 4) && (EqQ((4 * ~m - ~n) + 4, 0) && EqQ(~c * ~e + ~a * ~h, 0))))
    @apply_utils Antiderivative(((~d * ~x) ^ ~(m') * (~e + ~(f') * (~x) ^ ~(q') + ~(g') * (~x) ^ ~(r') + ~(h') * (~x) ^ ~(n'))) / (~a + ~(c') * (~x) ^ ~(n')) ^ (3 / 2), ~x) => ((~d * ~x) ^ ~m / (~x) ^ ~m) * Antiderivative(((~x) ^ ~m * (~e + ~f * (~x) ^ (~n / 4) + ~g * (~x) ^ ((3 * ~n) / 4) + ~h * (~x) ^ ~n)) / (~a + ~c * (~x) ^ ~n) ^ (3 / 2), ~x) <-- FreeQ([~a, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n], ~x) && (EqQ((4 * ~m - ~n) + 4, 0) && (EqQ(~q, ~n / 4) && (EqQ(~r, (3 * ~n) / 4) && EqQ(~c * ~e + ~a * ~h, 0))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~m * ~Pq * (~a + ~(b') * ~x) ^ ~p, ~x) => With([n = Denominator(~p)], (n / ~b) * Subst(Antiderivative((~x) ^ ((n * ~p + n) - 1) * ((-(~a) * ~c) / ~b + (~c * (~x) ^ n) / ~b) ^ ~m * ReplaceAll(~Pq, (~x->begin
                                            #= none:1 =#
                                            #= none:3 =#
                                            -(~a) / ~b + (~x) ^ n / ~b
                                        end)), ~x), ~x, (~a + ~b * ~x) ^ (1 / n))) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (PolyQ(~Pq, ~x) && (FractionQ(~p) && ILtQ(~m, -1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~(p'), ~x) => (1 / (~m + 1)) * Subst(Antiderivative(SubstFor((~x) ^ (~m + 1), ~Pq, ~x) * (~a + ~b * (~x) ^ Simplify(~n / (~m + 1))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (NeQ(~m, -1) && (IGtQ(Simplify(~n / (~m + 1)), 0) && PolyQ(~Pq, (~x) ^ (~m + 1))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~c * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~p, 0) || EqQ(~n, 1)))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~(p'), ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * SubstFor((~x) ^ ~n, ~Pq, ~x) * (~a + ~b * ~x) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && IntegerQ(Simplify((~m + 1) / ~n)))
    @apply_utils Antiderivative((~c * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~(p'), ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && IntegerQ(Simplify((~m + 1) / ~n)))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (~Pq * (~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * ~n * (~p + 1)) - (1 / (~b * ~n * (~p + 1))) * Antiderivative(D(~Pq, ~x) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~m, ~n], ~x) && (PolyQ(~Pq, ~x) && (EqQ((~m - ~n) + 1, 0) && LtQ(~p, -1)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (1 / ~d) * Antiderivative((~d * ~x) ^ (~m + 1) * PolynomialQuotient(~Pq, ~x, ~x) * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~d, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, ~x) && EqQ(Coeff(~Pq, ~x, 0), 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => Module([u = IntHide((~x) ^ ~m * ~Pq, ~x)], u * (~a + ~b * (~x) ^ ~n) ^ ~p - ~b * ~n * ~p * Antiderivative((~x) ^ (~m + ~n) * (~a + ~b * (~x) ^ ~n) ^ (~p - 1) * ExpandToSum(u / (~x) ^ (~m + 1), ~x), ~x)) <-- FreeQ([~a, ~b], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~n, 0) && (GtQ(~p, 0) && LtQ(~m + Expon(~Pq, ~x) + 1, 0))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => Module([q = Expon(~Pq, ~x), i], (~c * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ ~p * Sum((Coeff(~Pq, ~x, i) * (~x) ^ (i + 1)) / (~m + ~n * ~p + i + 1), [i, 0, q]) + ~a * ~n * ~p * Antiderivative((~c * ~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p - 1) * Sum((Coeff(~Pq, ~x, i) * (~x) ^ i) / (~m + ~n * ~p + i + 1), [i, 0, q]), ~x)) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (PolyQ(~Pq, ~x) && (IGtQ((~n - 1) / 2, 0) && GtQ(~p, 0)))
    @apply_utils Antiderivative(((~x) ^ 2 * ~P4) / (~a + ~(b') * (~x) ^ 4) ^ (3 / 2), ~x) => With([e = Coeff(~P4, ~x, 0), f = Coeff(~P4, ~x, 1), h = Coeff(~P4, ~x, 4)], -((f - 2 * h * (~x) ^ 3)) / (2 * ~b * sqrt(~a + ~b * (~x) ^ 4)) <-- EqQ(~b * e - 3 * ~a * h, 0)) <-- FreeQ([~a, ~b], ~x) && (PolyQ(~P4, ~x, 4) && (EqQ(Coeff(~P4, ~x, 2), 0) && EqQ(Coeff(~P4, ~x, 3), 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => With([q = ~m + Expon(~Pq, ~x)], Module([Q = PolynomialQuotient((~b) ^ (Floor((q - 1) / ~n) + 1) * (~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ ~n, ~x), R = PolynomialRemainder((~b) ^ (Floor((q - 1) / ~n) + 1) * (~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ ~n, ~x)], (-(~x) * R * (~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * ~n * (~p + 1) * (~b) ^ (Floor((q - 1) / ~n) + 1)) + (1 / (~a * ~n * (~p + 1) * (~b) ^ (Floor((q - 1) / ~n) + 1))) * Antiderivative((~a + ~b * (~x) ^ ~n) ^ (~p + 1) * ExpandToSum(~a * ~n * (~p + 1) * Q + ~n * (~p + 1) * R + D(~x * R, ~x), ~x), ~x)) <-- GeQ(q, ~n)) <-- FreeQ([~a, ~b], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~n, 0) && (LtQ(~p, -1) && IGtQ(~m, 0))))
    @apply_utils Antiderivative((~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], Module([Q = PolynomialQuotient(~a * (~b) ^ (Floor((q - 1) / ~n) + 1) * (~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ ~n, ~x), R = PolynomialRemainder(~a * (~b) ^ (Floor((q - 1) / ~n) + 1) * (~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ ~n, ~x), i], (-(~x) * R * (~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / ((~a) ^ 2 * ~n * (~p + 1) * (~b) ^ (Floor((q - 1) / ~n) + 1)) + (1 / (~a * ~n * (~p + 1) * (~b) ^ (Floor((q - 1) / ~n) + 1))) * Antiderivative((~x) ^ ~m * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * ExpandToSum(~n * (~p + 1) * (~x) ^ -(~m) * Q + Sum(((~n * (~p + 1) + i + 1) / ~a) * Coeff(R, ~x, i) * (~x) ^ (i - ~m), [i, 0, ~n - 1]), ~x), ~x))) <-- FreeQ([~a, ~b], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~n, 0) && (LtQ(~p, -1) && ILtQ(~m, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([g = GCD(~m + 1, ~n)], (1 / g) * Subst(Antiderivative((~x) ^ ((~m + 1) / g - 1) * ReplaceAll(~Pq, (~x->begin
                                                #= none:1 =#
                                                #= none:3 =#
                                                (~x) ^ (1 / g)
                                            end)) * (~a + ~b * (~x) ^ (~n / g)) ^ ~p, ~x), ~x, (~x) ^ g) <-- g != 1) <-- FreeQ([~a, ~b, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (IGtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative(((~(c') * ~x) ^ ~(m') * ~Pq) / (~a + ~(b') * (~x) ^ ~n), ~x) => With([v = Sum(((~c * ~x) ^ (~m + ii) * (Coeff(~Pq, ~x, ii) + Coeff(~Pq, ~x, ~n / 2 + ii) * (~x) ^ (~n / 2))) / ((~c) ^ ii * (~a + ~b * (~x) ^ ~n)), [ii, 0, ~n / 2 - 1])], Antiderivative(v, ~x) <-- SumQ(v)) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~n / 2, 0) && Expon(~Pq, ~x) < ~n))
    @apply_utils Antiderivative(~Pq / (~x * sqrt(~a + ~(b') * (~x) ^ ~n)), ~x) => Coeff(~Pq, ~x, 0) * Antiderivative(1 / (~x * sqrt(~a + ~b * (~x) ^ ~n)), ~x) + Antiderivative(ExpandToSum((~Pq - Coeff(~Pq, ~x, 0)) / ~x, ~x) / sqrt(~a + ~b * (~x) ^ ~n), ~x) <-- FreeQ([~a, ~b], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~n, 0) && NeQ(Coeff(~Pq, ~x, 0), 0)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => Module([q = Expon(~Pq, ~x), j, k], Antiderivative(Sum(((~c * ~x) ^ (~m + j) / (~c) ^ j) * Sum(Coeff(~Pq, ~x, j + (k * ~n) / 2) * (~x) ^ ((k * ~n) / 2), [k, 0, (2 * (q - j)) / ~n + 1]) * (~a + ~b * (~x) ^ ~n) ^ ~p, [j, 0, ~n / 2 - 1]), ~x)) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~n / 2, 0) && Not(PolyQ(~Pq, (~x) ^ (~n / 2)))))
    @apply_utils Antiderivative(((~(c') * ~x) ^ ~(m') * ~Pq) / (~a + ~(b') * (~x) ^ ~n), ~x) => Antiderivative(ExpandIntegrand(((~c * ~x) ^ ~m * ~Pq) / (~a + ~b * (~x) ^ ~n), ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (PolyQ(~Pq, ~x) && (IntegerQ(~n) && Not(IGtQ(~m, 0))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([Pq0 = Coeff(~Pq, ~x, 0)], (Pq0 * (~c * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * ~c * (~m + 1)) + (1 / (2 * ~a * ~c * (~m + 1))) * Antiderivative((~c * ~x) ^ (~m + 1) * ExpandToSum((2 * ~a * (~m + 1) * (~Pq - Pq0)) / ~x - 2 * ~b * Pq0 * (~m + ~n * (~p + 1) + 1) * (~x) ^ (~n - 1), ~x) * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- NeQ(Pq0, 0)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && (IGtQ(~n, 0) && (LtQ(~m, -1) && LeQ(~n - 1, Expon(~Pq, ~x)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], With([Pqq = Coeff(~Pq, ~x, q)], (Pqq * (~c * ~x) ^ (((~m + q) - ~n) + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * (~c) ^ ((q - ~n) + 1) * (~m + q + ~n * ~p + 1)) + (1 / (~b * (~m + q + ~n * ~p + 1))) * Antiderivative((~c * ~x) ^ ~m * ExpandToSum(~b * (~m + q + ~n * ~p + 1) * (~Pq - Pqq * (~x) ^ q) - ~a * Pqq * (((~m + q) - ~n) + 1) * (~x) ^ (q - ~n), ~x) * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x)) <-- NeQ(~m + q + ~n * ~p + 1, 0) && (q - ~n >= 0 && (IntegerQ(2 * ~p) || IntegerQ(~p + (q + 1) / (2 * ~n))))) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], -(Subst(Antiderivative((ExpandToSum((~x) ^ q * ReplaceAll(~Pq, (~x->begin
                                                        #= none:1 =#
                                                        #= none:3 =#
                                                        (~x) ^ -1
                                                    end)), ~x) * (~a + ~b * (~x) ^ -(~n)) ^ ~p) / (~x) ^ (~m + q + 2), ~x), ~x, 1 / ~x))) <-- FreeQ([~a, ~b, ~p], ~x) && (PolyQ(~Pq, ~x) && (ILtQ(~n, 0) && IntegerQ(~m)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([g = Denominator(~m), q = Expon(~Pq, ~x)], (-g / ~c) * Subst(Antiderivative((ExpandToSum((~x) ^ (g * q) * ReplaceAll(~Pq, (~x->begin
                                                        #= none:1 =#
                                                        #= none:3 =#
                                                        (~c) ^ -1 * (~x) ^ -g
                                                    end)), ~x) * (~a + ~b * (~c) ^ -(~n) * (~x) ^ (-g * ~n)) ^ ~p) / (~x) ^ (g * (~m + q + 1) + 1), ~x), ~x, 1 / (~c * ~x) ^ (1 / g))) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && (ILtQ(~n, 0) && FractionQ(~m)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([q = Expon(~Pq, ~x)], -((~c * ~x) ^ ~m) * ((~x) ^ -1) ^ ~m * Subst(Antiderivative((ExpandToSum((~x) ^ q * ReplaceAll(~Pq, (~x->begin
                                                        #= none:1 =#
                                                        #= none:3 =#
                                                        (~x) ^ -1
                                                    end)), ~x) * (~a + ~b * (~x) ^ -(~n)) ^ ~p) / (~x) ^ (~m + q + 2), ~x), ~x, 1 / ~x)) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (ILtQ(~n, 0) && Not(RationalQ(~m))))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => With([g = Denominator(~n)], g * Subst(Antiderivative((~x) ^ (g * (~m + 1) - 1) * ReplaceAll(~Pq, (~x->begin
                                            #= none:1 =#
                                            #= none:3 =#
                                            (~x) ^ g
                                        end)) * (~a + ~b * (~x) ^ (g * ~n)) ^ ~p, ~x), ~x, (~x) ^ (1 / g))) <-- FreeQ([~a, ~b, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && FractionQ(~n))
    @apply_utils Antiderivative((~c * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && FractionQ(~n))
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (1 / (~m + 1)) * Subst(Antiderivative(ReplaceAll(SubstFor((~x) ^ ~n, ~Pq, ~x), (~x->begin
                                        #= none:1 =#
                                        #= none:3 =#
                                        (~x) ^ Simplify(~n / (~m + 1))
                                    end)) * (~a + ~b * (~x) ^ Simplify(~n / (~m + 1))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n))))
    @apply_utils Antiderivative((~c * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, (~x) ^ ~n) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ ~n) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~c * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~n, ~p], ~x) && ((PolyQ(~Pq, ~x) || PolyQ(~Pq, (~x) ^ ~n)) && Not(IGtQ(~m, 0)))
    @apply_utils Antiderivative((~u) ^ ~(m') * ~Pq * (~a + ~(b') * (~v) ^ ~(n')) ^ ~p, ~x) => ((~u) ^ ~m / (Coeff(~v, ~x, 1) * (~v) ^ ~m)) * Subst(Antiderivative((~x) ^ ~m * SubstFor(~v, ~Pq, ~x) * (~a + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x, ~v) <-- FreeQ([~a, ~b, ~m, ~n, ~p], ~x) && (LinearPairQ(~u, ~v, ~x) && PolyQ(~Pq, (~v) ^ ~n))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a1 + ~(b1') * (~x) ^ ~(n')) ^ ~(p') * (~a2 + ~(b2') * (~x) ^ ~(n')) ^ ~(p'), ~x) => Antiderivative((~c * ~x) ^ ~m * ~Pq * (~a1 * ~a2 + ~b1 * ~b2 * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a1, ~b1, ~a2, ~b2, ~c, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~a2 * ~b1 + ~a1 * ~b2, 0) && (IntegerQ(~p) || GtQ(~a1, 0) && GtQ(~a2, 0))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a1 + ~(b1') * (~x) ^ ~(n')) ^ ~(p') * (~a2 + ~(b2') * (~x) ^ ~(n')) ^ ~(p'), ~x) => (((~a1 + ~b1 * (~x) ^ ~n) ^ FracPart(~p) * (~a2 + ~b2 * (~x) ^ ~n) ^ FracPart(~p)) / (~a1 * ~a2 + ~b1 * ~b2 * (~x) ^ (2 * ~n)) ^ FracPart(~p)) * Antiderivative((~c * ~x) ^ ~m * ~Pq * (~a1 * ~a2 + ~b1 * ~b2 * (~x) ^ (2 * ~n)) ^ ~p, ~x) <-- FreeQ([~a1, ~b1, ~a2, ~b2, ~c, ~m, ~n, ~p], ~x) && (PolyQ(~Pq, ~x) && (EqQ(~a2 * ~b1 + ~a1 * ~b2, 0) && Not(EqQ(~n, 1) && LinearQ(~Pq, ~x))))
    @apply_utils Antiderivative((~(h') * ~x) ^ ~(m') * (~e + ~(f') * (~x) ^ ~(n') + ~(g') * (~x) ^ ~(n2')) * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(p') * (~c + ~(d') * (~x) ^ ~(n')) ^ ~(p'), ~x) => (~e * (~h * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~p + 1)) / (~a * ~c * ~h * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~f, ~g, ~h, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~a * ~c * ~f * (~m + 1) - ~e * (~b * ~c + ~a * ~d) * (~m + ~n * (~p + 1) + 1), 0) && (EqQ(~a * ~c * ~g * (~m + 1) - ~b * ~d * ~e * (~m + 2 * ~n * (~p + 1) + 1), 0) && NeQ(~m, -1))))
    @apply_utils Antiderivative((~(h') * ~x) ^ ~(m') * (~e + ~(g') * (~x) ^ ~(n2')) * (~a + ~(b') * (~x) ^ ~(n')) ^ ~(p') * (~c + ~(d') * (~x) ^ ~(n')) ^ ~(p'), ~x) => (~e * (~h * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ ~n) ^ (~p + 1) * (~c + ~d * (~x) ^ ~n) ^ (~p + 1)) / (~a * ~c * ~h * (~m + 1)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~g, ~h, ~m, ~n, ~p], ~x) && (EqQ(~n2, 2 * ~n) && (EqQ(~m + ~n * (~p + 1) + 1, 0) && (EqQ(~a * ~c * ~g * (~m + 1) - ~b * ~d * ~e * (~m + 2 * ~n * (~p + 1) + 1), 0) && NeQ(~m, -1))))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.3 General/1.1.3.8 P(x) (c x)^m (a+b x^n)^p.jl


text:
# This example mirrors rest.jl, but via the CRUD API.

module CRUD

using Discord

const GUILD = 494962347434049536  # Update this to a guild you have access to.

function main()
    c = Client(ENV["DISCORD_TOKEN"])

    # Send a request.
    future = retrieve(c, Guild, GUILD)
    # Then await it.
    resp = fetch(future)
    # Check for success.
    if !resp.ok
        # You can access the HTTP response if it went through.
        if resp.http_response !== nothing
            println("HTTP status: ", resp.http_response.status)
        end
        # You can also access caught exceptions.
        if resp.exception !== nothing
            println("Request error: ", sprint(showerror, resp.exception))
        end
        error("Retrieving guild failed")
    end

    # Get the actual Guild.
    guild = resp.val
    # If we're really confident, we can skip the error checks.
    guild = fetchval(future)

    println("Guild name: $(guild.name)")

    # Create a channel.
    channel = fetchval(create(c, DiscordChannel, guild; name="foo"))

    # Send a message to the channel.
    message = fetchval(create(c, Message, channel; content="Hello, world!"))
    # React to it.
    create(c, Reaction, message, '👍')
    # Edit it
    update(c, message; content="Goodbye, world!")
    # Delete it.
    delete(c, message)

    # Delete the channel.
    delete(c, channel)
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    CRUD.main()
end


text:
@testset "IntDomain.jl" begin
    @testset "isempty()" begin
        trailer = SeaPearl.Trailer()
        domNotEmpty = SeaPearl.IntDomain(trailer, 20, 10)

        @test !isempty(domNotEmpty)

        emptyDom = SeaPearl.IntDomain(trailer, 0, 1)
        @test isempty(emptyDom)
    end

    @testset "length()" begin
        trailer = SeaPearl.Trailer()
        dom20 = SeaPearl.IntDomain(trailer, 20, 10)
        @test length(dom20) == 20

        dom2 = SeaPearl.IntDomain(trailer, 2, 0)

        @test length(dom2) == 2
    end

    @testset "in()" begin
        trailer = SeaPearl.Trailer()
        dom20 = SeaPearl.IntDomain(trailer, 20, 10)

        @test 11 in dom20
        @test !(10 in dom20)
        @test 21 in dom20
        @test !(1 in dom20)
        @test !(31 in dom20)
    end

    @testset "exchangePositions!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntDomain(trailer, 5, 10)

        @test dom.values == [1, 2, 3, 4, 5]
        @test dom.indexes == [1, 2, 3, 4, 5]

        SeaPearl.exchangePositions!(dom, 2, 5)

        @test dom.values == [1, 5, 3, 4, 2]
        @test dom.indexes == [1, 5, 3, 4, 2]

        SeaPearl.exchangePositions!(dom, 2, 1)

        @test dom.values == [2, 5, 3, 4, 1]
        @test dom.indexes == [5, 1, 3, 4, 2]
    end

    @testset "remove!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntDomain(trailer, 5, 10)

        SeaPearl.remove!(dom, 11)

        @test !(11 in dom)
        @test length(dom) == 4
        @test dom.min.value == 12
    end

    @testset "assign!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntDomain(trailer, 5, 10)

        removed = SeaPearl.assign!(dom, 14)

        @test !(12 in dom)
        @test 14 in dom
        @test length(dom) == 1
        @test sort(removed) == [11, 12, 13, 15]
        @test dom.max.value == 14 && dom.min.value == 14
    end

    @testset "iterate()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntDomain(trailer, 3, 10)

        j = 11
        for i in dom
            @test j == i
            j += 1
        end

        @test j == 14

    end

    @testset "removeAll!()" begin
        trailer = SeaPearl.Trailer()
        dom = SeaPearl.IntDomain(trailer, 3, 10)

        removed = SeaPearl.removeAll!(dom)

        @test isempty(dom)
        @test removed == [11, 12, 13]
    end

    @testset "updateMinFromRemovedVal!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)

        @test x.domain.min.value == 5

        SeaPearl.remove!(x.domain, 5)
        SeaPearl.updateMinFromRemovedVal!(x.domain, 5)

        @test x.domain.min.value == 6
    end

    @testset "updateMaxFromRemovedVal!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)

        @test x.domain.max.value == 10

        SeaPearl.remove!(x.domain, 10)
        SeaPearl.updateMaxFromRemovedVal!(x.domain, 10)

        @test x.domain.max.value == 9
    end

    @testset "updateBoundsFromRemovedVal!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)

        SeaPearl.remove!(x.domain, 5)
        SeaPearl.remove!(x.domain, 10)
        SeaPearl.updateBoundsFromRemovedVal!(x.domain, 5)
        SeaPearl.updateBoundsFromRemovedVal!(x.domain, 10)

        @test x.domain.min.value == 6
        @test x.domain.max.value == 9
    end

    @testset "removeAbove!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)

        @test sort(SeaPearl.removeAbove!(x.domain, 7)) == [8, 9, 10]

        @test length(x.domain) == 3
        @test 7 in x.domain
        @test 6 in x.domain
        @test !(8 in x.domain)

        @test sort(SeaPearl.removeAbove!(x.domain, 4)) == [5, 6, 7]

        @test isempty(x.domain)
    end

    @testset "removeBelow!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)

        @test sort(SeaPearl.removeBelow!(x.domain, 7)) == [5, 6]

        @test length(x.domain) == 4
        @test 7 in x.domain
        @test 8 in x.domain
        @test !(6 in x.domain)

        @test sort(SeaPearl.removeBelow!(x.domain, 11)) == [7, 8, 9, 10]

        @test isempty(x.domain)
    end

    @testset "minimum()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)

        @test SeaPearl.minimum(x.domain) == 5
    end

    @testset "maximum()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)

        @test SeaPearl.maximum(x.domain) == 10
    end

    @testset "reset_domain!()" begin
        trailer = SeaPearl.Trailer()
        x = SeaPearl.IntVar(5, 10, "x", trailer)
        SeaPearl.assign!(x, 10)

        @test x.domain.values == [6, 2, 3, 4, 5, 1]
        @test x.domain.indexes == [6, 2, 3, 4, 5, 1]
        @test SeaPearl.length(x.domain) == 1
        SeaPearl.reset_domain!(x.domain)
        @test x.domain.values == [1, 2, 3, 4, 5, 6]
        @test x.domain.indexes == [1, 2, 3, 4, 5, 6]
        @test SeaPearl.length(x.domain) == 6
    end
end


text:
# using PkgTemplates

# t = Template(;
#     user="vanhp",
#     authors="Vanh Phom",
#     julia=v"1.5",
#     plugins=[

#     ],
# )
# generate("MicroServiceJulia",t)

max_stars_repo_name: vanhp/MicroServiceJulia

max_issues_repo_name: vanhp/MicroServiceJulia

max_forks_repo_name: vanhp/MicroServiceJulia


text:
#######################################################################
# Coverage.jl
# Input: Code coverage and memory allocations
# Output: Useful things
# https://github.com/JuliaCI/Coverage.jl
#######################################################################
module Coverage
    using CoverageTools
    using LibGit2

    export FileCoverage
    export LCOV
    export analyze_malloc
    export amend_coverage_from_src!
    export clean_file
    export clean_folder
    export get_summary
    export merge_coverage_counts
    export process_cov
    export process_file
    export process_folder

    const CovCount = CoverageTools.CovCount
    const FileCoverage = CoverageTools.FileCoverage
    const amend_coverage_from_src! = CoverageTools.amend_coverage_from_src!
    const clean_file = CoverageTools.clean_file
    const clean_folder = CoverageTools.clean_folder
    const get_summary = CoverageTools.get_summary
    const iscovfile = CoverageTools.iscovfile
    const merge_coverage_counts = CoverageTools.merge_coverage_counts
    const process_cov = CoverageTools.process_cov
    const process_file = CoverageTools.process_file
    const process_folder = CoverageTools.process_folder

    include("coveralls.jl")
    include("codecovio.jl")
    include("lcov.jl")
    include("memalloc.jl")
    include("parser.jl")
end


text:
# Group format description https://www.c3d.org/HTML/groupformat1.htm
struct Group
    pos::Int
    nl::Int8 # Number of characters in group name
    isLocked::Bool # Locked if nl < 0
    gid::Int8 # Group ID
    name::String # Character set should be A-Z, 0-9, and _ (lowercase is ok)
    symname::Symbol
    np::Int16 # Pointer in bytes to the start of next group/parameter (officially supposed to be signed)
    dl::UInt8 # Number of characters in group description (nominally should be between 1 and 127)
    desc::String # Character set should be A-Z, 0-9, and _ (lowercase is ok)
    params::Dict{Symbol,AbstractParameter}
end

function Base.getproperty(g::Group, name::Symbol)
    if name ∈ fieldnames(typeof(g))
        return getfield(g, name)
    else
        return getfield(g, :params)[name].data
    end
end
Base.show(io::IO, g::Group) = show(io, keys(g.params))

function readgroup(f::IOStream, FEND::Endian, FType::Type{T}) where T <: Union{Float32,VaxFloatF}
    pos = position(f)
    nl = read(f, Int8)
    @assert nl != 0
    isLocked = nl < 0 ? true : false
    gid = read(f, Int8)
    @assert gid != 0
    name = transcode(String, read(f, abs(nl)))
    @assert any(!iscntrl, name)
    symname = Symbol(replace(strip(name), r"[^a-zA-Z0-9_]" => '_'))

    if occursin(r"[^a-zA-Z0-9_ ]", name)
        @debug "Group $name at $pos has unofficially supported characters.
            Unexpected results may occur"
    end

    np = saferead(f, Int16, FEND)
    dl = read(f, UInt8)
    desc = transcode(String, read(f, dl))

    return Group(pos, nl, isLocked, gid, name, symname, np, dl, desc, Dict{Symbol,AbstractParameter}())
end


text:
# Represents an N-dimensional ellipsoid
struct Ellipsoid
    ctr::Vector{Float64}  # center coordinates
    cov::Array{Float64, 2}
    icov::Array{Float64, 2}  # inverse of cov
    vol::Float64
end

# Draw a random point from within a unit N-ball
function randnball(ndim)
    z = randn(ndim)
    r2 = 0.
    for i=1:ndim
        r2 += z[i]*z[i]
    end
    factor = rand()^(1. /ndim) / sqrt(r2)
    for i=1:ndim
        z[i] *= factor
    end
    return z
end

# proportionality constant depending on dimension
# for n even:      (2pi)^(n    /2) / (2 * 4 * ... * n)
# for n odd :  2 * (2pi)^((n-1)/2) / (1 * 3 * ... * n)
function nball_vol_factor(ndim::Int)
    if ndim % 2 == 0
        c = 1.
        for i=2:2:ndim
            c *= 2pi / i
        end
        return c
    else
        c = 2.
        for i = 3:2:ndim
            c *= 2pi / i
        end
        return c
    end
end

function ellipsoid_volume(scaled_cov::Matrix{Float64})
    ndim = size(scaled_cov, 1)
    return nball_vol_factor(ndim) * sqrt(det(scaled_cov))
end

# find the bounding ellipsoid of points x where
function bounding_ellipsoid(x::Matrix{Float64}, enlarge=1.0)

    ndim, npoints = size(x)

    ctr = mean(x, dims=2)[:, 1]
    delta = x .- ctr
    cov = unscaled_covzm(delta, 2)
    icov = inv(cov)

    # Calculate expansion factor necessary to bound each point.
    # This finds the maximum of (delta_i' * icov * delta_i)
    fmax = -Inf
    for k in 1:npoints
        f = 0.0
        for j=1:ndim
            for i=1:ndim
                f += icov[i, j] * delta[i, k] * delta[j, k]
            end
        end
        fmax = max(fmax, f)
    end

    fmax *= enlarge
    cov .*= fmax
    icov .*= 1. /fmax
    vol = ellipsoid_volume(cov)

    return Ellipsoid(ctr, cov, icov, vol)
end

function sample_ellipsoid(ell::Ellipsoid)
    ndim = length(ell.ctr)

    # Get scaled eigenvectors (in columns): vs[:,i] is the i-th eigenvector.
    f = eigen(ell.cov)
    v, w = f.vectors, f.values
    for j=1:ndim
        tmp = sqrt(abs(w[j]))
        for i=1:ndim
            v[i, j] *= tmp
        end
    end

    return v*randnball(ndim) + ell.ctr
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "HiGHS"
version = v"0.1.4"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/ERGO-Code/HiGHS.git", "0b4cb7f882a6c5d208fbb2edb4d8c416291ed4f5"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
if [[ "${target}" == *86*-linux-musl* ]]; then
    pushd /opt/${target}/lib/gcc/${target}/*/include
    # Fix bug in Musl C library, see
    # https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/387
    atomic_patch -p0 $WORKSPACE/srcdir/patches/mm_malloc.patch
    popd
fi
mkdir -p HiGHS/build
cd HiGHS/build
apk add --upgrade cmake --repository=http://dl-cdn.alpinelinux.org/alpine/edge/main
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTING=OFF ..
make -j${nproc} highs
make install highs
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_gfortran_versions(
    expand_cxxstring_abis(supported_platforms())
)
filter!(p -> !isa(p, Windows), platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libhighs", :libhighs),
    LibraryProduct("libipx", :libipx)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"4.9")


text:
using SpectralDistances, Zygote
using SpectralDistances:
    ngradient, nhessian, njacobian, polyconv, eigsort, hproots, rev, companion
using SpectralDistances: getARXregressor, getARregressor
cosdist(x1, x2) = x1'x2 / norm(x1) / norm(x2)
using SpectralDistances, Zygote
Random.seed!(1)

using FiniteDifferences

# using Zygote
# function jacobian(m,x)
#     y  = m(x)
#     k  = length(y)
#     n  = length(x)
#     J  = Matrix{eltype(x)}(undef,k,n)
#     for i = 1:k
#         g = Zygote.gradient(x->m(x)[i], x)[1]
#         J[i,:] .= g
#     end
#     J
# end
# using ForwardDiff, FiniteDifferences
#
# function hessian(f, a)
#     H = jacobian(f', a)
# end
#
#
# function hessian_fd_reverse(f,a)
#     fd = central_fdm(9,1)
#     FiniteDifferences.jacobian(fd, f', a)
# end
#
# function nhessian_reverse(f,a)
#     njacobian(f', a)
# end
#
# a = randn(11)
# a[1] = 1
# using DoubleFloats
# m = AR(ContinuousRoots([-Float64(0.1)+im, -Float64(0.1)-im]))
# m2 = AR(ContinuousRoots([-Float64(0.11)+im, -Float64(0.12)-im]))
# a = Vector(m.ac)
# dist = CoefficientDistance(domain=Discrete(),distance=CosineDist())
# dist = EuclideanRootDistance(domain=Continuous(), p=2, weight=residueweight)
# dist = EuclideanRootDistance(domain=Continuous(), p=2, weight=unitweight)
# dist = OptimalTransportRootDistance(domain=Continuous(), p=2, iters=500, β=0.5)
# function forward(b)
#     real(dist(a, b))
# end
# forward(a)
# forward'(a)
# Zygote.refresh()
# dist(m,m2) |> typeof
# g = Zygote.gradient(Float64 ∘ dist, m.pc, m2.pc)
#
# mdist = ModelDistance(TLS(na=4), dist)
# t = 0:0.01:100
# x1 = sin.(70 .*t) .+ sin.(150 .* t) .+ 0.1 .*randn.()
# # x1 = filt(ones(10),[10],randn(10000))
# x2 = randn(10000)
# mdist(x1,x2)
#
# using TotalLeastSquares
# function TotalLeastSquares.tls(A::AbstractArray,y::AbstractArray)
#     AA  = [A y]
#     U,S,V   = svd(AA)
#     n   = size(A,2)
#     V21 = V[1:n,n+1:end]
#     V22 = V[n+1:end,n+1:end]
#     x   = -V21/V22
# end
#
# shift1(x) = x .+ x[2]
# g1,g2 = Zygote.gradient(abs ∘ mdist, x1, x2)
# gradfun(xo,x2) = xo .= Zygote.gradient(abs ∘ mdist, x2, x1)[1]
# gradfun(x2) = Zygote.gradient(abs ∘ mdist, x2, x1)[1]
# gradfun(x2)
# # x1 = filt(ones(10),[10],randn(10000))
# # x2 = randn(10000)
# for i = 1:1000
#     if i % 40 == 0
#         P = welch_pgram(x1)
#         plot(P)
#         P = welch_pgram(x2)
#         plot!(P) |> display
#     end
#     @show mdist(x1,x2)
#     x2 .-= 30 .* gradfun(x2)
# end
# using Optim
# res = Optim.optimize(x2->mdist(x2,x1), gradfun, x2, LBFGS(), Optim.Options(store_trace=true, show_trace=true, show_every=1, iterations=100, allow_f_increases=true, time_limit=100, x_tol=0, f_tol=0, g_tol=1e-8, f_calls_limit=100, g_calls_limit=0))
# x2 = res.minimizer

# ngradient(Float64 ∘ dist, m, m2)

# Zygote.refresh()
#
# # NOTE: Gradient at a should be almost zero as it is the minimum, everything else is incorrect. Hessian should be pos.def.
# @time forward'(a)
# @time H = hessian_fd_reverse(forward, a)
# @time H = nhessian_reverse(forward, a)
# c∇ = central_fdm(3,1)
#
# grad(c∇, forward, (a))
# ngradient(forward, copy(a))
#
#
# @btime H = hessian_fd_reverse($(forward), $a)
#
# eigvals(H[2:end,2:end])
#
# hessian(forward, a)

# true && get(ENV, "TRAVIS_BRANCH", nothing) == nothing && @testset "gradients" begin

y = (randn(10))
u = (randn(10))
@test ngradient((y) -> sum(getARXregressor(y, u, 2, 2)[2]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 2, 2)[2]), y)[1]
@test ngradient((y) -> sum(getARXregressor(y, u, 2, 2)[1]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 2, 2)[1]), y)[1]
@test ngradient((y) -> sum(getARXregressor(y, u, 5, 2)[2]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 5, 2)[2]), y)[1]
@test ngradient((y) -> sum(getARXregressor(y, u, 5, 2)[1]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 5, 2)[1]), y)[1]

@test ngradient((y) -> sum(getARXregressor(y, u, 2, 1)[2]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 2, 1)[2]), y)[1]
@test ngradient((y) -> sum(getARXregressor(y, u, 2, 1)[1]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 2, 1)[1]), y)[1]
@test ngradient((y) -> sum(getARXregressor(y, u, 5, 1)[2]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 5, 1)[2]), y)[1]
@test ngradient((y) -> sum(getARXregressor(y, u, 5, 1)[1]), y) ==
      Zygote.gradient((y) -> sum(getARXregressor(y, u, 5, 1)[1]), y)[1]

y = (randn(10))
@test ngradient((y) -> sum(getARregressor(y, 2)[2]), y) ==
      Zygote.gradient((y) -> sum(getARregressor(y, 2)[2]), y)[1]
@test ngradient((y) -> sum(getARregressor(y, 2)[1]), y) ==
      Zygote.gradient((y) -> sum(getARregressor(y, 2)[1]), y)[1]
@test ngradient((y) -> sum(getARregressor(y, 5)[2]), y) ==
      Zygote.gradient((y) -> sum(getARregressor(y, 5)[2]), y)[1]
@test ngradient((y) -> sum(getARregressor(y, 5)[1]), y) ==
      Zygote.gradient((y) -> sum(getARregressor(y, 5)[1]), y)[1]

a, b = randn(3), randn(4)
@test ngradient(a -> sum(SpectralDistances.polyconv(a, b)), a) ≈
      Zygote.gradient(a -> sum(SpectralDistances.polyconv(a, b)), a)[1] rtol = 1e-3
@test ngradient(b -> sum(SpectralDistances.polyconv(a, b)), b) ≈
      Zygote.gradient(b -> sum(SpectralDistances.polyconv(a, b)), b)[1] rtol = 1e-3

@test ngradient(a -> sum(SpectralDistances.polyvec(a)), a) ≈
      Zygote.gradient(a -> sum(SpectralDistances.polyvec(a)), a)[1] rtol = 1e-3

@test_skip let Gc = tf(1, [1, 1, 1, 1])
    w = c2d(Gc, 1).matrix[1] |> ControlSystems.denvec
    @test d2c(w) ≈ pole(Gc)
end

# y = randn(5000)
# fm = PLR(na=40, nc=2)
# @test_skip Zygote.gradient(y->sum(fm(y)[1]), y)[1] ≈ ForwardDiff.gradient(y->sum(fm(y)[1]), y)
# @test_skip Zygote.gradient(y->sum(fm(y)[2]), y)[1] ≈ ForwardDiff.gradient(y->sum(fm(y)[2]), y)

p = [1.0, 1, 1]
# @btime riroots(p)
fd = central_fdm(3, 1)
@test Zygote.gradient(p) do p
    r = roots(p)
    sum(abs2, r)
end[1][1:end-1] ≈ FiniteDifferences.grad(fd, p -> begin
    r = roots(p)
    sum(abs2, r)
end, p)[1][1:end-1]

fm = TLS(na = 4)
y = randn(50)
sum(abs2, fm(y).pc)
@test_broken Zygote.gradient(y) do y
    sum(abs2, fitmodel(fm, y, true).pc)
end

fdm = central_fdm(5, 1)

a = randn(5)
a[end] = 1
f = x -> sum(abs2(r) for r in (hproots(x)))

G = Zygote.gradient(f, a)[1]
@test G[1:end-1] ≈ FiniteDifferences.grad(fdm, f, a)[1][1:end-1]

fd = central_fdm(5, 1)
a = randn(30)
a[1] = 1
r = roots(reverse(a)) |> ContinuousRoots
@inferred residues(a, 1, r)
f = a -> sum(abs2, residues(a, 1, r))
g = a -> real(residues(complex.(a), 1, r)[2])
@test sum(abs, f'(complex.(a))[2:end] - grad(fd, f, a)[1][2:end]) < sqrt(eps())
@test sum(abs, g'((a))[2:end] - grad(fd, g, a)[1][2:end]) < sqrt(eps())

fd = central_fdm(7, 1)
a = randn(30)
a[1] = 1
residues(a, 1)
f = a -> sum(abs2, residues(a, 1))
g = a -> real(residues((a), 1)[2])
@test sum(abs, f'((a))[2:end] - grad(fd, f, a)[1][2:end]) < 1e-5
@test_skip sum(abs, g'((a))[2:end] - grad(fd, g, a)[1][2:end]) < sqrt(eps()) # Not robust

@testset "Numerical curvature" begin
    @info "Testing Numerical curvature"
    a = randn(11)
    a[1] = 1
    using DoubleFloats
    m = AR(ContinuousRoots([-Double64(0.1) + im, -Double64(0.1) - im]))
    a = Vector(m.ac)
    for dist in [
        EuclideanRootDistance(domain = Continuous(), p = 1, weight = residueweight),
        OptimalTransportRootDistance(domain = Continuous(), p = 2, β = 0.01),
    ]
        H = SpectralDistances.curvature(dist, Float64.(a))
        @test all(>(0) ∘ real, eigvals(H[2:end, 2:end]))
        @test all(==(0) ∘ imag, eigvals(H[2:end, 2:end]))
    end

end

# Zygote segfaults here
# @testset "Sinkhorn zygote" begin
#     @info "testing sinkhorn zygote"
#     function sinkdist(D,a,b)
#         ai = s1(a)
#         bi = s1(a+b)
#         P,u,v = sinkhorn(D,ai,bi, iters=1000, β=0.1)
#         dot(P, D)
#     end
#     a,b = abs.(randn(6)),abs.(randn(6))
#     D = SpectralDistances.distmat_euclidean(1:length(a), 1:length(a))
#     dD, da, db = Zygote.gradient(sinkdist, D,a,b)
#     @test n1(ForwardDiff.gradient(a->sinkdist(D,a,b), a))'n1(da) > 0.9
# end

a = -5:-1
@test SpectralDistances.roots2poly(a) ≈ SpectralDistances.roots2poly_zygote(a)

##

@testset "More roots diff tests" begin
    @info "Testing More roots diff tests"

    rootfun = eigsort ∘ hproots ∘ rev
    rootfunnr = eigsort ∘ hproots
    # Testing roots carefully
    a = polyconv([1.0, 1], [1.1, 1])
    b = polyconv([1.0, 1], [2, 1])

    @test SpectralDistances.roots(a) ≈ [-1.1, -1]
    @test SpectralDistances.roots(b) ≈ [-2, -1]

    @test all(hproots(a) .∈ Ref(Complex.([-1.1, -1])))
    @test all(hproots(b) .∈ Ref(Complex.([-2, -1])))

    g1 = Zygote.gradient(a -> sum(abs, rootfunnr(a)), a)[1]
    g2 = ngradient(a -> sum(abs, rootfunnr(a)), a)
    @test real(g1)[1:end-1] ≈ real(g2)[1:end-1]
    g1 = Zygote.gradient(b -> sum(abs, rootfunnr(b)), b)[1]
    g2 = ngradient(b -> sum(abs, rootfunnr(b)), b)
    @test real(g1)[1:end-1] ≈ real(g2)[1:end-1]

    a = [1, 0.1, 1]
    g1 = Zygote.gradient(a -> sum(abs, rootfunnr(a)), a)[1]
    g2 = ngradient(a -> sum(abs, rootfunnr(a)), a)
    @test real(g1)[1:end-1] ≈ real(g2)[1:end-1]

    a = [1, 0.1, 1]
    g1 = Zygote.gradient(a -> sum(abs2, rootfunnr(a)), a)[1]
    g2 = ngradient(a -> sum(abs2, rootfunnr(a)), a)
    @test real(g1)[1:end-1] ≈ real(g2)[1:end-1]

    a = [1, 0.1, 1]
    g1 = Zygote.gradient(a -> real(rootfunnr(a)[2]), a)[1]
    g2 = ngradient(a -> real(rootfunnr(a)[2]), a)
    @test real(g1)[1:end-1] ≈ real(g2)[1:end-1]

    g1 = Zygote.gradient(a -> imag(rootfunnr(a)[2]), a)[1]
    g2 = ngradient(a -> imag(rootfunnr(a)[2]), a)
    @test real(g1)[1:end-1] ≈ real(g2)[1:end-1]

    function getpoly(n, γ = 0.1)
        skew(x) = x - x' + 0.1 * (x'x)
        r = eigsort(eigvals(skew(γ .* randn(n, n))))
        a = Vector(SpectralDistances.roots2poly(r))
        a, r
    end
    a, r = getpoly(9)

    r1 = @inferred rootfun(a)
    r2 = @inferred eigsort(eigen(companion(rev(a))).values)
    @test cosdist(real(r1), real(r2)) > 0.99
    @test cosdist(real(r1), real(r)) > 0.99

    ##
    a = [1.0, 1]
    for i = 2:10
        isinteractive() && @show i
        # global a
        a = SpectralDistances.polyconv(a, [1, -i])
        isinteractive() && @show a
        # a[1] = 1
        # a[end] = 1

        g1 = Zygote.gradient(a -> sum(abs2, rootfun(a)), a)[1]
        g2 = ngradient(a -> sum(abs2, rootfun(a)), a)
        @test real(g1)[2:end] ≈ real(g2)[2:end] atol=1e-3

        g1 = Zygote.gradient(a -> real(rootfun(a)[2]), a)[1]
        g2 = ngradient(a -> real(rootfun(a)[2]), a)
        @test (g1)[2:end] ≈ (g2)[2:end] atol = 1e-4 atol=1e-3
        #
        g1 = Zygote.gradient(a -> imag(rootfun(a)[2]), a)[1]
        g2 = ngradient(a -> imag(rootfun(a)[2]), a)
        @test real(g1)[2:end] ≈ real(g2)[2:end] atol=1e-3
    end

    # NOTE: eigsort is improtant
    g1 = Zygote.gradient(a -> sum(abs2, rootfun((a))), (a))[1]
    g2 = ngradient(a -> sum(abs2, rootfun((a))), a)
    @test cosdist(real(g1)[2:end], real(g2)[2:end]) > 0.9
    @test real(g1)[2:end] ≈ real(g2)[2:end] rtol = 0.3

    g1 = Zygote.gradient(a -> real(rootfun(a)[2]), a)[1]
    g2 = ngradient(a -> real(rootfun(a)[2]), a)
    @test cosdist(g1[2:end], g2[2:end]) > 0.9
    @test g1[2:end] ≈ g2[2:end] rtol = 1e-3

    g1 = Zygote.gradient(a -> imag(rootfun(a)[2]), a)[1]
    g2 = ngradient(a -> imag(rootfun(a)[2]), a)
    @test real(g1)[2:end] ≈ real(g2)[2:end]

    ##
    # using FiniteDifferences
    # fdm = central_fdm(5, 1)
    t = 0:100
    x1 = sin.(0.70 .* t) .+ sin.(0.50 .* t) .+ 0.1 .* randn.()
    x2 = sin.(0.50 .* t) .+ sin.(0.30 .* t) .+ 0.1 .* randn.()
    fm = LS(na = 10)
    m1, m2 = fm.((x1, x2))

    testfun1(x1) = sum(abs2, fm(x1).a)
    @test Zygote.gradient(testfun1, x1)[1] ≈ grad(fdm, testfun1, x1)[1] rtol = 1e-3
    testfun2(x1) = sum(abs2, (fm(x1).p))
    g1 = Zygote.gradient(testfun2, x1)[1]
    g2 = ngradient(testfun2, x1, δ = 1e-4)
    # g3 = grad(fdm, testfun2, x1)[1]
    @test cosdist(g1, g2) > 0.95
    # cosdist(g3, g2)

    testfun3(x1) = sum(abs2, fm(x1).pc)
    g1 = Zygote.gradient(testfun3, x1)[1]
    g2 = ngradient(testfun3, x1)
    # g3 = grad(fdm, testfun3, x1)[1]
    @test cosdist(g1, g2) > 0.95
    # cosdist(g3, g2)

    dist = ModelDistance(fm, OptimalTransportRootDistance(domain = Continuous(), β = 0.1))
    Zygote.@nograd rand
    # Zygote.refresh()

    begin # This test is broken due to Zygote being bad at Complex: (MethodError: no method matching isless(::ComplexF64, ::Float64)) appears when it shouldn't
        # df = x -> real(evaluate(dist, real(x), x2))
        # df(x1)
        # g = Zygote.gradient(df, x1)[1]
        # gn = ngradient(x -> evaluate(dist, x, x2, iters = 50, solver = sinkhorn_log), x1)

        # @test cosdist(g, gn) > 0.9
        # isinteractive() && plot([gn g], layout = 2)
    end

end
# @btime Zygote.gradient(x->real(evaluate(dist,x,$x2)), $x1)
# @profiler Zygote.gradient(x->real(evaluate(dist,x,x2)), x1)
# using DSP
# function callback(x)
#     dump(x)
#     plot(welch_pgram(x))
#     display(plot!(welch_pgram(x2)))
#     false
# end
#
# using Optim
# res = Optim.optimize(
#     df,
#     df',
#     x1,
#     GradientDescent(),
#     Optim.Options(
#         store_trace = true,
#         show_trace = true,
#         show_every = 1,
#         iterations = 40,
#         allow_f_increases = true,
#         time_limit = 100,
#         x_tol = 0,
#         f_tol = 0,
#         g_tol = 1e-8,
#         f_calls_limit = 0,
#         g_calls_limit = 0,
#     ),
#     inplace = false,
# )
#
#
#
# # # Optimize using ADAM
# using Flux
# using Flux: Params, params#, gradient
# X = filt(ones(10), [10], randn(500)) |> v1
# Xp = randn(500) |> v1
#
# opt = ADAM(0.3)
# losses = Float64[]
# P = welch_pgram(X)
# shift1(x) = x .+ x[2]
# # plot3d(fill(0,length(P.power)), log10.(shift1(P.freq)), log10.(P.power), colorbar=false)
#
# lossobj = (
#     ModelDistance(
#         LS(na = 10),
#         OptimalTransportRootDistance(
#             domain = Continuous(),
#             β = 0.1,
#             p = 1,
#             weight = simplex_residueweight,
#         ),
#         # EuclideanRootDistance(domain = Continuous()),
#     ),
#     EnergyDistance(),
# )
#
# using Distances
# # lossobj = (
# #     ModelDistance(
# #         LS(na = 10),
# #         CoefficientDistance(domain = Continuous(), distance=Euclidean()),
# #     ),
# #     EnergyDistance(),
# # )
#
# # lossobj =
# #     ModelDistance(LS(na = 10), OptimalTransportRootDistance(domain = Continuous(), β = 0.1))
#
# loss = (X, Xp) -> evaluate(lossobj, X, Xp)
# gs = gradient(Xp -> real(loss(X, Xp)), Xp)
# loss(X, Xp)
# lossobj[1](X, Xp)
# lossobj[2](X, Xp)
# f1 = plot(P)
# iters = 200
# for i = 1:iters
#     gs = gradient(Xp -> real(loss(X, Xp)), Xp)[1]
#     l = loss(X, Xp)
#     push!(losses, l)
#     @show l
#     Flux.Optimise.update!(opt, Xp, gs)
#     if i % 20 == 0
#         P = welch_pgram(Xp)
#         plot!(f1, P, lab = "", line_z = i) |> display
#     end
# end
# display(current())
# #-
# plot(losses, yscale = :log10)

# gradcheck(f, xs...) =
#   all(isapprox.(ngradient(f, xs...),
#                 Zygote.gradient(f, xs...), rtol = 1e-5, atol = 1e-5))
#
# gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
# gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)
# i = 1

# correct = [
#       [2,3,1],
#       [1, 2, 3],
#       [1,2,3],
#       [2,1,3]
# ]
# for i = 1:3
#       @show i
#       @test Zygote.gradient(v->sort(v)[i], [3.,1,2])[1][correct[1][i]] == 1
#       @test Zygote.gradient(v->sort(v)[i], [1.,2,3])[1][correct[2][i]] == 1
#       @test Zygote.gradient(v->sort(v,by=x->x%10)[i], [11,2,99])[1][correct[3][i]] == 1
#       @test Zygote.gradient(v->sort(v,by=x->x%10)[i], [2,11,99])[1][correct[4][i]] == 1
# end

# a = randn(5)
# p = sortperm(a)
# as = a[p]
# a2 = as[invperm(p)]
# @test a == a2


text:
using Luxor, Test

# test empty drawings

Drawing()
[Drawing() for x in 1:10]
[Drawing(x, x) for x in 1:10]
[Drawing(x, x, :image) for x in 1:10]

p = Point(1, 2)
x, y = p

@test x == 1.0
@test y == 2.0

bb = BoundingBox(p, 2p)

rp = rand(bb)

@test rp in bb

println("...finished basic-test")


text:
using Test
using Checkpointing
using LinearAlgebra
# All tested AD tools
using ForwardDiff, ReverseDiff, Zygote, Enzyme
# Include all the AD tool interfaces through `jacobian()`
include("../examples/adtools.jl")

@testset "Testing Checkpointing.jl" begin
    @testset "Testing Enzyme..." begin
        include("speelpenning.jl")
        errf, errg = main()
        @test isapprox(errf, 0.0; atol = 1e-15)
        @test isapprox(errg, 0.0; atol = 1e-15)
    end
    @testset "Testing Revolve..." begin
        global steps = 50
        global checkpoints = 7
        global verbose = 0
        include("../examples/printaction.jl")
        revolve = main(steps, checkpoints)

        @test revolve.numfwd == 105
        @test revolve.numstore == 28
        @test revolve.numinv == 177
    end

    @testset "Testing Jacobian interface..." begin
        include("jacobian.jl")
    end
    @testset "Test optcontrol" begin
        include("../examples/deprecated/optcontrol.jl")
        @testset "AD Tool $adtool" for adtool in [EnzymeADTool(), ForwardDiffADTool(), ReverseDiffADTool(), ZygoteADTool()]
            @testset "Testing Revolve..." begin
                # Enzyme segfaults if the garbage collector is enabled in this immutable code
                if isa(adtool, EnzymeADTool)
                    GC.gc()
                    GC.enable(false)
                end
                global steps = 100
                global snaps = 3
                global info = 0

                function store(F_H, F_C,t, i)
                    F_C[1,i] = F_H[1]
                    F_C[2,i] = F_H[2]
                    F_C[3,i] = t
                    return
                end

                function restore(F_C, i)
                    F_H = [F_C[1,i], F_C[2,i]]
                    t = F_C[3,i]
                    return F_H, t
                end
                revolve = Revolve(steps, snaps, store, restore; verbose=info)
                F_opt, F_final, L_opt, L = optcontrol(revolve, steps, adtool)
                @test isapprox(F_opt, F_final, rtol=1e-4)
                @test isapprox(L_opt, L, rtol=1e-4)
            end

            @testset "Testing Periodic..." begin
                global steps = 100
                global snaps = 4
                global info = 1

                function store(F_H, F_C,t, i)
                    F_C[1,i] = F_H[1]
                    F_C[2,i] = F_H[2]
                    F_C[3,i] = t
                    return
                end

                function restore(F_C, i)
                    F_H = [F_C[1,i], F_C[2,i]]
                    t = F_C[3,i]
                    return F_H, t
                end
                periodic = Periodic(steps, snaps, store, restore; verbose=info)
                F_opt, F_final, L_opt, L = optcontrol(periodic, steps, adtool)
                @test isapprox(F_opt, F_final, rtol=1e-4)
                @test isapprox(L_opt, L, rtol=1e-4)
            end
        end
    end

    @testset "Test mutable optcontrol" begin
        include("../examples/optcontrol.jl")
        @testset "Testing Revolve..." begin
            global steps = 100
            global snaps = 3
            global info = 0

            revolve = Revolve(steps, snaps; verbose=info)
            F, L, F_opt, L_opt = muoptcontrol(revolve, steps)
            @test isapprox(F_opt, F, rtol=1e-4)
            @test isapprox(L_opt, L, rtol=1e-4)
        end

        @testset "Testing Periodic..." begin
            global steps = 100
            global snaps = 4
            global info = 0

            periodic = Periodic(steps, snaps; verbose=info)
            F, L, F_opt, L_opt = muoptcontrol(periodic, steps)
            @test isapprox(F_opt, F, rtol=1e-4)
            @test isapprox(L_opt, L, rtol=1e-4)
        end
    end
end


text:
module TimeSeriesCausality

# Imports
using Einsum
using FFTW: fft
using Statistics: mean, std, cov as covariance
using LinearAlgebra: det, I as identity

# Includes
include("utils.jl")
include("phase_slope_index.jl")
include("granger_estimation.jl")

# Exports
export psi_est
export granger_est, granger_aic, granger_bic

end # module


text:
# Run tests on functions in src/by_reference_doc/artist.jl

using Spotify.Artists

@testset verbose = true "GET-request endpoints for artists" begin

    # Input arguments from composite types defined in src/types.jl
    artist_id = SpArtistId()

    # Cycle through different input keywords for testing
    countries = ["US", "NL", "DE"]
    offsets = [25, 42]

    @test ~isempty(artist_get(artist_id)[1])

    @testset "For country = $(country_id)" for country_id in countries

        @testset "For offset = $(num_offset)" for num_offset in offsets

            @test ~isempty(artist_get_albums(artist_id, country = country_id, offset = num_offset)[1])

        end

        @test ~isempty(artist_top_tracks(artist_id, country = country_id)[1])

    end

    @test ~isempty(artist_get_related_artists(artist_id)[1])

end


text:
include("preprocessing/preprocessing.jl")
include("postprocessing/postprocessing.jl")


text:
using InstrumentConfig
using Test

@testset "InstrumentConfig.jl" begin
    # Write your tests here.
end


text:
using CCBlade
using Test

@testset "normal operation" begin

# --- verification using example from "Wind Turbine Blade Analysis using the Blade Element Momentum Method"
# --- by Grant Ingram: https://community.dur.ac.uk/g.l.ingram/download/wind_turbine_design.pdf

# Note: There were various problems with the pdf data. Fortunately the excel
# spreadsheet was provided: http://community.dur.ac.uk/g.l.ingram/download.php

# - They didn't actually use the NACA0012 data.  According to the spreadsheet they just used cl = 0.084*alpha with alpha in degrees.
# - There is an error in the spreadsheet where CL at the root is always 0.7.  This is because the classical approach doesn't converge properly at that station.
# - the values for gamma and chord were rounded in the pdf.  The spreadsheet has more precise values.
# - the tip is not actually converged (see significant error for Delta A).  I converged it further using their method.

# --- rotor definition ---
turbine = true
Rhub = 0.01
Rtip = 5.0
Rtip_eff = 5.0*100  # to eliminate tip effects as consistent with their study.
B = 3  # number of blades

# --- section definitions ---

r = [0.2, 1, 2, 3, 4, 5]
gamma = [61.0, 74.31002131, 84.89805553, 89.07195504, 91.25038415, 92.58003871]
theta = (90.0 .- gamma)*pi/180
chord = [0.7, 0.706025153, 0.436187551, 0.304517933, 0.232257636, 0.187279622]

function affunc(alpha, Re, M)

    cl = 0.084*alpha*180/pi

    return cl, 0.0
end

n = length(r)
airfoils = fill(affunc, n)

rotor = Rotor(Rhub, Rtip_eff, B, turbine)
sections = Section.(r, chord, theta, airfoils)

# --- inflow definitions ---

Vinf = 7.0
tsr = 8
Omega = tsr*Vinf/Rtip
rho = 1.0

ops = simple_op.(Vinf, Omega, r, rho)

# --- evaluate ---

out = solve.(Ref(rotor), sections, ops)

ivec = out.phi*180/pi .- theta*180/pi
betavec = 90 .- out.phi*180/pi

# outputs[1] is uncomparable because the classical method fails at the root so they fixed cl

@test isapprox(out.a[2], 0.2443, atol=1e-4)
@test isapprox(out.ap[2], 0.0676, atol=1e-4)
@test isapprox(out.a[3], 0.2497, atol=1e-4)
@test isapprox(out.ap[3], 0.0180, atol=1e-4)
@test isapprox(out.a[4], 0.2533, atol=1e-4)
@test isapprox(out.ap[4], 0.0081, atol=1e-4)
@test isapprox(out.a[5], 0.2556, atol=1e-4)
@test isapprox(out.ap[5], 0.0046, atol=1e-4)
@test isapprox(out.a[6], 0.25725, atol=1e-4)  # note that their spreadsheet is not converged so I ran their method longer.
@test isapprox(out.ap[6], 0.0030, atol=1e-4)

@test isapprox(betavec[2], 66.1354, atol=1e-3)
@test isapprox(betavec[3], 77.0298, atol=1e-3)
@test isapprox(betavec[4], 81.2283, atol=1e-3)
@test isapprox(betavec[5], 83.3961, atol=1e-3)
@test isapprox(betavec[6], 84.7113, atol=1e-3)  # using my more converged solution

#
#

# idx = 6

# aguess = avec[idx]
# apguess = apvec[idx]

# # aguess = 0.2557
# # apguess = 0.0046

# for i = 1:100

#     sigmap = B*chord[idx]/(2*pi*r[idx])
#     lambdar = inflows[idx].Vy/inflows[idx].Vx
#     beta2 = atan(lambdar*(1 + apguess)/(1 - aguess))
#     inc2 = gamma[idx]*pi/180 - beta2
#     cl2, cd2 = affunc(inc2, 1.0, 1.0)
#     global aguess = 1.0/(1 + 4*cos(beta2)^2/(sigmap*cl2*sin(beta2)))
#     global apguess = sigmap*cl2/(4*lambdar*cos(beta2))*(1 - a2)
#     # global apguess = 1.0 / (4*sin(beta2)/(sigmap*cl2) - 1)

# end

# aguess
# apguess

# --------------------------------------------------------------

# -------- verification: propellers.  using script at http://www.aerodynamics4students.com/propulsion/blade-element-propeller-theory.php ------
# I increased their tolerance to 1e-6

# inputs
chord = 0.10
D = 1.6
RPM = 2100
rho = 1.225
pitch = 1.0  # pitch distance in meters.

# --- rotor definition ---
turbine = false
Rhub = 0.0
Rtip = D/2
Rhub_eff = 1e-6  # something small to eliminate hub effects
Rtip_eff = 100.0  # something large to eliminate tip effects
B = 2  # number of blades

# --- section definitions ---

R = D/2.0
r = range(R/10, stop=R, length=11)
theta = atan.(pitch./(2*pi*r))

function affunc(alpha, Re, M)

    cl = 6.2*alpha
    cd = 0.008 - 0.003*cl + 0.01*cl*cl

    return cl, cd
end

n = length(r)
airfoils = fill(affunc, n)

# chord = chord*ones(n)
# rho = rho

rotor_no_F = Rotor(Rhub_eff, Rtip_eff, B, turbine)
rotor = Rotor(Rhub, Rtip, B, turbine)
sections = Section.(r, chord, theta, airfoils)

# --- inflow definitions ---

tsim = 1e3*[1.045361193032356, 1.025630300048415, 1.005234466788998, 0.984163367036026, 0.962407923825140, 0.939960208707079, 0.916813564966455, 0.892962691000145, 0.868403981825492, 0.843134981103815, 0.817154838249790, 0.790463442573673, 0.763063053839278, 0.734956576558370, 0.706148261507327, 0.676643975451150, 0.646450304160057, 0.615575090105131, 0.584027074365864, 0.551815917391907, 0.518952127358381, 0.485446691671386, 0.451311288662196, 0.416557935286392, 0.381199277009438, 0.345247916141561, 0.308716772800348, 0.271618894441869, 0.233967425339051, 0.195775319296371, 0.157055230270717, 0.117820154495231, 0.078082266879117, 0.037854059197644, -0.002852754149850, -0.044026182837742, -0.085655305814570, -0.127728999394140, -0.170237722799272, -0.213169213043848, -0.256515079286031, -0.300266519551194, -0.344414094748869, -0.388949215983616, -0.433863576642539, -0.479150401337354, -0.524801553114807, -0.570810405128802, -0.617169893200684, -0.663873474163182, -0.710915862524620, -0.758291877949762, -0.805995685105502, -0.854022273120508, -0.902366919041604, -0.951025170820984, -0.999992624287163, -1.049265666456123, -1.098840222937414, -1.148712509929845]
qsim = 1e2*[0.803638686218187, 0.806984572453978, 0.809709290183008, 0.811743686838315, 0.813015017103876, 0.813446921530685, 0.812959654049620, 0.811470393912576, 0.808893852696513, 0.805141916379142, 0.800124489784850, 0.793748780791057, 0.785921727832179, 0.776548246109426, 0.765532528164390, 0.752778882688809, 0.738190986274448, 0.721673076180745, 0.703129918771009, 0.682467282681955, 0.659592296506578, 0.634413303042323, 0.606840565246423, 0.576786093366321, 0.544164450503912, 0.508891967461804, 0.470887571011192, 0.430072787279711, 0.386371788290446, 0.339711042057184, 0.290019539402947, 0.237229503458026, 0.181274942660876, 0.122093307308376, 0.059623821454727, -0.006190834182631, -0.075406684829235, -0.148076528546541, -0.224253047813501, -0.303980950928302, -0.387309291734422, -0.474283793689904, -0.564946107631716, -0.659336973911858, -0.757495165410553, -0.859460291551374, -0.965266648683888, -1.074949504731187, -1.188540970723477, -1.306072104649531, -1.427575034895290, -1.553080300508925, -1.682614871422754, -1.816205997296014, -1.953879956474228, -2.095662107769925, -2.241576439746701, -2.391647474158875, -2.545897099743367, -2.704346566395035]

for i = 1:60

    Vinf = float(i)
    Omega = RPM * pi/30

    ops = simple_op.(Vinf, Omega, r, rho)

    # --- evaluate ---

    out = solve.(Ref(rotor_no_F), sections, ops)

    # Np, Tp = loads(outputs)
    # T, Q = thrusttorque(r[1], r, r[end], Np, Tp, B)
    # their spreadsheet did not use trapzezoidal rule, so just do a rect sum.
    T = sum(out.Np*(r[2]-r[1]))*B
    Q = sum(r.*out.Tp*(r[2]-r[1]))*B

    @test isapprox(T, tsim[i], atol=1e-2)  # 2 decimal places
    @test isapprox(Q, qsim[i], atol=2e-3)

end

Vinf = 20.0
Omega = RPM * pi/30

op = simple_op.(Vinf, Omega, r, rho)

out = solve.(Ref(rotor_no_F), sections, op)

T = sum(out.Np*(r[2]-r[1]))*B
Q = sum(r.*out.Tp*(r[2]-r[1]))*B
eff, CT, CQ = nondim(T, Q, Vinf, Omega, rho, rotor)

@test isapprox(CT, 0.056110238632657, atol=1e-7)
@test isapprox(CQ, 0.004337202960642, atol=1e-8)
@test isapprox(eff, 0.735350632777002, atol=1e-6)

# ---------------------------------------------

# ------ camber ----------

alpha0 = -3*pi/180

function affunc(alpha, Re, M)

    cl = 6.2*(alpha - alpha0)
    cd = 0.008 - 0.003*cl + 0.01*cl*cl

    return cl, cd
end

n = length(r)
airfoils = fill(affunc, n)

# theta .-= 3*pi/180

sections = Section.(r, chord, theta, airfoils)

Vinf = 5.0
Omega = RPM * pi/30
ops = simple_op.(Vinf, Omega, r, rho)

out = solve.(Ref(rotor_no_F), sections, ops)

T = sum(out.Np*(r[2]-r[1]))*B
Q = sum(r.*out.Tp*(r[2]-r[1]))*B

@test isapprox(T, 1223.0506862888788, atol=1e-8)
@test isapprox(Q, 113.79919472569034, atol=1e-8)

@test isapprox(out[4].Np, 427.3902632382494, atol=1e-8)
@test isapprox(out[4].Tp, 122.38414345762305, atol=1e-8)
@test isapprox(out[4].a, 2.2845512476210943, atol=1e-8)
@test isapprox(out[4].ap, 0.05024950801920044, atol=1e-8)
@test isapprox(out[4].u, 11.422756238105471, atol=1e-8)
@test isapprox(out[4].v, 3.2709314141649575, atol=1e-8)
@test isapprox(out[4].phi, 0.2596455971546484, atol=1e-8)
@test isapprox(out[4].alpha, 0.23369406105568025, atol=1e-8)
@test isapprox(out[4].W, 63.96697566502531, atol=1e-8)
@test isapprox(out[4].cl, 1.773534419416163, atol=1e-8)
@test isapprox(out[4].cd, 0.03413364011028978, atol=1e-8)
@test isapprox(out[4].cn, 1.7053239640124302, atol=1e-8)
@test isapprox(out[4].ct, 0.48832327407767123, atol=1e-8)
@test isapprox(out[4].F, 1.0, atol=1e-8)
@test isapprox(out[4].G, 1.0, atol=1e-8)

theta = atan.(pitch./(2*pi*r)) .- 3*pi/180

sections = Section.(r, chord, theta, airfoils)
out = solve.(Ref(rotor_no_F), sections, ops)
T = sum(out.Np*(r[2]-r[1]))*B
Q = sum(r.*out.Tp*(r[2]-r[1]))*B

@test isapprox(T, 1e3*0.962407923825140, atol=1e-3)
@test isapprox(Q, 1e2*0.813015017103876, atol=1e-4)

# -----------------------------------------------

end

@testset "qualitative check" begin

# --------- qualitative --------------
# The following tests are only qualitative.  They are not comparisons to known
# empirical data.  Rather they are from the figures used in the documentaiton.  Qualitatively
# the output is about right.  Minor changes in, for example, the airfoil interpolation method
# coudl slightly change the outputs.  The main purpose of these tests is to alert us if something
# significant changes.

Rhub = 1.5
Rtip = 63.0
B = 3
turbine = true
pitch = 0.0
precone = 2.5*pi/180

rotor = Rotor(Rhub, Rtip, B, turbine, pitch, precone)

r = [2.8667, 5.6000, 8.3333, 11.7500, 15.8500, 19.9500, 24.0500,
    28.1500, 32.2500, 36.3500, 40.4500, 44.5500, 48.6500, 52.7500,
    56.1667, 58.9000, 61.6333]
chord = [3.542, 3.854, 4.167, 4.557, 4.652, 4.458, 4.249, 4.007, 3.748,
    3.502, 3.256, 3.010, 2.764, 2.518, 2.313, 2.086, 1.419]
theta = pi/180*[13.308, 13.308, 13.308, 13.308, 11.480, 10.162, 9.011, 7.795,
    6.544, 5.361, 4.188, 3.125, 2.319, 1.526, 0.863, 0.370, 0.106]

# Define airfoils.  In this case we have 8 different airfoils that we load into an array.
# These airfoils are defined in files.
aftypes = Array{Any}(undef, 8)
aftypes[1] = af_from_files("airfoils/Cylinder1.dat")
aftypes[2] = af_from_files("airfoils/Cylinder2.dat")
aftypes[3] = af_from_files("airfoils/DU40_A17.dat")
aftypes[4] = af_from_files("airfoils/DU35_A17.dat")
aftypes[5] = af_from_files("airfoils/DU30_A17.dat")
aftypes[6] = af_from_files("airfoils/DU25_A17.dat")
aftypes[7] = af_from_files("airfoils/DU21_A17.dat")
aftypes[8] = af_from_files("airfoils/NACA64_A17.dat")

# indices correspond to which airfoil is used at which station
af_idx = [1, 1, 2, 3, 4, 4, 5, 6, 6, 7, 7, 8, 8, 8, 8, 8, 8]

# create airfoil array
airfoils = aftypes[af_idx]

sections = Section.(r, chord, theta, airfoils)

# operating point for the turbine
yaw = 0.0*pi/180
tilt = 5.0*pi/180
hubHt = 90.0
shearExp = 0.2

Vinf = 10.0
tsr = 7.55
rotorR = Rtip*cos(precone)
Omega = Vinf*tsr/rotorR
azimuth = 0.0*pi/180
rho = 1.225

op = windturbine_op.(Vinf, Omega, r, precone, yaw, tilt, azimuth, hubHt, shearExp, rho)

out = solve.(Ref(rotor), sections, op)

# plot distributed loads
# figure()
# plot(r/Rtip, out.Np/1e3)
# plot(r/Rtip, out.Tp/1e3)
# xlabel("r/Rtip")
# ylabel("distributed loads (kN/m)")
# legend(["flapwise", "lead-lag"])

nr = length(r)
Npnorm = [0.09718339956327719, 0.13149361678303992, 0.12220741751313423, 1.1634860128761517, 1.7001259694801125, 2.0716257635881257, 2.5120015027019678, 3.1336171133495045, 3.6916824972696465, 4.388661772599469, 5.068896486058948, 5.465165634664408, 6.035059239683594, 6.539134070994739, 6.831387531628286, 6.692665814418597, 4.851568452578296]
Tpnorm = [-0.03321034106737285, -0.08727189682145081, -0.12001897678344217, 0.4696423333976085 , 0.6226256283641799 , 0.6322961942049257 , 0.6474145670774534 , 0.6825021056687035 , 0.6999861694557595 , 0.7218774840801262 , 0.7365515905555542 , 0.7493905698765747 , 0.7529143446199785 , 0.7392483947274653, 0.6981206044592225, 0.614524256128813, 0.40353047553570615]
for i = 1:nr
    @test isapprox.(out.Np[i]/1e3, Npnorm[i], atol=1e-3)
    @test isapprox.(out.Tp[i]/1e3, Tpnorm[i], atol=1e-3)
end

# T, Q = thrusttorque(rotor, sections, out)

# azangles = pi/180*[0.0, 90.0, 180.0, 270.0]
# ops = windturbine_op.(Vinf, Omega, r, precone, yaw, tilt, azangles', hubHt, shearExp, rho)
# outs = solve.(Ref(rotor), sections, ops)

# T, Q = thrusttorque(rotor, sections, outs)

# CP, CT, CQ = nondim(T, Q, Vinf, Omega, rho, rotor)

ntsr = 20  # number of tip-speed ratios
tsrvec = range(2, stop=15, length=ntsr)
cpvec = zeros(ntsr)  # initialize arrays
ctvec = zeros(ntsr)

azangles = pi/180*[0.0, 90.0, 180.0, 270.0]

for i = 1:ntsr
    Omega = Vinf*tsrvec[i]/rotorR

    ops = windturbine_op.(Vinf, Omega, r, precone, yaw, tilt, azangles', hubHt, shearExp, rho)
    outs = solve.(Ref(rotor), sections, ops)
    T, Q = thrusttorque(rotor, sections, outs)

    cpvec[i], ctvec[i], _ = nondim(T, Q, Vinf, Omega, rho, rotor)
end

cpvec_test = [0.02350364982213745, 0.07009444382725231, 0.13891408294580307, 0.21795999362154536, 0.30793657066245234, 0.39220453169574504, 0.43584242088313013, 0.4590001212905916, 0.4695346630018769, 0.4680559467091491, 0.45853627726342805, 0.44346560401350865, 0.4247660950655427, 0.4031517121734798, 0.37850875685858926, 0.3506421939620393, 0.31947670368255254, 0.28492212475211376, 0.24684802187254656, 0.2053516821754716]
ctvec_test = [0.1234606206655405, 0.19135957414241966, 0.2753772801035205, 0.3637320188268945, 0.46087282130376256, 0.5702751788685602, 0.6510892321304749, 0.7117745993566226, 0.7600809643657047, 0.7994008741637466, 0.8332719078069747, 0.8638031244208144, 0.8924359419919573, 0.9199569340127343, 0.9465496984524919, 0.9723679911393518, 0.9975143642860824, 1.0220510627832238, 1.0458478359218666, 1.0685364647098319]

for i = 1:ntsr
    @test isapprox(cpvec[i], cpvec_test[i], atol=1e-3)
    @test isapprox(ctvec[i], ctvec_test[i], atol=1e-3)
end

# rotor definition
Rhub = 0.0254*.5
Rtip = 0.0254*3.0
B = 2  # number of blades
turbine = false

r = .0254*[0.7526, 0.7928, 0.8329, 0.8731, 0.9132, 0.9586, 1.0332,
     1.1128, 1.1925, 1.2722, 1.3519, 1.4316, 1.5114, 1.5911,
     1.6708, 1.7505, 1.8302, 1.9099, 1.9896, 2.0693, 2.1490, 2.2287,
     2.3084, 2.3881, 2.4678, 2.5475, 2.6273, 2.7070, 2.7867, 2.8661, 2.9410]
chord = .0254*[0.6270, 0.6255, 0.6231, 0.6199, 0.6165, 0.6125, 0.6054, 0.5973, 0.5887,
          0.5794, 0.5695, 0.5590, 0.5479, 0.5362, 0.5240, 0.5111, 0.4977,
          0.4836, 0.4689, 0.4537, 0.4379, 0.4214, 0.4044, 0.3867, 0.3685,
          0.3497, 0.3303, 0.3103, 0.2897, 0.2618, 0.1920]
theta = pi/180.0*[40.2273, 38.7657, 37.3913, 36.0981, 34.8803, 33.5899, 31.6400,
                   29.7730, 28.0952, 26.5833, 25.2155, 23.9736, 22.8421, 21.8075,
                   20.8586, 19.9855, 19.1800, 18.4347, 17.7434, 17.1005, 16.5013,
                   15.9417, 15.4179, 14.9266, 14.4650, 14.0306, 13.6210, 13.2343,
                   12.8685, 12.5233, 12.2138]

af = af_from_files("airfoils/NACA64_A17.dat")
airfoils = fill(af, length(r))

rotor = Rotor(Rhub, Rtip, B, turbine)
sections = Section.(r, chord, theta, airfoils)

rho = 1.225
Vinf = 10.0
Omega = 8000.0*pi/30.0

op = simple_op.(Vinf, Omega, r, rho)

outputs = solve.(Ref(rotor), sections, op)

Nptest = [1.8660880922356378, 2.113489633244873, 2.35855792055661, 2.60301402945597, 2.844874233881403, 3.1180230827072126, 3.560077224628854, 4.024057801497014, 4.480574891998562, 4.9279550384928275, 5.366395080074933, 5.79550918136406, 6.21594163851808, 6.622960527391846, 7.017012349498324, 7.3936834781240774, 7.751945902955048, 8.086176029603802, 8.393537672577372, 8.67090062789216, 8.912426896510306, 9.111379449026037, 9.264491105426602, 9.361598738055728, 9.397710628068818, 9.360730779314666, 9.236967116872792, 9.002418776792911, 8.617229305924996, 7.854554211296309, 5.839491141636506]
Tptest = [1.481919153409856, 1.5816880353415623, 1.6702432911163534, 1.7502397903925069, 1.822089134395204, 1.8965254874252537, 2.0022647148294554, 2.097706171361262, 2.178824887386094, 2.2475057498944886, 2.3058616094094666, 2.355253913018444, 2.3970643308370168, 2.4307239254050717, 2.4574034513165794, 2.4763893383410522, 2.488405728268889, 2.492461784055084, 2.4887264544021237, 2.4772963155708783, 2.457435891854637, 2.4282986089025607, 2.3902927838322237, 2.3418848562229155, 2.283388513786012, 2.2134191689454954, 2.130781255778788, 2.0328865955896, 1.9153448642630952, 1.7308451522888118, 1.2736544011110416]
unormtest = [0.1138639166930624, 0.1215785884908478, 0.12836706426605704, 0.13442694075150818, 0.13980334658952898, 0.14527249541450538, 0.15287706861957473, 0.15952130275430465, 0.16497198426904225, 0.16942902209255595, 0.17308482185019125, 0.17607059901193356, 0.17850357997819633, 0.1803781237713692, 0.18177831882512596, 0.18267665167815783, 0.18311924527883217, 0.18305367052760668, 0.182487470134022, 0.1814205780851561, 0.17980424363698078, 0.1775794222894007, 0.1747493664535781, 0.1712044765873724, 0.1669243629724566, 0.16177907188793106, 0.155616714947619, 0.14815634397029886, 0.13888287431637295, 0.12464247057085576, 0.09292645450646951]
vnormtest = [0.09042291182918308, 0.090986677078917, 0.09090479653774426, 0.09038728871285233, 0.08954144817337718, 0.08836143378908702, 0.08598138211326231, 0.08315706129440237, 0.08022297890584436, 0.07727195122065926, 0.07437202068064472, 0.07155384528141737, 0.06883664444997291, 0.0662014244622726, 0.06365995181515205, 0.061184457505937574, 0.05878201223442723, 0.056423985398129595, 0.05410845965501752, 0.05183227774671869, 0.04957767474023305, 0.04732717658478895, 0.04508635659098153, 0.04282827989707323, 0.04055808783300249, 0.03825394697198818, 0.0358976247398962, 0.033456013675480394, 0.030869388594900515, 0.027466462150913407, 0.020268236545135546]
nr = length(r)
for i = 1:nr
    @test isapprox(outputs.Np[i], Nptest[i], atol=1e-3)
    @test isapprox(outputs.Tp[i], Tptest[i], atol=1e-3)
    @test isapprox(outputs.u[i]/Vinf, unormtest[i], atol=1e-3)
    @test isapprox(outputs.v[i]/Vinf, vnormtest[i], atol=1e-3)
end

nJ = 20  # number of advance ratios

J = range(0.1, 0.9, length=nJ)  # advance ratio

Omega = 8000.0*pi/30
n = Omega/(2*pi)
D = 2*Rtip

eff = zeros(nJ)
CT = zeros(nJ)
CQ = zeros(nJ)

for i = 1:nJ
    Vinf = J[i] * D * n

    op = simple_op.(Vinf, Omega, r, rho)
    outputs = solve.(Ref(rotor), sections, op)
    T, Q = thrusttorque(rotor, sections, outputs)
    eff[i], CT[i], CQ[i] = nondim(T, Q, Vinf, Omega, rho, rotor)

end

CTtest = [0.12007272369269596, 0.11685168049381163, 0.113039580103141, 0.1085198115279716, 0.10319003920478352, 0.09725943970861005, 0.09086388629270399, 0.0839976424307642, 0.07671734636145504, 0.06920018426158829, 0.06146448769092457, 0.05351803580676398, 0.04536760940952122, 0.037021384017275026, 0.028444490366496933, 0.019460037067085302, 0.0101482712198453, 0.0009431056296630114, -0.008450506440868786, -0.018336716057292365]
CPtest = [0.04881364095025621 , 0.0495814889974431, 0.05010784517185014, 0.0503081363797726, 0.05016808708430908, 0.04959928236064606, 0.04857913603557237, 0.0470530147170992, 0.045034787733797786, 0.042554367414565766, 0.03958444136790737, 0.036104210007643946, 0.032085430933805684, 0.02750011065614528, 0.022307951002430132, 0.016419276860939747, 0.009846368564156775, 0.002833135813003493, -0.004876593891333574, -0.013591625085158215]
etatest = [0.24598190455626265, 0.3349080300487075, 0.4155652767326253, 0.48818637673414306, 0.5521115225679999, 0.6089123481436948, 0.6595727776885079, 0.7046724703349897, 0.7441662053086512, 0.7788447616541276, 0.8090611349633181, 0.8347808848055981, 0.8558196582739432, 0.8715046719672315, 0.8791362131978436, 0.8670633642311274, 0.7974063895510229, 0.2715632768892098, 0.0, 0.0]

for i = 1:nJ
    @test isapprox(CT[i], CTtest[i], atol=1e-3)
    @test isapprox(CQ[i]*2*pi, CPtest[i], atol=1e-3)
    @test isapprox(eff[i], etatest[i], atol=1e-3)
end

# ------ hover verification --------

# https://rotorcraft.arc.nasa.gov/Publications/files/RamasamyGB_ERF10_836.pdf

chord = 0.060
theta = 0.0
Rtip = 0.656
Rhub = 0.19*Rtip

turbine = false
B = 3

r = range(Rhub + 0.01*Rtip, Rtip - 0.01*Rtip, length=30)

af = af_from_files("naca0012v2.txt")
function af2(alpha, Re, M)
    cl, cd = af(alpha, Re, M)
    return cl, cd+0.014
end
sections = Section.(r, chord, theta, af2)

rho = 1.225
Omega = 800*pi/30
Vinf = 0.0
op = simple_op.(Vinf, Omega, r, rho)

nP = 40
pitch = range(1e-4, 20*pi/180, length=nP)

CT = zeros(nP)
CQ = zeros(nP)

for i = 1:nP
    rotor = Rotor(Rhub, Rtip, B, turbine, pitch[i])
    outputs = solve.(Ref(rotor), sections, op)
    T, Q = thrusttorque(rotor, sections, outputs)

    A = pi*Rtip^2
    CT[i] = T / (rho * A * (Omega * Rtip)^2)
    CQ[i] = Q / (rho * A * (Omega * Rtip)^2 * Rtip)

end

# these are not directly from the experimental data, but have been compared to the experimental data and compare favorably.
# this is more of a regression test on the Vx=0 case
CTcomp = [9.452864991304056e-9, 6.569947366946672e-5, 0.00022338783939012262 , 0.0004420355541809959, 0.0007048495858030926, 0.0010022162314665929, 0.0013268531109981317, 0.0016736995380106938, 0.0020399354072946035, 0.0024223576277264307, 0.0028189858460418893, 0.0032281290309981213, 0.003649357660426685, 0.004081628946214875, 0.004526034348853718, 0.004982651929181267, 0.0054553705714941, 0.005942700094508395, 0.006447634897014323, 0.006963626871239654, 0.007492654931894796, 0.00803866268066438, 0.008597914974368199, 0.009163315934297088, 0.00973817187875574, 0.010309276997090536, 0.010827599471613264, 0.011322361524464346, 0.01180210507896255, 0.012276543435307877, 0.012749323136224754, 0.013223371028562213, 0.013697833731701945, 0.01417556699620018, 0.014646124777465859, 0.015112116772851365, 0.015576452747370885, 0.01602507607909594, 0.016461827164870473, 0.016880126012974343]
CQcomp = [0.000226663607327854, 0.0002270862930229147, 0.0002292742856722754, 0.00023412703235791698 , 0.00024192624628054639 , 0.0002525855612031453, 0.00026638347417704255 , 0.00028314784456601373 , 0.00030299181501156373 , 0.0003259970210015136, 0.00035194661281707764 , 0.00038102864688744595 , 0.0004132249034847219, 0.00044859355432807347 , 0.0004873204055790553, 0.0005293656187218555, 0.0005753409000182888, 0.0006250099998058788, 0.0006788861946930185, 0.0007361096750412038, 0.0007970800153713466, 0.0008624036743669367, 0.0009315051772818803, 0.0010035766105979213, 0.0010791941808362153, 0.0011566643573792704, 0.001229236439467123, 0.0013007334425769355, 0.001372124993921022, 0.0014449961686871802, 0.0015197156782734364, 0.0015967388663224156, 0.0016761210460920718, 0.0017578748614666766, 0.0018409716992061841, 0.0019248522013432586, 0.0020103360819251357, 0.002096387027559033, 0.002182833604491109, 0.0022686470790128036]

for i = 1:nP
    @test isapprox(CT[i], CTcomp[i], atol=1e-4)
    @test isapprox(CQ[i], CQcomp[i], atol=1e-4)
end

end


text:
using Distributions
import Distributions.logpdf
using Parameters

# D is type type of the base distribution
# T is the type of Parameters
# X is the type of observations
export For
struct For # <: Distribution{Multivariate,S} where {T, X, D <: Distribution{V,X} where V <: VariateForm, S <: ValueSupport} # where {A, D <: Distribution{V,A} where V, T, X}
    f   # f(θ) returns a distribution of type D
    θs :: AbstractArray
end

# function For(f, θs)
#     T = eltype(θs)
#     d = f(θs[1])
#     D = typeof(d)
#     X = eltype(d)
#     For(f,θs)
# end

export rand

Base.rand(dist::For) = map(rand, map(dist.f,dist.θs))

# Distributions.logpdf(dist::For, xs) = logpdf.(map(dist.f, dist.θs), xs) |> sum

@inline function Distributions.logpdf(d::For,x::AbstractArray)
    @unpack (f,θs) = d

    s = 0.0
    @inbounds @simd for j in eachindex(x)
        θ = θs[j]
        s += logpdf(f(θ), x[j])
    end
    s
end

# function For(f, js; dist=nothing)
#     @match dist begin
#         fam::ExponentialFamily => ForEF(f, js, fam)
#         x => For(f, js)
#     end
# end


text:
@testset "Normalization" begin
    k = 2
    @test βp(1; k=k) ≈ βp(-1; k=k) atol = 1e-14
    @test βp(4; k=k) ≈ βp(-4; k=k) atol = 1e-14
end

@testset "Circular waves" begin
    k = 2
    @test b̃p(2; k=k)(1,π/2) ≈ +conj(b̃p(-2; k=k)(1,π/2)) atol = 1e-14
    @test b̃p(3; k=k)(1,π/2) ≈ -conj(b̃p(-3; k=k)(1,π/2)) atol = 1e-14
    @test bp(2; k=k)(1,π/2) ≈ +conj(bp(-2; k=k)(1,π/2)) atol = 1e-14
    @test bp(3; k=k)(1,π/2) ≈ -conj(bp(-3; k=k)(1,π/2)) atol = 1e-14
end

max_stars_repo_name: EmileParolin/STrAW

max_issues_repo_name: EmileParolin/STrAW

max_forks_repo_name: EmileParolin/STrAW


text:
using BHTsne
using Base.Test

@test 1 == 1