text:
function blocking_test_prep(obs=ℝ.([1.0, 1.2, 0.8, 1.3, 2.0]),
                            tt=[0.0, 1.0, 1.5, 2.3, 4.0],
                            knots=collect(1:length(obs)-2)[1:1:end],
                            change_pt_buffer=100)
    θ₀ = [10.0, -8.0, 25.0, 0.0, 3.0]
    P˟ = BSI.FitzhughDiffusion(θ₀...)
    P̃ = [BSI.FitzhughDiffusionAux(θ₀..., t₀, u[1], T, v[1]) for (t₀,T,u,v)
         in zip(tt[1:end-1], tt[2:end], obs[1:end-1], obs[2:end])]
    L = @SMatrix [1. 0.]
    Σdiagel = 10^(-10)
    Σ = @SMatrix [Σdiagel]

    Ls = [L for _ in P̃]
    Σs = [Σ for _ in P̃]
    τ(t₀,T) = (x) ->  t₀ + (x-t₀) * (2-(x-t₀)/(T-t₀))
    m = length(obs) - 1
    P = Array{BSI.ContinuousTimeProcess,1}(undef,m)
    dt = 1/50
    for i in m:-1:1
        num_pts = Int64(ceil((tt[i+1]-tt[i])/dt))+1
        t = τ(tt[i], tt[i+1]).( range(tt[i], stop=tt[i+1], length=num_pts) )
        P[i] = ( (i==m) ? BSI.GuidPropBridge(Float64, t, P˟, P̃[i], Ls[i], obs[i+1], Σs[i];
                                         change_pt=BSI.NoChangePt(change_pt_buffer),
                                         solver=BSI.Vern7()) :
                          BSI.GuidPropBridge(Float64, t, P˟, P̃[i], Ls[i], obs[i+1], Σs[i],
                                         P[i+1].H[1], P[i+1].Hν[1], P[i+1].c[1];
                                         change_pt=BSI.NoChangePt(change_pt_buffer),
                                         solver=BSI.Vern7()) )
    end
    create_blocks(BSI.ChequeredBlocking(), P,
                  (knots, 10^(-7), BSI.SimpleChangePt(change_pt_buffer)))
end

@testset "blocking object" begin
    obs = ℝ.([1.0, 1.2, 0.8, 1.3, 2.0])
    tt = [0.0, 1.0, 1.5, 2.3, 4.0]
    L = @SMatrix [1. 0.]
    ϵ = 10^(-7)
    Σ = @SMatrix [10^(-10)]

    blocks = blocking_test_prep(obs, tt)

    @testset "validity of initial set-up" begin
        @test blocks[1].knots == [1, 3]
        @test blocks[2].knots == [2]
        @test blocks[1].blocks == [[1], [2, 3], [4]]
        @test blocks[2].blocks == [[1, 2], [3, 4]]
        @test blocks[1].change_pts == [BSI.SimpleChangePt(100), BSI.NoChangePt(100), BSI.SimpleChangePt(100), BSI.NoChangePt(100)]
        @test blocks[2].change_pts == [BSI.NoChangePt(100), BSI.SimpleChangePt(100), BSI.NoChangePt(100), BSI.NoChangePt(100)]
        @test blocks[1].vs == blocks[2].vs == obs[2:end]
        @test blocks[1].Ls == [I, L, I, L]
        @test blocks[2].Ls == [L, I, L, L]
        @test blocks[1].Σs == [I*ϵ, Σ, I*ϵ, Σ]
        @test blocks[2].Σs == [Σ, I*ϵ, Σ, Σ]
        @test length(blocks[1]) == 3
        @test length(blocks[2]) == 2
        @test length(NoBlocking()) == 1
    end
end


text:
# Testing out LowRankModels in order to test

TrueA = repmat(randn(5,100),20,1); 			# Create a 100x100 matrix of rank 5, then add some noise

A = TrueA + 0.2 * randn(100,100);

u1,s1,v1 = svd(A)		# find svd of matrix and take the top 5 values. u1 * diagm(s1_trunkated) * v1'
						# should be close to A

SVDError = maximum(abs(u1 * diagm([s1[1:5],zeros(95)]) * v1' - TrueA))

using LowRankModels
m,n,k = 100,100,100
losses = fill(quadratic(),n)
rx = onesparse() # each row is assigned to exactly one cluster
ry = zeroreg() # no regularization on the cluster centroids
glrm = GLRM(A,losses,rx,ry,k)

X,Y,ch = fit!(glrm)

ext: jl

max_stars_repo_name: Darf-Ferrara/bayes-scripts

max_issues_repo_name: Darf-Ferrara/bayes-scripts

max_forks_repo_name: Darf-Ferrara/bayes-scripts


text:
using Compat

av_version(v) = VersionNumber(v>>16,(v>>8)&0xff,v&0xff)

have_avcodec()    = Libdl.dlopen_e(libavcodec)    != C_NULL
have_avformat()   = Libdl.dlopen_e(libavformat)   != C_NULL
have_avutil()     = Libdl.dlopen_e(libavutil)     != C_NULL
have_swscale()    = Libdl.dlopen_e(libswscale)    != C_NULL
have_avdevice()   = Libdl.dlopen_e(libavdevice)   != C_NULL
have_avfilter()   = Libdl.dlopen_e(libavfilter)   != C_NULL
#have_avresample() = isdefined(:libavresample) && Libdl.dlopen_e(libavresample) != C_NULL
#have_swresample() = isdefined(:libswresample) && Libdl.dlopen_e(libswresample) != C_NULL

_avcodec_version()    = have_avcodec()    ? av_version(ccall((:avcodec_version,    libavcodec),    UInt32, ())) : v"0"
_avformat_version()   = have_avformat()   ? av_version(ccall((:avformat_version,   libavformat),   UInt32, ())) : v"0"
_avutil_version()     = have_avutil()     ? av_version(ccall((:avutil_version,     libavutil),     UInt32, ())) : v"0"
_swscale_version()    = have_swscale()    ? av_version(ccall((:swscale_version,    libswscale),    UInt32, ())) : v"0"
_avdevice_version()   = have_avdevice()   ? av_version(ccall((:avdevice_version,   libavdevice),   UInt32, ())) : v"0"
_avfilter_version()   = have_avfilter()   ? av_version(ccall((:avfilter_version,   libavfilter),   UInt32, ())) : v"0"
#_avresample_version() = have_avresample() ? av_version(ccall((:avresample_version, libavresample), UInt32, ())) : v"0"
#_swresample_version() = have_swresample() ? av_version(ccall((:swresample_version, libswresample), UInt32, ())) : v"0"

ffmpeg_or_libav = _avutil_version().patch >= 100 ? "ffmpeg" : "libav"

avcodec_dir    = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "AVCodecs",   "v$(_avcodec_version().major)")
avformat_dir   = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "AVFormat",   "v$(_avformat_version().major)")
avutil_dir     = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "AVUtil",     "v$(_avutil_version().major)")
swscale_dir    = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "SWScale",    "v$(_swscale_version().major)")
avdevice_dir   = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "AVDevice",   "v$(_avdevice_version().major)")
avfilter_dir   = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "AVFilters",  "v$(_avfilter_version().major)")
#avresample_dir = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "AVResample", "v$(_avresample_version().major)")
#swresample_dir = joinpath(dirname(@__FILE__), ffmpeg_or_libav, "SWResample", "v$(_swresample_version().major)")

function versioninfo()
    println("Using $ffmpeg_or_libav")
    println("AVCodecs version $(_avcodec_version())")
    println("AVFormat version $(_avformat_version())")
    println("AVUtil version $(_avutil_version())")
    println("SWScale version $(_swscale_version())")
    println("AVDevice version $(_avdevice_version())")
    println("AVFilters version $(_avfilter_version())")
    #println("AVResample version $(_avresample_version())")
    #println("SWResample version $(_swresample_version())")
end


text:
# The MIT License (MIT)
#
# Copyright (c) 2015 Eric Tingey, Brigham Young University
#
# Modified 2018 Kevin Moore, Taylor McDonald BYU
# Refactored 2020 by Tyler Critchfield and Eduardo Alvarez.

module BPM

# Subroutine to use the BPM equations for turbine acoustics
export turbinepos,turbinepos_VAWT

# Default one-third octave band frequencies (Hz)
const default_f = [100.0,125.0,160.0,200.0,250.0,315.0,400.0,500.0,
                    630.0,800.0,1000.0,1250.0,1600.0,2000.0,2500.0,3150.0,
                    4000.0,5000.0,6300.0,8000.0,10000.0,12500.0,16000.0,
                    20000.0 ,25000.0,31500.0,40000.0]

    # A-weighting curve (dBA) for sound perception correction
const default_AdB = [-19.145,-16.190,-13.244,-10.847,-8.675,-6.644,
                        -4.774,-3.248,-1.908,-0.795,0.0,0.576,0.993,1.202,
                        1.271,1.202,0.964,0.556,-0.114,-1.144,-2.488,-4.250,
                        -6.701,-9.341 ,-12.322,-15.694,-19.402]

# cubic spline interpolation setup (for Tip Vortex Noise)
function splineint(n, x, y, xval)
    yval = 0.0
    # assuming the values of x are in accending order
    for i = 1:n
        if (xval < x[i])
            if (i == 2)
                x1 = x[1]
                x2 = x[2]
                x3 = x[3]
                y1 = y[1]
                y2 = y[2]
                y3 = y[3]
                yval = cubspline(x1,x2,x3,y1,y2,y3,xval)
            elseif (i == n)
                x1 = x[n-2]
                x2 = x[n-1]
                x3 = x[n]
                y1 = y[n-2]
                y2 = y[n-1]
                y3 = y[n]
                yval = cubspline(x1,x2,x3,y1,y2,y3,xval)
            else
                if (xval <= (x[i]+x[i-1])/2.0)
                    x1 = x[i-2]
                    x2 = x[i-1]
                    x3 = x[i]
                    y1 = y[i-2]
                    y2 = y[i-1]
                    y3 = y[i]
                    yval = cubspline(x1,x2,x3,y1,y2,y3,xval)
                else
                    x1 = x[i-1]
                    x2 = x[i]
                    x3 = x[i+1]
                    y1 = y[i-1]
                    y2 = y[i]
                    y3 = y[i+1]
                    yval = cubspline(x1,x2,x3,y1,y2,y3,xval)
                end
            end
            break
        elseif (xval == x[i])
            yval = y[i]
            break
        end

    end
    return yval
end

function cubspline(x1, x2, x3, y1, y2, y3, xval)

    a11 = 2.0/(x2-x1)
    a12 = 1.0/(x2-x1)
    a13 = 0.0
    a21 = 1.0/(x2-x1)
    a22 = 2.0*((1.0/(x2-x1))+(1.0/(x3-x2)))
    a23 = 1.0/(x3-x2)
    a31 = 0.0
    a32 = 1.0/(x3-x2)
    a33 = 2.0/(x3-x2)
    b1 = 3.0*(y2-y1)/(x2-x1)^2
    b2 = 3.0*(((y2-y1)/(x2-x1)^2)+((y3-y2)/(x3-x2)^2))
    b3 = 3.0*(y3-y2)/(x3-x2)^2

    bot = a11*a22*a33+a12*a23*a31+a13*a21*a32-a13*a22*a31-a12*a21*a33-a11*a23*a32
    if (xval < x2)
        xtop = b1*a22*a33+a12*a23*b3+a13*b2*a32-a13*a22*b3-a12*b2*a33-b1*a23*a32
        ytop = a11*b2*a33+b1*a23*a31+a13*a21*b3-a13*b2*a31-b1*a21*a33-a11*a23*b3

        k1 = xtop/bot
        k2 = ytop/bot

        a = k1*(x2-x1)-(y2-y1)
        b = -k2*(x2-x1)+(y2-y1)
        t = (xval-x1)/(x2-x1)

        yval = (1.0-t)*y1+t*y2+t*(1.0-t)*(a*(1.0-t)+b*t)
    else
        ytop = a11*b2*a33+b1*a23*a31+a13*a21*b3-a13*b2*a31-b1*a21*a33-a11*a23*b3
        ztop = a11*a22*b3+a12*b2*a31+b1*a21*a32-b1*a22*a31-a12*a21*b3-a11*b2*a32

        k2 = ytop/bot
        k3 = ztop/bot

        a = k2*(x3-x2)-(y3-y2)
        b = -k3*(x3-x2)+(y3-y2)
        t = (xval-x2)/(x3-x2)

        yval = (1.0-t)*y2+t*y3+t*(1.0-t)*(a*(1.0-t)+b*t)
    end
    return yval
end

# Function to compute directivity angles and distance
# Based on work by Luis Vargas (Wind Turbine Noise Prediction)
function direct(n, xt, yt, zt, c, c1, d, Hub, beta)

    theta_e = zeros(n)
    phi_e = zeros(n)
    c2 = zeros(n)
    r = zeros(n)

    # distance from pitch-axis to trailing edge
    c2[1:n] = c[1:n]-c1[1:n]

    # Calculating observer location from hub
    xo = xt # lateral direction
    yo = yt # downstream direction
    zo = zt-Hub # height direction

    for i=1:n
        # Calculating trailing edge position from hub
        xs = sin(beta)*d[i]-cos(beta)*c2[i]
        zs = cos(beta)*d[i]+sin(beta)*c2[i]

        # Calculating observer position from trailing edge
        xe_d = xo-xs
        ze_d = zo-zs

        # Rotating observer position with repsect to beta
        theta = pi-beta
        xe = cos(theta)*xe_d+sin(theta)*ze_d
        ze = -sin(theta)*xe_d+cos(theta)*ze_d

        # Calculating observer distance and directivity angles
        r[i] = sqrt(yo^2+xe^2+ze^2)
        theta_e[i] = atan(sqrt(yo^2+ze^2),xe)
        phi_e[i] = atan(yo,ze)

        # Quadratic smoothing when phi_e is close to 0 or 180 degrees
        if (abs(phi_e[i]) < 5.0*pi/180.0)
            if (phi_e[i] >= 0.0)
                sign = 1
            else
                sign = -1
            end
            phi_er = abs(phi_e[i])*180.0/pi
            phi_er = 0.1*phi_er^2+2.5
            phi_e[i] = sign*phi_er*pi/180.0
        elseif (abs(phi_e[i]) > 175.0*pi/180.0)
            if (phi_e[i] >= 0.0)
                sign = 1
            else
                sign = -1
            end
            phi_er = abs(phi_e[i])*180.0/pi
            phi_er = -0.1*(phi_er-180.0)^2+177.5
            phi_e[i] = sign*phi_er*pi/180.0
        end

    end
    return r,theta_e,phi_e
end #direct

#VAWT Directivty function
function directVAWT(n, xt, yt, zt, c, c1, ht, rad, Hub, rotdir, beta)

    theta_e = zeros(n)
    phi_e = zeros(n)
    c2 = zeros(n)
    r = zeros(n)

    # distance from pitch-axis to trailing edge
    c2[1:n] = c[1:n]-c1[1:n]

    # Calculating observer location from hub
    xo = xt # lateral direction
    yo = yt # downstream direction
    zo = zt-Hub # height direction

    for i=1:n
        # Calculating trailing edge position from hub
        if rotdir >= 0.0
            xs=-rad*cos(beta)-c2[i]*sin(beta)
            ys=-rad*sin(beta)+c2[i]*cos(beta)
            zs=ht[i]
        else
            xs=-rad*cos(beta)+c2[i]*sin(beta)
            ys=-rad*sin(beta)-c2[i]*cos(beta)
            zs=ht[i]
        end

        # Calculating observer position from trailing edge
        xe_d = xo-xs
        ye_d = yo-ys
        ze = zo-zs

        # Rotating observer position with repsect to beta
        theta = rotdir*(pi/2.0)+beta
        xe = cos(theta)*xe_d+sin(theta)*ye_d
        ye = -sin(theta)*xe_d+cos(theta)*ye_d

        # Calculating observer distance and directivity angles
        r[i] = sqrt(xe^2+ye^2+ze^2)
        theta_e[i] = atan(sqrt(ye^2+ze^2),xe)
        phi_e[i] = atan(ye,ze)
        # Quadratic smoothing when theta_e or phi_e are close to 0, +/-180 degrees
        if (abs(theta_e[i])< (5.0*pi/180.0))
            if (theta_e[i] >= 0.0)
                sign = 1.0
            else
                sign = -1.0
            end
            theta_er=abs(theta_e[i])*180.0/pi
            theta_er=0.1*(theta_er)^2+2.5
            theta_e[i] = sign*theta_er*pi/180.0
        elseif (abs(theta_e[i])> (175.0*pi/180))
            if (theta_e[i] >=0 )
                sign =1.0
            else
                sign=-1.0
            end
            theta_er = abs(theta_e[i])*180.0/pi
            theta_er = -0.1*(theta_er-180.0)^2+177.5
            theta_e[i] = sign*theta_er*pi/180.0
        end

        if (abs(phi_e[i]) < 5.0*pi/180.0)
            if (phi_e[i] >= 0.0)
                sign = 1.0
            else
                sign = -1.0
            end
            phi_er = abs(phi_e[i])*180.0/pi
            phi_er = 0.1*phi_er^2+2.5
            phi_e[i] = sign*phi_er*pi/180.0
        elseif (abs(phi_e[i]) > 175.0*pi/180.0)
            if (phi_e[i] >= 0.0)
                sign = 1
            else
                sign = -1
            end
            phi_er = abs(phi_e[i])*180.0/pi
            phi_er = -0.1*(phi_er-180.0)^2+177.5
            phi_e[i] = sign*phi_er*pi/180.0
        end
    end
    return r,theta_e,phi_e
end #directVAWT

# Directivity function for high-frequency noise
# not for high-angle separation; becomes inaccurate for theta_e approaching 180 deg
function Dhfunc(theta_e, phi_e,M)
    conv = 0.8 # convection factor for speed
    Mc = M*conv

    Dh = (2.0*(sin(theta_e/2.0))^2*(sin(phi_e))^2)/((1.0+M*cos(theta_e))*(1.0+(M-Mc)*cos(theta_e))^2)
    return Dh
end

# Directivity function for low-frequency noise
function Dlfunc(theta_e, phi_e,M)
    Dl = ((sin(theta_e))^2*(sin(phi_e))^2)/(1.0+M*cos(theta_e))^4
    return Dl
end #Dlfunc

# Spectral Function A
function Afunc(ain, Re)
    a = abs(log10(ain))

    # Calculating Amin
    if (a < 0.204)
        Amin = sqrt(67.552-886.788*a^2)-8.219
    elseif (a >= 0.204 && a <= 0.244)
        Amin = -32.665*a+3.981
    else
        Amin = -142.795*a^3+103.656*a^2-57.757*a+6.006
    end

    # Calculating Amax
    if (a < 0.13)
        Amax = sqrt(67.552-886.788*a^2)-8.219
    elseif (a >= 0.13 && a <= 0.321)
        Amax = -15.901*a+1.098
    else
        Amax = -4.669*a^3+3.491*a^2-16.699*a+1.149
    end

    # Calculating a0
    if (Re < 9.52e4)
        a0 = 0.57
    elseif (Re >= 9.52e4 && Re <= 8.57e5)
        a0 = -9.57e-13*(Re-8.57e5)^2+1.13
    else
        a0 = 1.13
    end

    # Calculating Amin(a0)
    if (a0 < 0.204)
        Amin0 = sqrt(67.552-886.788*a0^2)-8.219
    elseif (a0 >= 0.204 && a0 <= 0.244)
        Amin0 = -32.665*a0+3.981
    else
        Amin0 = -142.795*a0^3+103.656*a0^2-57.757*a0+6.006
    end

    # Calculating Amax(a0)
    if (a0 < 0.13)
        Amax0 = sqrt(67.552-886.788*a0^2)-8.219
    elseif (a0 >= 0.13 && a0 <= 0.321)
        Amax0 = -15.901*a0+1.098
    else
        Amax0 = -4.669*a0^3+3.491*a0^2-16.699*a0+1.149
    end

    AR0 = (-20.0-Amin0)/(Amax0-Amin0)

    Aspec = Amin+AR0*(Amax-Amin)
    return Aspec
end #Afunc

# Spectral Function B
function Bfunc(bin, Re)

    b = abs(log10(bin))

    # Calculating Bmin
    if (b < 0.13)
        Bmin = sqrt(16.888-886.788*b^2)-4.109
    elseif (b >= 0.13 && b <= 0.145)
        Bmin = -83.607*b+8.138
    else
        Bmin = -817.810*b^3+355.210*b^2-135.024*b+10.619
    end

    # Calculating Bmax
    if (b < 0.10)
        Bmax = sqrt(16.888-886.788*b^2)-4.109
    elseif (b >= 0.10 && b <= 0.187)
        Bmax = -31.330*b+1.854
    else
        Bmax = -80.541*b^3+44.174*b^2-39.381*b+2.344
    end

    # Calculating b0
    if (Re < 9.52e4)
        b0 = 0.30
    elseif (Re >= 9.52e4 && Re <= 8.57e5)
        b0 = -4.48e-13*(Re-8.57e5)^2+0.56
    else
        b0 = 0.56
    end

    # Calculating Bmin(b0)
    if (b0 < 0.13)
        Bmin0 = sqrt(16.888-886.788*b0^2)-4.109
    elseif (b0 >= 0.13 && b0 <= 0.145)
        Bmin0 = -83.607*b0+8.138
    else
        Bmin0 = -817.810*b0^3+355.210*b0^2-135.024*b0+10.619
    end

    # Calculating Bmax(b0)
    if (b0 < 0.10)
        Bmax0 = sqrt(16.888-886.788*b0^2)-4.109
    elseif (b0 >= 0.10 && b0 <= 0.187)
        Bmax0 = -31.330*b0+1.854
    else
        Bmax0 = -80.541*b0^3+44.174*b0^2-39.381*b0+2.344
    end

    BR0 = (-20.0-Bmin0)/(Bmax0-Bmin0)

    Bspec =  Bmin+BR0*(Bmax-Bmin)
    return Bspec
end #Bfunc

function G1func(e)
    if (e <= 0.5974)
        G1 = 39.8*log10(e)-11.12
    elseif (e <= 0.8545 && e > 0.5974)
        G1 = 98.409*log10(e)+2.0
    elseif (e <= 1.17 && e > 0.8545)
        G1 = sqrt(2.484-506.25*(log10(e))^2)-5.076
    elseif (e <= 1.674 && e > 1.17)
        G1 = -98.409*log10(e)+2.0
    else
        G1 = -39.8*log10(e)-11.12
    end
    return G1
end #G1func

function G2func(d)
    if (d <= 0.3237)
        G2 = 77.852*log10(d)+15.328
    elseif (d <= 0.5689 && d > 0.3237)
        G2 = 65.188*log10(d)+9.125
    elseif (d <= 1.7579 && d > 0.5689)
        G2 = -114.052*(log10(d))^2
    elseif (d <= 3.0889 && d > 1.7579)
        G2 = -65.188*log10(d)+9.125
    else
        G2 = -77.852*log10(d)+15.328
    end
    return G2
end #G2func

function G3func(alpha)
    G3 = 171.04-3.03*alpha
    return G3
end #G3func

function G4func(hdav, psi)
    if (hdav <= 5.0)
        G4 = 17.5*log10(hdav)+157.5-1.114*psi
    else
        G4 = 169.7-1.114*psi
    end
    return G4
end #G4func

function G5func(hdav, psi, StSt_peak)
    # finding G5 at phi = 14 deg
    eta = log10(StSt_peak)

    if (hdav < 0.25)
        mu = 0.1221
    elseif (hdav < 0.62 && hdav >= 0.25)
        mu = -0.2175*hdav+0.1755
    elseif (hdav < 1.15 && hdav >= 0.62)
        mu = -0.0308*hdav+0.0596
    else
        mu = 0.0242
    end

    if (hdav <= 0.02)
        m = 0.0
    elseif (hdav <= 0.5 && hdav > 0.02)
        m = 68.724*(hdav)-1.35
    elseif (hdav <= 0.62 && hdav > 0.5)
        m = 308.475*hdav-121.23
    elseif (hdav <= 1.15 && hdav > 0.62)
        m = 224.811*hdav-69.35
    elseif (hdav <= 1.2 && hdav > 1.15)
        m = 1583.28*hdav-1631.59
    else
        m = 268.344
    end

    eta_0 = -sqrt((m^2*mu^4)/(6.25+m^2*mu^2))
    k = 2.5*sqrt(1.0-(eta_0/mu)^2)-2.5-m*eta_0

    if (eta < eta_0)
        G14 = m*eta+k
    elseif (eta < 0.0 && eta >= eta_0)
        G14 = 2.5*sqrt(1.0-(eta/mu)^2)-2.5
    elseif (eta < 0.03616 && eta >= 0.0)
        G14 = sqrt(1.5625-1194.99*eta^2)-1.25
    else
        G14 = -155.543*eta+4.375
    end

    # finding G5 at psi = 0 deg
    hdav_prime = 6.724*hdav^2-4.019*hdav+1.107

    if (hdav_prime < 0.25)
        mu0 = 0.1221
    elseif (hdav_prime < 0.62 && hdav_prime >= 0.25)
        mu0 = -0.2175*hdav_prime+0.1755
    elseif (hdav_prime < 1.15 && hdav_prime >= 0.62)
        mu0 = -0.0308*hdav_prime+0.0596
    else
        mu0 = 0.0242
    end

    if (hdav_prime <= 0.02)
        m0 = 0.0
    elseif (hdav_prime <= 0.5 && hdav_prime > 0.02)
        m0 = 68.724*hdav_prime-1.35
    elseif (hdav_prime <= 0.62 && hdav_prime > 0.5)
        m0 = 308.475*hdav_prime-121.23
    elseif (hdav_prime <= 1.15 && hdav_prime > 0.62)
        m0 = 224.811*hdav_prime-69.35
    elseif (hdav_prime <= 1.2 && hdav_prime > 1.15)
        m0 = 1583.28*hdav_prime-1631.59
    else
        m0 = 268.344
    end

    eta_00 = -sqrt((m0^2*mu0^4)/(6.25+m0^2*mu0^2))
    k0 = 2.5*sqrt(1.0-(eta_00/mu0)^2)-2.5-m0*eta_00

    if (eta < eta_00)
        G0 = m0*eta+k0
    elseif (eta < 0.0 && eta >= eta_00)
        G0 = 2.5*sqrt(1.0-(eta/mu0)^2)-2.5
    elseif (eta < 0.03616 && eta >= 0.0)
        G0 = sqrt(1.5625-1194.99*eta^2)-1.25
    else
        G0 = -155.543*eta+4.375
    end

    G5 = G0+0.0714*psi*(G14-G0)
    return G5
end #G5func

# Turbulent Boundary Layer Trailing Edge Noise
function TBLTEfunc(f, V, L, c, r, theta_e, phi_e, alpha, nu, c0, trip)
    # constants
    M = V/c0
    Re = (V*c)/nu

    if (trip == false)
        # UNTRIPPED boundary layer at 0 deg- thickness, displacement thickness
        d0 = c*(10.0^(1.6569-0.9045*log10(Re)+0.0596*(log10(Re))^2))
        d0_d = c*(10.0^(3.0187-1.5397*log10(Re)+0.1059*(log10(Re))^2))
    else
        # TRIPPED boundary layer at 0 deg- thickness, displacement thickness
        d0 = c*(10.0^(1.892-0.9045*log10(Re)+0.0596*(log10(Re))^2))
        if (Re <= 0.3e6)
            d0_d = c*0.0601*Re^(-0.114)
        else
            d0_d = c*(10.0^(3.411-1.5397*log10(Re)+0.1059*(log10(Re))^2))
        end
    end

    # boundary layer on pressure side- thickness, displacement thickness
    dpr = d0*(10.0^(-0.04175*alpha+0.00106*alpha^2))
    dp_d = d0_d*(10.0^(-0.0432*alpha+0.00113*alpha^2))

    if (trip == false)
        # UNTRIPPED boundary layer on suction side- displacement thickness
        if (alpha <= 7.5 && alpha >= 0.0)
            ds_d = d0_d*10.0^(0.0679*alpha)
        elseif (alpha <= 12.5 && alpha > 7.5)
            ds_d = d0_d*0.0162*10.0^(0.3066*alpha)
            #elseif (alpha <= 25.0 && alpha > 12.5)
        else
            ds_d = d0_d*52.42*10.0^(0.0258*alpha)
        end
    else
        # TRIPPED boundary layer on suction side- displacement thickness
        if (alpha <= 5.0 && alpha >= 0.0)
            ds_d = d0_d*10.0^(0.0679*alpha)
        elseif (alpha <= 12.5 && alpha > 5.0)
            ds_d = d0_d*0.381*10.0^(0.1516*alpha)
            #elseif (alpha <= 25.0 && alpha > 12.5)
        else
            ds_d = d0_d*14.296*10.0^(0.0258*alpha)
        end
    end

    Dh = Dhfunc(theta_e,phi_e,M)
    Dl = Dlfunc(theta_e,phi_e,M)

    Stp = (f*dp_d)/V
    Sts = (f*ds_d)/V

    St1 = 0.02*M^(-0.6)

    if (alpha < 1.33)
        St2 = St1*1.0
    elseif (alpha <= 12.5 && alpha >= 1.33)
        St2 = St1*10.0^(0.0054*(alpha-1.33)^2)
    else
        St2 = St1*4.72
    end

    St_bar = (St1+St2)/2.0

    St_peak = max(St1,St2,St_bar)

    apre = Stp/St1
    asuc = Sts/St1
    bang = Sts/St2

    gamma = 27.094*M+3.31
    gamma0 = 23.43*M+4.651
    beta = 72.65*M+10.74
    beta0 = -34.19*M-13.82

    if (Re < 2.47e5)
        K1 = -4.31*log10(Re)+156.3
    elseif (Re >= 2.47e5 && Re <= 8.0e5)
        K1 = -9.0*log10(Re)+181.6
    else
        K1 = 128.5
    end

    if (alpha < (gamma0-gamma))
        K2 = K1-1000.0
    elseif (alpha >= (gamma0-gamma) && alpha <= (gamma0+gamma))
        K2 = K1+sqrt(beta^2-(beta/gamma)^2*(alpha-gamma0)^2)+beta0
    else
        K2 = K1-12.0
    end

    Re_dp = (V*dp_d)/nu

    if (Re_dp <= 5000.0)
        DeltaK1 = alpha*(1.43*log10(Re_dp)-5.29)
    else
        DeltaK1 = 0.0
    end

    # Keeping observer distance from getting too close to the turbine
    if (r < 1e-8)
        rc = 1e-8
    else
        rc = r
    end

    if (alpha > 12.5 || alpha > gamma0)
        # Turbulent Boundary Layer Separation Stall Noise (TBLSS); this is were the airfoil is stalling and stall noise dominates
        # SPLp = -infinity; 10^(SPLp/10) = 0
        # SPLs = -infinity; 10^(SPLs/10) = 0

        A = Afunc(bang,3.0*Re)

        SPLa = 10.0*log10((ds_d*M^5*L*Dl)/rc^2)+A+K2

        TBLTE = 10.0*log10(10.0^(SPLa/10.0))

    else
        Ap = Afunc(apre,Re)
        As = Afunc(asuc,Re)
        B = Bfunc(bang,Re)

        SPLp = 10.0*log10((dp_d*M^5*L*Dh)/rc^2)+Ap+(K1-3.0)+DeltaK1
        SPLs = 10.0*log10((ds_d*M^5*L*Dh)/rc^2)+As+(K1-3.0)
        SPLa = 10.0*log10((ds_d*M^5*L*Dh)/rc^2)+B+K2

        TBLTE =  10.0*log10(10.0^(SPLp/10.0)+10.0^(SPLs/10.0)+10.0^(SPLa/10.0))

    end
    return TBLTE
end #TBLTEfunc

# Turbulent Boundary Layer Tip Vortex Noise
function TBLTVfunc(f, V, c, r, theta_e, phi_e, atip, c0, tipflat, AR)

    # constants
    M = V/c0
    Mmax = M*(1.0+0.036*atip)

    Dh = Dhfunc(theta_e,phi_e,M)

    # Tip vortex noise correction based on data from "Airfoil Tip Vortex Formation Noise"
    AR_data = [2.0,2.67,4.0,6.0,12.0,24.0]
    atipcorr_data = [0.54,0.62,0.71,0.79,0.89,0.95]

    if ((AR >= 2.0) && (AR <= 24.0))
        atipcorr = splineint(6,AR_data,atipcorr_data,AR)
    elseif (AR > 24.0)
        atipcorr = 1.0
    else
        atipcorr = 0.5
    end

    atip_d = atip*atipcorr

    if (tipflat == false)
        # rounded tip
        l = 0.008*c*atip_d
    else
        # flat tip
        if (atip_d <= 2.0 && atip_d >= 0.0)
            l = c*(0.0230+0.0169*atip_d)
        else
            l = c*(0.0378+0.0095*atip_d)
        end
    end

    St = (f*l)/(V*(1.0+0.036*atip_d))

    # Keeping observer distance from getting too close to the turbine
    if (r < 1e-8)
        rc = 1e-8
    else
        rc = r
    end

    TBLTV = 10.0*log10((M^2*Mmax^3*l^2*Dh)/rc^2)-30.5*(log10(St)+0.3)^2+126.0
    return TBLTV
end #TBLTVfunc

# Laminar Boundary Layer Vortex Shedding
function LBLVSfunc(f, V, L, c, r, theta_e, phi_e, alpha, nu, c0, trip)

    # constants
    M = V/c0
    Re = (V*c)/nu

    if (trip == false)
        # UNTRIPPED boundary layer at 0 deg- thickness
        d0 = c*(10.0^(1.6569-0.9045*log10(Re)+0.0596*(log10(Re))^2))
    else
        # TRIPPED boundary layer at 0 deg- thickness
        d0 = c*(10.0^(1.892-0.9045*log10(Re)+0.0596*(log10(Re))^2))
    end
    # boundary layer on pressure side- thickness
    dpr = d0*(10.0^(-0.04175*alpha+0.00106*alpha^2))

    St = (f*dpr)/V

    Dh = Dhfunc(theta_e,phi_e,M)

    if (Re <= 1.3e5)
        St1 = 0.18
    elseif (Re <= 4.0e5 && Re > 1.3e5)
        St1 = 0.001756*Re^0.3931
    else
        St1 = 0.28
    end

    St_peak = St1*10.0^(-0.04*alpha)

    e = St/St_peak

    G1 = G1func(e)

    if (alpha <= 3.0)
        Re0 = 10.0^(0.215*alpha+4.978)
    else
        Re0 = 10.0^(0.12*alpha+5.263)
    end

    d = Re/Re0

    G2 = G2func(d)
    G3 = G3func(alpha)

    # Keeping observer distance from getting too close to the turbine
    if (r < 1e-8)
        rc = 1e-8
    else
        rc = r
    end

    LBLVS = 10.0*log10((dpr*M^5*L*Dh)/rc^2)+G1+G2+G3
    return LBLVS
end #LBLVSfunc

# Trailing Edge Bluntness Vortex Shedding Noise
function TEBVSfunc(f, V, L, c, h, r, theta_e, phi_e, alpha, nu, c0, psi, trip)

    # constants
    M = V/c0
    Re = (V*c)/nu

    if (trip == false)
        # UNTRIPPED boundary layer at 0 deg- thickness, displacement thickness
        d0 = c*(10.0^(1.6569-0.9045*log10(Re)+0.0596*(log10(Re))^2))
        d0_d = c*(10.0^(3.0187-1.5397*log10(Re)+0.1059*(log10(Re))^2))
    else
        # TRIPPED boundary layer at 0 deg- thickness, displacement thickness
        d0 = c*(10.0^(1.892-0.9045*log10(Re)+0.0596*(log10(Re))^2))
        if (Re <= 0.3e6)
            d0_d = c*0.0601*Re^(-0.114)
        else
            d0_d = c*(10.0^(3.411-1.5397*log10(Re)+0.1059*(log10(Re))^2))
        end
    end

    # boundary layer on pressure side- thickness, displacement thickness
    dpr = d0*(10.0^(-0.04175*alpha+0.00106*alpha^2))
    dp_d = d0_d*(10.0^(-0.0432*alpha+0.00113*alpha^2))

    if (trip == false)
        # UNTRIPPED boundary layer on suction side- displacement thickness
        if (alpha <= 7.5 && alpha >= 0.0)
            ds_d = d0_d*10.0^(0.0679*alpha)
        elseif (alpha <= 12.5 && alpha > 7.5)
            ds_d = d0_d*0.0162*10.0^(0.3066*alpha)
            #elseif (alpha <= 25 && alpha > 12.5)
        else
            ds_d = d0_d* 52.42* 10.0^(0.0258*alpha)
        end
    else
        # TRIPPED boundary layer on suction side- displacement thickness
        if (alpha <= 5.0 && alpha >= 0.0)
            ds_d = d0_d*10.0^(0.0679*alpha)
        elseif (alpha <= 12.5 && alpha > 5.0)
            ds_d = d0_d*0.381*10.0^(0.1516*alpha)
            #elseif (alpha <= 25.0 && alpha > 12.5)
        else
            ds_d = d0_d*14.296*10.0^(0.0258*alpha)
        end
    end

    Dh = Dhfunc(theta_e,phi_e,M)
    St = (f*h)/V
    dav = (dp_d+ds_d)/2.0

    hdav = h/dav

    if (hdav >= 0.2)
        St_peak = (0.212-0.0045*psi)/(1.0+0.235*(1.0/hdav)-0.0132*(1.0/hdav)^2)
    else
        St_peak = 0.1*(hdav)+0.095-0.00243*psi
    end

    StSt_peak = St/St_peak

    G4 = G4func(hdav,psi)
    G5 = G5func(hdav,psi,StSt_peak)

    # Keeping observer distance from getting too close to the turbine
    if (r < 1e-8)
        rc = 1e-8
    else
        rc = r
    end

    TEBVS = 10.0*log10((h*M^(5.5)*L*Dh)/rc^2)+G4+G5
    return TEBVS
end #TEBVSfunc

# Computing the overall sound pressure level (OASPL) of a turbine defined below (in dB)
function OASPL(ox, oy, oz, windvel, rpm, B, Hub, rad, c, c1, alpha, nu, c0, psi, AR; f=default_f, AdB=default_AdB)
    # constants

    nf = length(f)
    bf = 3

    n = length(rad)

    L = zeros(n-1)
    d = zeros(n-1)
    V = zeros(n-1)
    h = zeros(n-1)
    TV_t = zeros(B)
    TE_t = zeros((n-1)*B)
    BLVS_t = zeros((n-1)*B)
    BVS_t = zeros((n-1)*B)

    TE = zeros(nf)
    TV = zeros(nf)
    BLVS = zeros(nf)
    BVS = zeros(nf)
    SPLf = zeros(bf, nf)
    SPLfA = zeros(bf, nf)
    SPLo_d = zeros(bf)
    SPLoa_d = zeros(bf)

    # Using untripped or tripped boundary layer specficiation
    # trip = false # untripped
    trip = true # tripped

    # Tip specfication
    tipflat = false # round
    # tipflat = true # flat

    # Parameters of the wind turbine
    omega = (rpm*2.0*pi)/60.0  # angular velocity (rad/sec)

    for i = 1:n-1
        L[i] = rad[i+1]-rad[i] # length of each radial section (m)
        d[i] = rad[i] # radial section to be used in directivity calculations (m)
        V[i] = sqrt((omega*rad[i])^2+windvel^2) # wind speed over the blade (m/s)
    end

    h[1:n-1] = 0.01*c[1:n-1]  # trailing edge thickness; 1% of chord length (m)
    atip = alpha[n-1]  # angle of attack of the tip region (deg)

    # Blade rotation increments to rotate around (45 deg from Vargas paper)
    # beta = [0.0,0.25*pi,0.5*pi,0.75*pi,pi,1.25*pi,1.5*pi,1.75*pi] # 8 increments
    beta = [0.0,2.0*pi/9.0,4.0*pi/9.0] # 3 increments (equivalent of 9 for 3 blades)
    # beta = [0.0,pi] # 2 increments
    # beta = [0.0] # 1 increment (top blade facing straight up)

    B_int = 2.0*pi/B # Intervals between blades (from the first blade at 0 deg)

    for di=1:bf # for each rotation increment
        for j=1:nf # for each frequency
            for bi=1:B # for each blade
                # Calcuating observer distances and directivty angles for the given blade orientation
                r,theta_e,phi_e = direct(n-1,ox,oy,oz,c,c1,d,Hub,beta[di]+(bi-1)*B_int)

                TBLTV = TBLTVfunc(f[j],V[n-1],c[n-1],r[n-1],theta_e[n-1],phi_e[n-1],atip,c0,
                tipflat,AR)
                TV_t[bi] = TBLTV
                for k=1:n-1
                    # Calculating sound pressure level (dB) for each noise source at each radial position
                    TBLTE = TBLTEfunc(f[j],V[k],L[k],c[k],r[k],theta_e[k],phi_e[k],alpha[k],
                    nu,c0,trip)

                    if (trip == false)
                        LBLVS = LBLVSfunc(f[j],V[k],L[k],c[k],r[k],theta_e[k],phi_e[k],alpha[k],
                        nu,c0,trip)
                    else
                        LBLVS = 0.0
                    end
                    TEBVS = TEBVSfunc(f[j],V[k],L[k],c[k],h[k],r[k],theta_e[k],phi_e[k],
                    alpha[k],nu,c0,psi,trip)

                    # Assigning noise to blade segment
                    TE_t[k.+[n-1]*(bi-1)] .= TBLTE
                    BLVS_t[k.+[n-1]*(bi-1)] .= LBLVS
                    BVS_t[k.+[n-1]*(bi-1)] .= TEBVS
                end
            end

            # Adding sound pressure levels (dB)
            TE[j] = 10.0*log10(sum(10.0.^(TE_t/10.0)))
            TV[j] = 10.0*log10(sum(10.0.^(TV_t/10.0)))
            BLVS[j] = 10.0*log10(sum(10.0.^(BLVS_t/10.0)))
            BVS[j] = 10.0*log10(sum(10.0.^(BVS_t/10.0)))

            # Combining noise sources into overall SPL
            SPLf[di,j] = 10.0*log10(10.0^(TE[j]/10.0)+10.0^(TV[j]/10.0)+10.0^(BLVS[j]/10.0)+10.0^(BVS[j]/10.0))
        end

        # Correcting with A-weighting
        # SPLf[1:nf] = SPLf[1:nf]+AdB[1:nf]
        SPLfA[di,1:nf] = SPLf[di,1:nf]+AdB[1:nf]

        # Adding SPLs for each rotation increment
        SPLo_d[di] = 10.0*log10(sum(10.0.^(SPLf[di,:]/10.0))) #non-weighted OASPL per increment
        SPLoa_d[di] = 10.0*log10(sum(10.0.^(SPLfA[di,:]/10.0))) #A-weighted OASPL per increment

        # Protecting total calcuation from negative SPL values
        if (SPLoa_d[di] < 0.0)
            SPLoa_d[di] = 0.0
        end
        if (SPLo_d[di] < 0.0)
            SPLo_d[di] = 0.0
        end
    end

    # Performing root mean square calculation of SPLs at rotation increments for final value
    SPLo = sqrt(sum(SPLo_d.^2)/bf)
    SPLoa = sqrt(sum(SPLoa_d.^2)/bf)
    SPLf_overall = sqrt.(sum(SPLf.^2; dims=1)/bf)
    SPLfA_overall = sqrt.(sum(SPLfA.^2; dims=1)/bf)

    return SPLo, SPLoa, SPLf_overall, SPLfA_overall #OASPL, A-OASPL, SPL spectrum, A-SPL spectrum
end #OASPL

#Computing the overall sound pressure level (OASPL) of a turbine defined below (in dB)
function OASPLVAWT(p, ox, oy, oz, B, Hub, high, rad, c, c1, alpha, nu,
                    c0, psi, rot, Vinf, wakex, wakey, AR)

   # constants
   nf = 27
   bf = 8

   n = length(high)
   L = zeros(n-1)
   d = zeros(n-1)
   V = zeros(n-1)
   h = zeros(n-1)
   TV_t = zeros(2*B)
   TE_t = zeros((n-1)*B)
   BLVS_t = zeros((n-1)*B)
   BVS_t = zeros((n-1)*B)

   TE = zeros(nf)
   TV = zeros(nf)
   BLVS = zeros(nf)
   BVS = zeros(nf)
   SPLf = zeros(nf)
   SPLoa_d = zeros(bf)
   theta_vel = zeros(p)
   highmid = zeros(n-1)

   # Using untripped or tripped boundary layer specficiation
   trip = false # untripped
   # trip = true # tripped

   # Tip specfication
   tipflat = false # round
   # tipflat = true # flat
   for i = 1:n-1
       L[i] = high[i+1]-high[i] # length of each height section (m)
       highmid[i] = (high[i+1]+high[i])/2.0
   end
   h[1:n-1] = 0.01*c[1:n-1]
   atip1 = alpha[1] #angle of attack of the tip region on bottom (deg)
   atip2= alpha[n-1] #angle of attack of the tip region on top (deg)
   if (rot >= 0.0)
       rotdir = 1.0
   else
       rotdir = -1.0
   end

   # Blade rotation increments to rotate around (45 deg from Vargas paper)
   beta = [0.0,0.25*pi,0.5*pi,0.75*pi,pi,1.25*pi,1.5*pi,1.75*pi] # 8 increments
   #beta = [0.0,2.0*pi/9.0,4.0*pi/9.0] # 3 increments (equivalent of 9 for 3 blades)
   # beta = [0.0,pi] # 2 increments
   # beta = [0.0] # 1 increment (top blade facing straight up)

   B_int = 2.0*pi/B # Intervals betwqeen blades (from the first blade at 0 deg)

   # One-third octave band frequencies (Hz)
   f = [100.0,125.0,160.0,200.0,250.0,315.0,400.0,500.0,
   630.0,800.0,1000.0,1250.0,1600.0,2000.0,2500.0,3150.0,
   4000.0,5000.0,6300.0,8000.0,10000.0,12500.0,16000.0,
   20000.0] #,25000.0,31500.0,40000.0]

   # A-weighting curve (dBA) for sound perception correction
   AdB = [-19.145,-16.190,-13.244,-10.847,-8.675,-6.644,
   -4.774,-3.248,-1.908,-0.795,0.0,0.576,0.993,1.202,
   1.271,1.202,0.964,0.556,-0.114,-1.144,-2.488,-4.250,
   -6.701,-9.341] #,-12.322,-15.694,-19.402]

   for i=1:p
       theta_vel[i] = (2.0*pi/p)*i-(2.0*pi/p)/2.0
   end

   for di=1:bf # for each rotation increment
       for j=1:nf # for each frequency
           for bi=1:B # for each blade
               # Calcuating observer distances and directivty angles for the given blade orientation
               theta = beta[di]+(bi-1)*B_int
               r,theta_e,phi_e = directVAWT(n-1,ox,oy,oz,c,c1,highmid,rad,Hub,rotdir,theta)
               if ((theta >= theta_vel[1]) & (theta <= theta_vel[p]))
                   velwx=splineint(p,theta_vel,wakex,theta)
                   velwy=splineint(p,theta_vel,wakey,theta)
               else
                   velwx = (wakex[1]+wakex[p])/2.0
                   velwy = (wakey[1]+wakey[p])/2.0
               end

               Vx=rot*rad*cos(theta)+Vinf+velwx
               Vy=rot*rad*sin(theta) + velwy
               V=sqrt(Vx^2+Vy^2)
               TBLTV = TBLTVfunc(f[j],V[1],c[1],r[1],theta_e[1],phi_e[1],atip1,c0,
               tipflat,AR)
               TV_t[2*(bi-1)+1] = TBLTV
               TBLTV= TBLTVfunc(f[j],V,c[n-1],r[n-1],theta_e[n-1],phi_e[n-1],atip2,c0,
               tipflat,AR)
               TV_t[2*(bi-1)+2] = TBLTV

               for k=1:n-1
                   # Calculating sound pressure level (dB) for each noise source at each radial position
                   TBLTE = TBLTEfunc(f[j],V,L[k],c[k],r[k],theta_e[k],phi_e[k],alpha[k],
                   nu,c0,trip)

                   if (trip == false)
                       LBLVS = LBLVSfunc(f[j],V,L[k],c[k],r[k],theta_e[k],phi_e[k],alpha[k],
                       nu,c0,trip)
                   else
                       LBLVS = 0.0
                   end
                   TEBVS = TEBVSfunc(f[j],V,L[k],c[k],h[k],r[k],theta_e[k],phi_e[k],
                   alpha[k],nu,c0,psi,trip)
                   # Assigning noise to blade segment
                   TE_t[k.+[n-1]*(bi-1)] .= TBLTE
                   BLVS_t[k.+[n-1]*(bi-1)] .= LBLVS
                   BVS_t[k.+[n-1]*(bi-1)] .= TEBVS
               end
           end
           # Adding sound pressure levels (dB)
           TE[j] = 10.0*log10(sum(10.0.^(TE_t/10.0)))
           TV[j] = 10.0*log10(sum(10.0.^(TV_t/10.0)))
           BLVS[j] = 10.0*log10(sum(10.0.^(BLVS_t/10.0)))
           BVS[j] = 10.0*log10(sum(10.0.^(BVS_t/10.0)))

           # Combining noise sources into overall SPL
           SPLf[j] = 10.0*log10(10.0^(TE[j]/10.0)+10.0^(TV[j]/10.0)+10.0^(BLVS[j]/10.0)+10.0^(BVS[j]/10.0))
       end

       # Correcting with A-weighting
       SPLf[1:nf] = SPLf[1:nf]+AdB[1:nf]

       # Adding SPLs for each rotation increment
       SPLoa_d[di] = 10.0*log10(sum(10.0.^(SPLf/10.0)))

       # Protecting total calcuation from negative SPL values
       if (SPLoa_d[di] < 0.0)
           SPLoa_d[di] = 0.0
       end
   end

   # Performing root mean square calculation of SPLs at rotation increments for final value
   SPLoa = sqrt(sum(SPLoa_d.^2)/bf)
   return SPLoa
end #OASPLVAWT

# Placing a turbine in a specified location and finding the OASPL of the turbine with reference to an observer
"""

turbinepos(x,y,obs,winddir,windvel,rpm,B,Hub,rad,c,c1,alpha,nu,c0,psi,AR,noise_corr)

Calculating the sound pressure level for a HAWT

# Parameters
----------
- `turbx::array`:  x-positions of all the turbines heard by an observer (east to west, meter)
- `turby::array`:  y-positions of all the turbines heard by an observer (north to south, meter)
- `obs::array`:  x-, y-, and z-position of a specified observer (E-W, N-S, height; meter)
- `winddir::float`:  direction the wind blows from (180=N, 270=E, 0=S, 90=W; degree)
- `windvel::array`:  wind velocity at each turbine in the specified wind direction (m/s)
- `rpm::array`:  rotation rate of each turbine (RPM)
- `B::float`:  number of blades on a turbine
- `h::float`:  height of a turbine (meter)
- `rad::array`:  radial positions of the blade geometry (meter)
- `c::array`:  chord length at each radial segment (meter)
- `c1::array`:  distance from the pitch axis to leading edge at each radial segment (meter)
- `alpha::array`:  angle of attack of each radial segment (degree)
- `nu::float`:  kinematic viscosity of the air (m^2/s)
- `c0::float`:  speed of sound of the air (m/s)
- `psi::float`:  solid angle of turbine blades between upper and lower sides of trailing edge (degree)
- `AR::float`:  aspect ratio of turbine blades
- `noise_corr::float`:  correction factor for SPL calculations (1=none, use if calculations differ from expected)

# Returns
----------
- `SPL_HAWT::float`:  sound pressure level calculated at observer location (dB)
"""
function turbinepos(x, y, obs, winddir, windvel, rpm,
                    B, Hub, rad, c, c1, alpha, nu, c0,
                    psi, AR, noise_corr)

    nturb = length(x)

    tSPL = zeros(nturb)
    tSPLA = zeros(nturb)
    SPLf = zeros(nturb, 27) #based on nf from OASPL function
    SPLfA = zeros(nturb, 27)
    windrad = (winddir+180.0)*pi/180.0

    for i = 1:nturb # for each turbine
        # Centering the turbine at (0,0) with repect to the observer location
        ox = obs[1]-x[i]
        oy = obs[2]-y[i]
        oz = obs[3]

        # Adjusting the coordinates to turbine reference frame (wind moving in y-direction)
        rxy = sqrt(ox^2+oy^2)
        ang = atan(oy,ox)+windrad

        ox = rxy*cos(ang)
        oy = rxy*sin(ang)

        # Calculating the overall SPL of each of the turbines at the observer location
        tSPL[i], tSPLA[i], SPLf[i,:], SPLfA[i,:] = OASPL(ox, oy, oz, windvel[i], rpm[i], B, Hub, rad, c, c1, alpha, nu, c0, psi, AR)
    end

    # Combining the SPLs from each turbine and correcting the value based on the wind farm
    SPL_obs = (10.0*log10(sum(10.0.^(tSPL/10.0)))) * noise_corr
    SPLA_obs = (10.0*log10(sum(10.0.^(tSPLA/10.0)))) * noise_corr
    SPLf = (10.0*log10(sum(10.0.^(SPLf/10.0), 1))) * noise_corr
    SPLfA = (10.0*log10(sum(10.0.^(SPLfA/10.0), 1))) * noise_corr

    return SPL_obs, SPLA_obs, SPLf, SPLfA
end #turbinepos

function turbinepos_multi(x::Array{<:Real}, y::Array{<:Real}, obs::Array{<:Real},
                    winddir::Array{<:Real}, windvel::Array{<:Real},
                    rpm::Array{<:Real}, B::Array{<:Real}, Hub::Array{<:Real},
                    rad::Array{<:AbstractArray}, c::Array{<:AbstractArray},
                    c1::Array{<:AbstractArray}, alpha::Array{<:AbstractArray},
                    nu::Real, c0::Real,
                    psi::Array{<:Real}, AR::Array{<:Real}, noise_corr::Real; f=default_f, AdB=default_AdB)

    nf = length(f)
    nturb = length(x)
    tSPL = zeros(nturb)
    tSPLA = zeros(nturb)
    SPLf = zeros(nturb, nf)
    SPLfA = zeros(nturb, nf)
    windrad = (winddir .+ 180.0)*pi/180.0

    for i = 1:nturb # for each turbine
        # Centering the turbine at (0,0) with repect to the observer location
        ox = obs[1]-x[i]
        oy = obs[2]-y[i]
        oz = obs[3]

        # Adjusting the coordinates to turbine reference frame (wind moving in y-direction)
        rxy = sqrt(ox^2+oy^2)
        ang = atan(oy,ox)+windrad[i]

        ox = rxy*cos(ang)
        oy = rxy*sin(ang)

        # Calculating the overall SPL of each of the turbines at the observer location
        tSPL[i], tSPLA[i], SPLf[i,:], SPLfA[i,:] = OASPL(ox,oy,oz,windvel[i],rpm[i],B[i],Hub[i],rad[i],c[i],c1[i],alpha[i],nu,c0,psi[i],AR[i]; f=f, AdB=AdB)
    end

    # Combining the SPLs from each turbine and correcting the value based on the wind farm
    SPL_obs = (10.0*log10.(sum(10.0.^(tSPL/10.0)))) * noise_corr
    SPLA_obs = (10.0*log10(sum(10.0.^(tSPLA/10.0)))) * noise_corr
    SPLf = (10.0*log10.(sum(10.0.^(SPLf/10.0); dims=1))) * noise_corr
    SPLfA = (10.0*log10.(sum(10.0.^(SPLfA/10.0); dims=1))) * noise_corr

    return SPL_obs, SPLA_obs, SPLf, SPLfA
end #turbinepos

"""

turbinepos_VAWT(p,x,y,obs,winddir,B,Hub,high,rad,c,c1,alpha,nu,c0,psi,AR,noise_corr,rot,Vinf,wakex,wakey)

Calculating the sound pressure level for a HAWT

# Parameters
----------
- `p::integer`: number of points along blade flight path to calculate velocities
- `turbx::array`:  x-positions of all the turbines heard by an observer (east to west, meter)
- `turby::array`:  y-positions of all the turbines heard by an observer (north to south, meter)
- `obs::array`:  x-, y-, and z-position of a specified observer (E-W, N-S, height; meter)
- `winddir::float`:  direction the wind blows from (180=N, 270=E, 0=S, 90=W; degree)
- `B::float`:  number of blades on a turbine
- `Hub::float`: hub height of a turbine (meter)
- `high::array`: height positions along the turbine blade (meter)
- `rad::array`:  turbine radius (meter)
- `c::array`:  chord length at each radial segment (meter)
- `c1::array`:  distance from the pitch axis to leading edge at each radial segment (meter)
- `alpha::array`:  angle of attack of each radial segment (degree)
- `nu::float`:  kinematic viscosity of the air (m^2/s)
- `c0::float`:  speed of sound of the air (m/s)
- `psi::float`:  solid angle of turbine blades between upper and lower sides of trailing edge (degree)
- `AR::float`:  aspect ratio of turbine blades
- `noise_corr::float`:  correction factor for SPL calculations (1=none, use if calculations differ from expected)
- `rot::array`: rotation rate of each turbine (rad/s)
- `Vinf::float`: free stream wind speed (m/s)
- `wakex::array`: the wake influenced x-velcoity of the turbine at each point along the blade flight path (m/s)
- `wakey::array`: the wake influenced y-velcoity of the turbine at each point along the blade flight path (m/s)
# Returns
----------
- `SPL_VAWT::float`:  sound pressure level calculated at observer location (dB)
"""
function turbinepos_VAWT(p, x, y, obs, winddir, B, Hub, high,
                            rad, c, c1, alpha, nu, c0, psi, AR,
                            noise_corr, rot, Vinf, wakex, wakey)

    nturb = length(x)
    tSPL = zeros(nturb)
    # SPLf = zeros(nturb)
    # SPLfA = zeros(nturb)

    wakexd = zeros(p)
    wakeyd = zeros(p)
    windrad = (winddir+180.0)*pi/180.0
    for i = 1:nturb # for each turbine
        # Centering the turbine at (0,0) with repect to the observer location
        ox = obs[1]-x[i]
        oy = obs[2]-y[i]
        oz = obs[3]

        # Adjusting the coordinates to turbine reference frame (wind moving in y-direction)
        rxy = sqrt(ox^2+oy^2)
        ang = atan(oy,ox)+windrad

        ox = rxy*cos(ang)
        oy = rxy*sin(ang)

        for j=1:p
            k=p*(i-1)+j
            wakexd[j] = wakex[k]
            wakeyd[j] = wakey[k]
        end
        # Calculating the overall SPL of each of the turbines at the observer location
        tSPL[i] = OASPLVAWT(p,ox,oy,oz,B,Hub,high,rad,c,c1,alpha,nu,c0,psi,rot[i],Vinf,wakexd,wakeyd,AR)
    end

    # Combining the SPLs from each turbine and correcting the value based on the wind farm
    SPL_obs = (10.0*log10.(sum(10.0.^(tSPL/10.0))))*noise_corr

    return SPL_obs
end #turbinepos_VAWT

end # module


text:
#=
utilities for real polynomial interpolation
=#

"""
$(SIGNATURES)

Compute interpolation data for a real weighted sum-of-squares conic constraint
on a domain.
"""
function interpolate(
    dom::Domain{T},
    d::Int;
    constraint_degs::Vector{Int} = Int[],
    calc_V::Bool = false,
    get_quadr::Bool = false,
    sample = nothing,
    sample_factor::Int = 0,
    ) where {T <: Real}
    @assert d >= 1
    @assert all(constraint_degs .>= 0)
    @assert all(2d - degree(dom) .>= constraint_degs)
    n = dimension(dom)
    U = binomial(n + 2d, n)

    if isnothing(sample)
        sample = (!isa(dom, BoxDomain) || n >= 7 || prod_consec(n, d) > 35_000)
    end

    if sample
        if sample_factor <= 0
            if U <= 12_000
                sample_factor = 10
            elseif U <= 15_000
                sample_factor = 5
            elseif U <= 22_000
                sample_factor = 2
            else
                sample_factor = 1
                if U > 35_000 && get_quadr
                    error("dimensions are too large to compute quadrature weights")
                end
            end
        end
        return interp_sample(dom, d, constraint_degs, get_quadr, sample_factor)
    else
        return interp_box(dom, n, d, constraint_degs, calc_V, get_quadr)
    end
end

get_L(n::Int, d::Int) = binomial(n + d, n)
get_U(n::Int, d::Int) = binomial(n + 2d, n)

# sampling-based point selection for general domains
function interp_sample(
    dom::Domain{T},
    d::Int,
    constraint_degs::Vector{Int},
    get_quadr::Bool,
    sample_factor::Int,
    ) where {T <: Real}
    n = dimension(dom)
    U = get_U(n, d)

    cand_pts = sample(dom, U * sample_factor)
    (pts, P0, P0sub, V, w) = make_wsos_arrays(dom, cand_pts, d, get_quadr)

    sqrtg = [sqrt.(gi) for gi in weights(dom, pts)]
    PWts = Matrix{T}[gi .* P0sub for gi in sqrtg]
    Ps = Matrix{T}[P0, PWts...]

    deg_weights = Dict(0 => Ps)
    constraint_Ps = similar(constraint_degs, Vector{Matrix{T}})
    for (i, deg) in enumerate(constraint_degs)
        if !haskey(deg_weights, deg)
            Lsub = get_L(n, div(2d - degree(dom) - deg, 2))
            deg_weights[deg] = Matrix{T}[view(P0, :, 1:get_L(n, div(2d - deg, 2))),
                                         [gi .* view(P0, :, 1:Lsub) for gi in sqrtg]...]
        end
        constraint_Ps[i] = deg_weights[deg]
    end

    return (U = U, pts = pts, Ps = Ps, constraint_Ps = constraint_Ps, V = V, w = w)
end

# slow but high-quality hyperrectangle/box point selections
# unlike sampling functions, P0 is always formed using points in [-1, 1]
function interp_box(
    dom::FreeDomain{T},
    n::Int,
    d::Int,
    constraint_degs::Vector{Int},
    calc_V::Bool,
    get_quadr::Bool,
    ) where {T <: Real}
    @assert all(2d .>=- constraint_degs)
    (U, pts, P0, P0sub, V, w) = interp_box(T, n, d, calc_V, get_quadr)
    return (U = U, pts = pts, Ps = Matrix{T}[P0,],
        constraint_Ps = [Matrix{T}[view(P0, :, 1:get_L(n, div(2d - deg, 2))),]
                         for deg in constraint_degs],
        V = V, w = w)
end

function interp_box(
    dom::BoxDomain{T},
    n::Int,
    d::Int,
    constraint_degs::Vector{Int},
    calc_V::Bool,
    get_quadr::Bool,
    ) where {T <: Real}
    (U, pts, P0, P0sub, V, w) = interp_box(T, n, d, calc_V, get_quadr)

    # TODO refactor/cleanup below
    # scale and shift points, get WSOS matrices
    pscale = [T(0.5) * (dom.u[mod(j - 1, n) + 1] -
        dom.l[mod(j - 1, dimension(dom)) + 1]) for j in 1:n]
    pshift = [T(0.5) * (dom.u[mod(j - 1, n) + 1] +
        dom.l[mod(j - 1, dimension(dom)) + 1]) for j in 1:n]
    @views Wtsfun = (j -> sqrt.(1 .- abs2.(pts[:, j])) * pscale[j])
    PWts = Matrix{T}[Wtsfun(j) .* P0sub for j in 1:n]
    Ps = Matrix{T}[P0, PWts...]

    deg_weights = Dict(0 => Ps)
    constraint_Ps = similar(constraint_degs, Vector{Matrix{T}})
    for (i, deg) in enumerate(constraint_degs)
        if !haskey(deg_weights, deg)
            Lsub = get_L(n, div(2d - degree(dom) - deg, 2))
            deg_weights[deg] = Matrix{T}[view(P0, :, 1:get_L(n, div(2d - deg, 2))),
                                         [Wtsfun(j) .* view(P0, :, 1:Lsub) for j in 1:n]...]
        end
        constraint_Ps[i] = deg_weights[deg]
    end

    trpts = pts .* pscale' .+ pshift'
    return (U = U, pts = trpts, Ps = Ps, constraint_Ps = constraint_Ps, V = V, w = w)
end

function interp_box(
    T::Type{<:Real},
    n::Int,
    d::Int,
    calc_V::Bool,
    get_quadr::Bool,
    )
    if n == 1
        return cheb2_data(T, d, calc_V, get_quadr)
    elseif n == 2
        return padua_data(T, d, calc_V, get_quadr) # or approxfekete_data(n, d)
    elseif n > 2
        return approxfekete_data(T, n, d, get_quadr)
    end
end

# return k Chebyshev points of the second kind
cheb2_pts(T::Type{<:Real}, k::Int) = [-cospi(T(j) / T(k - 1)) for j in 0:(k - 1)]

function calc_univariate_chebyshev(
    pts_i::AbstractVector{T},
    d::Int;
    calc_gradient::Bool = false,
    calc_hessian::Bool = false,
    ) where {T <: Real}
    @assert d > 0
    u = zeros(T, length(pts_i), d + 1)
    @. @views u[:, 1] = 1
    @. @views u[:, 2] = pts_i
    for t in 3:(d + 1)
        @. @views u[:, t] = 2 * pts_i * u[:, t - 1] - u[:, t - 2]
    end

    if !calc_gradient && !calc_hessian
        return u
    end
    @assert calc_gradient

    # calculate gradient
    ug = zero(u)
    @. @views ug[:, 1] = 0
    @. @views ug[:, 2] = 1
    for t in 3:(d + 1)
        @. @views ug[:, t] = 2 * (u[:, t - 1] + pts_i *
            ug[:, t - 1]) - ug[:, t - 2]
    end

    if !calc_hessian
        return (u, ug)
    end
    @assert d > 1

    # calculate hessian
    uh = zero(u)
    @. @views uh[:, 1:2] = 0
    for t in 3:(d + 1)
        @. @views uh[:, t] = 2 * (2 * ug[:, t - 1] + pts_i *
            uh[:, t - 1]) - uh[:, t - 2]
    end

    return (u, ug, uh)
end

function cheb2_data(
    T::Type{<:Real},
    d::Int,
    calc_V::Bool,
    get_quadr::Bool,
    )
    @assert d > 0
    U = get_U(1, d)
    L = get_L(1, d)

    # Chebyshev points for degree 2d
    pts = reshape(cheb2_pts(T, U), :, 1)

    # evaluations
    if calc_V
        V = make_chebyshev_vandermonde(pts, 2d)
        P0 = V[:, 1:L]
    else
        V = zeros(T, 0, 0)
        P0 = make_chebyshev_vandermonde(pts, d)
    end
    P0sub = view(P0, :, 1:get_L(1, d - 1))

    # weights for Clenshaw-Curtis quadrature at pts
    if get_quadr
        wa = T[2 / T(1 - j^2) for j in 0:2:(U - 1)]
        @views append!(wa, wa[div(U, 2):-1:2])
        tempconst = pi / T(length(wa)) * 2 * im
        # inverse FFT
        w = [abs(sum(wa[j] * exp(tempconst * (i - 1) * j) for
            j in eachindex(wa)) / length(wa)) for i in eachindex(wa)]
        w[1] /= 2
        push!(w, w[1])
    else
        w = T[]
    end

    return (U, pts, P0, P0sub, V, w)
end

function padua_data(
    T::Type{<:Real},
    d::Int,
    calc_V::Bool,
    get_quadr::Bool,
    )
    @assert d > 0
    U = get_U(2, d)
    L = get_L(2, d)

    # Padua points for degree 2d
    cheba = cheb2_pts(T, 2d + 1)
    chebb = cheb2_pts(T, 2d + 2)
    pts = zeros(T, U, 2)
    j = 1
    for a in 0:2d, b in 0:(2d + 1)
        if iseven(a + b)
            pts[j, 1] = -cheba[a + 1]
            pts[(U + 1 - j), 2] = -chebb[2d + 2 - b]
            j += 1
        end
    end

    # evaluations
    if calc_V
        V = make_chebyshev_vandermonde(pts, 2d)
        P0 = V[:, 1:L]
    else
        V = zeros(T, 0, 0)
        P0 = make_chebyshev_vandermonde(pts, d)
    end
    P0sub = view(P0, :, 1:get_L(2, d - 1))

    # cubature weights at Padua points
    # even-degree Chebyshev polynomials on the subgrids
    if get_quadr
        te1 = [cospi(T(i * j) / T(2d)) for i in 0:2:2d, j in 0:2:2d]
        to1 = [cospi(T(i * j) / T(2d)) for i in 0:2:2d, j in 1:2:2d]
        te2 = [cospi(T(i * j) / T(2d + 1)) for i in 0:2:2d, j in 0:2:(2d + 1)]
        to2 = [cospi(T(i * j) / T(2d + 1)) for i in 0:2:2d, j in 1:2:(2d + 1)]
        @views te1[2:(d + 1), :] .*= sqrt(T(2))
        @views to1[2:(d + 1), :] .*= sqrt(T(2))
        @views te2[2:(d + 1), :] .*= sqrt(T(2))
        @views to2[2:(d + 1), :] .*= sqrt(T(2))
        # even, even moments matrix
        mom = T(2) * sqrt(T(2)) ./ [T(1 - i^2) for i in 0:2:2d]
        mom[1] = 2
        Mmom = zeros(T, d + 1, d + 1)
        f = inv(T(d * (2d + 1)))
        for j in 1:(d + 1), i in 1:(d + 2 - j)
            Mmom[i, j] = mom[i] * mom[j] * f
        end
        Mmom[1, d + 1] /= 2
        # cubature weights as matrices on the subgrids
        W = zeros(T, d + 1, 2d + 1)
        @views W[:, 1:2:(2d + 1)] .= to2' * Mmom * te1
        @views W[:, 2:2:(2d + 1)] .= te2' * Mmom * to1
        @views W[:, [1, (2d + 1)]] ./= 2
        @views W[1, 2:2:(2d + 1)] ./= 2
        @views W[d + 1, 1:2:(2d + 1)] ./= 2
        w = vec(W)
    else
        w = T[]
    end

    return (U, pts, P0, P0sub, V, w)
end

prod_consec(n::Int, d::Int, j::Int = 0) = prod(big(2d + 1 + j):big(2d + n))

function approxfekete_data(
    T::Type{<:Real},
    n::Int,
    d::Int,
    get_quadr::Bool,
    )
    @assert d > 0
    @assert n > 1

    # points in the initial interpolation grid
    npts = prod_consec(n, d)
    cand_pts = zeros(T, npts, n)

    for j in 1:n
        ig = prod_consec(n, d, j)
        cs = cheb2_pts(T, 2d + j)
        i = 1
        l = 1
        while true
            @views cand_pts[i:(i + ig - 1), j] .= cs[l]
            i += ig
            l += 1
            if l >= 2d + 1 + j
                if i >= npts
                    break
                end
                l = 1
            end
        end
    end

    dom = BoxDomain{T}(-ones(T, n), ones(T, n))
    (pts, P0, P0sub, V, w) = make_wsos_arrays(dom, cand_pts, d, get_quadr)

    return (size(pts, 1), pts, P0, P0sub, V, w)
end

# TODO could merge this function and choose_interp_pts
function make_wsos_arrays(
    dom::Domain{T},
    cand_pts::Matrix{T},
    d::Int,
    get_quadr::Bool,
    ) where {T <: Real}
    n = size(cand_pts, 2)
    (V, keep_pts, w) = choose_interp_pts(cand_pts, d, get_quadr)
    pts = cand_pts[keep_pts, :]
    P0 = V[:, 1:get_L(n, d)] # subset of poly evaluations up to total degree d
    Lsub = get_L(n, div(2d - degree(dom), 2))
    P0sub = view(P0, :, 1:Lsub)
    return (pts, P0, P0sub, V, w)
end

n_deg_exponents(n::Int, deg::Int) = [xp for t in 0:deg for
    xp in Combinatorics.multiexponents(n, t)]

# indices of points to keep and quadrature weights at those points
function choose_interp_pts(
    cand_pts::Matrix{T},
    d::Int,
    get_quadr::Bool,
    ) where {T <: Real}
    n = size(cand_pts, 2)
    U = get_U(n, d)

    V = make_chebyshev_vandermonde(cand_pts, 2d)

    if !get_quadr && size(cand_pts, 1) == U && U > 35_000
        # large matrix and don't need to perform QR procedure, so don't
        # TODO this is hacky; later the interpolate function and functions it
        # calls should take options (and have better defaults) for whether to
        # perform the QR or not
        return (V, 1:U, T[])
    end

    F = qr!(Array(V'), ColumnNorm())
    keep_pts = F.p[1:U]
    V = V[keep_pts, :]

    if get_quadr
        n > 32 && @warn("quadrature weights not numerically stable for large n")
        m = zeros(T, U)
        for (col, xp) in enumerate(n_deg_exponents(n, 2d))
            if all(iseven, xp)
                @inbounds m[col] = prod(2 / (1 - xp[j]^2) for j in 1:n)
            end
        end
        Qtm = F.Q' * m
        w = UpperTriangular(F.R[:, 1:U]) \ Qtm
        return (V, keep_pts, w)
    end

    return (V, keep_pts, T[])
end

# make Vandermonde with rows corresponding to points and columns to Chebyshev polys
function make_chebyshev_vandermonde(pts::Matrix{T}, deg::Int) where {T <: Real}
    n = size(pts, 2)
    expos = n_deg_exponents(n, deg)
    univ_chebs = [calc_univariate_chebyshev(view(pts, :, i), deg) for i in 1:n]
    return make_product_vandermonde(univ_chebs, expos)
end

function make_product_vandermonde(
    u::Vector{Matrix{T}},
    expos::Vector,
    ) where {T <: Real}
    npts = size(u[1], 1)
    n = length(u)
    V = zeros(T, npts, length(expos))

    for (col, xp) in enumerate(expos)
        @inbounds @. @views V[:, col] = u[1][:, xp[1] + 1]
        @inbounds for j in 2:n
            @. @views V[:, col] *= u[j][:, xp[j] + 1]
        end
    end

    return V
end


text:
export Model, set_input!, get_output, reward!

struct Model
    cfg::Dict
    cont::Controller
    cells::Array{Cell}
    inputs::Array{Ref{Cell}}
    outputs::Array{Ref{Cell}}
    rewards::Array{Ref{Cell}}
    state::Array{Int64}
end

function Model(cfg::Dict, cont::Controller)
    cells = Array{Cell}(undef, 0)
    inputs = Array{Ref{Cell}}(undef, 0)
    outputs = Array{Ref{Cell}}(undef, 0)
    rewards = Array{Ref{Cell}}(undef, 0)
    state = zeros(1)
    Model(cfg, cont, cells, inputs, outputs, rewards, state)
end

function set_input!(m::Model, inputs::Array{Float64})
    for i in eachindex(inputs)
        m.inputs[i][].input = inputs[i]
    end
    nothing
end

function get_output(m::Model)
    map(i->m.outputs[i][].output, eachindex(m.outputs))
end

function reward!(m::Model, reward::Array{Float64})
    for i in eachindex(reward)
        m.rewards[i][].input = reward[i]
    end
    nothing
end

function get_conns(m::Model, cell::Cell)
    conns = Array{Conn}(undef, 0)
    for c in m.cells
        for conn in c.conns
            if conn.dest[] == cell
                push!(conns, conn)
            end
        end
    end
    conns
end


text:
#=
# Optional log writing
=#
export printerr, init_log, set_log_level, log_entry, log_dump_config, log_dump_prob, close_log

# verbosity levels
# 0 = minimal, not very useful
# 1 = basic messages about steps to show progress
# 2 = default, details about each step
# 3 = debug, lots of info for debugging that is hard to read
log_level = 2;

function set_log_level(lev)
    global log_level = lev;
end

function printerr(msg)
    log_entry("Error: "*msg, 0);
    println("Error: "*msg);
end

function init_log(name, dir, level=2)
    global log_file = dir*"/"*name*".txt";
    global use_log = true;
    global log_level = level;
    file = open(log_file, "w");
    println(file, "######################################");
    println(file, "# Femshop Log for: "*project_name);
    println(file, "######################################");
    println(file, "(verbosity = "*string(log_level)*")");
    close(file)
end

function log_entry(text, level=2)
    global log_line_index;
    global use_log;
    if use_log && (level <= log_level)
        file = open(log_file, "a");
        println(file, string(log_line_index)*".\t"*text);
        log_line_index += 1;
        close(file);
    end
end

function log_dump_config(c=config)
    global log_line_index;
    global use_log;
    if use_log
        file = open(log_file, "a");
        println(file, string(log_line_index)*".\tDumping configuration:");
        log_line_index += 1;
        for f in fieldnames(Femshop_config)
            println(file, string(log_line_index)*".\t\t"*string(f)*" = "*string(getfield(c, f)));
            log_line_index += 1;
        end
        close(file);
    end
end

function log_dump_prob(p = prob)
    global log_line_index;
    global use_log;
    if use_log
        file = open(log_file, "a");
        println(file, string(log_line_index)*".\tDumping problem details:");
        log_line_index += 1;
        for f in fieldnames(Femshop_prob)
            println(file, string(log_line_index)*".\t\t"*string(f)*" = "*string(getfield(p, f)));
            log_line_index += 1;
        end
        close(file);
    end
end

function close_log()
    global use_log;
    if use_log
        # log_entry("Finalizing. Dumping state.", 3);
        # if log_level >= 3
        #     log_dump_config();
        #     log_dump_prob();
        # end

        log_entry("Completed. Closing Log.", 0);
        use_log = false;
    end
end


text:
#=
 █████╗  ██████╗██████╗ ██╗
██╔══██╗██╔════╝██╔══██╗██║
███████║██║     ██████╔╝██║
██╔══██║██║     ██╔══██╗██║
██║  ██║╚██████╗██║  ██║███████╗
╚═╝  ╚═╝ ╚═════╝╚═╝  ╚═╝╚══════╝
File:       rendezvous_review.jl
Author:     Gabriel Barsi Haberfeld, 2020. gbh2@illinois.edu
Function:   This program simulates all results in the paper "Geometry-Informed
            Minimum Risk Rendezvous Algorithm for Heterogeneous Agents in Urban
            Environments"
Instructions:   Run this file in juno with Julia 1.2.0 or later.
Requirements:   JuMP, Ipopt, Plots, LinearAlgebra, BenchmarkTools.
=#

using Plots, LinearAlgebra, Random, Statistics, ColorSchemes, LazySets
using JuMP, Ipopt, Measures, Printf, LaTeXStrings, Distributions
using BenchmarkTools, StatsPlots

mass = [3, 1]
m = mass
alpha = 20
vmax = 20
vo = 50
ho = 20
pv1 = [0 800; 600 800; 600 900; 1300 900]
pv2 = [0 800; 1000 800; 1000 1200]
Mp = 1000
MaxDriverSpeed = 8
tb = Mp / MaxDriverSpeed

include("Map.jl")
include("GPTools.jl")
include("OptimizationTools.jl")
include("PlotTools.jl")
include("DriverFunctions.jl")

function rowNormSquared(mat)
    A = zeros(Float64, size(mat, 2))
    for j = 1:size(mat, 2)
        for i = 1:size(mat, 1)
            A[j] += mat[i, j]^2
        end
        #A[j] = sqrt(A[j])
    end
    A
end

function rowNorm(mat)
    A = zeros(Float64, size(mat, 2))
    for j = 1:size(mat, 2)
        for i = 1:size(mat, 1)
            A[j] += mat[i, j]^2
        end
        A[j] = sqrt(A[j])
    end
    A
end

function rankMultiPath(
    TimeSamples,
    UASPos,
    Ns,
    p = [1, 2],
    gp = nothing,
    ti = 0.0,
    DriverPos = 0.0,
    rmethod = "BestFirst",
    Er = Inf,
    sol = nothing,
    riska = false,
)
    np = size(TimeSamples, 2)
    if np != length(p)
        error("Invalid #map v #smpl")
    end
    lt = size(TimeSamples, 1)
    x = hcat(UASPos[1] * ones(lt), UASPos[2] * ones(lt))
    #x contains a copy-vector of x-y coordinates of the UAS
    E = zeros(lt, np)
    PosSamples = DriverPosition(ti, TimeSamples, gp) .+ DriverPos
    if riska
        USamples = DriverUncertainty(ti, TimeSamples, gp)
        pU = PosSamples .+ USamples #Position + σ
        nU = PosSamples .- USamples #Position - σ
    end
    Threads.@threads for i = 1:np #for each path
        EuclideanPositions = pathSample2Array(path.(PosSamples[:, i], i))
        EuclideanDistances = EuclideanPositions .- x #mean distance
        if riska
            pUEuclidean = pathSample2Array(path.(pU[:, i], i)) .- x
            nUEuclidean = pathSample2Array(path.(nU[:, i], i)) .- x
            mRadii =
                max.(
                    rowNorm(pUEuclidean'), #r_j^+
                    rowNorm(nUEuclidean'), #r_j^-
                    rowNorm(EuclideanDistances'), #r_j
                )
        else
            mRadii = rowNorm(EuclideanDistances')
        end

        v = mRadii ./ TimeSamples[:, i]

        if sol == nothing
            E[:, i] =
                m[1] .* v .^ 2 .* TimeSamples[:, i] .+
                m[1] * alpha * TimeSamples[:, i]
        else #we have a solution available
            vs = sol[1]
            ts = sol[2]
            LPos = sol[3]
            tSol = ts[3] #just use the time provided by the solver

            EuclideanDistances = LPos' .- EuclideanPositions #mean distance
            if riska
                pUEuclidean = LPos' .- pathSample2Array(path.(pU[:, i], i))
                nUEuclidean = LPos' .- pathSample2Array(path.(nU[:, i], i))
                mRadii =
                    max.(
                        rowNorm(pUEuclidean'), #r_j^+
                        rowNorm(nUEuclidean'), #r_j^-
                        rowNorm(EuclideanDistances'), #r_j
                    )
            else
                mRadii = rowNorm(EuclideanDistances')
            end

            LPos = hcat(LPos[1] * ones(lt), LPos[2] * ones(lt))
            vNew = mRadii ./ tSol
            E[:, i] =
                m[1] .* v .^ 2 .* TimeSamples[:, i] .+
                m[1] * alpha * TimeSamples[:, i] .+
                1 .*
                (m[2] .* rowNormSquared(vNew') .* tSol .+ m[2] * alpha * tSol)
        end
    end
    # E now contains energies for each sample group. Next choose a rmethod
    replace!(E, NaN => Inf)
    if Ns == 1
        return [argmin(E[:, 1]) argmin(E[:, 2])]
    end
    if rmethod == "BestFirst"
        Emins = zeros(np)
        elites = zeros(Int16, Ns, np)
        Threads.@threads for i = 1:np
            elites[:, i] = partialsortperm(E[:, i], 1:min(Ns, lt))
            idx = elites[1, i]
            Emins[i] = E[idx, i]
        end
        #to select the path using this policy
        #get the best from each, and output the best
        ptgt = argmin(Emins)
        return elites, ptgt
    elseif rmethod == "WorstFirst"
        Emins = zeros(np)
        elites = zeros(Int16, Ns, np)
        Threads.@threads for i = 1:np
            elites[:, i] = partialsortperm(E[:, i], 1:min(Ns, lt))
            idx = elites[1, i]
            Emins[i] = E[idx, i]
        end
        #to select the path using this policy
        #get the best from each, and output the best
        ptgt = argmax(Emins)
        return elites, ptgt
    else
        error("Invalid rmethod Argument")
    end
    Es = E[:, 1] .+ E[:, 2]
    SumElites = partialsortperm(Es, 1:min(n, lt))
end

function pathSample2Array(pathSample)
    return hcat([y[1] for y in pathSample], [y[2] for y in pathSample])
end

function SampleTime(μ, Σ, N = 1, t0 = 0.0, tf = 75.0)
    if size(μ, 1) != size(Σ, 1) || size(μ, 2) != 1 || size(Σ, 2) != 1
        error("Wrong distribution dimensions!")
    end
    ndist = size(μ, 1)
    tf = tb
    TimeSamples = μ' .+ randn(N, ndist) .* Σ'
    TimeSamples[TimeSamples.>=tf] .= tf
    TimeSamples[TimeSamples.<=t0] .= t0
    return TimeSamples
end

function TestSampling(μ, Σ, N)
    TimeSamples = SampleTime(μ, Σ, N)
    plotTimeSamples(TimeSamples, 1)
end

function CEM(μ, Σ, np, elites, OptTimeSample, TimeSamples, Ns)
    β = 0.9
    γ = 1.0
    elite = elites[1, :]
    Threads.@threads for j = 1:length(np)
        μn = sum(TimeSamples[elites[:, j], j]) / Ns
        μ[j] = γ * μn + (1 - γ) * μ[j]
        Σ[j] =
            β * (sum((TimeSamples[elites[:, j], j] .- μ[j]) .^ 2) ./ Ns + 1.0) +
            (1 - β) * Σ[j]
        Σ[j] = min(Σ[j], 100)
        OptTimeSample[j] = TimeSamples[elite[j], j]
    end
    return μ, Σ, OptTimeSample
end

function iterateCEM(μ, Σ, np)
    OptTimeSample = zeros(length(p))
    for i = 1:1000
        TimeSamples = SampleTime(μ, Σ, N)
        elites = rankMultiPath(TimeSamples, UASPos, Ns, np)
        CEM(μ, Σ, np, elites, OptTimeSample, TimeSamples)
    end
    return μ, Σ, OptTimeSample
end

function animateGpCem(tt = 20, ti = 5, dt = 1)
    #first build an initial dataset and gp
    t = range(0.0, stop = ti, step = dt)
    li = length(t)
    D = zeros(li, 2)
    D[:, 1] = VelocityPrior.(t)
    D[:, 2] = Deviation.(VelocityPrior.(t)) + NoiseStd .* randn(li)
    gp = LearnDeviationFunction(D, true, "full")
    #now we move the driver, collect new data, update gp and sample on the gp
    t = range(0.0, stop = tt, step = dt)
    l = length(t)
    D = [D; zeros(l, 2)]
    DriverPos = 0.0
    DriverPosVec = zeros(l)
    DriverPosEstimateVec = zeros(l)
    DriverPosErrorVec = zeros(l)
    PosSamples = zeros(N, length(p))
    OptTimeSample = zeros(length(p))
    anim = @animate for i = 1:l
        #sample driver velocity
        DriverVelSample = DriverVelocity(t[i]) + NoiseStd * randn(1)[1]
        DriverDeviationSample = DriverVelSample - VelocityPrior(t[i])
        #add to dataset
        D[li+i, 1] = VelocityPrior(t[i])
        D[li+i, 2] = DriverDeviationSample
        #re-train gp
        gp = LearnDeviationFunction(D[1:(li+i), :], true, "DTC")
        #sample rendezvous candidates
        TimeSamples = SampleTime(μ, Σ, N)
        elites =
            rankMultiPath(TimeSamples, UASPos, Ns, np, gp, t[i], DriverPos)
        CEM(μ, Σ, np, elites, OptTimeSample, TimeSamples)
        PosSamples = DriverPosition(t[i], TimeSamples, gp) .+ DriverPos
        l = @layout [a; b]
        pp = plot()
        pp = plotpath!(np)
        pp = plotTimeSamples!(TimeSamples, np, t[i], DriverPos, gp)
        pp = drawMultiConvexHull(PosSamples, p)
        po = plot(gp, xlims = (6, 10), ylim = (-2, 2))
        pp = plot(pp, po, layout = l)
        DriverPosVec[i] = DriverPos
        DriverPosEstimateVec[i] = DriverPosition(0.0, t[i])
        DriverPosErrorVec[i] = DriverPosVec[i] - DriverPosEstimateVec[i]
        DriverPos = DriverPos + DriverVelocity(t[i]) * dt
        @show μ, Σ, i
    end
    gif(anim, "GpCem.gif", fps = 15)
end

#=
Closed loop mission using BestFirst. Order of things:
-Collect sensor data
-Update driver model with gp
-Update distribution with CEM
-Compute trajectory using MPC
-Send control iputs to UAS
-Update driver position
-Update UAS position
=#

function mission(Er = 18000, rmethod = "WorstFirst", tt = 80, ti = 5, dt = 1)
    Random.seed!(125)
    UASPos = [800, 450]
    LPos = [1000, 600]
    #LPos = [800, 450]
    ts = 0
    PNRStat = false
    N = 100
    p = [1, 2]
    ptgt = 2
    dt = 1
    PrevPNR = [0.0, 0.0]
    μ = 60 * ones(2)
    Σ = 10 * ones(2)
    Ns = 5
    UASPosVec = zeros(N, 2)
    #TimeSamples = ts:120
    t = range(0.0, stop = ti, step = dt)
    li = length(t)
    D = zeros(li, 2)
    D[:, 1] = VelocityPrior.(t)
    D[:, 2] = Deviation.(VelocityPrior.(t)) + NoiseStd .* randn(li)
    gp = LearnDeviationFunction(D, true, "DTC")
    #now we move the driver, collect new data, update gp and sample on the gp
    tv = range(0.0, stop = tt, step = dt)
    l = length(tv)
    D = [D; zeros(l, 2)]
    DriverPos = 0.0
    DriverPosVec = zeros(l)
    DriverPosEstimateVec = zeros(l)
    DriverPosErrorVec = zeros(l)
    PosSamples = zeros(N, length(p))
    OptTimeSample = zeros(length(p))
    sol = nothing
    solr = 0.0
    riska = true
    anim = @animate for i = 1:l
        #sample driver velocity
        DriverVelSample = DriverVelocity(tv[i]) + NoiseStd * randn(1)[1]
        DriverDeviationSample = DriverVelSample - VelocityPrior(tv[i])
        #add to dataset
        D[li+i, 1] = VelocityPrior(tv[i])
        D[li+i, 2] = DriverDeviationSample
        #re-train gp
        gp = LearnDeviationFunction(D[1:(li+i), :], true, "DTC")
        #sample rendezvous candidates
        UASPosVec[i, :] = UASPos
        TimeSamples = SampleTime(μ, Σ, N)
        elites, ptgt = rankMultiPath(
            TimeSamples,
            UASPos,
            Ns,
            p,
            gp,
            tv[i],
            DriverPos,
            rmethod,
            Er,
            sol,
            riska,
        )
        CEM(μ, Σ, p, elites, OptTimeSample, TimeSamples, Ns)
        PosSamples = DriverPosition(tv[i], TimeSamples, gp) .+ DriverPos
        #MPC goes here
        v, t = RendezvousPlanner(
            UASPos,
            LPos,
            OptTimeSample,
            Er,
            tv[i],
            PrevPNR,
            false,
            ptgt,
            gp,
            DriverPos,
        )
        sol = (v, t, LPos, UASPos, DriverPos, tv[i], OptTimeSample, ptgt)
        solr = sol
        pp = plot()
        pp = plot!(UASPosVec[1:i, 1], UASPosVec[1:i, 2], legend = false)
        pp = drawMultiConvexHull(PosSamples, p)
        pp = plotPlan!(
            UASPos,
            LPos,
            v,
            t,
            TimeSamples,
            length(p),
            tv[i],
            DriverPos,
            gp,
        )
        pp = plotTimeSamples!(TimeSamples, p, tv[i], DriverPos, gp)
        pp = scatter!(
            path(DriverPos, ptgt),
            markersize = 6.0,
            label = "",
            markercolor = :green,
        )
        PrevPNR = v[:, 1] .* t[1] .+ UASPos
        Ed =
            m[1] * v[1, 1]^2 * t[1] +
            m[1] * v[1, 2]^2 * t[2] +
            m[2] * v[1, 3]^2 * t[3] +
            m[1] * v[2, 1]^2 * t[1] +
            m[1] * v[2, 2]^2 * t[2] +
            m[2] * v[2, 3]^2 * t[3] +
            m[1] * alpha * t[1] +
            m[1] * alpha * t[2] +
            m[2] * alpha * t[3]
        x, y, Er = uav_dynamics(
            UASPos[1],
            UASPos[2],
            v[1, 1],
            v[2, 1],
            dt,
            Er,
            vmax,
            PNRStat ? m[2] : m[1],
        )
        UASPos = [x, y]
        if ptgt == 2
            pv = pv2
        else
            pv = pv1
        end
        # if (t[1] <= 1.0 ||
        #    EuclideanDistance(
        #     path(DriverPosition(tv[i], OptTimeSample[ptgt], gp) + DriverPos),
        #     pv[2, :],
        # ) <= 10) && tv[i] > 10
        #     break
        # end
        if (t[1] <= 1.0 && tv[i] > 10)
            break
        end
        #@show t, i, Σ, μ, Er, Ed
        DriverPosVec[i] = DriverPos
        DriverPosEstimateVec[i] = DriverPosition(0.0, tv[i])
        DriverPosErrorVec[i] = DriverPosVec[i] - DriverPosEstimateVec[i]
        DriverPos = DriverPos + DriverVelocity(tv[i]) * dt
        vsq = rowNorm(sol[1])
        @show i, t[1], ptgt, vsq
    end
    #risk1 = GPCVaR(gp, solr, 1)
    #risk2 = GPCVaR(gp, solr, 2)
    #if (risk1 > 200 || risk2 > 200)
    #    println("Risk too high, aborting!")
    #end
    gif(anim, "RDV_Anim_MP.gif", fps = 15)
    return solr
end

function GPCVaR(
    gp,
    sol, #v, t, LPos, UASPos, DriverPos, ti, OptTimeSample, ptgt
    ptgt = 1,
    n = 10000,
    c_α = 0.01,
)
    vs = sol[1]
    ts = sol[2]
    LPos = sol[3]
    UASPos = sol[4]
    DriverPos = sol[5]
    ti = sol[6]

    OptTimeSample = sol[7][ptgt]

    μ = DriverPosition(ti, OptTimeSample, gp) + DriverPos #mean
    Σ = DriverUncertainty(ti, OptTimeSample, gp) #variance
    RDVPos = path(μ, ptgt, true)
    gpStd = sqrt(Σ + 1.0)
    gp_distr = Normal(μ, gpStd)
    x = rand(gp_distr, n)
    UASPosv = hcat(UASPos[1] * ones(n), UASPos[2] * ones(n))
    LPosv = hcat(LPos[1] * ones(n), LPos[2] * ones(n))
    EuclideanPositions = pathSample2Array(path.(x, ptgt))
    EuclideanDistancesUAS = EuclideanPositions .- UASPosv
    EuclideanDistancesL = EuclideanPositions .- LPosv

    EDS = rowNorm(EuclideanDistancesUAS')
    EDL = rowNorm(EuclideanDistancesL')

    t1 = (OptTimeSample - ti)
    #nominal
    sEDS = norm(UASPos - RDVPos)
    sEDL = norm(LPos - RDVPos)
    sVS = sEDS ./ t1
    Em, tm, vm = minReturnEnergy(RDVPos, LPos)

    #compute new velocities
    VS = EDS ./ t1
    VL = EDL ./ tm

    #compute nominal energy
    Γ1 = sVS^2 * m[1] / 2 * t1 + m[1] * alpha * t1 + Em

    #compute new energies
    E =
        VS .^ 2 .* m[1] ./ 2 .* t1 + .+VL .^ 2 .* m[2] ./ 2 .* tm .+
        m[1] * alpha * t1 .+ m[2] * alpha * tm
    xd = Γ1 .- E
    # To find our elusive distribution, we will sample from gp and fit.
    ddistr = fit(Normal, xd)

    #now its easy to compute CVaR of a normal distribution
    x_α = VaR(ddistr, c_α)
    CVaRα = CVaR(ddistr, x_α, c_α)
    gainmain = -CVaRα #extra distance we have to deal with on the main path

    #now redo everything for the second path, but subs t3 with new one
    maintgt = ptgt
    ptgt = ptgt == 1 ? 2 : 1

    OptTimeSample = sol[7][ptgt]

    μ = DriverPosition(ti, OptTimeSample, gp) + DriverPos #mean
    Σ = DriverUncertainty(ti, OptTimeSample, gp) #variance
    RDVPos = path(μ, ptgt, true)
    gpStd = sqrt(Σ + 1.0)
    gp_distr = Normal(μ, gpStd)
    x = rand(gp_distr, n)
    UASPosv = hcat(UASPos[1] * ones(n), UASPos[2] * ones(n))
    LPosv = hcat(LPos[1] * ones(n), LPos[2] * ones(n))
    EuclideanPositions = pathSample2Array(path.(x, ptgt))
    EuclideanDistancesUAS = EuclideanPositions .- UASPosv
    EuclideanDistancesL = EuclideanPositions .- LPosv

    EDS = rowNorm(EuclideanDistancesUAS')
    EDL = rowNorm(EuclideanDistancesL')

    t1 = (OptTimeSample - ti)
    #nominal
    sEDS = norm(UASPos - RDVPos)
    sEDL = norm(LPos - RDVPos)
    sVS = sEDS ./ t1
    Em, tm, vm = minReturnEnergy(RDVPos, LPos)

    #compute new velocities
    VS = EDS ./ t1
    VL = EDL ./ tm

    #compute nominal energy
    Γ1 = sVS^2 * m[1] / 2 * t1 + m[1] * alpha * t1 + Em

    #compute new energies
    E =
        VS .^ 2 .* m[1] ./ 2 .* t1 + .+VL .^ 2 .* m[2] ./ 2 .* tm .+
        m[1] * alpha * t1 .+ m[2] * alpha * tm
    xd = Γ1 .- E
    # To find our elusive distribution, we will sample from gp and fit.
    ddistr = fit(Normal, xd)

    #now its easy to compute CVaR of a normal distribution
    x_α = VaR(ddistr, c_α)
    CVaRα = CVaR(ddistr, x_α, c_α)
    gainalt = -CVaRα

    return gainmain, gainalt, maintgt

end

function VaR(d, α::Real)
    #truncate:
    q1 = 0.001
    q2 = 1 - q1
    @assert(α > q1)
    @assert(α < q2)
    q1 = quantile(d, q1)
    q2 = quantile(d, q2)
    x = range(q1, stop = q2, length = 10^3)
    Fx = cdf(d, x)
    i = findfirst(x -> x >= α, Fx)
    return x[i]
end

function CVaR(d, x_α::Real, α::Real)
    x_i = quantile(d, 1e-6)
    CVaRα, err = quadgk(x -> x * pdf(d, x), x_i, x_α, 1e-9)
    CVaRα = 1 / α * CVaRα
end

max_stars_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_issues_repo_name: gbarsih/Multi-Path-Safe-Rendezvous

max_forks_repo_name: gbarsih/Multi-Path-Safe-Rendezvous


text:
module TestData
  using Test
  using EcologicalNetworks

  all_wol = web_of_life()
  for wol in all_wol
    N = web_of_life(wol.ID)
    @test richness(N) == wol.Species
    @test links(N) == wol.Interactions
  end

  for n in nz_stream_foodweb()
    @test typeof(n) <: BinaryNetwork
    @test typeof(n) <: AbstractUnipartiteNetwork
  end

end