text:
abstract type AbstractWakeDeficitModel end

"""
    JensenTopHat(alpha)

Container for parameters related to the Jensen Top Hat deficit model

# Arguments
- `alpha::Float`: parameter controlling the wake spreading rate and deficit decay. Default value is 0.1
"""
struct JensenTopHat{TF} <: AbstractWakeDeficitModel
    alpha::TF
end
JensenTopHat() = JensenTopHat(0.2)

"""
    JensenCosine(alpha)

Container for parameters related to the Jensen Cosine deficit model

# Arguments
- `alpha::Float`: parameter controlling the wake deficit decay rate. Default value is 0.1
- `beta::Float`: parameter controlling the width of the cosine function. Default value is 20.0 deg., given in radians.
"""
struct JensenCosine{TF,ATF} <: AbstractWakeDeficitModel
    alpha::TF
    beta::TF
    wec_factor::ATF
end
JensenCosine() = JensenCosine(0.1, 20.0*pi/180.0, [1.0])
JensenCosine(x) = JensenCosine(x, 20.0*pi/180.0, [1.0])
JensenCosine(x, y) = JensenCosine(x, y, [1.0])

"""
    Multizone(me, ke, MU, aU, bU)

Container for parameters related to the Multizone deficit model

# Arguments
- `me::Float`: parameter controlling general wake expansion. Default value is 0.065
- `ke::Array{Float}(3)`: parameters controlling the wake expansion of each zone respectively. Default values are [-0.5 0.22 1.0].
- `MU::Array{Float}(3)`: parameters controlling the wake deficit decay of each zone respectively. Default values are [0.5 1.0 5.5].
- `aU::Float`: parameter impacting the wake deficit decay for a constant wake deflection. Default value is 5.0.
- `bU::Float`: parameter changing the wake deficit decay under yawed conditions. Default value is 1.66.
"""
struct Multizone{ATF, TF} <: AbstractWakeDeficitModel
    me::ATF
    ke::TF
    MU::ATF
    aU::TF
    bU::TF
end
Multizone() = Multizone([-0.5 0.22 1.0], 0.065, [0.5 1.0 5.5], 5.0, 1.66)

"""
    GaussOriginal(k_star)

Container for parameters related to the origina Gaussian deficit model presented by Bastankhah and Porte-Agel 2014

# Arguments
- `k_star::Float`: parameter controlling the wake spreading rate and deficit decay. Default value is 0.075
"""
struct GaussOriginal{TF} <: AbstractWakeDeficitModel
    k_star::TF
end
GaussianOriginal() = GaussOriginal(0.075)

"""
    GaussYaw(turbulence_intensity, horizontal_spread_rate, vertical_spread_rate, alpha_star, beta_star)

Container for parameters related to the Gaussian deficit model with yaw presented by Bastankhah and Porte-Agel 2016

# Arguments
- `horizontal_spread_rate::Float`: parameter controlling the horizontal spread of the deficit model. Default value is 0.022.
- `vertical_spread_rate::Float`: parameter controlling the vertical spread of the deficit model. Default value is 0.022.
- `alpha_star::Float`: parameter controlling the impact of turbulence intensity on the length of the near wake. Default value is 2.32.
- `beta_star::Float`: parameter controlling the impact of the thrust coefficient on the length of the near wake. Default value is 0.154.
"""
struct GaussYaw{TF, ATF} <: AbstractWakeDeficitModel
    horizontal_spread_rate::TF
    vertical_spread_rate::TF
    alpha_star::TF
    beta_star::TF
    wec_factor::ATF
end
GaussYaw() = GaussYaw(0.022, 0.022, 2.32, 0.154, [1.0])
GaussYaw(a, b, c, d) = GaussYaw(a, b, c, d, [1.0])

"""
    GaussYawVariableSpread(turbulence_intensity, horizontal_spread_rate, vertical_spread_rate, alpha_star, beta_star)

Container for parameters related to the Gaussian deficit model with yaw presented by Bastankhah and Porte-Agel 2016

# Arguments
- `alpha_star::Float`: parameter controlling the impact of turbulence intensity on the length of the near wake. Default value is 2.32.
- `beta_star::Float`: parameter controlling the impact of the thrust coefficient on the length of the near wake. Default value is 0.154.
"""
struct GaussYawVariableSpread{TF, ATF} <: AbstractWakeDeficitModel
    alpha_star::TF
    beta_star::TF
    k1::TF
    k2::TF
    wec_factor::ATF
end
GaussYawVariableSpread() = GaussYawVariableSpread(2.32, 0.154, 0.3837, 0.003678, [1.0])
GaussYawVariableSpread(x, y, z) = GaussYawVariableSpread(x, y, 0.3837, 0.003678, z)
GaussYawVariableSpread(x, y) = GaussYawVariableSpread(x, y, 0.3837, 0.003678, [1.0])

"""
    GaussSimple(k, wec_factor)

Container for parameters related to the Gaussian deficit model with yaw presented by Bastankhah and Porte-Agel 2016

# Arguments
- `k::Float`: parameter controlling the spread of the wake
- `wec_factor::Array{Float}`: parameter artificial wake spreading for wake expansion continuation (WEC) optimization
"""
struct GaussSimple{TF, ATF} <: AbstractWakeDeficitModel
    k::TF
    wec_factor::ATF
end
GaussSimple(k) = GaussSimple(k, [1.0])

"""
    wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::JensenTopHat)

Computes the wake deficit according to the original Jensen top hat wake model, from the paper:
"A Note on Wind Generator Interaction" by N.O. Jensen (1983)

# Arguments
- `locx::Float`: x coordinate where wind speed is calculated
- `locy::Float`: y coordinate where wind speed is calculated
- `locz::Float`: z coordinate where wind speed is calculated
- `turbine_x::Array(Float)`: vector containing x coordinates for all turbines in farm
- `turbine_y::Array(Float)`: vector containing y coordinates for all turbines in farm
- `turbine_z::Array(Float)`: vector containing z coordinates for all turbines in farm
- `deflection_y::Float`: deflection in the y direction of downstream wake
- `deflection_z::Float`: deflection in the z direction of downstream wake
- `upstream_turbine_id::Int`: index of the upstream wind turbine creating the wake
- `downstream_turbine_id::Int`: index of the downstream turbine feeling the wake (if not referencing a turbine set to zero)
- `hub_height::Array(Float)`: vector containing hub heights for all turbines in farm
- `rotor_diameter::Array(Float)`: vector containing rotor diameters for all turbines in farm
- `turbine_ai::Array(Float)`: vector containing initial velocity deficits for all turbines in farm
- `turbine_local_ti::Array(Float)`: vector containing local turbulence intensities for all turbines in farm
- `turbine_ct::Array(Float)`: vector containing thrust coefficients for all turbines in farm
- `turbine_yaw::Array(Float)`: vector containing the yaw angle? for all turbines in farm
- `model::JensenTopHat`: indicates the wake model in use

"""
function wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::JensenTopHat)

    r0 = rotor_diameter[upstream_turbine_id]/2.0 #turbine rotor radius

    if downstream_turbine_id == 0

        # find delta x, y, and z. dx is the downstream distance from the turbine to
        # the point of interest. dy and dz are the distances from the point of interest
        # and the wake center (in y and z)
        dx = locx-turbine_x[upstream_turbine_id]
        dy = locy-(turbine_y[upstream_turbine_id]+deflection_y)
        dz = locz-(turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

        del = sqrt(dy^2+dz^2) #distance from wake center to the point of interest
        r = model.alpha*dx + r0 #figure (1) from the paper

        if (dx < 0.0) || (del > r)
            totloss = 0.0
        else
            loss = 2.0*turbine_ai[upstream_turbine_id]*(r0/(r))^2 #equation (2) from the paper
            totloss = loss
        end
    else
        # find delta x, y, and z. dx is the downstream distance from the upstream turbine hub to
        # the downstream turbine hub of interest. dy and dz are the distances from the point of interest
        # and the wake center (in y and z)
        dxt = turbine_x[downstream_turbine_id]-turbine_x[upstream_turbine_id]
        dyt = turbine_y[downstream_turbine_id]-(turbine_y[upstream_turbine_id]+deflection_y)
        dzt = turbine_z[downstream_turbine_id]-(turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

        r = model.alpha*dxt + r0 #figure (1) from the paper

        posz = turbine_z[downstream_turbine_id] + hub_height[downstream_turbine_id]

        # finds the wakediameter
        wake_center_y = (turbine_y[upstream_turbine_id]+deflection_y)
        wake_center_z = (turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)
        wake_diameter = 2*r

        loss = 2.0*turbine_ai[upstream_turbine_id]*(r0/(r))^2 #equation (2) from the paper

        overlap = overlap_area_func(turbine_y[downstream_turbine_id], posz,
        rotor_diameter[downstream_turbine_id], wake_center_y, wake_center_z, wake_diameter; tol=1E-6)

        ovlpercentage = overlap/(pi*r0^2)
        totloss = loss*ovlpercentage

    end

    return totloss
end

"""
    wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::JensenCosine)

Computes the wake deficit according to the original Jensen cosine wake model, from the paper:
"A Note on Wind Generator Interaction" by N.O. Jensen (1983)

# Arguments
- `locx::Float`: x coordinate where wind speed is calculated
- `locy::Float`: y coordinate where wind speed is calculated
- `locz::Float`: z coordinate where wind speed is calculated
- `turbine_x::Array(Float)`: vector containing x coordinates for all turbines in farm
- `turbine_y::Array(Float)`: vector containing y coordinates for all turbines in farm
- `turbine_z::Array(Float)`: vector containing z coordinates for all turbines in farm
- `deflection_y::Float`: deflection in the y direction of downstream wake
- `deflection_z::Float`: deflection in the z direction of downstream wake
- `upstream_turbine_id::Int`: index of the upstream wind turbine creating the wake
- `downstream_turbine_id::Int`: index of the downstream turbine feeling the wake (if not referencing a turbine set to zero)
- `hub_height::Array(Float)`: vector containing hub heights for all turbines in farm
- `rotor_diameter::Array(Float)`: vector containing rotor diameters for all turbines in farm
- `turbine_ai::Array(Float)`: vector containing initial velocity deficits for all turbines in farm
- `turbine_local_ti::Array(Float)`: vector containing local turbulence intensities for all turbines in farm
- `turbine_ct::Array(Float)`: vector containing thrust coefficients for all turbines in farm
- `turbine_yaw::Array(Float)`: vector containing the yaw angle? for all turbines in farm
- `model::JensenCosine`: indicates the wake model in use

"""
function wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::JensenCosine)
    """the original Jensen cosine wake model, from the paper: "A Note on Wind
    Generator Interaction" by N.O. Jensen (1983)"""

    # get wec factor (See Thomas and Ning 2018)
    wec_factor = model.wec_factor[1]

    # find delta x, y, and z. dx is the downstream distance from the turbine to
    # the point of interest. dy and dz are the distances from the point of interest
    # and the wake center (in y and z)
    dx = locx-turbine_x[upstream_turbine_id]
    dy = locy-(turbine_y[upstream_turbine_id]+deflection_y)
    dz = locz-(turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

    r0 = rotor_diameter[upstream_turbine_id]/2.0 #turbine rotor radius
    del = sqrt(dy^2+dz^2) #distance from wake center to the point of interest

    if dx < 0.
        loss = 0.0 # no loss outside the wake
    else
        d = wec_factor*r0/tan(model.beta) # distance from fulcrum of wake cone to wind turbine hub
        theta = atan(dy/(dx+d)) # angle from center of wake to point of interest
        if theta > model.beta # if you're outside the wake
            loss = 0.0
        else # if you're inside the wake
            n = pi / model.beta # see Jensen 1983 eq. 3. Value used for n in paper was 9, corresponding to beta = 20.0 deg.
            ftheta = (1.0 + cos(n * theta)) / 2.0 # cosine term to be applied to loss equation as per Jensen 1983
            loss = 2.0*turbine_ai[upstream_turbine_id]*(ftheta*r0/(r0+model.alpha*dx))^2 #equation (2) from the paper
        end
    end

    return loss
end

"""
    wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::Multizone)

Computes the wake deficit at a given location using the original multizone "FLORIS" wake model, from the paper:
"Wind plant power optimization through yaw control using a parametric model for wake effects—a CFD simulation study" by Gebraad et al. (2014)

# Arguments
- `locx::Float`: x coordinate where wind speed is calculated
- `locy::Float`: y coordinate where wind speed is calculated
- `locz::Float`: z coordinate where wind speed is calculated
- `turbine_x::Array(Float)`: vector containing x coordinates for all turbines in farm
- `turbine_y::Array(Float)`: vector containing y coordinates for all turbines in farm
- `turbine_z::Array(Float)`: vector containing z coordinates for all turbines in farm
- `deflection_y::Float`: deflection in the y direction of downstream wake
- `deflection_z::Float`: deflection in the z direction of downstream wake
- `upstream_turbine_id::Int`: index of the upstream wind turbine creating the wake
- `downstream_turbine_id::Int`: index of the downstream turbine feeling the wake (if not referencing a turbine set to zero)
- `hub_height::Array(Float)`: vector containing hub heights for all turbines in farm
- `rotor_diameter::Array(Float)`: vector containing rotor diameters for all turbines in farm
- `turbine_ai::Array(Float)`: vector containing initial velocity deficits for all turbines in farm
- `turbine_local_ti::Array(Float)`: vector containing local turbulence intensities for all turbines in farm
- `turbine_ct::Array(Float)`: vector containing thrust coefficients for all turbines in farm
- `turbine_yaw::Array(Float)`: vector containing the yaw angle? for all turbines in farm
- `model::Multizone`: indicates the wake model in use

"""
function wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::Multizone)

    dt = rotor_diameter[upstream_turbine_id]

    # extract model parameters
    ke = model.ke
    me = model.me
    MU = model.MU
    aU = model.aU
    bU = model.bU

    # aU is given in degrees when it should be radians
    aU = deg2rad(aU)

    if downstream_turbine_id == 0             # Finds the velocity at that point if not looking for a turbine

        # find delta x, y, and z. dx is the downstream distance from the turbine to
        # the point of interest. dy and dz are the distances from the point of interest
        # and the wake center (in y and z)

        dx = locx-turbine_x[upstream_turbine_id]
        dy = locy-(turbine_y[upstream_turbine_id]+deflection_y)
        dz = locz-(turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

        if dx < 0.
            c = 0.0
        else
            del = sqrt(dy^2+dz^2) #distance from wake center to the point of interest

            # calculate the diameter of the wake in each of the three zones (at the specified dx)
            Dw = zeros(3)
            for i = 1:3
                Dw[i] = max(dt+2*ke*me[i]*dx,0) # equation (13) from the paper
            end
            Rw = Dw./2 # radius of the wake zones

            # calculate the coefficient c of the appropriate wake zone:
            if del > Rw[3]
                c = 0.0
            else
                # equations (15, 16, and 17) from the paper
                if del < Rw[1]
                    MUi = MU[1]
                elseif del < Rw[2]
                    MUi = MU[2]
                elseif del < Rw[3]
                    MUi = MU[3]
                end

                mU = MUi/(cos(aU+bU*turbine_yaw[upstream_turbine_id]))
                c = (dt/(dt+2.0*ke*mU*dx))^2
            end

        end

        loss = 2.0*turbine_ai[upstream_turbine_id]*c # calculate the wake loss. Equation (14) fro mthe paper

    else            # this section will calculate the wind effect on an entire turbine.

        # find delta x, y, and z. dx is the downstream distance from the turbine to
        # the point of interest. dy and dz are the distances from the point of interest
        # and the wake center (in y and z)
        dxt = turbine_x[downstream_turbine_id]-turbine_x[upstream_turbine_id]
        dyt = turbine_y[downstream_turbine_id]-(turbine_y[upstream_turbine_id]+deflection_y)
        dzt = (turbine_z[downstream_turbine_id]+hub_height[downstream_turbine_id])-(turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

        posz = turbine_z[downstream_turbine_id] + hub_height[downstream_turbine_id] # finds the z coordinate of the turbine hub

        if dxt < 0.
            c1,c2,c3 = 0,0,0
        else
            del = sqrt(dxt^2+dzt^2) #distance from wake center to the point of interest

            wake_center_y = (turbine_y[upstream_turbine_id]+deflection_y)
            wake_center_z = (turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

            # calculate the diameter of the wake in each of the three zones (at the specified dx)
            # after calculating the size of the wake the model will then find its overlap with the turbine

            Dw = zeros(3)
            ovlp = zeros(3)
            for i = 1:3
                Dw[i] = max(dt+2*ke*me[i]*dxt,0) # equation (13) from the paper
                overlap = overlap_area_func(turbine_y[downstream_turbine_id], posz,
                rotor_diameter[downstream_turbine_id], wake_center_y, wake_center_z, Dw[i]; tol=1E-6)
                area = pi*(rotor_diameter[downstream_turbine_id]^2)/4
                ovlp[i] = overlap/area
            end

            Rw = Dw./2 # radius of the wake zones

            # equations (15, 16, and 17) from the paper calculated for all 3 zones
            # losses computed for each zone and the subsequently added together

            mU1 = MU[1]/(cos(aU+bU*turbine_yaw[upstream_turbine_id]))
            c1 = (dt/(dt+2.0*ke*mU1*dxt))^2

            mU2 = MU[2]/(cos(aU+bU*turbine_yaw[upstream_turbine_id]))
            c2 = (dt/(dt+2.0*ke*mU2*dxt))^2

            mU3 = MU[3]/(cos(aU+bU*turbine_yaw[upstream_turbine_id]))
            c3 = (dt/(dt+2.0*ke*mU3*dxt))^2

        end

        # losses calculated for all 3 zones with equation [14] from the paper
        loss1 = 2.0*turbine_ai[upstream_turbine_id]*c1*ovlp[1]
        loss2 = 2.0*turbine_ai[upstream_turbine_id]*c2*ovlp[2]
        loss3 = 2.0*turbine_ai[upstream_turbine_id]*c3*ovlp[3]

        # losses of 3 zones are summed to get total loss
        loss = loss1 + loss2 + loss3
    end

    return loss
end

"""
    wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussOriginal)

Computes the wake deficit at a given location using the Gaussian wake model presented by Bastankhah and Porte-Agel in the paper: "A new analytical model for wind-turbine wakes" (2014)

# Arguments
- `locx::Float`: x coordinate where wind speed is calculated
- `locy::Float`: y coordinate where wind speed is calculated
- `locz::Float`: z coordinate where wind speed is calculated
- `turbine_x::Array(Float)`: vector containing x coordinates for all turbines in farm
- `turbine_y::Array(Float)`: vector containing y coordinates for all turbines in farm
- `turbine_z::Array(Float)`: vector containing z coordinates for all turbines in farm
- `deflection_y::Float`: deflection in the y direction of downstream wake
- `deflection_z::Float`: deflection in the z direction of downstream wake
- `upstream_turbine_id::Int`: index of the upstream wind turbine creating the wake
- `downstream_turbine_id::Int`: index of the downstream turbine feeling the wake (if not referencing a turbine set to zero)
- `hub_height::Array(Float)`: vector containing hub heights for all turbines in farm
- `rotor_diameter::Array(Float)`: vector containing rotor diameters for all turbines in farm
- `turbine_ai::Array(Float)`: vector containing initial velocity deficits for all turbines in farm
- `turbine_local_ti::Array(Float)`: vector containing local turbulence intensities for all turbines in farm
- `turbine_ct::Array(Float)`: vector containing thrust coefficients for all turbines in farm
- `turbine_yaw::Array(Float)`: vector containing the yaw angle? for all turbines in farm
- `model::GaussOriginal`: indicates the wake model in use

"""
function wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussOriginal)

    dx = locx-turbine_x[upstream_turbine_id]
    dy = locy-(turbine_y[upstream_turbine_id]+deflection_y)
    dz = locz-(turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

    # extract turbine properties
    dt = rotor_diameter[upstream_turbine_id]
    yaw = turbine_yaw[upstream_turbine_id]
    ct = turbine_ct[upstream_turbine_id]

    as = model.alpha_star
    bs = model.beta_star
    ti = model.turbulence_intensity

    # extract model parameters
    ks = model.k_star       # wake spread rate (k* in 2014 paper)
    ky = model.horizontal_spread_rate
    kz = model.vertical_spread_rate
    as = model.alpha_star
    bs = model.beta_star

    ti = turbine_local_ti[upstream_turbine_id]

    # calculate beta (paper eq: 6)
    beta = 0.5*(1.0+sqrt(1.0-ct))/sqrt(1.0-ct)

    # calculate the length of the potential core (2016 paper eq: 7.3)
    x0 = ((1.0+sqrt(1.0+ct)))/(sqrt(2.0)*as*ti+bs*(1.0-sqrt(1.0-ct)))

    # calculate loss (paper eq: 23)
    enum =((dz/dt)^2+(dy/dt)^2)
    if dx > x0
        denom = (ks*dx/dt+0.2*sqrt(beta))^2
    else
        denom = (ks*x0/dt+0.2*sqrt(beta))^2
    end

    loss = (1.0 - sqrt(1.0-ct/(8.0*denom)))*exp(-enum/(2.0*denom))

end

"""
    _gauss_yaw_potential_core(dt, yaw, ct, as, ti, bs)

Helper function for wake_deficit_model when using the GaussYaw model. Computes the length of the near wake potential core.
"""
function _gauss_yaw_potential_core(d, yaw, ct, as, ti, bs)
    # from Bastankhah and Porte-Agel 2016 eqn 7.3
    x0 = d*(cos(yaw)*(1.0+sqrt(1.0-ct)))/(sqrt(2.0)*(as*ti+bs*(1.0-sqrt(1.0-ct))))
    return x0
end

"""
    _gauss_yaw_spread(dt, k, dx, x0, yaw)

Helper function for wake_deficit_model when using the GaussYaw model. Computes the standard deviation of the wake.
"""
function _gauss_yaw_spread(dt, k, dx, x0, yaw)
    # from Bastankhah and Porte-Agel 2016 eqn 7.2
    sigma = dt*(k*(dx-x0)/dt+cos(yaw)/sqrt(8.0))

    return sigma

end

function _gauss_yaw_model_deficit(dx, dy, dz, dt, yaw, ct, ti, as, bs, ky, kz, wf)

    if dx > 0.0 # loss in the wake

        # calculate the length of the potential core (paper eq: 7.3)
        x0 = _gauss_yaw_potential_core(dt, yaw, ct, as, ti, bs)

        # calculate wake spread
        if dx > x0
            # calculate horizontal wake spread (paper eq: 7.2)
            sigma_y = _gauss_yaw_spread(dt, ky, dx, x0, yaw)

            # calculate vertical wake spread (paper eq: 7.2)
            sigma_z = _gauss_yaw_spread(dt, kz, dx, x0, 0.0)

        else # linear interpolation in the near wakes
            # calculate horizontal wake spread
            sigma_y = _gauss_yaw_spread(dt, ky, x0, x0, yaw)

            # calculate vertical wake spread
            sigma_z = _gauss_yaw_spread(dt, kz, x0, x0, 0.0)
        end

        # calculate velocity deficit
        ey = exp(-0.5*(dy/(wf*sigma_y))^2)
        ez = exp(-0.5*(dz/(wf*sigma_z))^2)

        if (1.0-ct*cos(yaw)/(8.0*(sigma_y*sigma_z/dt^2))) >= 1e-8
            loss = (1.0-sqrt(1.0-ct*cos(yaw)/(8.0*(sigma_y*sigma_z/dt^2))))*ey*ez
        else
            loss = ey*ez
        end
    else # loss upstream of turbine
        loss = 0.0
    end

    return loss

end

"""
    wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussYaw)

Computes the wake deficit at a given location using the The Gaussian wake model presented by Bastankhah and Porte-Agel in the paper: "Experimental and theoretical study of wind turbine wakes in yawed conditions" (2016)

# Arguments
- `locx::Float`: x coordinate where wind speed is calculated
- `locy::Float`: y coordinate where wind speed is calculated
- `locz::Float`: z coordinate where wind speed is calculated
- `turbine_x::Array(Float)`: vector containing x coordinates for all turbines in farm
- `turbine_y::Array(Float)`: vector containing y coordinates for all turbines in farm
- `turbine_z::Array(Float)`: vector containing z coordinates for all turbines in farm
- `deflection_y::Float`: deflection in the y direction of downstream wake
- `deflection_z::Float`: deflection in the z direction of downstream wake
- `upstream_turbine_id::Int`: index of the upstream wind turbine creating the wake
- `downstream_turbine_id::Int`: index of the downstream turbine feeling the wake (if not referencing a turbine set to zero)
- `hub_height::Array(Float)`: vector containing hub heights for all turbines in farm
- `rotor_diameter::Array(Float)`: vector containing rotor diameters for all turbines in farm
- `turbine_ai::Array(Float)`: vector containing initial velocity deficits for all turbines in farm
- `turbine_local_ti::Array(Float)`: vector containing local turbulence intensities for all turbines in farm
- `turbine_ct::Array(Float)`: vector containing thrust coefficients for all turbines in farm
- `turbine_yaw::Array(Float)`: vector containing the yaw angle? for all turbines in farm
- `model::GaussYaw`: indicates the wake model in use

"""
function wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussYaw)

    dx = locx-turbine_x[upstream_turbine_id]
    dy = locy-(turbine_y[upstream_turbine_id]+deflection_y)
    dz = locz-(turbine_z[upstream_turbine_id]+hub_height[upstream_turbine_id]+deflection_z)

    # extract turbine properties
    dt = rotor_diameter[upstream_turbine_id]
    yaw = turbine_yaw[upstream_turbine_id]
    ct = turbine_ct[upstream_turbine_id]

    # extract model parameters
    # ks = model.k_star       # wake spread rate (k* in 2014 paper)
    ti = turbine_local_ti[upstream_turbine_id]
    ky = model.horizontal_spread_rate
    kz = model.vertical_spread_rate
    as = model.alpha_star
    bs = model.beta_star
    wec_factor = model.wec_factor[1]

    loss = _gauss_yaw_model_deficit(dx, dy, dz, dt, yaw, ct, ti, as, bs, ky, kz, wec_factor)

    return loss

end

"""
    wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussYawVariableSpread)

Computes the wake deficit at a given location using the The Gaussian wake model presented by Bastankhah and Porte-Agel in the paper: "Experimental and theoretical study of wind turbine wakes in yawed conditions" (2016)
The spread rate is adjusted based on local turbulence intensity as in Niayifar and Porte-Agel 2016

# Arguments
- `locx::Float`: x coordinate where wind speed is calculated
- `locy::Float`: y coordinate where wind speed is calculated
- `locz::Float`: z coordinate where wind speed is calculated
- `turbine_x::Array(Float)`: vector containing x coordinates for all turbines in farm
- `turbine_y::Array(Float)`: vector containing y coordinates for all turbines in farm
- `turbine_z::Array(Float)`: vector containing z coordinates for all turbines in farm
- `deflection_y::Float`: deflection in the y direction of downstream wake
- `deflection_z::Float`: deflection in the z direction of downstream wake
- `upstream_turbine_id::Int`: index of the upstream wind turbine creating the wake
- `downstream_turbine_id::Int`: index of the downstream turbine feeling the wake (if not referencing a turbine set to zero)
- `hub_height::Array(Float)`: vector containing hub heights for all turbines in farm
- `rotor_diameter::Array(Float)`: vector containing rotor diameters for all turbines in farm
- `turbine_ai::Array(Float)`: vector containing initial velocity deficits for all turbines in farm
- `turbine_local_ti::Array(Float)`: vector containing local turbulence intensities for all turbines in farm
- `turbine_ct::Array(Float)`: vector containing thrust coefficients for all turbines in farm
- `turbine_yaw::Array(Float)`: vector containing the yaw angle? for all turbines in farm
- `model::GaussYawVariableSpread`: indicates the wake model in use

"""
function wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussYawVariableSpread)

    dx = locx-turbine_x[upstream_turbine_id]
    dy = locy-(turbine_y[upstream_turbine_id]+deflection_y)
    dz = locz-(turbine_z[upstream_turbine_id]+hub_height[1]+deflection_z)

    # extract turbine properties
    dt = rotor_diameter[upstream_turbine_id]
    yaw = turbine_yaw[upstream_turbine_id]
    ct = turbine_ct[upstream_turbine_id]

    # extract model parameters
    # ks = model.k_star       # wake spread rate (k* in 2014 paper)
    ti = turbine_local_ti[upstream_turbine_id]
    ky = kz = _k_star_func(ti, model.k1, model.k2)

    as = model.alpha_star
    bs = model.beta_star
    wec_factor = model.wec_factor[1]

    loss = _gauss_yaw_model_deficit(dx, dy, dz, dt, yaw, ct, ti, as, bs, ky, kz, wec_factor)

    return loss
end

"""
    wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussSimple)

Computes the wake deficit at a given location using the Gaussian wake model presented by Bastankhah and Porte-Agel in the paper: "A new analytical model for wind-turbine wakes" (2014)
    as modified for IEA Task 37 Case Studies 3 and 4

# Arguments
- `locx::Float`: x coordinate where wind speed is calculated
- `locy::Float`: y coordinate where wind speed is calculated
- `locz::Float`: z coordinate where wind speed is calculated
- `turbine_x::Array(Float)`: vector containing x coordinates for all turbines in farm
- `turbine_y::Array(Float)`: vector containing y coordinates for all turbines in farm
- `turbine_z::Array(Float)`: vector containing z coordinates for all turbines in farm
- `deflection_y::Float`: deflection in the y direction of downstream wake
- `deflection_z::Float`: deflection in the z direction of downstream wake
- `upstream_turbine_id::Int`: index of the upstream wind turbine creating the wake
- `downstream_turbine_id::Int`: index of the downstream turbine feeling the wake (if not referencing a turbine set to zero)
- `hub_height::Array(Float)`: vector containing hub heights for all turbines in farm
- `rotor_diameter::Array(Float)`: vector containing rotor diameters for all turbines in farm
- `turbine_ai::Array(Float)`: vector containing initial velocity deficits for all turbines in farm
- `turbine_local_ti::Array(Float)`: vector containing local turbulence intensities for all turbines in farm
- `turbine_ct::Array(Float)`: vector containing thrust coefficients for all turbines in farm
- `turbine_yaw::Array(Float)`: vector containing the yaw angle? for all turbines in farm
- `model::GaussSimple`: indicates the wake model in use

"""
function wake_deficit_model(locx, locy, locz, turbine_x, turbine_y, turbine_z, deflection_y, deflection_z, upstream_turbine_id, downstream_turbine_id, hub_height, rotor_diameter, turbine_ai, turbine_local_ti, turbine_ct, turbine_yaw, model::GaussSimple)

    dx = locx-turbine_x[upstream_turbine_id]
    dy = locy-(turbine_y[upstream_turbine_id]+deflection_y)

    # extract turbine properties
    dt = rotor_diameter[upstream_turbine_id]
    ct = turbine_ct[upstream_turbine_id]

    # extract model properties
    k = model.k
    wf = model.wec_factor[1]

    # calculate loss
    sigmay = k*dx+dt/sqrt(8.0)
    radical = 1.0-ct/(8.0*(sigmay^2)/(dt^2))
    exponent = -0.5*(dy/(wf*sigmay))^2
    loss = (1.0 - sqrt(radical))*exp(exponent)

    return loss

end


text:
"""
    estimating_function_template(nobs::Function,
                                 ef_contribution::Function)

[Composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) for defining an `estimating_function_template`.

Arguments
===
+ `nobs`: a function of `data` that computes the number of observations of the particular data type,
+ `ef_contribution`: a function of the parameters `theta`, the `data` and the observation index `i` that returns a vector of length `length(theta)`.

Result
===
An `estimating_function_template` object with fields `nobs` and `obj_contributions`.
"""
struct estimating_function_template
    nobs::Function
    ef_contribution::Function
end

"""
    estimating_function(theta::Vector,
                        data::Any,
                        template::estimating_function_template,
                        br::Bool = false,
                        concentrate::Vector{Int64} = Vector{Int64}())

Evaluate a vector of estimating functions at `theta` by adding up all contributions in `data`, according to an [`estimating_function_template`](@ref).

Arguments
===
+ `theta`: a `Vector` of parameter values at which to evaluate the estimating functions
+ `data`: typically an object of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) with all the data required to compute the `estimating_function`.
+ `template`: an [`estimating_function_template`](@ref) object.
+ `br`: a `Bool`. If `false` (default), the estimating functions is constructed by adding up all contributions in
`data`, according to [`estimating_function_template`](@ref), before it is evaluated at `theta`. If `true` then the empirical bias-reducing adjustments in [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786) are computed and added to the estimating functions.
+ `concentrate`: a `Vector{Int64}`; if specified, empirical bias-reducing adjustments are added only to the subset of estimating functions indexed by `concentrate`. The default is to add empirical bias-reducing adjustments to all estimating functions.

Result
===
A `Vector`.

Details
===
`data` can be used to pass additional constants other than the actual data to the objective.
"""
function estimating_function(theta::Vector,
                             data::Any,
                             template::estimating_function_template,
                             br::Bool = false,
                             concentrate::Vector{Int64} = Vector{Int64}())
    p = length(theta)
    n_obs = template.nobs(data)
    ef = zeros(p)
    for i in 1:n_obs
        ef += template.ef_contribution(theta, data, i)
    end
    if (br)
        quants = ef_quantities(theta, data, template, br, concentrate)
        ef + quants[1]
    else
        ef
    end
end

"""
    get_estimating_function(data::Any,
                            template::estimating_function_template,
                            br::Bool = false,
                            concentrate::Vector{Int64} = Vector{Int64}(),
                            regularizer::Any = Vector{Int64}())

Construct the estimating functions by adding up all contributions in the `data` according to [`estimating_function_template`](@ref).

Arguments
===
+ `data`: typically an object of [composite type](https://docs.julialang.org/en/v1/manual/types/#Composite-Types-1) with all the data required to compute the `estimating_function`.
+ `template`: an [`estimating_function_template`](@ref) object.
+ `br`: a `Bool`. If `false` (default), the estimating functions is constructed by adding up all contributions in
`data`, according to [`estimating_function_template`](@ref), before it is evaluated at `theta`. If `true` then the empirical bias-reducing adjustments in [Kosmidis & Lunardon, 2020](http://arxiv.org/abs/2001.03786) are computed and added to the estimating functions.
+ `concentrate`: a `Vector{Int64}`; if specified, empirical bias-reducing adjustments are added only to the subset of estimating functions indexed by `concentrate`. The default is to add empirical bias-reducing adjustments to all estimating functions.
+ `regularizer`: a function of `theta` and `data` returning a `Vector` of dimension equal to the number of the estimating functions, which is added to the (bias-reducing) estimating function; the default value will result in no regularization.

Result
===
An in-place function that stores the value of the estimating functions inferred from `template`, in a preallocated vector passed as its first argument, ready to be used withing `NLsolve.nlsolve`. This is the in-place version of [`estimating_function`](@ref) with the extra `regularizer` argument.
"""
function get_estimating_function(data::Any,
                                 template::estimating_function_template,
                                 br::Bool = false,
                                 concentrate::Vector{Int64} = Vector{Int64}(),
                                 regularizer::Any = Vector{Int64}())
    ## regulizer here has different type and default than fit, because
    ## the dimension of theta cannot be inferred
    has_regularizer = typeof(regularizer) <: Function
    function (F, theta::Vector)
        if has_regularizer
            out = estimating_function(theta, data, template, br, concentrate) + regularizer(theta, data)
        else
             out = estimating_function(theta, data, template, br, concentrate)
        end
        for i in 1:length(out)
            F[i] = out[i]
        end
    end
end

function ef_quantities(theta::Vector,
                       data::Any,
                       template::estimating_function_template,
                       adjustment::Bool = false,
                       concentrate::Vector{Int64} = Vector{Int64}())
    estfun_i = (pars, i) -> template.ef_contribution(pars, data, i)
    function ja_i(eta::Vector, i::Int)
        out = similar(eta, p, p)
        ForwardDiff.jacobian!(out, pars -> estfun_i(pars, i), eta)
    end
    p = length(theta)
    n_obs = template.nobs(data)
    psi = zeros(p)
    emat = zeros(p, p)
    jmat = zeros(p, p)

    if adjustment
        function u(eta::Vector, i::Int)
            out = similar(eta, p * p, p)
            ForwardDiff.jacobian!(out, beta -> ja_i(beta, i), eta)
        end
        psi2 = Vector(undef, p)
        for j in 1:p
            psi2[j] = zeros(p, p)
        end
        # psi2 = zeros(p, p, p)
        umat = zeros(p * p, p)
        for i in 1:n_obs
            cpsi = estfun_i(theta, i)
            psi += cpsi
            emat += cpsi * cpsi'
            jaco = ja_i(theta, i)
            jmat += -jaco
            umat += u(theta, i)
            for j in 1:p
                psi2[j] += jaco[j, :] * cpsi'
            end
        end
    else
        for i in 1:n_obs
            cpsi = estfun_i(theta, i)
            psi += cpsi
            emat += cpsi * cpsi'
            jmat += - ja_i(theta, i)
        end
    end

    jmat_inv = try
        inv(jmat)
    catch
        fill(NaN, p, p)
    end
    vcov = jmat_inv * (emat * jmat_inv')

    if adjustment
        Afun(j::Int64) = -tr(jmat_inv * psi2[j] + vcov * umat[j:p:(p * p - p + j), :] / 2)
        A = map(Afun, 1:p)
        ## if concentrate then redefine A
        if length(concentrate) > 0
            if any((concentrate .> p) .| (concentrate .< 1))
                error(concentrate, " should be a vector of integers between ", 1, " and ", p)
            else
                ist = concentrate
                nce = deleteat!(collect(1:p), concentrate)
                A_ist = A[ist]
                A_nce = A[nce]
                A = vcat(A_ist + inv(jmat_inv[ist, ist]) * jmat_inv[ist, nce] * A_nce,
                         zeros(length(nce)))
            end
        end
        [A, jmat_inv, emat, psi]
    else
        [vcov, jmat_inv, emat, psi]
    end
end

ext: jl


text:
"""
    TurbulenceClosures

Functions for turbulence, sub-grid scale modelling. These include
viscosity terms, diffusivity and stress tensors.

- [`ConstantViscosity`](@ref)
- [`ViscousSponge`](@ref)
- [`SmagorinskyLilly`](@ref)
- [`Vreman`](@ref)
- [`AnisoMinDiss`](@ref)

"""
module TurbulenceClosures

# ## Turbulence Closures
# In `turbulence.jl` we specify turbulence closures. Currently,
# pointwise models of the eddy viscosity/eddy diffusivity type are
# supported for turbulent shear and tracer diffusivity. Methods currently supported
# are:\
# [`ConstantViscosity`](@ref constant-viscosity)\
# [`ViscousSponge`](@ref viscous-sponge)\
# [`SmagorinskyLilly`](@ref smagorinsky-lilly)\
# [`Vreman`](@ref vreman)\
# [`AnisoMinDiss`](@ref aniso-min-diss)\

#md # !!! note
#md #     Usage: This is a quick-ref guide to using turbulence models as a subcomponent
#md #     of `BalanceLaw` \
#md #     $\nu$ is the kinematic viscosity, $C_smag$ is the Smagorinsky Model coefficient,
#md #     `turbulence=ConstantDynamicViscosity(ρν)`\
#md #     `turbulence=ConstantKinematicViscosity(ν)`\
#md #     `turbulence=ViscousSponge(ν, z_max, z_sponge, α, γ)`\
#md #     `turbulence=SmagorinskyLilly(C_smag)`\
#md #     `turbulence=Vreman(C_smag)`\
#md #     `turbulence=AnisoMinDiss(C_poincare)`

using ClimateMachine
using DocStringExtensions
using LinearAlgebra
using StaticArrays

import CLIMAParameters: AbstractParameterSet
import ClimateMachine.Mesh.Geometry:
    LocalGeometry, resolutionmetric, lengthscale

using ClimateMachine.Orientations
using ClimateMachine.VariableTemplates
using ClimateMachine.BalanceLaws
using CLIMAParameters.Atmos.SubgridScale: inv_Pr_turb

import ClimateMachine.BalanceLaws:
    vars_state,
    flux_first_order!,
    flux_second_order!,
    source!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    transform_post_gradient_laplacian!,
    init_state_prognostic!,
    update_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

export TurbulenceClosureModel,
    ConstantViscosity,
    ConstantDynamicViscosity,
    ConstantKinematicViscosity,
    SmagorinskyLilly,
    Vreman,
    AnisoMinDiss,
    HyperDiffusion,
    NoHyperDiffusion,
    DryBiharmonic,
    EquilMoistBiharmonic,
    NoViscousSponge,
    UpperAtmosSponge,
    turbulence_tensors,
    init_aux_turbulence!,
    init_aux_hyperdiffusion!,
    turbulence_nodal_update_auxiliary_state!,
    sponge_viscosity_modifier

# ### Abstract Type
# We define a `TurbulenceClosureModel` abstract type and
# default functions for the generic turbulence closure
# which will be overloaded with model specific functions.

"""
    Abstract type with default do-nothing behaviour for
arbitrary turbulence closure models.
"""
abstract type TurbulenceClosureModel end

vars_state(::TurbulenceClosureModel, ::AbstractStateType, FT) = @vars()

"""
    ConstantViscosity <: TurbulenceClosureModel
Abstract type for constant viscosity models
"""
abstract type ConstantViscosity <: TurbulenceClosureModel end

"""
    Abstract type for Hyperdiffusion models
"""
abstract type HyperDiffusion end

"""
    Abstract type for viscous sponge layers.
Modifier for viscosity computed from existing turbulence closures.
"""
abstract type ViscousSponge end

"""
    init_aux_turbulence!
Initialise auxiliary variables for turbulence models.
Overload for specific turbulence closure type.
"""
function init_aux_turbulence!(
    ::TurbulenceClosureModel,
    ::BalanceLaw,
    aux::Vars,
    geom::LocalGeometry,
) end

"""
    turbulence_nodal_update_auxiliary_state!
Update auxiliary variables for turbulence models.
Overload for specific turbulence closure type.
"""
function turbulence_nodal_update_auxiliary_state!(
    ::TurbulenceClosureModel,
    ::BalanceLaw,
    state::Vars,
    aux::Vars,
    t::Real,
) end

"""
    compute_gradient_argument!
Assign pre-gradient-transform variables specific to turbulence models.
"""
function compute_gradient_argument!(
    ::TurbulenceClosureModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
"""
    compute_gradient_flux!(::TurbulenceClosureModel, _...)
Post-gradient-transformed variables specific to turbulence models.
"""
function compute_gradient_flux!(
    ::TurbulenceClosureModel,
    ::Orientation,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end

# Fallback functions for hyperdiffusion model
vars_state(::HyperDiffusion, ::AbstractStateType, FT) = @vars()

function init_aux_hyperdiffusion!(
    ::HyperDiffusion,
    ::BalanceLaw,
    aux::Vars,
    geom::LocalGeometry,
) end
function hyperdiffusion_nodal_update_auxiliary_state!(
    ::HyperDiffusion,
    ::BalanceLaw,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_argument!(
    ::HyperDiffusion,
    ::BalanceLaw,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function transform_post_gradient_laplacian!(
    h::HyperDiffusion,
    bl::BalanceLaw,
    hyperdiffusive::Vars,
    gradvars::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_second_order!(
    h::HyperDiffusion,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    h::HyperDiffusion,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end

function turbulence_tensors end
function hyperviscosity_tensors end

"""
    ν, D_t, τ = turbulence_tensors(
                    ::TurbulenceClosureModel,
                    orientation::Orientation,
                    param_set::AbstractParameterSet,
                    state::Vars,
                    diffusive::Vars,
                    aux::Vars,
                    t::Real
                )
Compute the kinematic viscosity (`ν`), the diffusivity (`D_t`) and SGS momentum flux tensor (`τ`)
for a given turbulence closure. Each closure overloads this method with the appropriate calculations
for the returned quantities.

# Arguments

- `::TurbulenceClosureModel` = Struct identifier for turbulence closure model
- `orientation` = `BalanceLaw.orientation`
- `param_set` parameter set
- `state` = Array of prognostic (state) variables. See `vars_state` in `BalanceLaw`
- `diffusive` = Array of diffusive variables
- `aux` = Array of auxiliary variables
- `t` = time
"""
function turbulence_tensors(
    m::TurbulenceClosureModel,
    bl::BalanceLaw,
    state,
    diffusive,
    aux,
    t,
)
    ν, D_t, τ = turbulence_tensors(
        m,
        bl.orientation,
        bl.param_set,
        state,
        diffusive,
        aux,
        t,
    )
    ν, D_t, τ = sponge_viscosity_modifier(bl, bl.viscoussponge, ν, D_t, τ, aux)
    return (ν, D_t, τ)
end

# We also provide generic math functions for use within the turbulence closures,
# commonly used quantities such as the [principal tensor invariants](@ref tensor-invariants), handling of
# [symmetric tensors](@ref symmetric-tensors) and [tensor norms](@ref tensor-norms)are addressed.

# ### [Pricipal Invariants](@id tensor-invariants)
# ```math
# \textit{I}_{1} = \mathrm{tr(X)} \\
# \textit{I}_{2} = (\mathrm{tr(X)}^2 - \mathrm{tr(X^2)}) / 2 \\
# \textit{I}_{3} = \mathrm{det(X)} \\
# ```
"""
    principal_invariants(X)

Calculates principal invariants of a tensor `X`. Returns 3 element tuple containing the invariants.
"""
function principal_invariants(X)
    first = tr(X)
    second = (first^2 - tr(X .^ 2)) / 2
    third = det(X)
    return (first, second, third)
end

# ### [Symmetrize](@id symmetric-tensors)
# ```math
# \frac{\mathrm{X} + \mathrm{X}^{T}}{2} \\
# ```
"""
    symmetrize(X)

Given a (3,3) second rank tensor X, compute `(X + X')/2`, returning a
symmetric `SHermitianCompact` object.
"""
function symmetrize(X::StaticArray{Tuple{3, 3}})
    SHermitianCompact(SVector(
        X[1, 1],
        (X[2, 1] + X[1, 2]) / 2,
        (X[3, 1] + X[1, 3]) / 2,
        X[2, 2],
        (X[3, 2] + X[2, 3]) / 2,
        X[3, 3],
    ))
end

# ### [2-Norm](@id tensor-norms)
# Given a tensor X, return the tensor dot product
# ```math
# \sum_{i,j} S_{ij}^2
# ```
"""
    norm2(X)
Given a tensor `X`, computes X:X.
"""
function norm2(X::SMatrix{3, 3, FT}) where {FT}
    abs2(X[1, 1]) +
    abs2(X[2, 1]) +
    abs2(X[3, 1]) +
    abs2(X[1, 2]) +
    abs2(X[2, 2]) +
    abs2(X[3, 2]) +
    abs2(X[1, 3]) +
    abs2(X[2, 3]) +
    abs2(X[3, 3])
end
function norm2(X::SHermitianCompact{3, FT, 6}) where {FT}
    abs2(X[1, 1]) +
    2 * abs2(X[2, 1]) +
    2 * abs2(X[3, 1]) +
    abs2(X[2, 2]) +
    2 * abs2(X[3, 2]) +
    abs2(X[3, 3])
end

# ### [Strain-rate Magnitude](@id strain-rate-magnitude)
# By definition, the strain-rate magnitude, as defined in
# standard turbulence modelling is computed such that
# ```math
# |\mathrm{S}| = \sqrt{2 \sum_{i,j} \mathrm{S}_{ij}^2}
# ```
# where
# ```math
# \vec{S}(\vec{u}) = \frac{1}{2}  \left(\nabla\vec{u} +  \left( \nabla\vec{u} \right)^T \right)
# ```
# \mathrm{S} is the rate-of-strain tensor. (Symmetric component of the velocity gradient). Note that the
# skew symmetric component (rate-of-rotation) is not currently computed.
"""
    strain_rate_magnitude(S)
Given the rate-of-strain tensor `S`, computes its magnitude.
"""
function strain_rate_magnitude(S::SHermitianCompact{3, FT, 6}) where {FT}
    return sqrt(2 * norm2(S))
end

"""
    WithDivergence
A divergence type which includes the divergence term in the momentum flux tensor
"""
struct WithDivergence end
export WithDivergence
"""
    WithoutDivergence
A divergence type which does not include the divergence term in the momentum flux tensor
"""
struct WithoutDivergence end
export WithoutDivergence

# ### [Constant Viscosity Model](@id constant-viscosity)
# `ConstantViscosity` requires a user to specify the constant viscosity (dynamic or kinematic)
# and appropriately computes the turbulent stress tensor based on this term. Diffusivity can be
# computed using the turbulent Prandtl number for the appropriate problem regime.
# ```math
# \tau =
#     \begin{cases}
#     - 2 \nu \mathrm{S} & \mathrm{WithoutDivergence},\\
#     - 2 \nu \mathrm{S} + \frac{2}{3} \nu \mathrm{tr(S)} I_3 & \mathrm{WithDivergence}.
#     \end{cases}
# ```

"""
    ConstantDynamicViscosity <: ConstantViscosity

Turbulence with constant dynamic viscosity (`ρν`).
Divergence terms are included in the momentum flux tensor if divergence_type is WithDivergence.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct ConstantDynamicViscosity{FT, DT} <: ConstantViscosity
    "Dynamic Viscosity [kg/m/s]"
    ρν::FT
    divergence_type::DT
    function ConstantDynamicViscosity(
        ρν::FT,
        divergence_type::Union{WithDivergence, WithoutDivergence} = WithoutDivergence(),
    ) where {FT}
        return new{FT, typeof(divergence_type)}(ρν, divergence_type)
    end
end

"""
    ConstantKinematicViscosity <: ConstantViscosity

Turbulence with constant kinematic viscosity (`ν`).
Divergence terms are included in the momentum flux tensor if divergence_type is WithDivergence.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct ConstantKinematicViscosity{FT, DT} <: ConstantViscosity
    "Kinematic Viscosity [m2/s]"
    ν::FT
    divergence_type::DT
    function ConstantKinematicViscosity(
        ν::FT,
        divergence_type::Union{WithDivergence, WithoutDivergence} = WithoutDivergence(),
    ) where {FT}
        return new{FT, typeof(divergence_type)}(ν, divergence_type)
    end
end

vars_state(::ConstantViscosity, ::Gradient, FT) = @vars()
vars_state(::ConstantViscosity, ::GradientFlux, FT) =
    @vars(S::SHermitianCompact{3, FT, 6})

function compute_gradient_flux!(
    ::ConstantViscosity,
    ::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
end

compute_stress(div_type::WithoutDivergence, ν, S) = (-2 * ν) * S
compute_stress(div_type::WithDivergence, ν, S) =
    (-2 * ν) * S + (2 * ν / 3) * tr(S) * I

function turbulence_tensors(
    m::ConstantDynamicViscosity,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    S = diffusive.turbulence.S
    ν = m.ρν / state.ρ
    D_t = ν * _inv_Pr_turb
    τ = compute_stress(m.divergence_type, ν, S)
    return ν, D_t, τ
end

function turbulence_tensors(
    m::ConstantKinematicViscosity,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    S = diffusive.turbulence.S
    ν = m.ν
    D_t = ν * _inv_Pr_turb
    τ = compute_stress(m.divergence_type, ν, S)
    return ν, D_t, τ
end

# ### [Smagorinsky-Lilly](@id smagorinsky-lilly)
# The Smagorinsky turbulence model, with Lilly's correction to
# stratified atmospheric flows, is included in ClimateMachine.
# The input parameter to this model is the Smagorinsky coefficient.
# For atmospheric flows, the coefficient `C_smag` typically takes values between
# 0.15 and 0.23. Flow dependent `C_smag` are currently not supported (e.g. Germano's
# extension). The Smagorinsky-Lilly model does not contain explicit filtered terms.
# #### Equations
# ```math
# \nu = (C_{s} \mathrm{f}_{b} \Delta)^2 \sqrt{|\mathrm{S}|}
# ```
# with the stratification correction term
# ```math
# f_{b} =
#    \begin{cases}
#    1 & \mathrm{Ri} \leq 0 ,\\
#    \max(0, 1 - \mathrm{Ri} / \mathrm{Pr}_{t})^{1/4} & \mathrm{Ri} > 0 .
#    \end{cases}
# ```
# ```math
# \mathrm{Ri} =  \frac{N^2}{{|S|}^2}
# ```
# ```math
# N = \left( \frac{g}{\theta_v} \frac{\partial \theta_v}{\partial z}\right)^{1/2}
# ```
# Here, $\mathrm{Ri}$ and $\mathrm{Pr}_{t}$ are the Richardson and
# turbulent Prandtl numbers respectively.  $\Delta$ is the mixing length in the
# relevant coordinate direction. We use the DG metric terms to determine the
# local effective resolution (see `src/Mesh/Geometry.jl`), and modify the vertical lengthscale by the
# stratification correction factor $\mathrm{f}_{b}$ so that $\Delta_{vert} = \Delta z f_b$.

"""
    SmagorinskyLilly <: TurbulenceClosureModel

# Fields

$(DocStringExtensions.FIELDS)

# Smagorinsky Model Reference

See [Smagorinsky1963](@cite)

# Lilly Model Reference

See [Lilly1962](@cite)

# Brunt-Väisälä Frequency Reference

See [Durran1982](@cite)

"""
struct SmagorinskyLilly{FT} <: TurbulenceClosureModel
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end

vars_state(::SmagorinskyLilly, ::Auxiliary, FT) = @vars(Δ::FT)
vars_state(::SmagorinskyLilly, ::Gradient, FT) = @vars(θ_v::FT)
vars_state(::SmagorinskyLilly, ::GradientFlux, FT) =
    @vars(S::SHermitianCompact{3, FT, 6}, N²::FT)

function init_aux_turbulence!(
    ::SmagorinskyLilly,
    ::BalanceLaw,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end

function compute_gradient_argument!(
    m::SmagorinskyLilly,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end

function compute_gradient_flux!(
    ::SmagorinskyLilly,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::SmagorinskyLilly,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    S = diffusive.turbulence.S
    normS = strain_rate_magnitude(S)
    k̂ = vertical_unit_vector(orientation, param_set, aux)

    # squared buoyancy correction
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(FT(1) - Richardson * _inv_Pr_turb, FT(0), FT(1)))
    ν₀ = normS * (m.C_smag * aux.turbulence.Δ)^2 + FT(1e-5)
    ν = SVector{3, FT}(ν₀, ν₀, ν₀)
    ν_v = k̂ .* dot(ν, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

# ### [Vreman Model](@id vreman)
# Vreman's turbulence model for anisotropic flows, which provides a
# less dissipative solution (specifically in the near-wall and transitional regions)
# than the Smagorinsky-Lilly method. This model
# relies of first derivatives of the velocity vector (i.e., the gradient tensor).
# By design, the Vreman model handles transitional as well as fully turbulent flows adequately.
# The input parameter to this model is the Smagorinsky coefficient - the coefficient is modified
# within the model functions to account for differences in model construction.
# #### Equations
# ```math
# \nu_{t} = 2.5 C_{s}^2 \sqrt{\frac{B_{\beta}}{u_{i,j}u_{i,j}}},
# ```
# where ($i,j, m = (1,2,3)$)
# ```math
# \begin{align}
# B_{\beta} &= \beta_{11}\beta_{22} + \beta_{11}\beta_{33} + \beta_{22}\beta_{33} - (\beta_{13}^2 + \beta_{12}^2 + \beta_{23}^2) \\
# \beta_{ij} &= \Delta_{m}^2 u_{i, m} u_{j, m} \\
# u_{i,j} &= \frac{\partial u_{i}}{\partial x_{j}}.
# \end{align}
# ```

"""
    Vreman{FT} <: TurbulenceClosureModel

Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Smagorinsky coefficient
is specified and used to compute the equivalent Vreman coefficient.

1) ν_e = √(Bᵦ/(αᵢⱼαᵢⱼ)) where αᵢⱼ = ∂uⱼ∂uᵢ with uᵢ the resolved scale velocity component.
2) βij = Δ²αₘᵢαₘⱼ
3) Bᵦ = β₁₁β₂₂ + β₂₂β₃₃ + β₁₁β₃₃ - β₁₂² - β₁₃² - β₂₃²
βᵢⱼ is symmetric, positive-definite.
If Δᵢ = Δ, then β = Δ²αᵀα

# Fields

$(DocStringExtensions.FIELDS)

# Reference

 - [Vreman2004](@cite)
"""
struct Vreman{FT} <: TurbulenceClosureModel
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end
vars_state(::Vreman, ::Auxiliary, FT) = @vars(Δ::FT)
vars_state(::Vreman, ::Gradient, FT) = @vars(θ_v::FT)
vars_state(::Vreman, ::GradientFlux, FT) =
    @vars(∇u::SMatrix{3, 3, FT, 9}, N²::FT)

function init_aux_turbulence!(
    ::Vreman,
    ::BalanceLaw,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end
function compute_gradient_argument!(
    m::Vreman,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end
function compute_gradient_flux!(
    ::Vreman,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.turbulence.∇u = ∇transform.u
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::Vreman,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    α = diffusive.turbulence.∇u
    S = symmetrize(α)
    k̂ = vertical_unit_vector(orientation, param_set, aux)

    normS = strain_rate_magnitude(S)
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * _inv_Pr_turb, 0, 1))

    β = (aux.turbulence.Δ)^2 * (α' * α)
    Bβ = principal_invariants(β)[2]

    ν₀ = m.C_smag^2 * FT(2.5) * sqrt(abs(Bβ / (norm2(α) + eps(FT))))

    ν = SVector{3, FT}(ν₀, ν₀, ν₀)
    ν_v = k̂ .* dot(ν, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

# ### [Anisotropic Minimum Dissipation](@id aniso-min-diss)
# This method is based Vreugdenhil and Taylor's minimum-dissipation eddy-viscosity model.
# The principles of the Rayleigh quotient minimizer are applied to the energy dissipation terms in the
# conservation equations, resulting in a maximum dissipation bound, and a model for
# eddy viscosity and eddy diffusivity.
# ```math
# \nu_e = (\mathrm{C}\delta)^2  \mathrm{max}\left[0, - \frac{\hat{\partial}_k \hat{u}_{i} \hat{\partial}_k \hat{u}_{j} \mathrm{\hat{S}}_{ij}}{\hat{\partial}_p \hat{u}_{q}} \right]
# ```
"""
    AnisoMinDiss{FT} <: TurbulenceClosureModel

Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Poincare coefficient
is specified and used to compute the equivalent AnisoMinDiss coefficient (computed as the solution to the
eigenvalue problem for the Laplacian operator).

# Fields
$(DocStringExtensions.FIELDS)

# Reference

See [Vreugdenhil2018](@cite)

"""
struct AnisoMinDiss{FT} <: TurbulenceClosureModel
    C_poincare::FT
end
vars_state(::AnisoMinDiss, ::Auxiliary, FT) = @vars(Δ::FT)
vars_state(::AnisoMinDiss, ::Gradient, FT) = @vars(θ_v::FT)
vars_state(::AnisoMinDiss, ::GradientFlux, FT) =
    @vars(∇u::SMatrix{3, 3, FT, 9}, N²::FT)
function init_aux_turbulence!(
    ::AnisoMinDiss,
    ::BalanceLaw,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end
function compute_gradient_argument!(
    m::AnisoMinDiss,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end
function compute_gradient_flux!(
    ::AnisoMinDiss,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.∇u = ∇transform.u
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end
function turbulence_tensors(
    m::AnisoMinDiss,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    k̂ = vertical_unit_vector(orientation, param_set, aux)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)

    ∇u = diffusive.turbulence.∇u
    S = symmetrize(∇u)
    normS = strain_rate_magnitude(S)

    δ = aux.turbulence.Δ
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * _inv_Pr_turb, 0, 1))

    δ_vec = SVector(δ, δ, δ)
    δ_m = δ_vec ./ transpose(δ_vec)
    ∇û = ∇u .* δ_m
    Ŝ = symmetrize(∇û)
    ν₀ =
        (m.C_poincare .* δ_vec) .^ 2 * max(
            FT(1e-5),
            -dot(transpose(∇û) * (∇û), Ŝ) / (dot(∇û, ∇û) .+ eps(normS)),
        )

    ν_v = k̂ .* dot(ν₀, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

"""
  NoHyperDiffusion <: HyperDiffusion
Defines a default hyperdiffusion model with zero hyperdiffusive fluxes.
"""
struct NoHyperDiffusion <: HyperDiffusion end

hyperviscosity_tensors(m::HyperDiffusion, bl::BalanceLaw, args...) =
    hyperviscosity_tensors(m, bl.orientation, bl.param_set, args...)

"""
  EquilMoistBiharmonic{FT} <: HyperDiffusion

Assumes equilibrium thermodynamics in compressible flow.
Horizontal hyperdiffusion methods for application in GCM and LES settings
Timescales are prescribed by the user while the diffusion coefficient is
computed as a function of the grid lengthscale.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct EquilMoistBiharmonic{FT} <: HyperDiffusion
    τ_timescale::FT
    τ_timescale_q_tot::FT
end

EquilMoistBiharmonic(τ_timescale::FT) where {FT} =
    EquilMoistBiharmonic(τ_timescale, τ_timescale)

vars_state(::EquilMoistBiharmonic, ::Auxiliary, FT) = @vars(Δ::FT)
vars_state(::EquilMoistBiharmonic, ::Gradient, FT) =
    @vars(u_h::SVector{3, FT}, h_tot::FT, q_tot::FT)
vars_state(::EquilMoistBiharmonic, ::GradientLaplacian, FT) =
    @vars(u_h::SVector{3, FT}, h_tot::FT, q_tot::FT)
vars_state(::EquilMoistBiharmonic, ::Hyperdiffusive, FT) = @vars(
    ν∇³u_h::SMatrix{3, 3, FT, 9},
    ν∇³h_tot::SVector{3, FT},
    ν∇³q_tot::SVector{3, FT}
)

function init_aux_hyperdiffusion!(
    ::EquilMoistBiharmonic,
    ::BalanceLaw,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.hyperdiffusion.Δ = lengthscale(geom)
end

function compute_gradient_argument!(
    h::EquilMoistBiharmonic,
    bl::BalanceLaw,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    u = state.ρu * ρinv
    k̂ = vertical_unit_vector(bl, aux)
    u_h = (SDiagonal(1, 1, 1) - k̂ * k̂') * u
    transform.hyperdiffusion.u_h = u_h
    transform.hyperdiffusion.h_tot = transform.h_tot
    transform.hyperdiffusion.q_tot = state.moisture.ρq_tot * ρinv
end

function transform_post_gradient_laplacian!(
    h::EquilMoistBiharmonic,
    bl::BalanceLaw,
    hyperdiffusive::Vars,
    hypertransform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    _inv_Pr_turb = eltype(state)(inv_Pr_turb(bl.param_set))
    ∇Δu_h = hypertransform.hyperdiffusion.u_h
    ∇Δh_tot = hypertransform.hyperdiffusion.h_tot
    ∇Δq_tot = hypertransform.hyperdiffusion.q_tot
    # Unpack
    τ_timescale = h.τ_timescale
    τ_timescale_q_tot = h.τ_timescale_q_tot
    # Compute hyperviscosity coefficient
    ν₄ = (aux.hyperdiffusion.Δ / 2)^4 / 2 / τ_timescale
    ν₄_q_tot = (aux.hyperdiffusion.Δ / 2)^4 / 2 / τ_timescale_q_tot
    hyperdiffusive.hyperdiffusion.ν∇³u_h = ν₄ * ∇Δu_h
    hyperdiffusive.hyperdiffusion.ν∇³h_tot = ν₄ * ∇Δh_tot
    hyperdiffusive.hyperdiffusion.ν∇³q_tot = ν₄_q_tot * ∇Δq_tot
end

function flux_second_order!(
    h::EquilMoistBiharmonic,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρu += state.ρ * hyperdiffusive.hyperdiffusion.ν∇³u_h
    flux.ρe += hyperdiffusive.hyperdiffusion.ν∇³u_h * state.ρu
    flux.ρe += hyperdiffusive.hyperdiffusion.ν∇³h_tot * state.ρ
    flux.moisture.ρq_tot += hyperdiffusive.hyperdiffusion.ν∇³q_tot * state.ρ
end

"""
  DryBiharmonic{FT} <: HyperDiffusion

Assumes dry compressible flow.
Horizontal hyperdiffusion methods for application in GCM and LES settings
Timescales are prescribed by the user while the diffusion coefficient is
computed as a function of the grid lengthscale.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct DryBiharmonic{FT} <: HyperDiffusion
    τ_timescale::FT
end
vars_state(::DryBiharmonic, ::Auxiliary, FT) = @vars(Δ::FT)
vars_state(::DryBiharmonic, ::Gradient, FT) =
    @vars(u_h::SVector{3, FT}, h_tot::FT)
vars_state(::DryBiharmonic, ::GradientLaplacian, FT) =
    @vars(u_h::SVector{3, FT}, h_tot::FT)
vars_state(::DryBiharmonic, ::Hyperdiffusive, FT) =
    @vars(ν∇³u_h::SMatrix{3, 3, FT, 9}, ν∇³h_tot::SVector{3, FT})

function init_aux_hyperdiffusion!(
    ::DryBiharmonic,
    ::BalanceLaw,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.hyperdiffusion.Δ = lengthscale(geom)
end

function compute_gradient_argument!(
    h::DryBiharmonic,
    bl::BalanceLaw,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    u = state.ρu * ρinv
    k̂ = vertical_unit_vector(bl, aux)
    u_h = (SDiagonal(1, 1, 1) - k̂ * k̂') * u
    transform.hyperdiffusion.u_h = u_h
    transform.hyperdiffusion.h_tot = transform.h_tot
end

function transform_post_gradient_laplacian!(
    h::DryBiharmonic,
    bl::BalanceLaw,
    hyperdiffusive::Vars,
    hypertransform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    _inv_Pr_turb = eltype(state)(inv_Pr_turb(bl.param_set))
    ∇Δu_h = hypertransform.hyperdiffusion.u_h
    ∇Δh_tot = hypertransform.hyperdiffusion.h_tot
    # Unpack
    τ_timescale = h.τ_timescale
    # Compute hyperviscosity coefficient
    ν₄ = (aux.hyperdiffusion.Δ / 2)^4 / 2 / τ_timescale
    hyperdiffusive.hyperdiffusion.ν∇³u_h = ν₄ * ∇Δu_h
    hyperdiffusive.hyperdiffusion.ν∇³h_tot = ν₄ * ∇Δh_tot
end

function flux_second_order!(
    h::DryBiharmonic,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρu += state.ρ * hyperdiffusive.hyperdiffusion.ν∇³u_h
    flux.ρe += hyperdiffusive.hyperdiffusion.ν∇³u_h * state.ρu
    flux.ρe += hyperdiffusive.hyperdiffusion.ν∇³h_tot * state.ρ
end

# ### [Viscous Sponge](@id viscous-sponge)
# `ViscousSponge` requires a user to specify a constant viscosity (kinematic),
# a sponge start height, the domain height, a sponge strength, and a sponge
# exponent.
# Given viscosity, diffusivity and stresses from arbitrary turbulence models,
# the viscous sponge enhances diffusive terms within a user-specified layer,
# typically used at the top of the domain to absorb waves. A smooth onset is
# ensured through a weight function that increases weight height from the sponge
# onset height.
# ```
"""
    NoViscousSponge
No modifiers applied to viscosity/diffusivity in sponge layer
# Fields
#
$(DocStringExtensions.FIELDS)
"""
struct NoViscousSponge <: ViscousSponge end
function sponge_viscosity_modifier(
    bl::BalanceLaw,
    m::NoViscousSponge,
    ν,
    D_t,
    τ,
    aux,
)
    return (ν, D_t, τ)
end

"""
    Upper domain viscous relaxation
Applies modifier to viscosity and diffusivity terms
in a user-specified upper domain sponge region
# Fields
#
$(DocStringExtensions.FIELDS)
"""
struct UpperAtmosSponge{FT} <: ViscousSponge
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Sponge exponent"
    γ::FT
end

function sponge_viscosity_modifier(
    bl::BalanceLaw,
    m::UpperAtmosSponge,
    ν,
    D_t,
    τ,
    aux::Vars,
)
    z = altitude(bl.orientation, bl.param_set, aux)
    if z >= m.sponge
        r = (z - m.z_sponge) / (m.z_max - m.z_sponge)
        β_sponge = m.α_max * sinpi(r / 2)^m.γ
        ν += β_sponge * ν
        D_t += β_sponge * D_t
        τ += β_sponge * τ
    end
    return (ν, D_t, τ)
end

end #module TurbulenceClosures.jl


text:
@testset "reverse" begin
    @testset "indices" begin
        inds = Indices([1,2,3])

        @test first(inds) == 1
        @test firstindex(inds) == 1
        @test last(inds) == 3
        @test lastindex(inds) == 3

        @test reverse(inds)::Indices == Indices([3, 2, 1])
        @test Iterators.reverse(inds)::ReverseIndices == Indices([3, 2, 1])

        @test reverse(convert(ArrayIndices, inds))::ArrayIndices == Indices([3, 2, 1])
        @test Iterators.reverse(convert(ArrayIndices, inds))::ReverseIndices == Indices([3, 2, 1])
    end

    @testset "dictionaries" begin
        d = Dictionary([:a, :b], [1, 2])
        rd = Dictionary([:b, :a], [2, 1])

        @test reverse(d)::Dictionary == rd
        @test Iterators.reverse(d)::ReverseDictionary == rd

        @test reverse(ArrayDictionary(d))::ArrayDictionary == rd
        @test Iterators.reverse(ArrayDictionary(d))::ReverseDictionary == rd

        @test first(d) == 1
        @test last(d) == 2
        @test firstindex(d) == :a
        @test lastindex(d) == :b
        #@test d[begin] == 1 # Parsing issues on earlier versions of Julia...
        @test d[end] == 2
    end
end


text:
## Faro ##

"""
    Faro{S}
    Weave{S}

Type with no fields representing the [Faro (weave)](https://en.wikipedia.org/wiki/Faro_shuffle)
card shuffle where `S` can be `:in` for an in-shuffle or `:out` for an out-shuffle.

See singleton instances [`infaro`](@ref), [`outfaro`](@ref), [`inweave`](@ref),
[`outweave`](@ref).

# Examples
```jldoctest
julia> shuffle([1, 2, 3, 4, 5, 6, 7, 8], Faro{:out}())
8-element Vector{Int64}:
 1
 5
 2
 6
 3
 7
 4
 8

julia> nshuffle(collect(1:52), 26, inweave) == collect(52:-1:1)
true
```
"""
struct Faro{S} <: AbstractDeterministicShuffle end

const Weave{S} = Faro{S}

"""
    infaro

The singleton instance of type [`Faro{:in}`](@ref)
"""
const infaro = Faro{:in}()

"""
    outfaro

The singleton instance of type [`Faro{:out}`](@ref)
"""
const outfaro = Faro{:out}()

"""
    inweave

Alias for [`infaro`](@ref).
"""
const inweave = Weave{:in}()

"""
    outweave

Alias for [`outfaro`](@ref).
"""
const outweave = Weave{:out}()

@inline _bo(::Array, i::Int, ::Faro{:in}) = 2i
@inline _bo(::Array, i::Int, ::Faro{:out}) = 2i - 1
@inline _bo(c::AbstractArray, i::Int, ::Faro{:in}) = 2i + firstindex(c) - 1
@inline _bo(c::AbstractArray, i::Int, ::Faro{:out}) = 2i + firstindex(c) - 2

@inline _eo(::Array, i::Int, ::Faro{:in}) = 2i - 1
@inline _eo(::Array, i::Int, ::Faro{:out}) = 2i
@inline _eo(c::AbstractArray, i::Int, ::Faro{:in}) = 2i + firstindex(c) - 2
@inline _eo(c::AbstractArray, i::Int, ::Faro{:out}) = 2i + firstindex(c) - 1

function shuffle!(c::AbstractArray, f::Faro)
    iseven(length(c)) || error("Faro (Weave) shuffling requires even length array")

    hlf = length(c) ÷ 2 + firstindex(c) - 1

    tmp::Array = c[firstindex(c):hlf]

    for i in 1:length(c) ÷ 2
        @inbounds begin
            c[_eo(c, i, f)] = c[hlf + i]
            c[_bo(c, i, f)] = tmp[i]
        end
    end

    return c
end

## Cut ##

"""
    Cut(n::Integer)

represents a [cut](https://en.wikipedia.org/wiki/Cut_(cards)) at `n`, meaning elements up to
and including `n` are moved to the bottom or end of a collection, while the rest shift up to
the beginning.

No in-place [`shuffle!`](@ref) or [`nshuffle!`](@ref) exist for this shuffling type.

# Examples

```jldoctest
julia> shuffle([1, 2, 3, 4, 5, 6, 7], Cut(3))
7-element Vector{Int64}:
 4
 5
 6
 7
 1
 2
 3
```
"""
struct Cut <: AbstractDeterministicShuffle
    n::Int
end

@inline _cuterr(A, n) = DomainError(
    n, "cannot cut array with indices $(eachindex(A)) at $n")
@inline _validcut(A, n) = firstindex(A) - 1 <= n <= lastindex(A)
@inline _checkcut(A, n) = _validcut(A, n) || throw(_cuterr(A, n))

shuffle(A::AbstractArray, c::Cut) = (_checkcut(A, c.n); circshift(A, -c.n))

nshuffle(A::AbstractArray, n::Integer, c::Cut) = (_checkcut(A, c.n); circshift(A, -c.n * n))


text:
export visualizeevents

function visualizeevents(Collection,Meta,Events)
    pygui(true)
    close("all")
    set_context("poster")
    set_context("poster")
    if length(Events)>100000
        nevents = length(Events)
        error("We can not visualize that many events ($nevents).")
    end

    nx,ny,x0,x1,y0,y1 = Meta.nx, Meta.ny, Meta.edgesx[1], Meta.edgesx[end], Meta.edgesy[1], Meta.edgesy[end]
    print(nx)
    fig, ax = plt.subplots(1,1,figsize=(6,6))
    cols = plt.cm.tab10(range(0,stop=1,length=1+nx*ny))
    cols[:,end] .= 0.3 # alpha
    edgesx = range(x0,stop=x1,length=nx+1)
    edgesy = range(y0,stop=y1,length=ny+1)
    for i=1:nx
        for j=1:ny
            #ax.axvline(edgesx[i],y0,y1,color="grey",lw=0.2,linestyle="--")
            #ax.axhline(edgesy[j],x0,x1,color="grey",lw=0.2,linestyle="--")
            for n=1:Collection[i,j].n
                x,y,r = Collection[i,j].x[n], Collection[i,j].y[n], Collection[i,j].r[n]
                ax.add_artist(plt.Circle((x,y), r,color=tuple([0,0,0]/255...)))#cols[j+(i-1)*ny,:]))
                xx = round(x,digits=2)
                yy = round(y,digits=2)
                #ax.text(x,y,"($i,$j,$n)",horizontalalignment="center")
            end
        end
    end

    ax.set_facecolor(tuple([0,150,130]/255...))

    #ax.axvline(edgesx[nx+1],y0,y1,color="grey",lw=1,linestyle="--")
    #ax.axhline(edgesy[ny+1],x0,x1,color="grey",lw=1,linestyle="--")
    ax.set_xlim([x0-2*Collection[1,1].r[1],x1+2*Collection[1,1].r[1]])
    ax.set_ylim([y0-2*Collection[1,1].r[1],y1+2*Collection[1,1].r[1]])

    #cols = plt.cm.Blues(range(0.2,stop=1,length=length(Events)))
    for i=1:length(Events)
        e = Events[i]
        if !e.crossinterfaceperiodic
            col = ifelse(e.collision,"red","blue")
            #col = cols[i,:]
            col = "white"
            col = ifelse(i==1,"white",col)
            ax.plot([e.oldx,e.newx],[e.oldy,e.newy],color=col,lw = 1.0,alpha=1)
        end
    end
    fig.tight_layout()
    ax.set_aspect("equal")
    ax.set_xlim([x0,x1])
    ax.set_ylim([y0,y1])
    ax.set_xticks([])
    ax.set_yticks([])
    plt.tight_layout()

    plt.savefig("Bounce.png",dpi=200)
    return fig, ax
end


text:
using Test, Random, StatsBase
using Zygote
using BetaVQE.VAN
using BetaVQE.VAN:bitarray, get_energy, free_energy, network

@testset "normalize" begin
    Random.seed!(2)
    nbits = 6
    model = PSAModel(nbits)
    configs = bitarray(collect(0:(1<<nbits-1)), nbits)
    logp = get_logp(model, configs)
    norm = sum(exp.(logp))

    @test isapprox(norm, 1.0, rtol=1e-5)
end

@testset "sample" begin
    Random.seed!(3)
    nbits = 4
    K = randn(nbits, nbits)
    K = (K+K')/2
    nsamples = 5000
    β = 1.0
    model = PSAModel(nbits)

    configs = bitarray(collect(0:(1<<nbits-1)), nbits)
    logp = get_logp(model, configs)
    f = sum(exp.(logp) .* (get_energy(K, configs) + get_logp(model, configs)))

    samples = gen_samples(model, nsamples)
    f_sample = free_energy(K, model, samples)
    @test isapprox(f, f_sample, rtol=1e-2)
end

@testset "gradscorefunction" begin
    Random.seed!(2)
    nbits = 6
    nbatchs = 10000
    model = PSAModel(nbits)
    samples = gen_samples(model, nbatchs)

    score(model, samples) = mean(get_logp(model, samples))
    grad = gradient(score, model, samples)[1]
    @show grad.w
    @test all(isapprox.(grad.w, 0.0, atol=1e-1))
end


text:
function load_genelist(genelist, description="entest")
    url = "http://maayanlab.cloud/Enrichr/addList"
    genestr = join(genelist, "\n")
    payload = Dict("list" => genestr, "description" => description)
    r = requests.post(url, files=payload)
    !r.ok && error("Error loading list")
    JSON.parse(r.text)
end

function view_genelist(config)
    url = string("http://maayanlab.cloud/Enrichr/view?userListId=", config["userListId"])
    r = requests.get(url)
    !r.ok && error("Error fetching results")
    JSON.parse(r.text)
end

function enrichment(config, geneset)
    url = "http://maayanlab.cloud/Enrichr/enrich"
    querystring = string("?userListId=", config["userListId"], "&backgroundType=", geneset)
    r = requests.get(string(url, querystring))
    !r.ok && error("Error fetching results")
    data = JSON.parse(r.text)

    df = DataFrame(GeneSet=String[], Rank=Int[], Term=String[], pvalue=Float64[], zscore=Float64[], Combined=Float64[], Genes=Vector{String}[], padj=Float64[])

    for record in data[geneset]
        push!(df, ([geneset ; record[1:end-2]]))
    end
    df
end


text:
### Export Data to HDF5

using HDF5

### Power System Data
function powesystem()
    reference = "https://labs.ece.uw.edu/pstca/pf14/pg_tca14bus.htm"
    grid = "transmission"
    name = "case14"

    basePower = 100

    bus = [1	3   0	    0      0	0	1	1.06	 0	    0	1	1.06	0.94
    	   2	2   21.7	12.7   0	0	1	1.045	-4.98	0	1	1.06	0.94
    	   3	2	94.2	19     0	0	1	1.01	-12.72	0	1	1.06	0.94
    	   4	1	47.8   -3.9    0	0	1	1.019	-10.33	0	1	1.06	0.94
    	   5	1	7.6	    1.6	   0	0	1	1.02	-8.78	0	1	1.06	0.94
    	   6	2	11.2	7.5    0	0	1	1.07	-14.22	0	1	1.06	0.94
    	   7	1	0	    0      0	0	1	1.062	-13.37	0	1	1.06	0.94
    	   8	2	0	    0      0	0	1	1.09	-13.36	0	1	1.06	0.94
    	   9	1	29.5	16.6   0	19	1	1.056	-14.94	0	1	1.06	0.94
    	  10	1	9	    5.8	   0	0	1	1.051	-15.1	0	1	1.06	0.94
    	  11	1	3.5	    1.8	   0	0	1	1.057	-14.79	0	1	1.06	0.94
    	  12	1	6.1	    1.6	   0	0	1	1.055	-15.07	0	1	1.06	0.94
    	  13	1	13.5	5.8	   0	0	1	1.05	-15.16	0	1	1.06	0.94
    	  14	1	14.9	5	   0	0	1	1.036	-16.04	0	1	1.06	0.94]

    branch = [1     1	   2	0.01938	 0.05917	0.0528	0	0	0	0	   0	1	-360	360
    	      2     1	   5	0.05403	 0.22304	0.0492	0	0	0	0	   0	1	-360	360
    	      3     2	   3	0.04699	 0.19797	0.0438	0	0	0	0	   0	1	-360	360
    	      4     2	   4	0.05811	 0.17632	0.034	0	0	0	0      0	1	-360	360
              5     2	   5	0.05695	 0.17388	0.0346	0	0	0	0	   0	1	-360	360
              6     3	   4	0.06701	 0.17103	0.0128	0	0	0	0	   0	1	-360	360
              7     4	   5	0.01335	 0.04211	0	    0	0	0	0	   0	1	-360	360
              8     4	   7	0	     0.20912	0	    0	0	0	0.978  0	1	-360	360
              9     4	   9	0	     0.55618	0	    0	0	0	0.969  0	1	-360	360
              10    5	   6	0	     0.25202	0	    0	0	0	0.932  0	1	-360	360
              11    6	  11	0.09498	 0.1989	    0	    0	0	0	0	   0	1	-360	360
              12    6	  12	0.12291	 0.25581	0	    0	0	0	0	   0	1	-360	360
              13    6	  13	0.06615	 0.13027	0	    0	0	0	0	   0	1	-360	360
              14    7	   8	0	     0.17615	0	    0	0	0	0	   0	1	-360	360
              15    7	   9	0	     0.11001	0	    0	0	0	0	   0	1	-360	360
              16    9	  10	0.03181	 0.0845	    0	    0	0	0	0	   0	1	-360	360
              17    9	  14	0.12711	 0.27038	0	    0	0	0	0	   0	1	-360	360
              18    10    11	0.08205	 0.19207	0	    0	0	0	0	   0	1	-360	360
              19    12    13	0.22092	 0.19988	0	    0	0	0	0	   0	1	-360	360
              20    13    14	0.17093	 0.34802	0	    0	0	0	0	   0	1	-360	360]

    generator = [1	232.4  -16.9	10	 0	 1.06	100	 1	332.4	0	0	0	0	0	0	0	0	0	0	0	0
              	 2	40	    42.4	50	-40	 1.045	100	 1	140	    0	0	0	0	0	0	0	0	0	0	0	0
              	 3	0	    23.4	40	 0	 1.01	100	 1	100	    0	0	0	0	0	0	0	0	0	0	0	0
              	 6	0	    12.2	24	-6	 1.07	100	 1	100	    0	0	0	0	0	0	0	0	0	0	0	0
              	 8	0	    17.4	24	-6	 1.09	100	 1 	100	    0	0	0	0	0	0	0	0	0	0	0	0]

    generatorcost = [2  0	0	3	0.0430292599	20	 0
             	     2	0	0	3	0.25	        20	 0
             	     2	0	0	3	0.01	        40	 0
             	     2	0	0	3	0.01	        40	 0
             	     2	0	0	3	0.01	        40	 0]

    return bus, branch, generator, generatorcost, basePower, reference, grid, name
end

### Save Data to HDF5
function savetoHDF5(ARGS...; reference = "", grid = "", path = "", name = "")
    busi = ["Bus" "Type" "Demand" "Demand" "Shunt Conductance" "Shunt Susceptance" "Area" "Voltage" "Voltage" "Base Voltage" "Loss Zone" "Maximum Voltage" "Minimum Voltage"]
    generatori = ["Bus" "Generation" "Generation" "Maximum Generation" "Minimum Generation" "Voltage" "Base" "Status" "Maximum Generation" "Minimum Generation" "Lower of PQ Curve" "Upper of PQ Curve" "Minimum at PC1" "Maximum at PC1" "Minimum at PC2" "Maximum at PC2" "Ramp Rate ACG" "Ramp Rate 10" "Ramp Rate 30" "Ramp Rate Q" "Area Factor"]
    branchi = ["Branch" "From Bus" "To Bus" "Series Parameter" "Series Parameter" "Charging Parameter" "Long Term Rate" "Short Term Rate" "Emergency Rate" "Transformer" "Transformer" "Status" "Minimum Voltage Difference" "Maximum Voltage Difference"]
    basei = ["Base Power"]

    busu = ["Integer" "PQ(1), PV(2), Slack(3)" "Active Power [MW]" "Reactive Power [MVAr]" "Active Power [MW]" "Reactive Power [MVAr]" "Integer" "Magnitude [p.u.]" "Angle [deg]" "Magnitude [kV]" "Integer" "Magnitude [p.u.]" "Magnitude [p.u.]"]
    generatoru = ["Integer" "Active Power [MW]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Magnitude [p.u.]" "Power [MVA]" "Integer" "Active Power [MW]" "Active Power [MW]" "Active Power [MW]" "Active Power [MW]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Reactive Power [MVAr]" "Active Power per Minut [MW/min]" "Active Power [MW]" "Active Power [MW]" "Reactive Power per Minut [MVAr/min]" "Integer"]
    branchu = ["Integer" "Integer" "Integer" "Resistance [p.u.]" "Reactance [p.u.]" "Susceptance [p.u.]" "Power [MVA]" "Power [MVA]" "Power [MVA]" "Turns Ratio" "Shift Angle [deg]" "Integer" "Angle [deg]" "Angle [deg]"]
    baseu = ["Power [MVA]"]

    generatorcosti = ["Cost Model" "Cost" "Cost" "Cost Model"]
    generatorcostu = ["Piecewise(1), Polynomial(2)" "Startup [currency]" "Shutdown [currency]" "Number of Data Points"]
    generatorcosta = ["Cost Model" for n = 1:generatorcost[1, 4]]
    generatorcostb = ["Coefficient c$(trunc(Int, (n-1)))" for n = generatorcost[1, 4]:-1:1]
    generatorcosti = [generatorcosti hcat(generatorcosta...)]
    generatorcostu = [generatorcostu hcat(generatorcostb...)]

    header1 = [busi, branchi, generatori, generatorcosti, basei]
    header2 = [busu, branchu, generatoru, generatorcostu, baseu]
    group = ["bus", "branch", "generator", "generatorcost", "basePower"]

    Nbus = size(bus, 1); Nbranch = size(branch, 1); Ngen = size(generator, 1)
    Ntrain = 0; Ntra = 0; pv = 0
    for i = 1:Nbranch
        if branch[i, 12] == 1 && (branch[i, 10] != 0 || branch[i, 11] != 0)
            Ntrain += 1
        end
        if branch[i, 10] != 0 || branch[i, 11] != 0
            Ntra += 1
        end
    end
    pv = length(unique(generator[:, 1]))

    info = [["Reference" string(reference) ""]
            ["Data" string(name) ""]
            ["Grid" string(grid) ""]
            ["" "" ""]
            ["Bus" string(Nbus) ""]
            ["PV bus" string(pv) ""]
            ["PQ bus" string(Nbus - pv - 1) ""]
            ["Shunt element" string(count(x->x != 0, abs.(bus[:, 5]) + abs.(bus[:, 6]))) ""]
            ["Generator" string(Ngen) "$(trunc(Int, sum(generator[:, 8]))) in-service"]
            ["Branch" string(Nbranch) "$(trunc(Int, sum(branch[:, 12]))) in-service"]
            ["Transformer" string(Ntra) "$Ntrain in-service"]]

    path = joinpath(path, string(name, ".h5"))
    h5open(path, "w") do file
        for i = 1:length(ARGS)
            write(file, group[i], ARGS[i])
            atr = Dict(string("row", k) => if !isempty(header2[i][k]) string(s, ": ", header2[i][k]) else s end for (k, s) in enumerate(header1[i]))
            h5writeattr(path, group[i], atr)
        end
        write(file, "info", info)
    end
end

### Run
pathtosave = "D:/Dropbox/"
bus, branch, generator, generatorcost, basePower, reference, grid, name = powesystem()
savetoHDF5(bus, branch, generator, generatorcost, basePower; reference = reference, grid = grid, name = name, path = pathtosave)

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text:
struct name
	field :: type_undefined
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Revtok
using Base.Test

# write your own tests here
@test begin
    text = replace(normalize_string(readstring(
        download("http://www.gutenberg.org/cache/epub/1661/pg1661.txt"))[4:end],
        newline2lf=true), r"\n+", s->length(s) == 1 ? " " : "\n" ^ (length(s) - 1));
    detokenize(tokenize(text)) == text
end