text:
module RI

include(joinpath(@__DIR__, "EiFunction.jl"))
import .EiFunction: expint

export lorentz_lorenz,
	   bruggeman,
       looyenga,
       maxwell,
       monecke,
       gedf,
       gem,
	   sellmeier,
	   cauchy_urbach,
	   drude_lorentz,
	   tauc_lorentz,
	   forouhi_bloomer,
	   lorentz_plasmon,
	   cody_lorentz,
       InfoF

function Info()
    tmp1 = "\n " *
	"\n " *
    "\n Available indices of refraction models:" *
	"\n " *
    "\n     bruggeman(p, n; df=1/3) -> EMA" *
    "\n     looyenga(p, n; df=1/3) -> EMA" *
    "\n     lorentz_lorenz(p, n) -> EMA" *
    "\n     maxwell(p, n; df=1/3) -> EMA" *
    "\n     monecke(p, n; dfm=:spheres, c=0.0) -> EMA" *
    "\n     gedf(x, n) -> EMA" *
    "\n     gem(x, n) -> EMA" *
    "\n     sellmeier(x, λ)" *
    "\n     cauchy_urbach(x, λ)" *
	"\n     drude_lorentz(x, ħω)" *
	"\n     tauc_lorentz(x, ħω)" *
	"\n     forouhi_bloomer(x, ħω)" *
	"\n     lorentz_plasmon(x, ħω)" *
	"\n     cody_lorentz(x, ħω; cme=:dipole)" *
	"\n " *
	"\n     To use any of these functions type: ?RI.function" *
	"\n " *
	"\n 	EMA: Effective medium approximation." *
	"\n "
	return println(tmp1)
end

## Function that checks the input for some EMA.
function _checkInputs_np(p, n)
	length(p) == size(n,2) || throw("The number of proportion values should be
    equal to the number of columns of n: (length(p) == size(n,2).")
    # sum(p) == 1.0 || throw("The sum over the p array should be equal to 1.")
end

function _checkInputs_ndf(n, df)
	length(df) == size(n,2) || throw("The length of depolarisation factors should be
    equal to the number of columns of n: (length(df) == size(n,2).")
end

"""

    Returns the effective index of refraction of a mixture of two liquids.

        neff = RI.lorentz_lorenz(p, n)

			p: proportion (real number, 0<p<1) of material n[:,1].
			n: column-array with the complex indices of refraction of the materials in the
 		       mixture. Each column represents the index of refraction of a material for a
 			   range of wavelengths.

            n: effective refractive index

		Source: https://en.wikipedia.org/wiki/Clausius–Mossotti_relation
				Thin Solid Films 519 (2011) 2994-2997

"""
function lorentz_lorenz(
    p::T0, n::Array{T1,2},
) where {T0<:Real, T1<:ComplexF64}
	ϵ = n.^2
	ϵeff = @. @views (2.0*p*ϵ[:,1] + ϵ[:,1]*ϵ[:,2] - 2.0*p*ϵ[:,2] + 2.0*ϵ[:,2])/(ϵ[:,1] + p*(ϵ[:,2] - ϵ[:,1]) + 2.0)
	return sqrt.(ϵeff)
end

"""

    Calculates the complex refractive index of a porous material using Bruggeman
    model, for nonmagnetic materials.

	Supports two components due to its complexity for a general solution,
    with different types of inclusions' geometries. To input the geometrical
    inclusions choosing the optional parameter, the depolarisation factor, df.
    By default takes df = 1/3.

        neff = RI.bruggeman(p, n; df=1/3)

			p: proportion parameter (real number, 0<p<1), with the proportion of material 1.
			n: column-array with the complex indices of refraction of the materials in the
			   mixture. The first column. n[:,1] has proportion p (inclusions), and n[:,2]
			   proportion 1-p (host material). Each column represents the index of refraction of
			   a material for a range of wavelengths.

               df: depolarisation factor (0.00001<=df<=0.99999, cannot be exactly zero or 1)
			   	  = 1/3 for spheres, default
				  = 1/2 for cylinders
				  = 0.99999 (close to one) for laminar structures with the fields directed parallel
				    to the direction along which layers alternate (stratification direction)
				  = 1/3<=df<=1 for prolate spheroidal shapes
                  = -1 for laminar structures for the fields directed perpendicularly
				   to the direction along which layers alternate (stratification direction)
            neff: effective index of refraction

    Source: Crystallography Reports, 2007, Vol. 52, No. 4, pp. 672-685
			Choy 2016, Effective medium theory: principles and applications
			Physical Review B, Vol 61, No 15, 2000-I, p 10437

"""
function bruggeman(
	p::T0, n::Array{T1,2};
	df=1.0/3.0,
) where {T0<:Real, T1<:ComplexF64}
    L = isequal(df, -1) ? :perpendicular : float(df)
    ϵeff = _bruggeman(float.([p, 1.0-p]), n.^2, L)
    return sqrt.(ϵeff)
end

function _bruggeman(
	p::Array{T0,1}, ϵ::Array{T1,2}, df::T0,
) where {T0<:Float64, T1<:ComplexF64}
	β = @. @views ϵ[:,1]*(p[1]/df - 1.0) + ϵ[:,2]*(p[2]/df - 1.0)
	ϵeff = @. @views 0.5*(β + sqrt(β^2 + 4*ϵ[:,1]*ϵ[:,2]*(1/df - 1.0)))/(1.0/df - 1.0)
	return ϵeff
end

function _bruggeman(
	p::Array{T0,1}, ϵ::Array{T1,2}, df::T2,
) where {T0<:Float64, T1<:ComplexF64, T2<:Symbol}
	ϵeff = @. @views ϵ[:,1]*p[1] + ϵ[:,2]*p[2]
	return ϵeff
end

"""

    Calculates the complex refractive index of the porous material using
    Looyenga-Landau-Lifshitz model, for nonmagnetic and isotropic materials.
    (See also Lichtenecker mixing formula: Stenzel, Optical characterization of
    thin solid films - 2018).

	Supports arbitrary number of components and the inclusions can be set
    through the polarisation factor df.

       neff = RI.looyenga(p, n; df=1.0/3.0)

	   	   p: array with proportion parameters (real number, 0<p<1). p[i] is the proportion of
		      material i respect to the total.
	   	   n: column-array with the complex indices of refraction of the materials in the
		      mixture. The first column. n[:,1] has proportion p[1], n[:,2] proportion p[2]
			  and so on. Each column represents the index of refraction of a material for a
			  range of wavelengths.
			  df: depolarization factors
   		   		 = 1/3 for spheres (default)
				 = 1/2 for cylinders
				 = 1 for laminar structures with the fields directed parallel to the direction
				   along which layers alternate (stratification direction)
				 = -1 for laminar structures with the fields directed perpendicular to the
				   direction along which layers alternate (stratification direction)
   			  The df entered as single values represent the same factors for all components, i.e.,
              there is only one geometrical type of inclusions for all of them.

           neff: effective index of refraction

    Sources: Physical Review B, Vol 61, No 15, 2000-I, p 10437
			 IEEE Transactions on Geoscience and Remote Sensing, Vol. 38, 2000.

"""
function looyenga(
	p::Array{T0,N0}, n::Array{T1,2};
	df=1.0/3.0,
) where {T0<:Real, N0, T1<:ComplexF64}
    p = float.(vec(p))
    _checkInputs_np(p, n)
    ϵeff = _looyenga(p, n.^2, float.(df))
    return sqrt.(ϵeff)
end

function _looyenga(
	p::Vector{T0}, ϵ::Array{T1,2}, df::T0,
) where {T0<:Float64, T1<:ComplexF64}
	ϵeff = zeros(ComplexF64, size(ϵ,1))
	for i in collect(1:size(ϵ,2))
		@. @views ϵeff += ϵ[:,i]^(1.0 - 2.0*df)*p[i]
	end
	@. ϵeff = ϵeff^(1.0/(1.0 - 2.0*df))
	return ϵeff
end

function _looyenga(
	p::Vector{T0}, ϵ::Array{T1,2}, df::Array{Array{T0,1},1},
) where {T0<:Float64, T1<:ComplexF64}
	_checkInputs_ndf(ϵ, df)
	ϵeff = zeros(ComplexF64, size(ϵ,1))
	for i in collect(1:size(ϵ,2)), j in collect(1:length(df))
        @. @views ϵeff += ϵ[:,i]^(1.0 - 2.0*df[i][j])*p[i]/length(df[i])
	end
	@. ϵeff = ϵeff^(1.0/(1.0 - 2.0*df))
	return ϵeff
end

"""

	Calculates the complex refractive index of a porous material using
    Maxwell-Garnett model, for nonmagnetic and isotropic materials.

		neff = RI.maxwell(p, n; df=1.0/3.0)

			p: array with proportion parameters (real number, 0<p<1). p[i] is the proportion of
			   material i respect to the total.
			n: column-array with the complex indices of refraction of the materials in the
			   mixture. The first column. n[:,1] has proportion p[1], and n[:,2] proportion p[2]
			   and so on. NOTICE that the host material is taken as the last column in the array
			   n. Each column represents the index of refraction of a material for a range of
			   wavelengths.
			   df: depolarization factors
				  = 1/3 for spheres (default)
				  = 1/2 for cylinders
				  = 1 for laminar structures with the fields directed parallel to the direction
					along which layers alternate (stratification direction)
				  = -1 for laminar structures with the fields directed perpendicular to the
					direction along which layers alternate (stratification direction)
			   The df factors are the same for all components, i.e., there is only one
			   geometrical type of inclusions for all of them.

			neff: effective index of refraction

	Source: IEEE Transactions on Geoscience and Remote Sensing, Vol. 38, 2000.
			Choy 2016, Effective medium theory: principles and applications
			Physical Review B, Vol 61, No 15, 2000-I, p 10437
			Crystallography Reports, 2007, Vol. 52, No. 4, pp. 672-685

"""
function maxwell(
	p::Array{T0,N0}, n::Array{T1,2};
	df::T0=1.0/3.0,
) where {T0<:Real, N0, T1<:ComplexF64}
	_checkInputs_np(p, n)
	ϵ = n.^2
	ϵeff = zeros(ComplexF64, size(n,1))
	for i in collect(1:size(n,2))
		@. @views ϵeff += p[i]*(ϵ[:,i] - ϵ[:,end])/(ϵ[:,end] + (ϵ[:,i] - ϵ[:,end])*df)
	end
	@. @views ϵeff = ϵ[:,end]*(ϵeff*(1.0 - df) + 1.0)/(1.0 - ϵeff*df)
    return sqrt.(ϵeff)
end

"""

    Calculates the complex refractive index of the porous material using
    Monecke model, for nonmagnetic and isotropic materials.

	Supports two components with different types of inclusions' geometries. To input the
	geometrical inclusions choosing the optional parameter dfm = :spheres or :cylinders.
	The last one are considered to be aligned infinite cylinders whose axes are perpendicular
	to the electric field. By default takes :spheres.

        neff = RI.monecke(p, n; dfm=:spheres, c=0.0)

			p: proportion parameter (real number, 0<p<1), with the proportion of n[:,1].
			n: column-array with the complex indices of refraction of the materials in the
			   mixture. The first column. n[:,1] has proportion p (inclusions), and n[:,2]
			   proportion 1-p (host material). Each column represents the index of refraction
			   of a material for a range of wavelengths
			   dfm: inclusions geometry
                   :spheres, default
                   :cylinders, aligned infinite cylinders whose axes are perpendicular to
				   the electric field
			   c: parameter to define the concentration of the host see source 2, to be used
			   	  with the :cylinders geometry, is = 0.0 by default

			neff: effective index of refraction

    Source: Physical Review B, Vol 61, No 15, 2000-I, p 10437
			Physical Review B, Vol 55, No 11, 1997-I, p 6739

"""
function monecke(
	p::T0, n::Array{T1,2};
	dfm::T2=:spheres, c::T0=0.0,
) where {T0<:Real, T1<:ComplexF64, T2<:Symbol}
    x = float.([p c])
    ϵ = n.^2
    if isequal(dfm, :spheres)
        ϵeff = _monecke_spheres(x, ϵ)
    elseif isequal(dfm, :cylinders)
        ϵeff = _monecke_cylinders(x, ϵ)
    else
        error("The accepted df values are: :spheres or :cylinders.")
    end
    return sqrt.(ϵeff)
end

function _monecke_spheres(
	x::Array{T0,2}, ϵ::Array{T1,2},
) where {T0<:Float64, T1<:ComplexF64}
	p = x[1]
	β = @. @views (ϵ[:,2] - ϵ[:,1])/(2.0*ϵ[:,1] + ϵ[:,2])
	ϵeff = @. @views (ϵ[:,1]*β + 4.0*(1.0 - p)*ϵ[:,1] + 2.0*(1.0 - p)^2*(ϵ[:,2] - ϵ[:,1]))/(β + 1.0 - p)
	return ϵeff
end

function _monecke_cylinders(
	x::Array{T0,2}, ϵ::Array{T1,2},
) where {T0<:Float64, T1<:ComplexF64}
	p, c = x # Physical Review B, Vol 55, No 11, 1997-I, p 6739
	s = @. @views ϵ[:,1]/(ϵ[:,1] - ϵ[:,2])
	s0 = @. 0.5 - (1.0 - p)*0.5 + 1.5*(1.0 - p)^2 - (1.0 - p)^3
	η = @. 0.5*s0*p*(1.0 - p)
	ϵeff = @. @views ϵ[:,1]*(1.0 - (c - η)/s - η/(s - s0))
	return ϵeff
end

"""

    Calculates the complex refractive index of a porous material using a general
    dielectric model, for nonmagnetic and isotropic materials. It considers
    n[:,1] hosted in another material n[:,2].

        neff = RI.gedf(x, n)

        	x: [p, β]
        		p: porosity parameter (real number, 0<p<1), proportion of n[:,1] over total
	            β: either 1 or 3
            n: column-array with the complex indices of refraction of the materials in the
			   mixture. The first column. n[:,1] has proportion p (inclusions), and n[:,2]
			   proportion 1-p (host material). Each column represents the index of refraction
			   of a material for a range of wavelengths.

            neff: effective refractive index

    Source: Appl. Phys. Lett. 108, 102902 (2016); doi: 10.1063/1.4943639.

"""
function gedf(x::AbstractArray{S}, n::Array{T,2}) where {S<:Float64, T<:ComplexF64}
	p, β = x
    ϵ = n.^2
	ϵeff = @. @views ϵ[:,2]*(ϵ[:,1] + (β*ϵ[:,2]) + (β*p*ϵ[:,1]) - (β*p*ϵ[:,2])) / (ϵ[:,1] + (β*ϵ[:,2]) - (p*ϵ[:,1]) + (p*ϵ[:,2]))
    return sqrt.(ϵeff)
end

"""

    Here we used equation (3) and the parameters ϕ, ϕc and tp as fitted. For a
    more detailed analysis and relation to the depolarization factors see the
    reference. It considers n[:,1] hosted in another material n[:,2].

        neff = RI.gem(x, n)

        	x: [ϕ, ϕc, tp]
        		ϕ: porosity parameter (real number, 0<ϕ<1), proportion of n1 over total
        		ϕc: critical volume fraction at which the high-index phase first percolates
        		tp: exponent related both to ϕc and to the shape of the grains
			n: column-array with the complex indices of refraction of the materials in the
			   mixture. The first column. n[:,1] has proportion p (inclusions), and n[:,2]
			   proportion 1-p (host material). Each column represents the index of refraction
			   of a material for a range of wavelengths.

            neff: effective refractive index

    Source: Carbon 40 (2002) 2801–2815.

"""
function gem(x::AbstractArray{S}, n::Array{T,2}) where {S<:Float64, T<:ComplexF64}
	ϕ, ϕc, tp = x
    α = (1-ϕc)/ϕc
    β = 1/tp
	ϵ = (n.^2).^β
	ϵeff = @. @views ((α*ϵ[:,1] - ϕ*ϵ[:,1] - α*ϕ*ϵ[:,1] - ϵ[:,2] + ϕ*ϵ[:,2] + α*ϕ*ϵ[:,2] + sqrt(4.0*α*ϵ[:,1]*ϵ[:,2] + ((α*(ϕ-1) + ϕ)*ϵ[:,1] + ϵ[:,2] - (1.0 + α)*ϕ*ϵ[:,2])^2))/(2.0*α))^tp
    return sqrt.(ϵeff)
end

"""

    Returns the index of refraction using the Sellmeier equation.

        n = RI.sellmeier(x, λ)

            x: Array containing the six values for the Sellmeier equation
               as described in https://en.wikipedia.org/wiki/Sellmeier_equation.
               (x = [B1, B2, B3, C1, C2, C3])
            λ: range of wavelengths in μm

            n: complex index of refraction

    Source: https://en.wikipedia.org/wiki/Sellmeier_equation

	Notice λ in μm.

"""
function sellmeier(x::AbstractArray{T}, λ::AbstractArray{T}) where {T<:Float64}
    λ = λ.^2
	n = @. sqrt(complex(1.0 + (x[1]*λ/(λ - x[4])) + (x[2]*λ/(λ - x[5])) + (x[3]*λ/(λ - x[6]))))
	return n
end

"""

    Returns the index of refraction using the Cauchy-Urbach model:
	n = Aₙ + Bₙ*1e4/λ^2 + Cₙ*1e9/λ^4 + im*α₀*exp((1240.0/λ - E₀)/Eᵤ).

        n = RI.cauchy_urbach(x, λ)

            x: Array containing the six values used in the model.
               (x = [
			   		 Aₙ, Bₙ, Cₙ, # parameters of the real part of the index of refraction
					 α₀, E₀, Eᵤ, # parameters of the imag part of the index of refraction
				]
			   )
            λ: range of wavelengths in nm

            n: complex index of refraction

    Source: E. Krous. Characterization of scandium oxide thin films for use in interference
	        coatings for high-power lasers operating in the near-infrared. Master’s thesis,
			Colorado State University, 2010. Page 60.

"""
function cauchy_urbach(x::AbstractArray{T}, λ::AbstractArray{T}) where {T<:Float64}
	Aₙ, Bₙ, Cₙ, α₀, E₀, Eᵤ = x
	# n = @. Aₙ + Bₙ*1e4/λ^2 + Cₙ*1e9/λ^4 + im*Aₖ*exp(Bₖ + (1240.0/λ/Cₖ))
	# n = @. Aₙ + Bₙ*1e4/λ^2 + Cₙ*1e9/λ^4 + im*Aₖ*exp(Bₖ*(1.0/λ - 1.0/Cₖ))
	n = @. Aₙ + Bₙ*1e4/λ^2 + Cₙ*1e9/λ^4 + im*α₀*exp((1240.0/λ - E₀)/Eᵤ)
	return n
end

"""

    Returns the index of refraction using the Drude-Lorentz oscillator model.

        n = RI.drude_lorentz(x, ħω)

            x: Array of array containing the values used in the model.
               (x = [
                      [ϵinf], # dielectric function at infinite frequency (offset)
					  [f1, ħωp1, ħω01, Γ1], # first oscillator
                      [f2, ħωp2, ħω02, Γ2], # second oscillator
                      [...],  # kth oscillator
                ])
            ħω: range of energies (1240/λ(nm))

            n: complex index of refraction

    Sources: F. Wooten, Optical properties of solids, Academic Press (1972)
			 http://www.horiba.com/fileadmin/uploads/Scientific/Downloads/OpticalSchool_CN/TN/ellipsometer/Lorentz_Dispersion_Model.pdf
             Stenzel, Optical characterization of thin solid films - 2018

"""
function drude_lorentz(
    x::Array{Array{T0,1},1}, ħω::AbstractArray{T1},
) where {T0<:Real, T1<:Real}
    ϵ = zeros(ComplexF64, length(vec(ħω)))
	ϵinf = x[1][1]
    for i = 2:length(x)
        f, ħωₚ, ħω₀, Γ = x[i]
        @. ϵ += (f*ħωₚ^2/(ħω₀^2 - ħω^2 - im*Γ*ħω))
    end
	@. ϵ += ϵinf
	return sqrt.(ϵ)
end

"""

    Returns the index of refraction using the Tauc-Lorentz dispersion model.

        n = RI.tauc_lorentz(x, ħω)

            x: Array of array containing the values used in the model.
               (x = [
			   		  [ϵinf, Eg],
                      [A1, E01, Γ1], # first oscillator
                      [A2, E02, Γ2], # second oscillator
                      [...],  # kth oscillator
                ])
               where A, E0 and Γ are in eV units.
            ħω: range of energies (1240/λ(nm))

            n: complex index of refraction

    Sources: Appl. Phys. Lett. 69 (3), 1996
    		 Appl. Phys. Lett. 69 (14), 1996
             J. Phys.: Condens. Matter 20 015216, 2008 (doi.org/10.1088/0953-8984/20/01/015216)

"""
function tauc_lorentz(
    x::Array{Array{T0,1},1}, ħω::AbstractArray{T1},
) where {T0<:Real, T1<:Real}
	ϵ = zeros(ComplexF64, length(vec(ħω)))
	ϵinf, Eg = x[1]
	for i = 2:length(x)
		A, E₀, C = x[i]
		ħω2 = ħω.^2
		E02 = E₀^2
		Eg2 = Eg^2
		C2 = C^2
		γ2 = sqrt(complex(E02 - C2/2))^2
		α = sqrt(complex(4.0*E02 - C2))
		ζ4 = @. (ħω2 - γ2)^2 + α^2*C2/4
		aL = @. (Eg2 - E02)*ħω2 + Eg2*C2 - E02*(E02 + 3.0*Eg2)
		aA = @. (ħω2 - E02)*(E02 + Eg2) + Eg2*C2
		# ϵf1 = @. A*C*aL/2.0/π/ζ4/α/E₀*log((E02 + Eg2 + α*Eg)/(E02 + Eg2 - α*Eg)) -
		# 		 A*aA/π/ζ4/E₀*(π - atan((2.0*Eg + α)/C) + atan((α - 2.0*Eg)/C)) +
		# 		 2.0*A*E₀*Eg/π/ζ4/α*(ħω2 - γ2)*(π + 2.0*atan(2.0/α/C*(γ2 - Eg2))) -
		# 		 A*E₀*C*(ħω2 + Eg2)/π/ζ4/ħω*log(abs(ħω - Eg)/(ħω + Eg)) +
		# 		 2.0*A*E₀*C*Eg/π/ζ4*log((abs(ħω - Eg)*(ħω + Eg))/sqrt((E02 - Eg2)^2 + Eg2*C2))
		ϵf1 = @. A*C*aL/2.0/π/ζ4/α/E₀*log(complex((E02 + Eg2 + α*Eg)/(E02 + Eg2 - α*Eg))) -
 				 A*aA/π/ζ4/E₀*(π - atan((2.0*Eg + α)/C) + atan((α - 2.0*Eg)/C)) +
 				 2.0*A*E₀*Eg/π/ζ4/α*(ħω2 - γ2)*(π + 2.0*atan(2.0/α/C*(γ2 - Eg2))) -
 				 A*E₀*C*(ħω2 + Eg2)/π/ζ4/ħω*log(abs(ħω - Eg)/(ħω + Eg)) +
 				 2.0*A*E₀*C*Eg/π/ζ4*log(abs(ħω - Eg)*(ħω + Eg)/sqrt(complex((E02 - Eg2)^2 + Eg2*C2)))
		ϵf2 = @. (A*E₀*C*(ħω - Eg)^2)/((ħω2 - E02)^2 + C2*ħω2)/ħω*(ħω > Eg)
		@. ϵ += ϵf1 + im*ϵf2
	end
	@. ϵ += ϵinf
	return sqrt.(ϵ)
end

"""

    Returns the index of refraction using the Forouhi-Bloomer model.
    Example marerials: Ta2O5, HfO2, Si3N4, Al2O3, a-Si, AlN.

        n = RI.forouhi_bloomer(x, ħω)

		x: Array of array containing the values used in the model.
		   (x = [
		   		  [ninf, Egap],
				  [A1, B1, C1], # first oscillator
				  [A2, B2, C2], # second oscillator
				  [...],  # kth oscillator
			])
            where B and Egap are in eV units, and C in eV^2.
            ħω: range of energies (1240/λ(nm))

            n: complex index of refraction

    Source: Physical Review B, 34, 7018 (1986)
			Physical Review B, 38, 1865 (1988)

"""
function forouhi_bloomer(
    x::Array{Array{T0,1},1}, ħω::AbstractArray{T1},
) where {T0<:Real, T1<:Real}
	N = zeros(ComplexF64, length(vec(ħω)))
	ninf, Egap = x[1]
	for i = 2:length(x)
		A, B, C = x[i]
		Q = real(0.5*sqrt(complex(4*C - B^2)))
		B₀ = A/Q*(-(0.5*B^2) + (Egap*B) - Egap^2 + C)
		C₀ = A/Q*(0.5*B*(Egap^2 + C) - (2.0*Egap*C))
		tmp1 = @. ħω^2 - B*ħω + C
		n = @. (B₀*ħω + C₀)/tmp1
		k = @. (A*(ħω - Egap)^2)/tmp1
		@. N += n + im*k
	end
	@. N += ninf
	return N
end

"""

    Returns the index of refraction using the Lorentz oscillator with plasmon damping.

        n = RI.lorentz_plasmon(x, ħω)

			x: Array of array containing the values used in the model.
		   		(x = [
						[ϵinf, ħωp, Γ],
				  		[ħωL1, ħωT1, γ1], # first oscillator
						[ħωL2, ħωT2, γ2], # second oscillator
						[...],  # kth oscillator
				])
               where all the parameters are in eV units, and ϵinf non-dimensional.
            ħω: range of energies (1240./λ(nm))

            n: complex index of refraction

    Sources: Physical Review B, 61, 10437 (2000-I)
			 J. Appl. Phys. 116, 233105 (2014)

"""
function lorentz_plasmon(
	x::Array{Array{T0,1},1}, ħω::AbstractArray{T1},
) where {T0<:Real, T1<:Real}
	ϵ = zeros(ComplexF64, length(vec(ħω)))
	ϵinf, ħωp, Γ = x[1]
	for i = 2:length(x)
		ħωL, ħωT, γ = x[i]
		@. ϵ += (ħω^2 - ħωL^2 + im*γ*ħω)/(ħω^2 - ħωT^2 + im*γ*ħω)
	end
	@. ϵ -= ħωp^2/(ħω^2 + im*Γ*ħω)
	@. ϵ *= ϵinf
	return sqrt.(ϵ)
end

"""

    Returns the index of refraction using the Cody-Lorentz dispersion model.

        n = RI.cody_lorentz(x, ħω; cme=:dipole)

            x: Array of array containing the values used in the model.
               (x = [
			   		  [ϵinf, Eg, Et, Ep, Eu],
                      [A1, E01, Γ1], # first oscillator
                      [A2, E02, Γ2], # second oscillator
                      [...],  # kth oscillator
                ])
               where all the parameters are in eV units.
            ħω: range of energies [eV] (1240./λ(nm))
			cme: There are two versions of Cody-Lorentz as defined by Ferlauto et al.
			 	depending on the constant matrix element (cme): :dipole (default) or
				:momentum. If you choose cme=:dipole, then the G function is calculated
				based on the Cody constant dipole matrix element. If you set cme=:momentum
				then G is calculated based on the Tauc theory with constant momentum matrix
				element.

            n: complex index of refraction

    Sources: Ferlauto et al. J. Appl. Phys., Vol. 92, No. 5, 2002
			 N. Malkova et al. Thin Solid Films 595 (2015) 32-35

	Notes:

	1. This function was not tested much since there is a scarce of biliography with data.

	2. This function uses the expint function (exponential integral function) from the
	   Bridge.jl package developed by @stevenjg (and modified by @mschauer) since it is not
	   yet included in the SpecialFunctions.jl module.
	   For more, check:
	   	https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb
	   	https://github.com/JuliaMath/SpecialFunctions.jl/issues/19

"""
# 1. The parameter Eu is calculated as described in Malkova et al., and the real part of the
# dielectric function is assumed to be the same as described in Ferlauto, considering that Eu is a contant.
function cody_lorentz(
    x::Array{Array{T0,1},1}, ħω::AbstractArray{T1};
	cme::T2=:dipole,
) where {T0<:Real, T1<:Real, T2<:Symbol}
	isequal(cme, :dipole) && return _cody_lorentz_dipole(float.(x), vec(float.(ħω)))
	isequal(cme, :momentum) && return _cody_lorentz_momentum(float.(x), vec(float.(ħω)))
end

# Lorentz oscillator function _L and variable band edge function _G
_G(ħω::T1, Eg::T1) where {T1<:Float64} = ((ħω - Eg)/ħω)^2
_G(ħω::T1, Eg::T1, Ep::T1) where {T1<:Float64} = (ħω - Eg)^2/((ħω - Eg)^2 + Ep^2)
_L(A::T1, E0::T1, Γ::T1, ħω::T1) where {T1<:Float64} = A*E0*Γ*ħω/((ħω^2 - E0^2)^2 + (Γ*ħω)^2)
_Ld(E0::T1, Γ::T1, ħω::T1) where {T1<:Float64} = (ħω^2 - E0^2)^2 + (Γ*ħω)^2

# Determination of IcL assuming constant dipole matrix element
function _cody_lorentz_dipole(x::Array{Array{T1,1},1}, ħω::Array{T1,1}) where {T1<:Float64}
	ϵinf, Eg, Et, Ep, Eu = x[1]
	G = _G.(ħω, Eg, Ep)
	ϵ = zeros(ComplexF64, length(vec(ħω)))
	Iu = similar(ϵ)
	for i = 2:length(x)
		A, E0, Γ = x[i]
		# Eg<=Et<=E0 || throw("The Cody-Lorentz model requires that Eg<=Et<=E0.")
		D = (E0^2 - Et^4)/((E0^2 - Et^2)^2 + (Γ*Et)^2) + Ep^2*Et/(Et - Eg)/((Et - Eg)^2 + Ep^2)
		Eu = 0.5*Et/D
		L = _L.(A, E0, Γ, ħω)
		E1 = Et*_G.(Et, Eg, Ep)*_L.(A, E0, Γ, Et)
		_urbach_tail!.(Iu, E1, ħω, Eu, Et)
		Icl = _matrix_element_dipole.(A, E0, Eg, Et, Eu, Ep, Γ, ħω, G)
		ϵf1 = @. real(Iu) + real(Icl)
		ϵf2 = @. G*L*(ħω>Et) + E1*exp((ħω - Et)/Eu)/ħω*(ħω<=Et)
		@. ϵf2 *= ϵf2>0.0
		@. ϵ += ϵf1 + im*ϵf2
	end
	@. ϵ += ϵinf
	return sqrt.(ϵ)
end

function _matrix_element_dipole(
	A::T1, E0::T1, Eg::T1, Et::T1, Eu::T1, Ep::T1, Γ::T1, ħω::T1, G::T1,
) where {T1<:Float64}
	I0c = (π/2.0 - atan((Et - Ep)/Ep))/Ep
	c0c = 0.5*ħω*G/_Ld(E0,Γ,ħω)
	d0c = - 0.5*ħω*(ħω + Eg)^2/((ħω + Eg)^2 + Ep^2)/_Ld(E0,Γ,ħω)
	Fsq = Ep^2 + Eg^2
	ζ = E0^2 - Γ^2/2.0
	χ = sqrt(complex(4.0*E0^2 - Γ^2))
	Ksq = 2.0*Fsq + 2.0*ζ - 4.0*Eg^2
	Y4 = E0^4 + Fsq*(Ksq - Fsq) - 4.0*Eg^2*Ksq
	b0c = Y4*Fsq*(_Ld(E0,Γ,ħω)*((c0c - d0c)/ħω + 2.0*Eg*Ksq/Y4*(c0c + d0c)) - 1.0)/
		  ((Ksq - Fsq)*Fsq*Y4 + E0^4*Y4 + 4.0*Eg^2*Fsq*Ksq^2)
	b1c = (2.0*Eg*Ksq*b0c - _Ld(E0,Γ,ħω)*(c0c + d0c))/Y4
	a3c = - (b1c + c0c + d0c)
	a2c = - (b0c + 2.0*Eg*b1c + ħω*(c0c - d0c))
	a1c = - (2.0*Eg*b0c - (Ksq - Fsq)*b1c + (ħω^2 - 2.0*ζ)*(c0c + d0c))
	a0c = 1.0 + (Ksq - Fsq)*b0c + 2.0*Eg*Ksq*b1c - ħω*(ħω^2 - 2.0*ζ)*(c0c - d0c)
	_Icl = 2.0*A*E0*Γ/π*(b1c*(Eg*I0c - sqrt(complex(log((Et - Eg)^2 + Ep^2)))) + b0c*I0c)
	I1t = 0.5/χ/Γ*(π - 2.0*atan(2.0*(Et^2 - ζ)/χ/Γ))
	I0at = 0.5/Γ*(π - atan((2.0*Et + χ)/Γ) + atan((-2.0*Et + χ)/Γ))
	I0bt = 0.25/χ*log((Et^2 + E0^2 + χ*Et)/(Et^2 + E0^2 - χ*Et))
	Itl = 2.0*A*E0*Γ/π*(a3c*(ζ*I1t - log(_Ld(E0,Γ,Et))^(0.25)) +
		  a2c*(I0at + I0bt) + a1c*I1t + a0c*((I0at - I0bt)/E0^2) -
		  c0c*log(abs(ħω - Et)) - d0c*log(ħω + Et))
	Icl = _Icl + Itl
	return Icl
end

# Determination of ItL assuming constant momentum matrix element
function _cody_lorentz_momentum(x::Array{Array{T1,1},1}, ħω::Array{T1,1}) where {T1<:Float64}
	ϵinf, Eg, Et, Ep, Eu = x[1]
	G = _G.(ħω, Eg)
	ϵ = zeros(ComplexF64, length(vec(ħω)))
	Iu = similar(ϵ)
	for i = 2:length(x)
		A, E0, Γ = x[i]
		# Eg<=Et<=E0 || throw("The Cody-Lorentz model requires that Eg<=Et<=E0.")
		D = (E0^2 - Et^4)/((E0^2 - Et^2)^2 + (Γ*Et)^2) + Ep^2*Et/(Et - Eg)/((Et - Eg)^2 + Ep^2)
		Eu = 0.5*Et/D
		L = _L.(A, E0, Γ, ħω)
		E1 = Et*_G.(Et, Eg)*_L.(A, E0, Γ, Et)
		_urbach_tail!.(Iu, E1, ħω, Eu, Et)
		Itl = _matrix_element_momentum.(A, E0, Eg, Et, Eu, Γ, ħω, G)
		ϵf1 = @. real(Iu) + real(Itl)
		ϵf2 = @. G*L*(ħω>Et) + E1*exp((ħω - Et)/Eu)/ħω*(ħω<=Et)
		@. ϵf2 *= ϵf2>0.0
		@. ϵ += ϵf1 + im*ϵf2
	end
	@. ϵ += ϵinf
	return sqrt.(ϵ)
end

function _matrix_element_momentum(
	A::T1, E0::T1, Eg::T1, Et::T1, Eu::T1, Γ::T1, ħω::T1, G::T1,
) where {T1<:Float64}
	ζ = E0^2 - Γ^2/2.0 # I define ζ as ζ squared by removing the sqrt
	χ = sqrt(complex(4.0*E0^2 - Γ^2))
	c0t = 0.5*ħω*G/_Ld(E0,Γ,ħω)
	d0t = - 0.5*(ħω + Eg)^2/ħω/_Ld(E0,Γ,ħω)
	a3t = - (c0t + d0t)
	a2t = - ħω*(c0t - d0t)
	a0t = 1.0 + (ħω^2 - 2.0*ζ)*a2t
	a1t = (ħω^2 - 2.0*ζ)*a3t
	I1t = 0.5/χ/Γ*(π - 2.0*atan(2.0*(Et^2 - ζ)/χ/Γ))
	I0at = 0.5/Γ*(π - atan((2.0*Et + χ)/Γ) + atan((-2.0*Et + χ)/Γ))
	I0bt = 0.25/χ*log((Et^2 + E0^2 + χ*Et)/(Et^2 + E0^2 - χ*Et))
	Itl = 2.0*A*E0*Γ/π*(a3t*(ζ*I1t - log(_Ld(E0,Γ,Et))^(0.25)) +
		  a2t*(I0at + I0bt) + a1t*I1t + a0t*((I0at - I0bt)/E0^2) -
		  c0t*log(abs(ħω - Et)) - d0t*log(ħω + Et))
	return Itl
end

function _urbach_tail!(
		Iu::T0, E1::T1, ħω::T1, Eu::T1, Et::T1,
) where {T0<:ComplexF64, T1<:Float64}
	Iu = E1/π/ħω*(
			exp((ħω - Et)/Eu)*(expint(complex((Et - ħω)/Eu)) - expint(complex(-ħω/Eu))) -
			exp(-(ħω + Et)/Eu)*(expint(complex((Et + ħω)/Eu)) - expint(complex(ħω/Eu)))
	)
	return nothing
end

end # module RI


text:
module Vectors3D
import Base.show, Base.convert, Base.abs, Base.dot
export Vector3D, getindex, phi
immutable Vector3D{T <: Real}
   x::T
   y::T
   z::T
end
+(u::Vector3D, v::Vector3D) = Vector3D(u.x + v.x, u.x + v.y, u.z + v.z)
-(u::Vector3D, v::Vector3D) = Vector3D(u.x - v.x, u.x - v.y, u.z - v.z)
*(a::Number, u::Vector3D) = Vector3D(a*u.x, a*u.y, a*u.z)
*(u::Vector3D, a::Number) = Vector3D(a*u.x, a*u.y, a*u.z)
/(a::Number, u::Vector3D) = Vector3D(u.x/a, u.y/a, u.z/a)
/(u::Vector3D, a::Number) = Vector3D(u.x/a, u.y/a, u.z/a)
show{T}(io::IO, u::Vector3D{T}) = print(io, "[$(u.x), $(u.y), $(u.z)]")

abs{T}(u::Vector3D{T}) = (u.x*u.x + u.y*u.y + u.z*u.z)^0.5
function phi(u::Vector3D, v::Vector3D)
  w0 = u.x*v.x + u.y*v.y + u.z*v.z
  w1 = abs(u) * abs(v)
  try
    return atan(w0 / w1)
  catch
    error("Division by zero")
  end
end
dot(u::Vector3D, v::Vector3D) = u.x*v.x + u.y*v.y + u.z+v.z
getindex(u::Vector3D, i) = (i == 1 ? u.x : (i==2 ? u.y : u.z))
convert(::Type{Vector3D}, u::Array{Float64,1}) =
                                    Vector3D(u[1],u[2],u[3])
end

ext: jl


text:
# ---
# title: 935. Knight Dialer
# id: problem935
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/knight-dialer/description/>
# hidden: true
# ---
#
# The chess knight has a **unique movement** , it may move two squares
# vertically and one square horizontally, or two squares horizontally and one
# square vertically (with both forming the shape of an **L** ). The possible
# movements of chess knight are shown in this diagaram:
#
# A chess knight can move as indicated in the chess diagram below:
#
# ![](https://assets.leetcode.com/uploads/2020/08/18/chess.jpg)
#
# We have a chess knight and a phone pad as shown below, the knight **can only
# stand on a numeric cell**  (i.e. blue cell).
#
# ![](https://assets.leetcode.com/uploads/2020/08/18/phone.jpg)
#
# Given an integer `n`, return how many distinct phone numbers of length `n` we
# can dial.
#
# You are allowed to place the knight **on any numeric cell** initially and then
# you should perform `n - 1` jumps to dial a number of length `n`. All jumps
# should be **valid** knight jumps.
#
# As the answer may be very large, **return the answer modulo** `109 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 1
#     Output: 10
#     Explanation: We need to dial a number of length 1, so placing the knight over any numeric cell of the 10 cells is sufficient.
#
#
# **Example 2:**
#
#
#
#     Input: n = 2
#     Output: 20
#     Explanation: All the valid number we can dial are [04, 06, 16, 18, 27, 29, 34, 38, 40, 43, 49, 60, 61, 67, 72, 76, 81, 83, 92, 94]
#
#
# **Example 3:**
#
#
#
#     Input: n = 3
#     Output: 46
#
#
# **Example 4:**
#
#
#
#     Input: n = 4
#     Output: 104
#
#
# **Example 5:**
#
#
#
#     Input: n = 3131
#     Output: 136006598
#     Explanation: Please take care of the mod.
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 5000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
if haskey(ENV, "CI") && ENV["CI"] == "true"
    ENV["PYTHON"] = ""
    using Pkg
    Pkg.build("PyCall")
end

using Test
using CUDA
using PyPlot
using LaTeXStrings

using Oceananigans
using Oceananigans.Advection

using ConvergenceTests
using ConvergenceTests.OneDimensionalGaussianAdvectionDiffusion: run_test
using ConvergenceTests.OneDimensionalUtils: unpack_errors, defaultcolors, removespines

""" Run advection test for all Nx in resolutions. """
function run_convergence_test(κ, U, resolutions, advection_scheme, arch)

    # Determine safe time-step
           Lx = 2.5
            h = Lx / maximum(resolutions)
           Δt = min(0.01 * h / U, 0.1 * h^2 / κ)
    stop_time = Δt

    # Run the tests
    results = [run_test(architecture=arch, Nx=Nx, Δt=Δt, advection=advection_scheme,
                        stop_iteration=1, U=U, κ=κ) for Nx in resolutions]

    return results
end

#####
##### Run tests
#####

arch = CUDA.has_cuda() ? GPU() : CPU()

advection_schemes = (CenteredSecondOrder(), CenteredFourthOrder(), UpwindBiasedThirdOrder(),
                     UpwindBiasedFifthOrder(), WENO5())

U = 1
κ = 1e-8
Nx = [8, 16, 32, 64, 96, 128, 192, 256, 384, 512]

results = Dict()
for scheme in advection_schemes
    t_scheme = typeof(scheme)
    results[t_scheme] = run_convergence_test(κ, U, Nx, scheme, arch)
end

rate_of_convergence(::CenteredSecondOrder) = 2
rate_of_convergence(::CenteredFourthOrder) = 4
rate_of_convergence(::UpwindBiasedThirdOrder) = 3
rate_of_convergence(::UpwindBiasedFifthOrder) = 5
rate_of_convergence(::WENO5) = 5
# rate_of_convergence(::WENO{K}) where K = 2K-1

test_resolution(::CenteredSecondOrder)    = 512
test_resolution(::CenteredFourthOrder)    = 512
test_resolution(::UpwindBiasedThirdOrder) = 512
test_resolution(::UpwindBiasedFifthOrder) = 512
test_resolution(::WENO5)                  = 512

tolerance(::CenteredSecondOrder)    = 0.02
tolerance(::CenteredFourthOrder)    = 0.06
tolerance(::UpwindBiasedThirdOrder) = 0.08
tolerance(::UpwindBiasedFifthOrder) = 0.2
tolerance(::WENO5)                  = 0.4

colors = ("xkcd:royal blue", "xkcd:light red")

for scheme in advection_schemes

    t_scheme = typeof(scheme)
    name = string(t_scheme)

    @testset "$name" begin

        fig, ax = subplots()

        roc = rate_of_convergence(scheme)
        atol = tolerance(scheme)
        Ntest = test_resolution(scheme)
        itest = searchsortedfirst(Nx, Ntest)

        (cx_L₁, cy_L₁, cz_L₁, uy_L₁, uz_L₁, vx_L₁, vz_L₁, wx_L₁, wy_L₁,
        cx_L∞, cy_L∞, cz_L∞, uy_L∞, uz_L∞, vx_L∞, vz_L∞, wx_L∞, wy_L∞) = unpack_errors(results[typeof(scheme)])

        test_rate_of_convergence(cx_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cx_L₁")
        test_rate_of_convergence(cy_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cy_L₁")
        test_rate_of_convergence(cz_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cz_L₁")

        test_rate_of_convergence(uy_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uy_L₁")
        test_rate_of_convergence(uz_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uz_L₁")

        test_rate_of_convergence(vx_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vx_L₁")
        test_rate_of_convergence(vz_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vz_L₁")

        test_rate_of_convergence(wx_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wx_L₁")
        test_rate_of_convergence(wy_L₁, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wy_L₁")

        test_rate_of_convergence(cx_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cx_L∞")
        test_rate_of_convergence(cy_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cy_L∞")
        test_rate_of_convergence(cz_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" cz_L∞")

        test_rate_of_convergence(uy_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uy_L∞")
        test_rate_of_convergence(uz_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" uz_L∞")

        test_rate_of_convergence(vx_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vx_L∞")
        test_rate_of_convergence(vz_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" vz_L∞")

        test_rate_of_convergence(wx_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wx_L∞")
        test_rate_of_convergence(wy_L∞, Nx, Ntest=Ntest, expected=-roc, atol=atol, name=name*" wy_L∞")

        @test cx_L₁ ≈ cy_L₁
        @test cx_L₁ ≈ cz_L₁
        @test uy_L₁ ≈ uz_L₁
        @test vx_L₁ ≈ vz_L₁
        @test wx_L₁ ≈ wy_L₁

        @test cx_L∞ ≈ cy_L∞
        @test cx_L∞ ≈ cz_L∞
        @test uy_L∞ ≈ uz_L∞
        @test vx_L∞ ≈ vz_L∞
        @test wx_L∞ ≈ wy_L∞

        common_kwargs = (linestyle="None", color=colors[1], mfc="None", alpha=0.8)

        loglog(Nx, cx_L₁; marker="*", label="\$L_1\$-norm, \$c(x)\$ $name", common_kwargs...)
        loglog(Nx, cy_L₁; marker="+", label="\$L_1\$-norm, \$c(y)\$ $name", common_kwargs...)
        loglog(Nx, cz_L₁; marker="_", label="\$L_1\$-norm, \$c(z)\$ $name", common_kwargs...)

        loglog(Nx, uy_L₁; marker="1", label="\$L_1\$-norm, \$u(y)\$ $name", common_kwargs...)
        loglog(Nx, uz_L₁; marker="^", label="\$L_1\$-norm, \$u(z)\$ $name", common_kwargs...)

        loglog(Nx, vx_L₁; marker="s", label="\$L_1\$-norm, \$v(x)\$ $name", common_kwargs...)
        loglog(Nx, vz_L₁; marker="v", label="\$L_1\$-norm, \$v(z)\$ $name", common_kwargs...)

        loglog(Nx, wx_L₁; marker="X", label="\$L_1\$-norm, \$w(x)\$ $name", common_kwargs...)
        loglog(Nx, wy_L₁; marker="D", label="\$L_1\$-norm, \$w(y)\$ $name", common_kwargs...)

        common_kwargs = (linestyle="None", color=colors[2], mfc="None", alpha=0.8)

        loglog(Nx, cx_L∞; marker="*", label="\$L_\\infty\$-norm, \$c(x)\$ $name", common_kwargs...)
        loglog(Nx, cy_L∞; marker="+", label="\$L_\\infty\$-norm, \$c(y)\$ $name", common_kwargs...)
        loglog(Nx, cz_L∞; marker="_", label="\$L_\\infty\$-norm, \$c(z)\$ $name", common_kwargs...)

        loglog(Nx, uy_L∞; marker="1", label="\$L_\\infty\$-norm, \$u(y)\$ $name", common_kwargs...)
        loglog(Nx, uz_L∞; marker="^", label="\$L_\\infty\$-norm, \$u(z)\$ $name", common_kwargs...)

        loglog(Nx, vx_L∞; marker="s", label="\$L_\\infty\$-norm, \$v(x)\$ $name", common_kwargs...)
        loglog(Nx, vz_L∞; marker="v", label="\$L_\\infty\$-norm, \$v(z)\$ $name", common_kwargs...)

        loglog(Nx, wx_L∞; marker="X", label="\$L_\\infty\$-norm, \$w(x)\$ $name", common_kwargs...)
        loglog(Nx, wy_L∞; marker="D", label="\$L_\\infty\$-norm, \$w(y)\$ $name", common_kwargs...)

        label = raw"\sim N_x^{-" * "$roc" * raw"}" |> latexstring

        loglog(Nx[itest-3:itest], uy_L₁[itest] .* (Nx[itest] ./ Nx[itest-3:itest]) .^ roc, color=colors[1], alpha=0.8, label=label)

        xscale("log", base=2)
        yscale("log", base=10)
        title("Convergence for $name advection scheme")
        xlabel(L"N_x")
        ylabel("\$L\$-norms of \$ | c_\\mathrm{sim} - c_\\mathrm{analytical} |\$")
        removespines("top", "right")
        lgd = legend(loc="upper right", bbox_to_anchor=(1.4, 1.0), prop=Dict(:size=>6))

        filename = "one_dimensional_convergence_$(name)_$(typeof(arch)).png"
        filepath = joinpath(@__DIR__, "figs", filename)
        mkpath(dirname(filepath))
        savefig(filepath, dpi=480, bbox_extra_artists=(lgd,), bbox_inches="tight")
    end
end


text:
using Zygote
import Zygote: accum

const dtol = 1e-12 # threshold value for precise recalculation of distances

function accum(D::Diagonal{T}, B::AbstractMatrix) where {T}
    A = Matrix{Union{T, Nothing}}(undef, size(D))
    A[diagind(A)] .= D.diag
    return accum(A, B)
end
accum(A::AbstractMatrix, D::Diagonal) = accum(D, A)
accum(A::Diagonal, B::Diagonal) = Diagonal(accum(diag(A), diag(B)))

#
# Some very specific broadcasting hacks while Zygote has crappy broadcasting.
#

import Base.Broadcast: broadcasted

function ChainRulesCore.rrule(::typeof(broadcasted), ::typeof(-), x::AbstractArray)
    function broadcasted_minus_pullback(Δ)
        return (NO_FIELDS, DoesNotExist(), .-Δ)
    end
    return .-x, broadcasted_minus_pullback
end

function ChainRulesCore.rrule(::typeof(broadcasted), ::typeof(exp), x::AbstractArray)
    y = exp.(x)
    return y, Δ->(NO_FIELDS, DoesNotExist(), Δ .* y)
end


text:
"""
    drop_missing!(table::CleanTable; missing_values::Vector=[])

Removes in-place from a `CleanTable` all rows where `missing` or `missing_values` have been
found.
"""
function drop_missing!(table::CleanTable; missing_values::Vector=[])
    columns = cols(table)
    nrows = length(columns[1])
    row_state = replace!(Vector{Bool}(undef, nrows), true => false)

    if Threads.nthreads() > 1 && length(columns) > 1 && length(columns[1]) >= 1_000_000
        Threads.@threads for col in columns
            i = 1

            while i <= nrows
                if !row_state[i]
                    if ismissing(col[i]) || in(col[i], missing_values)
                        row_state[i] = true
                    end
                end
                i += 1
            end
        end

        to_delete = findall(row_state)
        Threads.@threads for col in columns
            for num in 0:(length(to_delete) - 1)
                deleteat!(col, to_delete[num + 1] - num)
            end
        end
    else
        for col in columns
            i = 1

            while i <= nrows
                if !row_state[i]
                    if ismissing(col[i]) || in(col[i], missing_values)
                        row_state[i] = true
                    end
                end
                i += 1
            end
        end

        to_delete = findall(row_state)
        for col in columns
            for num in 0:(length(to_delete) - 1)
                deleteat!(col, to_delete[num + 1] - num)
            end
        end
    end

    return table
end

"""
    drop_missing(table; missing_values::Vector=[])

Creates a `CleanTable` with copied columns and removes from it all rows where `missing` or
`missing_values` have been found.
"""
function drop_missing(table; missing_values::Vector=[])
    return drop_missing!(CleanTable(table); missing_values=missing_values)
end

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
package java.security;

import java.security.Key;

public class Cipher {
    static int DECRYPT_MODE = 0;
    static int ENCRYPT_MODE = 1;
    static int PRIVATE_KEY = 2;
    static int PUBLIC_KEY = 3;
    static int SECRET_KEY = 4;
    static int UNWRAP_MODE = 5;
    static int WRAP_MODE = 6;

    private Cipher() {}

    public static Cipher getInstance(String s) {
	return new Cipher();
    }

    public void init(int mode, Key k) {}

    public byte[] doFinal(byte[] input) {
	return new byte[1];
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
function destroy(band::AbstractRasterBand)
    band.ptr = GDALRasterBand(C_NULL)
    return band
end

function destroy(band::IRasterBand)
    band.ptr = GDALRasterBand(C_NULL)
    band.ownedby = Dataset()
    return band
end

"""
    blocksize(band::AbstractRasterBand)

Fetch the "natural" block size of this band.

GDAL contains a concept of the natural block size of rasters so that
applications can organized data access efficiently for some file formats.
The natural block size is the block size that is most efficient for accessing
the format. For many formats this is simple a whole scanline in which case
`*pnXSize` is set to `GetXSize()`, and *pnYSize is set to 1.

However, for tiled images this will typically be the tile size.

Note that the X and Y block sizes don't have to divide the image size evenly,
meaning that right and bottom edge blocks may be incomplete. See `ReadBlock()`
for an example of code dealing with these issues.
"""
function blocksize(band::AbstractRasterBand)
    xy = Array{Cint}(undef, 2); x = pointer(xy); y = x + sizeof(Cint)
    GDAL.gdalgetblocksize(band.ptr, x, y)
    return xy
end

"""
    pixeltype(band::AbstractRasterBand)

Fetch the pixel data type for this band.
"""
pixeltype(band::AbstractRasterBand{T}) where T = T

"""
    width(band::AbstractRasterBand)

Fetch the width in pixels of this band.
"""
width(band::AbstractRasterBand) = GDAL.gdalgetrasterbandxsize(band.ptr)

"""
    height(band::AbstractRasterBand)

Fetch the height in pixels of this band.
"""
height(band::AbstractRasterBand) = GDAL.gdalgetrasterbandysize(band.ptr)

"""
    accessflag(band::AbstractRasterBand)

Return the access flag (e.g. `OF_ReadOnly` or `OF_Update`) for this band.
"""
accessflag(band::AbstractRasterBand) = GDAL.gdalgetrasteraccess(band.ptr)

"""
    indexof(band::AbstractRasterBand)

Fetch the band number (1+) within its dataset, or 0 if unknown.

This method may return a value of 0 to indicate overviews, or free-standing
`GDALRasterBand` objects without a relationship to a dataset.
"""
indexof(band::AbstractRasterBand) = GDAL.gdalgetbandnumber(band.ptr)

"""
    getdataset(band::AbstractRasterBand)

Fetch the handle to its dataset handle, or `NULL` if this cannot be determined.

Note that some `RasterBand`s are not considered to be a part of a dataset,
such as overviews or other "freestanding" bands.
"""
getdataset(band::AbstractRasterBand) =
    Dataset(GDAL.gdalgetbanddataset(band.ptr))
# ↑ GDAL wrapper checks null by default, but it is a valid result in this case

"""
    getunittype(band::AbstractRasterBand)

Return a name for the units of this raster's values. For instance, it might be
"m" for an elevation model in meters, or "ft" for feet.
"""
getunittype(band::AbstractRasterBand) = GDAL.gdalgetrasterunittype(band.ptr)

"""
    setunittype!(band::AbstractRasterBand, unitstring::AbstractString)

Set unit type of `band` to `unittype`.

Values should be one of \"\" (the default indicating it is unknown), \"m\"
indicating meters, or \"ft\" indicating feet, though other nonstandard values
are allowed.
"""
function setunittype!(band::AbstractRasterBand, unitstring::AbstractString)
    result = GDAL.gdalsetrasterunittype(band.ptr, unitstring)
    @cplerr result "Failed to set unit type"
    return band
end

"""
    getoffset(band::AbstractRasterBand)

Fetch the raster value offset.

This (in combination with `GetScale()`) is used to transform raw pixel values
into the units returned by `GetUnits()`. For e.g. this might be used to store
elevations in `GUInt16` bands with a precision of 0.1, starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically, a value of 0 is returned.
"""
getoffset(band::AbstractRasterBand) = GDAL.gdalgetrasteroffset(band.ptr, C_NULL)

"""
    setoffset!(band::AbstractRasterBand, value::Real)

Set scaling offset.
"""
function setoffset!(band::AbstractRasterBand, value::Real)
    result = GDAL.gdalsetrasteroffset(band.ptr, value)
    @cplerr result "Failed to set scaling offset."
    return band
end

"""
    getscale(band::AbstractRasterBand)

Fetch the raster value scale.

This value (in combination with the `GetOffset()` value) is used to transform
raw pixel values into the units returned by `GetUnits()`. For example this
might be used to store elevations in GUInt16 bands with a precision of 0.1,
and starting from -100.

    Units value = (raw value * scale) + offset

For file formats that don't know this intrinsically a value of one is returned.
"""
getscale(band::AbstractRasterBand) = GDAL.gdalgetrasterscale(band.ptr, C_NULL)

"""
    setscale!(band::AbstractRasterBand, ratio::Real)

Set scaling ratio.
"""
function setscale!(band::AbstractRasterBand, ratio::Real)
    result = GDAL.gdalsetrasterscale(band.ptr, ratio)
    @cplerr result "Failed to set scaling ratio"
    return band
end

"""
    getnodatavalue(band::AbstractRasterBand)

Fetch the no data value for this band.

If there is no out of data value, an out of range value will generally be
returned. The no data value for a band is generally a special marker value
used to mark pixels that are not valid data. Such pixels should generally
not be displayed, nor contribute to analysis operations.

### Returns
the nodata value for this band.
"""
function getnodatavalue(band::AbstractRasterBand)
    # ### Parameters
    # * `pbSuccess`   pointer to a boolean to use to indicate if a value is
    #     actually associated with this layer. May be `NULL` (default).
    return GDAL.gdalgetrasternodatavalue(band.ptr, C_NULL)
end

"""
    setnodatavalue!(band::AbstractRasterBand, value::Real)

Set the no data value for this band.
"""
function setnodatavalue!(band::AbstractRasterBand, value::Real)
    result = GDAL.gdalsetrasternodatavalue(band.ptr, value)
    @cplerr result "Could not set nodatavalue"
    return band
end

function deletenodatavalue!(band::AbstractRasterBand)
    result = GDAL.gdaldeleterasternodatavalue(band.ptr)
    @cplerr result "Could not delete nodatavalue"
    return band
end

"""
    setcategorynames!(band::AbstractRasterBand, names)

Set the category names for this band.
"""
function setcategorynames!(band::AbstractRasterBand, names)
    result = GDAL.gdalsetrastercategorynames(band.ptr, names)
    @cplerr result "Failed to set category names"
    return band
end

"""
    minimum(band::AbstractRasterBand)

Fetch the minimum value for this band.
"""
minimum(band::AbstractRasterBand) = GDAL.gdalgetrasterminimum(band.ptr, C_NULL)

"""
    maximum(band::AbstractRasterBand)

Fetch the maximum value for this band.
"""
maximum(band::AbstractRasterBand) = GDAL.gdalgetrastermaximum(band.ptr, C_NULL)

"""
    getdefaultRAT(band::AbstractRasterBand)

Fetch default Raster Attribute Table.
"""
getdefaultRAT(band::AbstractRasterBand) = GDAL.gdalgetdefaultrat(band.ptr)

"""
    setdefaultRAT!(band::AbstractRasterBand, rat::RasterAttrTable)

Set default Raster Attribute Table.

Associates a default RAT with the band. If not implemented for the format a
CPLE_NotSupported error will be issued. If successful a copy of the RAT is made,
the original remains owned by the caller.
"""
function setdefaultRAT!(band::AbstractRasterBand, rat::RasterAttrTable)
    result = GDAL.gdalsetdefaultrat(band.ptr, rat.ptr)
    @cplerr result "Failed to set default raster attribute table"
    return band
end

"""
    copywholeraster!( source::AbstractRasterBand, dest::AbstractRasterBand;
        [options, [progressdata, [progressfunc]]])

Copy all raster band raster data.

This function copies the complete raster contents of one band to another
similarly configured band. The source and destination bands must have the same
width and height. The bands do not have to have the same data type.

It implements efficient copying, in particular "chunking" the copy in
substantial blocks.

Currently the only `options` value supported is : "COMPRESSED=YES" to
force alignment on target dataset block sizes to achieve best compression.
More options may be supported in the future.

### Parameters
* `source`        the source band
* `dest`          the destination band
* `options`       transfer hints in "StringList" Name=Value format.
* `progressfunc`  progress reporting function.
* `progressdata`  callback data for progress function.
"""
function copywholeraster!(
        source::AbstractRasterBand,
        dest::AbstractRasterBand;
        options                 = StringList(C_NULL),
        progressdata            = C_NULL,
        progressfunc::Function  = GDAL.gdaldummyprogress
    )
    result = GDAL.gdalrasterbandcopywholeraster(source.ptr, dest.ptr, options,
        @cplprogress(progressfunc), progressdata)
    @cplerr result "Failed to copy whole raster"
    return source
end

"""
    noverview(band::AbstractRasterBand)

Return the number of overview layers available, zero if none.
"""
noverview(band::AbstractRasterBand) = GDAL.gdalgetoverviewcount(band.ptr)

"""
    getoverview(band::IRasterBand, i::Integer)

Fetch overview raster band object.
"""
getoverview(band::IRasterBand, i::Integer) =
    IRasterBand(GDAL.gdalgetoverview(band.ptr, i), ownedby = band.ownedby)

unsafe_getoverview(band::AbstractRasterBand, i::Integer) =
    RasterBand(GDAL.gdalgetoverview(band.ptr, i))

"""
    sampleoverview(band::IRasterBand, nsamples::Integer)

Fetch best overview satisfying `nsamples` number of samples.

Returns the most reduced overview of the given band that still satisfies the
desired number of samples `nsamples`. This function can be used with zero as the
number of desired samples to fetch the most reduced overview. The same band as
was passed in will be returned if it has not overviews, or if none of the
overviews have enough samples.
"""
function sampleoverview(band::IRasterBand, nsamples::Integer)
    return IRasterBand(
        GDAL.gdalgetrastersampleoverviewex(band.ptr, UInt64(nsamples)),
        ownedby = band.ownedby
    )
end

unsafe_sampleoverview(band::AbstractRasterBand, nsamples::Integer) =
    RasterBand(GDAL.gdalgetrastersampleoverviewex(band.ptr, UInt64(nsamples)))

"""
    getcolorinterp(band::AbstractRasterBand)

Color Interpretation value for band
"""
getcolorinterp(band::AbstractRasterBand) =
    GDAL.gdalgetrastercolorinterpretation(band.ptr)

"""
    setcolorinterp!(band::AbstractRasterBand, color::GDALColorInterp)

Set color interpretation of a band.
"""
function setcolorinterp!(band::AbstractRasterBand, color::GDALColorInterp)
    result = GDAL.gdalsetrastercolorinterpretation(band.ptr, color)
    @cplerr result "Failed to set color interpretation"
    return band
end

"""
    unsafe_getcolortable(band::AbstractRasterBand)

Returns a clone of the color table associated with the band.

(If there is no associated color table, the original result is `NULL`. The
original color table remains owned by the `GDALRasterBand`, and can't be
depended on for long, nor should it ever be modified by the caller.)
"""
function unsafe_getcolortable(band::AbstractRasterBand)
    result = ColorTable(GDALColorTable(GDAL.gdalgetrastercolortable(band.ptr)))
    if result.ptr == C_NULL
        return result
    else
        return unsafe_clone(result)
    end
end

"""
    setcolortable!(band::AbstractRasterBand, colortable::ColorTable)

Set the raster color table.

The driver will make a copy of all desired data in the colortable. It remains
owned by the caller after the call.

### Parameters
* `colortable` color table to apply (where supported).
"""
function setcolortable!(band::AbstractRasterBand, colortable::ColorTable)
    result = GDAL.gdalsetrastercolortable(band.ptr, colortable.ptr)
    @cplwarn result "CPLError $(result): action is unsupported by the driver"
    return band
end

function clearcolortable!(band::AbstractRasterBand)
    result = GDAL.gdalsetrastercolortable(band.ptr, GDALColorTable(C_NULL))
    @cplwarn result "CPLError $(result): action is unsupported by the driver"
    return band
end

"""
    regenerateoverviews!(band::AbstractRasterBand, overviewbands::Vector{<:AbstractRasterBand},
        resampling = "NEAREST")

Generate downsampled overviews.

This function will generate one or more overview images from a base image using
the requested downsampling algorithm. Its primary use is for generating
overviews via BuildOverviews(), but it can also be used to generate downsampled
images in one file from another outside the overview architecture.

### Parameters
* `band`              the source (base level) band.
* `overviewbands`   the list of downsampled bands to be generated.

### Keyword Arguments
* `resampling`      (optional) Resampling algorithm (eg. "AVERAGE"). default to
                    "NEAREST".
* `progressfunc`    (optional) progress report function.
* `progressdata`    (optional) progress function callback data.

### Additional Remarks
The output bands need to exist in advance.

This function will honour properly `NODATA_VALUES` tuples (special dataset
metadata) so that only a given RGB triplet (in case of a RGB image) will be
considered as the nodata value and not each value of the triplet independantly
per band.
"""
function regenerateoverviews!(
        band::AbstractRasterBand,
        overviewbands::Vector{<:AbstractRasterBand},
        resampling::AbstractString  = "NEAREST",
        # progressfunc::Function      = GDAL.gdaldummyprogress,
        progressdata                = C_NULL
    )
    cfunc = @cfunction(GDAL.gdaldummyprogress, Cint,
        (Cdouble, Cstring, Ptr{Cvoid}))
    result = GDAL.gdalregenerateoverviews(
        band.ptr,
        length(overviewbands),
        GDALRasterBand[band.ptr for band in overviewbands],
        resampling,
        cfunc,
        progressdata
    )
    @cplerr result "Failed to regenerate overviews"
    return band
end

# "Advise driver of upcoming read requests."
# _rasteradviseread(hRB::GDALRasterBandH,
#                   nDSXOff::Integer,
#                   nDSYOff::Integer,
#                   nDSXSize::Integer,
#                   nDSYSize::Integer,
#                   nBXSize::Integer,
#                   nBYSize::Integer,
#                   eBDataType::GDALDataType,
#                   papszOptions::Ptr{Cstring}) =
#     GDALRasterAdviseRead(hRB, nDSXOff, nDSYOff, nDSXSize, nDSYSize, nBXSize,
#                          nBYSize, eBDataType, papszOptions)::CPLErr

"""
    getcategorynames(band::AbstractRasterBand)

Fetch the list of category names for this raster.

The return list is a "StringList" in the sense of the CPL functions. That is a
NULL terminated array of strings. Raster values without associated names will
have an empty string in the returned list. The first entry in the list is for
raster values of zero, and so on.
"""
getcategorynames(band::AbstractRasterBand) =
    GDAL.gdalgetrastercategorynames(band.ptr)

"""
    setcategorynames!(band::AbstractRasterBand, names::Vector{String})

Set the category names for this band.
"""
function setcategorynames!(band::AbstractRasterBand, names::Vector{String})
    result = GDAL.gdalsetrastercategorynames(band.ptr, names)
    @cplerr result "Failed to set category names for this band"
    return band
end

# """
# Flush raster data cache.

# This call will recover memory used to cache data blocks for this raster band,
# and ensure that new requests are referred to the underlying driver.
# """
# function flushcache!(band::AbstractRasterBand)
#     result = GDAL.flushrastercache(band.ptr)
#     @cplerr result "Failed to flush raster data cache"
#     result
# end

"""
    fillraster!(band::AbstractRasterBand, realvalue::Real, imagvalue::Real = 0)

Fill this band with a constant value.

GDAL makes no guarantees about what values pixels in newly created files are
set to, so this method can be used to clear a band to a specified "default"
value. The fill value is passed in as a double but this will be converted to
the underlying type before writing to the file. An optional second argument
allows the imaginary component of a complex constant value to be specified.

### Parameters
* `realvalue`: Real component of fill value
* `imagvalue`: Imaginary component of fill value, defaults to zero
"""
function fillraster!(
        band::AbstractRasterBand,
        realvalue::Real,
        imagvalue::Real = 0
    )
    result = GDAL.gdalfillraster(band.ptr, realvalue, imagvalue)
    @cplerr result "Failed to fill raster band"
    return band
end

"""
    getmaskband(band::IRasterBand)

Return the mask band associated with the band.

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
GDALNoDataValuesMaskBand class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will
return `GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags
`GMF_PER_DATASET` and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new
`GDALAllValidRasterBand` class will be returned that has 255 values for all
pixels. The null flags will return `GMF_ALL_VALID`.

Note that the `GetMaskBand()` should always return a `GDALRasterBand` mask,
even if it is only an all 255 mask with the flags indicating `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
getmaskband(band::IRasterBand) =
    IRasterBand(GDAL.gdalgetmaskband(band.ptr), ownedby = band.ownedby)

unsafe_getmaskband(band::AbstractRasterBand) =
    RasterBand(GDAL.gdalgetmaskband(band.ptr))

"""
    maskflags(band::AbstractRasterBand)

Return the status flags of the mask band associated with the band.

The GetMaskFlags() method returns an bitwise OR-ed set of status flags with the
following available definitions that may be extended in the future:

* `GMF_ALL_VALID` (`0x01`):    There are no invalid pixels, all mask values
will be 255. When used this will normally be the only flag set.
* `GMF_PER_DATASET` (`0x02`):  The mask band is shared between all bands on
the dataset.
- `GMF_ALPHA` (`0x04`):        The mask band is actually an alpha band and may
have values other than 0 and 255.
- `GMF_NODATA` (`0x08`):       Indicates the mask is actually being generated
from nodata values. (mutually exclusive of `GMF_ALPHA`)

The `GDALRasterBand` class includes a default implementation of `GetMaskBand()`
that returns one of four default implementations:

- If a corresponding .msk file exists it will be used for the mask band.
- If the dataset has a `NODATA_VALUES` metadata item, an instance of the new
`GDALNoDataValuesMaskBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA | GMF_PER_DATASET`.
- If the band has a nodata value set, an instance of the new
`GDALNodataMaskRasterBand` class will be returned. `GetMaskFlags()` will return
`GMF_NODATA`.
- If there is no nodata value, but the dataset has an alpha band that seems to
apply to this band (specific rules yet to be determined) and that is of type
`GDT_Byte` then that alpha band will be returned, and the flags
`GMF_PER_DATASET` and `GMF_ALPHA` will be returned in the flags.
- If neither of the above apply, an instance of the new `GDALAllValidRasterBand`
class will be returned that has 255 values for all pixels. The null flags will
return `GMF_ALL_VALID`.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask

### Returns
a valid mask band.
"""
maskflags(band::AbstractRasterBand) = GDAL.gdalgetmaskflags(band.ptr)

"""
    createmaskband!(band::AbstractRasterBand, nflags::Integer)

Adds a mask band to the current band.

The default implementation of the `CreateMaskBand()` method is implemented
based on similar rules to the `.ovr` handling implemented using the
`GDALDefaultOverviews` object. A `TIFF` file with the extension `.msk` will be
created with the same basename as the original file, and it will have as many
bands as the original image (or just one for `GMF_PER_DATASET`). The mask
images will be deflate compressed tiled images with the same block size as the
original image if possible.

If you got a mask band with a previous call to `GetMaskBand()`, it might be
invalidated by `CreateMaskBand()`. So you have to call `GetMaskBand()` again.

See also: http://trac.osgeo.org/gdal/wiki/rfc15_nodatabitmask
"""
function createmaskband!(band::AbstractRasterBand, nflags::Integer)
    result = GDAL.gdalcreatemaskband(band.ptr, nflags)
    @cplerr result "Failed to create mask band"
    return band
end


text:
        # Exponential cone constrained problems

using Test, LinearAlgebra, SparseArrays
using COSMO

@testset "Exponential cone problems" begin

    @testset "Feasible" begin

        # solve the following exponential cone problem
        # max  x
        # s.t. y * exp(x / y) <= z
        #      y == 1, z == exp(5)

        P = spzeros(3, 3)
        q = [-1.; 0; 0]

        # y * exp( x / y) <= z
        A1 = sparse(Matrix(1.0I, 3, 3))
        b1 = zeros(3)

        # y == 1 and z == exp(5) = 148.4131591025766
        A2 = sparse([0 1. 0; 0 0 1])
        b2 = [-1.; -exp(5)]
        cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)
        cs2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)

        model = COSMO.Model()
        assemble!(model, P, q, [cs1; cs2])

        res = COSMO.optimize!(model)
        @test res.status == :Solved
        @test isapprox(res.obj_val, -5, atol=1e-3)
    end

    @testset "Primal infeasible 1" begin
        # min  x
        # s.t. y * exp( x / y) <= z
        #       y == 1
        #       z == -1

        P = spzeros(3, 3)
        q = [1.; 0; 0]

        # y * exp( x / y) <= z
        A1 = sparse(Matrix(1.0I, 3, 3))
        b1 = zeros(3)

        # y == 1
        A2 = [0 -1. 0]
        b2 = [-1.]

        # z == -1
        A3 = [0 0 -1.]
        b3 = [1]

        cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)
        cs2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)
        cs3 = COSMO.Constraint(A3, b3, COSMO.ZeroSet)

        model = COSMO.Model()
        assemble!(model, P, q, [cs1; cs2; cs3], settings = COSMO.Settings())

        res = COSMO.optimize!(model)
        @test res.status == :Primal_infeasible
    end

    @testset "Primal infeasible 2" begin
        # min  x
        # s.t. y * exp( x / y) <= z - 0.2
        #      y * exp( x / y) >= z + 0.2

        P = spzeros(3, 3)
        q = [1.; 0; 0]

        # y * exp( x / y) <= z - 0.2
        A1 = sparse(Matrix(1.0I, 3, 3))
        b1 = [0.; 0; -0.2]

        # y * exp( x / y) >= z + 0.3
        A2 = sparse(-Matrix(1.0I, 3, 3))
        b2 = [0.; 0; -0.3]

        cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)
        cs2 = COSMO.Constraint(A2, b2, COSMO.ExponentialCone)

        model = COSMO.Model()
        assemble!(model, P, q, [cs1; cs2], settings = COSMO.Settings())

        res = COSMO.optimize!(model)
        @test res.status == :Primal_infeasible
    end

    @testset "Dual infeasible" begin
        # max z
        # s.t. y * exp(x / y) <= z

        P = spzeros(3, 3)
        q = [0; 0; -1.0]

        # y * exp( x / y) <= z
        A1 = sparse(Matrix(1.0I, 3, 3))
        b1 = zeros(3)

        cs1 = COSMO.Constraint(A1, b1, COSMO.ExponentialCone)

        model = COSMO.Model()
        assemble!(model, P, q, cs1)

        res = COSMO.optimize!(model)
        @test res.status == :Dual_infeasible
    end

end

@testset "Dual Exponential cone problems" begin

    @testset "Feasible" begin

        # solve the following dual exponential cone problem
        # min  y
        # s.t. -x * exp(y / x) <= z
        #      x == -1, z == exp(5)

        P = spzeros(3, 3)
        q = [0; 1.; 0]

        # -x * exp( y / x) <= z
        A1 = sparse(Matrix(1.0I, 3, 3))
        b1 = zeros(3)
        cs1 = COSMO.Constraint(A1, b1, COSMO.DualExponentialCone)

        # y == 1 and z == exp(5)
        A2 = sparse([1. 0 0; 0 0 1])
        b2 = [1.; -exp(5)]
        cs2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)

        model = COSMO.Model()
        assemble!(model, P, q, [cs1; cs2])

        res = COSMO.optimize!(model)
        @test res.status == :Solved
        @test isapprox(res.obj_val, -6., atol=1e-3)
    end
end

nothing


text:
using FermiCG
using Test
using Random

Random.seed!(1234567)

@testset "FermiCG" begin
    include("test_hamiltonians.jl")
    include("test_cmf.jl")
    include("test_Clusters.jl")
    include("test_FCI.jl")
    include("test_davidson.jl")
    include("test_TDMs.jl")
    include("test_tpsci.jl")
    include("test_bounds.jl")
    include("test_s2.jl")
    include("test_bst.jl")
    include("test_schmidt.jl")
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG