text:
using Luna
import Luna.PhysData: wlfreq
import Logging
import FFTW
import NumericalIntegration: integrate, SimpsonEven
import LinearAlgebra: mul!, ldiv!
Logging.disable_logging(Logging.BelowMinLevel)

import PyPlot:pygui, plt

a = 225e-6
gas = :Ar
pres = 0.4

τfwhm = 30e-15
λ0 = 1800e-9
energy = 500e-6

modes = (Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0),
         Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=π/2))
nmodes = length(modes)

grid = Grid.RealGrid(250e-2, λ0, (200e-9, 3000e-9), 1e-12)

energyfun, energyfunω = Fields.energyfuncs(grid)

gausspulse(t) = sqrt.(Maths.gauss.(t; fwhm=τfwhm)) .* cos.(wlfreq(λ0).*t)

Etlin = gausspulse(grid.t)
cenergy = energyfun(Etlin)
Etlin = sqrt(energy)/sqrt(cenergy) .* Etlin

Et = [Etlin zero(grid.t)]
#Etd = size(Et)
#Etp = Array{ComplexF64,2}(undef, Etd[2], Etd[1])
#permutedims!(Etp, Et, [2, 1])

Ew = FFTW.rfft(Et, 1);
Ewd = size(Ew)
Ewp = Array{ComplexF64,2}(undef, Ewd[2], Ewd[1])
permutedims!(Ewp, Ew, [2, 1])

γ = deg2rad(10.0)
CP = Polarisation.rotate(Polarisation.WP(π/2), π/4 - γ)
Ewp .= CP * Ewp
permutedims!(Ew, Ewp, [2, 1])
#Et = FFTW.irfft(Ew, length(grid.t), 1);

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)

responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
             #Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

inputs = Fields.GaussField(λ0=λ0, τfwhm=τfwhm, energy=energy)

Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs, modes, :xy; full=false)

Eω .= Ew

statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
linop = LinearOps.make_const_linop(grid, modes, λ0)

Luna.run(Eω, grid, linop, transform, FT, output)

ω = grid.ω
t = grid.t

zout = output.data["z"]
Eout = output.data["Eω"]

Eoutd = size(Eout)
Ewp = Array{ComplexF64,2}(undef, Eoutd[2], Eoutd[1])
permutedims!(Ewp, Eout[:,:,end], [2, 1])

CP = Polarisation.rotate(Polarisation.WP(π/2), γ)
LP = Polarisation.rotate(Polarisation.LP(), π/2)
Ewp .= LP * CP * Ewp

Eoutp = Array{ComplexF64,2}(undef, Eoutd[1], Eoutd[2])
permutedims!(Eoutp, Ewp, [2, 1])

pygui(true)
plt.figure()
plt.plot(ω./2π.*1e-15, log10.(abs2.(Eout[:,:,end])))
plt.plot(ω./2π.*1e-15, log10.(abs2.(Eoutp)))
#plt.ylim(-6, 0)
plt.xlim(0,2.0)

##
plt.figure()
It = abs2.(Maths.hilbert(FFTW.irfft(Eout[:,:,end], length(grid.t), 1)))
plt.plot(t/1e-15, It[:,1]./maximum(It[:,1]))
plt.plot(t/1e-15, It[:,2]./maximum(It[:,2]))
Itp = abs2.(Maths.hilbert(FFTW.irfft(Eoutp, length(grid.t), 1)))
plt.plot(t/1e-15, Itp[:,1]./maximum(Itp[:,1]))
plt.plot(t/1e-15, Itp[:,2]./maximum(Itp[:,2]))
plt.xlim(-50.0,50.0)

println("$(sum(Itp)/sum(It))")
##
Etout = FFTW.irfft(Eout, length(grid.t), 1)
It = Maths.normbymax(abs2.(Maths.hilbert(Etout)))
Ilog = log10.(Maths.normbymax(abs2.(Eout)))
for i = 1:nmodes
    plt.figure()
    plt.subplot(121)
    plt.pcolormesh(ω./2π.*1e-15, zout, transpose(Ilog[:,i,:]))
    plt.clim(-6, 0)
    plt.xlim(0,2.0)
    plt.colorbar()
    plt.subplot(122)
    plt.pcolormesh(t.*1e15, zout, transpose(It[:,i,:]))
    plt.xlim(-30.0,100.0)
    plt.colorbar()
end

S = Array{Float64,3}(undef, 4, Eoutd[1], Eoutd[3])
elip = Array{Float64,2}(undef, Eoutd[1], Eoutd[3])
for i = 1:Eoutd[3]
    permutedims!(Ewp, Eout[:,:,i], [2, 1])
    for j = 1:Eoutd[1]
        S[:,j,i] .= Polarisation.Stokes(Ewp[:,j])
        elip[j,i] = Polarisation.ellipticity(S[:,j,i])
    end
end

plt.figure()
plt.pcolormesh(zout, ω./2π.*1e-15, S[3,:,:]./S[1,:,:])
plt.colorbar()

plt.figure()
plt.pcolormesh(zout, ω./2π.*1e-15, S[2,:,:]./S[1,:,:])
plt.colorbar()

plt.figure()
plt.pcolormesh(zout, ω./2π.*1e-15, elip)
plt.colorbar()

Eenvo = Maths.hilbert(Etout)
Etoutd = size(Eenvo)
Etp = Array{ComplexF64,2}(undef, Etoutd[2], Etoutd[1])

St = Array{Float64,3}(undef, 4, Etoutd[1], Etoutd[3])
elipt = Array{Float64,2}(undef, Etoutd[1], Etoutd[3])
Ss = Array{Float64,2}(undef, 4, Etoutd[3])
Elipt = Array{Float64,1}(undef, Etoutd[3])
CP = Polarisation.rotate(Polarisation.WP(π/2), γ)
for i = 1:Etoutd[3]
    permutedims!(Etp, Eenvo[:,:,i], [2, 1])
    for j = 1:Etoutd[1]
        St[:,j,i] .= Polarisation.Stokes(CP * Etp[:,j])
        elipt[j,i] = Polarisation.ellipticity(St[:,j,i])
    end
    k = argmax(St[1,:,i])
    Ss[:,i] .= St[:,k,i] ./ St[1,k,i]
    Elipt[i] = elipt[k,i]
end

plt.figure()
plt.plot(zout, Ss[1,:], label="S1")
plt.plot(zout, Ss[2,:], label="S2")
plt.plot(zout, Ss[3,:], label="S3")
plt.plot(zout, Ss[4,:], label="S4")
plt.plot(zout, Elipt, label="ε")
plt.legend()


text:
DATA_SIZE = 360
data_vector = randn(DATA_SIZE)
index_timetype = Date(2007, 1,1) + Day.(0:(DATA_SIZE - 1))
df_timetype_index = DataFrame(Index = index_timetype, data = data_vector)
ts_daily = TS(df_timetype_index, 1)
ts_daily_matrix = TS(DataFrames.innerjoin(df_timetype_index, df_timetype_index, on=:Index, makeunique=true))

# function apply(ts::TS, period::Union{T,Type{T}}, fun::V, index_at::Function=first) where {T<:Union{DatePeriod,TimePeriod}, V<:Function}

# test_apply_first
ts_monthly = apply(ts_daily, Dates.Month, first)
@test typeof(ts_monthly) == TSx.TS
@test DataFrames.nrow(ts_monthly.coredata) == 12 # 360 days
@test ts_monthly.coredata[1, :Index] == df_timetype_index[1, :Index]

# yearly
@test_broken ts_yearly = apply(ts_daily_matrix, Dates.Year, Statistics.mean)
@test typeof(ts_yearly) == TSx.TS
@test typeof(ts_yearly.coredata) == DataFrame
@test DataFrames.nrow(ts_yearly.coredata) == 1

ext: jl


text:
using Oceananigans.ImmersedBoundaries: ImmersedBoundaryGrid, GridFittedBoundary, RasterDepthMask
using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization

@testset "Immersed boundaries with hydrostatic free surface models" begin
    @info "Testing immersed boundaries with hydrostatic free surface models..."

    for arch in archs
        Nx = 60
        Ny = 60

        # A spherical domain
        underlying_grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1),
                                               longitude = (-30, 30),
                                               latitude = (15, 75),
                                               z = (-4000, 0))

        @inline raster_depth(i, j) = 30 < i < 35 && 42 < j < 48

        grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBoundary(raster_depth, mask_type=RasterDepthMask()))

        free_surface = ImplicitFreeSurface(gravitational_acceleration=0.1)
        coriolis = HydrostaticSphericalCoriolis(scheme = VectorInvariantEnstrophyConserving())

        surface_wind_stress_parameters = (τ₀ = 1e-4,
                                          Lφ = grid.Ly,
                                          φ₀ = 15)
        surface_wind_stress(λ, φ, t, p) = p.τ₀ * cos(2π * (φ - p.φ₀) / p.Lφ)

        surface_wind_stress_bc = FluxBoundaryCondition(surface_wind_stress,
                                               parameters = surface_wind_stress_parameters)

        μ = 1 / 60days

        @inline u_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.u[i, j, 1]
        @inline v_bottom_drag(i, j, grid, clock, fields, μ) = @inbounds - μ * fields.v[i, j, 1]

        u_bottom_drag_bc = FluxBoundaryCondition(u_bottom_drag,
                                         discrete_form = true,
                                         parameters = μ)

        v_bottom_drag_bc = FluxBoundaryCondition(v_bottom_drag,
                                         discrete_form = true,
                                         parameters = μ)

        u_bcs = FieldBoundaryConditions(top = surface_wind_stress_bc,
                                bottom = u_bottom_drag_bc)

        v_bcs = FieldBoundaryConditions(bottom = v_bottom_drag_bc)

        νh₀ = 5e3 * (60 / grid.Nx)^2
        constant_horizontal_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh₀)

        model = HydrostaticFreeSurfaceModel(grid = grid,
                                    architecture = CPU(),
                                    momentum_advection = VectorInvariant(),
                                    free_surface = free_surface,
                                    coriolis = coriolis,
                                    boundary_conditions = (u=u_bcs, v=v_bcs),
                                    closure = constant_horizontal_diffusivity,
                                    tracers = nothing,
                                    buoyancy = nothing)

        simulation = Simulation(model, Δt = 3600, stop_time = 3600)
        run!(simulation)

        # If reached here it didn't error, so pass for now!
        @test true
    end
end


text:
#=

# [Dispatch mechanism](@id Dispatch)

Some API functions cannot be called directly from the Vulkan library. In particular, extension functions *must* be called through a pointer obtained via API commands. In addition, most functions can be made faster by calling directly into their function pointer, instead of going through the loader trampoline which resolves the function pointer every time.

To circumvent that, we provide a handy way for retrieving and calling into function pointers, as well as a thread-safe global dispatch table that can automate this work for you.

## Retrieving function pointers

API Function pointers can be obtained with the [`function_pointer`](@ref) function, using the API function name.

=#

using SwiftShader_jll # hide
using Vulkan
@set_driver :SwiftShader # hide

const instance = Instance([], [])
const pdevice = first(unwrap(enumerate_physical_devices(instance)))
const device = Device(pdevice, [DeviceQueueCreateInfo(0, [1.0])], [], [])

#=

```@repl dispatch
function_pointer("vkCreateInstance")
function_pointer(instance, "vkDestroyInstance")
function_pointer(device, "vkCreateFence")
```

It is essentially a wrapper around `get_instance_proc_addr` and `get_device_proc_addr`, leveraging multiple dispatch to make it more intuitive to use.

## Providing function pointers

Every wrapper function or handle constructor has a signature which accepts the standard function arguments plus a function pointer to call into. For example, instead of doing

=#

foreach(println, unwrap(enumerate_instance_layer_properties()))

# you can do

fptr = function_pointer("vkEnumerateInstanceLayerProperties")
foreach(println, unwrap(enumerate_instance_layer_properties(fptr)))

# In the case of a handle constructor which calls both a creation and a destruction function, there is an argument for each corresponding function pointer:

fptr_create = function_pointer(device, "vkCreateFence")
fptr_destroy = function_pointer(device, "vkDestroyFence")
Fence(device, fptr_create, fptr_destroy)

#=

## Automatic dispatch

Querying and retrieving function pointers every time results in a lot of boilerplate code. To remedy this, we provide a concurrent global dispatch table which loads all available function pointers for loader, instance and device commands. It has one dispatch table per instance and per device to support using multiple instances and devices at the same time.

This feature can be disabled by setting the [preference](@ref Package-options) `USE_DISPATCH_TABLE` to `"false"`.

!!! warn
    All instance and device creations must be externally synchronized if the dispatch table is enabled. This is because all function pointers are retrieved and stored right after the creation of an instance or device handle.

=#

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
module CalibrationTests

using Reexport

@reexport using CalibrationErrors
using ConsistencyResampling: ConsistencyResampling
using HypothesisTests: HypothesisTests
@reexport using KernelFunctions
using LinearAlgebra: LinearAlgebra
using Random: Random
using Statistics: Statistics
using StatsFuns: StatsFuns
using StructArrays: StructArrays

using CalibrationErrors: CalibrationErrorEstimator, SKCE, unsafe_skce_eval

export ConsistencyTest
export DistributionFreeSKCETest, AsymptoticBlockSKCETest, AsymptoticSKCETest
export AsymptoticCMETest

# re-export
using HypothesisTests: pvalue, confint
export pvalue, confint

include("consistency.jl")

include("skce/asymptotic.jl")
include("skce/asymptotic_block.jl")
include("skce/distribution_free.jl")

include("cme.jl")

include("deprecated.jl")

end # module


text:
struct WindowManager{WM,W} <: AbstractWindowManager
    impl::WM
    active_keys::Dictionary{Symbol,EventDetails{KeyEvent{KeyPressed},W,Float64}}
    active_buttons::Dictionary{MouseButton,EventDetails{<:MouseEvent{ButtonPressed},W,Float64}}
    history::Dictionary{WindowAbstractions.Event,EventDetails}
end

function WindowManager(wm::XWindowManager)
    WindowManager{XWindowManager,XCBWindow}(
        wm,
        Dictionary{Symbol,EventDetails{KeyEvent{KeyPressed},XCBWindow,Float64}}(),
        Dictionary{MouseButton,EventDetails{<:MouseEvent{ButtonPressed},XCBWindow,Float64}}(),
        Dictionary{WindowAbstractions.Event,EventDetails}()
    )
end

@forward WindowManager.impl terminate_window!, get_window, get_window_symbol, callbacks, poll_for_event, wait_for_event, set_callbacks!

run(wm::WindowManager, mode::ExecutionMode, execute_callback; kwargs...) = run(wm.impl, mode, execute_callback; kwargs...)

function execute_callback(wm::WindowManager, ed::EventDetails)
    WindowAbstractions.execute_callback(wm.impl, ed)
    update!(wm, ed)
end

function update!(wm::WindowManager, ed::EventDetails{T}) where {T}
    if T <: KeyEvent{KeyPressed}
        set!(wm.active_keys, ed.data.key_name, ed)
    elseif T <: KeyEvent{KeyReleased}
        delete!(wm.active_keys, ed.data.key_name)
    elseif T <: MouseEvent{ButtonPressed}
        set!(wm.active_buttons, ed.data.button, ed)
    elseif T <: MouseEvent{ButtonReleased}
        delete!(wm.active_buttons, ed.data.button)
    end
    set!(wm.history, action(ed)(), ed)
end


text:
module RunBeast #module for running BEAST

export run_beast,
       check_beast,
       find_beast

const BEAST_JAR = "beast.jar"
const BEAST_HOME = "BEAST_HOME"

function find_beast(;beast_home::String="")
    if isempty(beast_home)
        beast_home = haskey(ENV, BEAST_HOME) ? ENV[BEAST_HOME] : pwd()
    end

    if basename(beast_home) == BEAST_JAR
        return abspath(beast_home)
    end

    path = joinpath(beast_home, BEAST_JAR)
    if isfile(path)
        beast_path = path
    else
        path = joinpath(beast_home, "build", "dist", BEAST_JAR)
        if isfile(path)
            beast_path = path
        else
            error("Could not find beast.jar file. Consider setting the " *
                  "$BEAST_HOME environment variable to the location of your " *
                  "BEAST installation or git repo. You may also specify the " *
                  "optional argument 'beast_jar=path/to/your/beast.jar'")
        end
    end

    return abspath(beast_path)
end

function check_beast(;beast_jar::String = find_beast())
    println("Checking Java installation...")
    run(`java -version`)

    println("\n\nChecking BEAST installation...")
    run(`java -jar $beast_jar -version`)

    println("\n\nJava and BEAST checks succeeded.")
end

function run_beast(xml_path::String;
                    seed::Int = -1, # don't set seed by default
                    overwrite::Bool = false, # don't overwrite log files by default
                    directory::String = pwd(), #don't change working directory by default
                    beast_jar::String = find_beast()
                    )

    old_directory = pwd()

    cd(directory)
    cmds = ["java", "-jar", beast_jar]

    if seed != -1
        push!(cmds, "-seed")
        push!(cmds, string(seed))
    end

    if overwrite
        push!(cmds, "-overwrite")
    end

    push!(cmds, "-fail_threads")

    push!(cmds, xml_path)

    out = run(Cmd(cmds))
    if out.exitcode != 0
        error("BEAST run threw an error. See console for error message.")
    end

    cd(old_directory)
end

end


text:
using TaskDispatcher
using Test

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


text:
# ---------------------------------------------------------------------------- #
#
#   egads.jl
#
#   Julia wrapper for egads (part of Engineering Sketch Pad)
#
#   egads.jl
#   Spring 2018
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

__precompile__()

"""
    sterno

Julia wrapper for egads (part of Engineering Sketch Pad)

Max Opgenoord

Spring 2018
"""
module egads

import Compat

# find egads library file
if Compat.Sys.isapple()
    ext = ".dylib"
elseif Compat.Sys.islinux()
    ext = ".so"
else
    error("Unsupported operating system")
end
libname      = string("libegads", ext)
commonlocesp = ["Code/ESP112/","code/","Code/"]
libdir       = "EngSketchPad/lib"
commonloc    = similar(commonlocesp)
for jj in 1:length(commonloc)
    commonloc[jj] = joinpath(homedir(),commonlocesp[jj],libdir)
end
libfind = Base.Libdl.find_library( libname, commonloc )
if isempty(libfind)
    error("EGADS library not found, update the location in egads.jl or update commonlocesp")
else
    const libegads_include = libfind
end

# include wrapper

include("wrapped_api/types_and_consts_egads.jl")
include("wrapped_api/include_egads.jl")

include("export.jl")

# include added functions
include("added_functionality.jl")

end


text:
# Autogenerated wrapper script for CUDAFastRegistrationLib_jll for x86_64-linux-gnu-cxx03
export libcuda_fast_registration

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("CUDAFastRegistrationLib")
JLLWrappers.@declare_library_product(libcuda_fast_registration, "libcuda_fast_registration.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcuda_fast_registration,
        "lib/libcuda_fast_registration.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "Monin_Obukhov_length" begin
    datetime, ustar, Tair, pressure, H = tha48[24,Cols(:datetime, :ustar, :Tair, :pressure, :H)]
    MOL24 = @inferred Monin_Obukhov_length(Tair, pressure, ustar, H)
    MOL24 = @inferred Monin_Obukhov_length(Tair, pressure, ustar, H; constants =BigleafConstants())
    @test ≈(MOL24, -104.3, rtol = 1/1000)
    #
    df = copy(tha48)
    @inferred Monin_Obukhov_length!(df)
    @test df.MOL[24] == MOL24
end

@testset "stability_parameter" begin
    df = copy(tha48)
    df24 = df[24,:]
    MOL24 = Monin_Obukhov_length(df24.Tair, df24.pressure, df24.ustar, df24.H)
    z=40;d=15
    zeta = @inferred stability_parameter(z,d,MOL24)
    @test ≈(zeta , -0.240, rtol = 1/100)
    zeta2 = @inferred stability_parameter(z,d,df24.Tair, df24.pressure, df24.ustar, df24.H)
    @test zeta2 == zeta
    #
    @inferred stability_parameter!(df;z,d)
    @test df.zeta[24] == zeta
    zeta_scalar = df.zeta
    #
    # also test zr as a vector
    df = copy(tha48)
    df[!,:zi] .= z
    df.zi[1] = z/2
    zetas = stability_parameter!(df; z=df.zi, d).zeta
    @test zetas[2:24] == zeta_scalar[2:24]
    @test zetas[1] != zeta_scalar[1]
    #
    # specifyngn MOL directly
    df4 = copy(tha48)
    Monin_Obukhov_length!(df4)
    zetas2 = stability_parameter!(df4; z, d, MOL = df4.MOL./2).zeta
    @test all(zetas2[2:end] .== zetas[2:end] .* 2)
end

@testset "stability_correction zeta" begin
    zeta = -2:0.5:0.5
    @inferred stability_correction(first(zeta))
    df2 = DataFrame(stability_correction.(zeta))
    @test all(isapprox.(df2.psi_h, SA[2.431,2.197,1.881,1.386,0,-2.5], rtol=1e-3))
    @test all(isapprox.(df2.psi_m, SA[2.275, 2.061, 1.772, 1.317, 0, -2.5], rtol=1e-3))
    df2 = DataFrame(stability_correction.(zeta; stab_formulation=Businger1971()))
    @test all(isapprox.(df2.psi_h, SA[2.085, 1.862, 1.564, 1.106,0, -3.9], rtol=1e-3))
    @test all(isapprox.(df2.psi_m, SA[2.418, 2.200, 1.904, 1.435,0, -3], rtol=1e-3))
    #
    # https://bitbucket.org/juergenknauer/bigleaf/issues/8/logical-inconsistency-possible-bug-with
    resm = @inferred stability_correction(missing)
    @test keys(resm) == (:psi_h, :psi_m)
    @test all(ismissing.(values(resm)))
    #
    resm = @inferred stability_correction(first(zeta);
        stab_formulation = NoStabilityCorrection())
    @test resm == (psi_h = 0.0, psi_m = 0.0)
end

@testset "stability_correction metvars" begin
    z=40;d=15
    df = DataFrame(Tair=25.0, pressure=100.0, ustar=0.2:0.1:1.0, H=40:20.0:200)
    df1 = df[1,:]
    zeta1 = stability_parameter(z,d,df1.Tair, df1.pressure, df1.ustar, df1.H)
    res1 = @inferred stability_correction(z,d, df1.Tair, df1.pressure, df1.ustar, df1.H)
    @test res1 == stability_correction(zeta1)
end

@testset "stability_correction DataFrame variant" begin
    z=40;d=15
    dfo = DataFrame(Tair=25.0, pressure=100.0, ustar=0.2:0.1:1.0, H=40:20.0:200)
    df = copy(dfo)
    res = stability_correction(df; z, d)
    # for type stability, use columntable(df)
    stability_correction!(df; z, d)
    propertynames(df)[(end-1):end] == SA[:psi_h, :psi_m]
    @test DataFrame(res) == df[!, (end-1):end]
    #
    dfm = allowmissing(dfo)
    dfm.ustar[1] = missing
    res = stability_correction(
        dfm; z, d, stab_formulation = NoStabilityCorrection())
    stability_correction!(
        dfm; z, d, stab_formulation = NoStabilityCorrection())
    propertynames(dfm)[(end-1):end] == SA[:psi_h, :psi_m]
    @test all(iszero.(dfm.psi_h))
    @test all(iszero.(dfm.psi_m))
    @test DataFrame(res) == dfm[!, (end-1):end]
    #
    # test specifying zeta instead of z and d
    df2 = copy(dfo)
    stability_parameter!(df2; z, d) # adds zeta
    stability_correction!(df2, zeta=df2.zeta)
    @test df2.psi_h == df.psi_h
    @test df2.psi_m == df.psi_m
    #
    # test specifying zr as a vector
    df3 = copy(dfo)
    df3[!,:zi] .= z
    df3.zi[1] = z/2
    res = stability_correction(df3; z=df3.zi, d)
    @inferred stability_correction!(df3; z=df3.zi, d)
    @test df3.psi_h[2:end] == df.psi_h[2:end]
    @test df3.psi_h[1] != df.psi_h[1]
    @test DataFrame(res) == df3[!, (end-1):end]
end

@testset "stability_correction DataFrame variant Float32" begin
    z=40;d=15
    dfo = DataFrame(Tair=25.0, pressure=100.0, ustar=0.2:0.1:1.0, H=40:20.0:200)
    float_cols = names(dfo, All())
    dfo = transform(dfo,  float_cols .=> ByRow(passmissing(Float32)) .=> float_cols)

    df = copy(dfo)
    res = stability_correction(df; z, d)
    @test eltype(res.psi_h) == Float32
    # for type stability, use columntable(df)
    stability_correction!(df; z, d)
    propertynames(df)[(end-1):end] == SA[:psi_h, :psi_m]
    @test eltype(df.psi_h) == Float32
    @test DataFrame(res) == df[!, (end-1):end]
    #
    dfm = allowmissing(dfo)
    dfm.ustar[1] = missing
    res = stability_correction(
        dfm; z, d, stab_formulation = NoStabilityCorrection())
    @test eltype(res.psi_h) == Float32
    stability_correction!(
        dfm; z, d, stab_formulation = NoStabilityCorrection())
    propertynames(dfm)[(end-1):end] == SA[:psi_h, :psi_m]
    @test all(iszero.(dfm.psi_h))
    @test all(iszero.(dfm.psi_m))
    #@test eltype(dfm.psi_h) == typeof(z) #Float32
    @test DataFrame(res) == dfm[!, (end-1):end]
    #
    # test specifying zeta instead of z and d
    df2 = copy(dfo)
    stability_parameter!(df2; z, d) # adds zeta
    stability_correction!(df2, zeta=df2.zeta)
    @test eltype(df2.psi_h) == Float32
    @test df2.psi_h == df.psi_h
    @test df2.psi_m == df.psi_m
    #
    # test specifying zr as a vector
    df3 = copy(dfo)
    df3[!,:zi] .= z
    df3.zi[1] = z/2
    df3_before = copy(df3)
    res = stability_correction(df3; z=df3.zi, d)
    @test df3 == df3_before
    @test eltype(res.psi_h) == Float32
    @inferred stability_correction!(df3; z=df3.zi, d)
    @test eltype(df3.psi_h) == Float32
    @test df3.psi_h[2:end] == df.psi_h[2:end]
    @test df3.psi_h[1] != df.psi_h[1]
    @test DataFrame(res) == df3[!, (end-1):end]
end


text:
"""
    Fermi.CoupledCluster

Module for running CoupledCluster computations in Fermi.
"""
module CoupledCluster
# Import Fermi basics
using Fermi
using Fermi.Options
using Fermi.Integrals
using Fermi.Orbitals

function cc_header()
    banner =
raw"""
================================================================================
//    _____                   _          _   _____ _           _              \\
//   /  __ \                 | |        | | /  __ \ |         | |             \\
//   | /  \/ ___  _   _ _ __ | | ___  __| | | /  \/ |_   _ ___| |_ ___ _ __   \\
//   | |    / _ \| | | | '_ \| |/ _ \/ _` | | |   | | | | / __| __/ _ \ '__|  \\
//   | \__/\ (_) | |_| | |_) | |  __/ (_| | | \__/\ | |_| \__ \ ||  __/ |     \\
//    \____/\___/ \__,_| .__/|_|\___|\__,_|  \____/_|\__,_|___/\__\___|_|     \\
//                    | |                                                     \\
//                    |_|                                                     \\
//                                                                            \\
//                 Module by G.J.R. Aroeira and M. M. Davis                   \\
================================================================================
"""
    output(banner)
end

"""
    Fermi.CoupledCluster.AbstractCCWavefunction

Fermi abstract type common to all Coupled Cluster wavefunctions

_struct tree:_

**AbstractCCWavefunction** <: AbstractWavefunction
"""
abstract type AbstractCCWavefunction <: Fermi.AbstractWavefunction end

#include("RCCD/RCCD.jl")
include("RCCSD/RCCSD.jl")
include("PerturbativeTriples/PerturbativeTriples.jl")

end #module CC

max_stars_repo_name: FermiQC/Fermi

max_issues_repo_name: FermiQC/Fermi

max_forks_repo_name: FermiQC/Fermi


text:
module SEAL

# SEAL_jll provides `libsealc`, which we will use in this package
using SEAL_jll

"""
    SEALObject

Abstract parent type for all types based on SEAL classes.
"""
abstract type SEALObject end

"""
    gethandle(object::SEALObject)

Return the raw C pointer to where `object` resides in memory.
"""
@inline gethandle(object::SEALObject) = object.handle

"""
    sethandle!(object::SEALObject, handle)

Set the underlying raw C pointer to where `object` resides in memory to `handle`.
"""
@inline sethandle!(object::SEALObject, handle) = object.handle = handle

"""
    destroy!(object::SEALObject)

Call the corresponding destruction function on `object` to free up memory and reset object handle to
a null pointer. If `object` is not allocated, `destroy!` will not do anything.
"""
function destroy!(object::SEALObject) end

"""
    isnull(object::SEALObject)

Return true if the object handle is a null pointer and false otherwise.
"""
@inline isnull(object::SEALObject) = gethandle(object) == C_NULL

"""
    isallocated(object::SEALObject)

Return true if the object is allocated, i.e., if it is not null.
"""
@inline isallocated(object::SEALObject) = !isnull(object)

export SEALObject, gethandle, sethandle!, destroy!, isnull, isallocated

Base.unsafe_convert(::Type{Ptr{Cvoid}}, object::SEALObject) = gethandle(object)

include("auxiliary.jl")
# Julia-only auxiliary methods -> no exports

include("version.jl")
export version_major, version_minor, version_patch, version

include("modulus.jl")
export Modulus, SecLevelType, bit_count, value, coeff_modulus_create, coeff_modulus_bfv_default

include("serialization.jl")
export ComprModeType, SEALHeader, load_header!

include("encryptionparams.jl")
export EncryptionParameters, SchemeType, poly_modulus_degree,
       set_poly_modulus_degree!, set_coeff_modulus!, coeff_modulus,
       scheme, plain_modulus, set_plain_modulus!, plain_modulus_batching, parms_id, save!,
       save_size, load!

include("context.jl")
export SEALContext, first_parms_id, last_parms_id, get_context_data, key_context_data,
       first_context_data, parameter_error_message, using_keyswitching
export ContextData, chain_index, parms, parms_id, total_coeff_modulus_bit_count, qualifiers,
       next_context_data
export EncryptionParameterQualifiers, using_batching

include("publickey.jl")
export PublicKey, parms_id

include("secretkey.jl")
export SecretKey, parms_id, save!, load!

include("galoiskeys.jl")
export GaloisKeys, parms_id

include("relinkeys.jl")
export RelinKeys, parms_id, save_size, save!, load!

include("keygenerator.jl")
export KeyGenerator, create_public_key!, create_public_key, secret_key, create_relin_keys!,
       create_relin_keys, create_galois_keys!

include("plaintext.jl")
export Plaintext, scale, scale!, parms_id, to_string, save_size, save!

include("ciphertext.jl")
export Ciphertext, scale, scale!, parms_id, size, length, save_size, save!, load!, reserve!

include("encryptor.jl")
export Encryptor, set_secret_key!, encrypt!, encrypt_symmetric, encrypt_symmetric!

include("evaluator.jl")
export Evaluator, square!, square_inplace!, relinearize!, relinearize_inplace!, rescale_to_next!,
       rescale_to_next_inplace!, multiply_plain!, multiply_plain_inplace!, multiply!,
       multiply_inplace!, mod_switch_to!, mod_switch_to_inplace!, mod_switch_to_next!,
       mod_switch_to_next_inplace!, add!, add_inplace!, add_plain!, add_plain_inplace!,
       rotate_vector!, rotate_vector_inplace!, rotate_rows!, rotate_rows_inplace!,
       rotate_columns!, rotate_columns_inplace!, complex_conjugate!, complex_conjugate_inplace!,
       negate!, negate_inplace!, sub!, sub_inplace!

include("decryptor.jl")
export Decryptor, decrypt!, invariant_noise_budget

include("ckks.jl")
export CKKSEncoder, slot_count, encode!, decode!

include("batchencoder.jl")
export BatchEncoder, slot_count, encode!, decode!

include("memorymanager.jl")
export MemoryPoolHandle, alloc_byte_count, memory_manager_get_pool

end


text:
@testset "Coordinate system transformations" begin
    # Global coordinates
    lla = LLA(lat=-24.007712762068806, lon=96.44002819812094, alt=1374.7804632852078)
    lla2 = LLA(lat=-24.1, lon=96.5, alt=1659.18948659927)

    ecef_lla = ECEFfromLLA(wgs84)
    lla_ecef = LLAfromECEF(wgs84)

    ecef = ecef_lla(lla)
    ecef2 = ecef_lla(lla2)

    @test ecef ≈ ECEF(-654007.2768949205,5.794061748224952e6,-2.5796231395720555e6)
    @test lla ≈ lla_ecef(ecef)

    lla_utm = LLAfromUTM(utm_zone(lla)..., wgs84)
    utm_lla = UTMfromLLA(utm_zone(lla)..., wgs84)

    utm = utm_lla(lla)
    utm2 = utm_lla(lla2)

    @test utm ≈ UTM(239579.67583179142,7.342551466042985e6,1374.7804632852078)
    @test lla ≈ lla_utm(utm)
    # Test also the poles
    @test UTMfromLLA(0,true,wgs84)(LLA(89.0, 89.0, 89.0)) ≈ UTM(2.111009610242531e6, 1.9980623200455606e6, 89.0)
    @test LLAfromUTM(0,true,wgs84)(UTM(2.111009610242531e6, 1.9980623200455606e6, 89.0)) ≈ LLA(89.0, 89.0, 89.0)

    lla_utmz = LLAfromUTMZ(wgs84)
    utmz_lla = UTMZfromLLA(wgs84)

    utmz = utmz_lla(lla)
    utmz2 = utmz_lla(lla2)
    # Test also the poles go to UPS
    @test utmz_lla(LLA(89.0, 89.0, 89.0)) ≈ UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true)
    @test lla_utmz(UTMZ(2.111009610242531e6, 1.9980623200455606e6, 89.0, 0, true)) ≈ LLA(89.0, 89.0, 89.0)

    @test utmz ≈ UTMZ(239579.67583179142,7.342551466042985e6,1374.7804632852078, UInt8(47), false)
    @test lla ≈ lla_utmz(utmz)

    utmz_utm = UTMZfromUTM(47, false, wgs84)
    utm_utmz = UTMfromUTMZ(47, false, wgs84)

    @test utmz == utmz_utm(utm)
    @test utm == utm_utmz(utmz)
    # Sometimes this has to do something non-trivial - make sure it does
    utm_utmz2 = UTMfromUTMZ(46, false, wgs84)
    @test utm_utmz2(utmz) ≈ UTM(850009.7418418773, 7.340641097689279e6, 1374.7804632852078)

    # ENU coordinates
    enu_ecef = ENUfromECEF(lla2, wgs84)
    ecef_enu = ECEFfromENU(lla2, wgs84)

    enu = enu_ecef(ecef)

    @test enu_ecef ≈ ENUfromECEF(ecef2, wgs84) # (Can supply origin an any format)
    @test enu_ecef ≈ ENUfromECEF(utm2, 47, false, wgs84)
    @test enu_ecef ≈ ENUfromECEF(utmz2, wgs84)

    @test enu ≈ ENU(-6103.186938282723, 10222.54767562731, -295.55857190545794)
    @test ecef_enu(ENU(-6103.186938282723, 10222.54767562731, -295.55857190545794)) ≈ ecef

    # https://github.com/JuliaGeo/Geodesy.jl/issues/49
    @test LLA(ECEF(1, 0, 0), wgs84) ≈ LLA(89.99866260445, 0.00000000000, -6356752.314234)

    # Inverses
    @test inv(ecef_lla) == lla_ecef
    @test inv(lla_ecef) == ecef_lla
    @test inv(utm_lla) == lla_utm
    @test inv(lla_utm) == utm_lla
    @test inv(utmz_lla) == lla_utmz
    @test inv(lla_utmz) == utmz_lla
    @test inv(utmz_utm) == utm_utmz
    @test inv(utm_utmz) == utmz_utm
    @test inv(enu_ecef) == ecef_enu
    @test inv(ecef_enu) == enu_ecef

    # Composed transformation functions
    @test LLAfromENU(ecef, wgs84) == LLAfromECEF(wgs84) ∘ ECEFfromENU(ecef, wgs84)
    @test ENUfromLLA(ecef, wgs84) == ENUfromECEF(ecef, wgs84) ∘ ECEFfromLLA(wgs84)

    @test UTMfromECEF(47, false, wgs84) == UTMfromLLA(47, false, wgs84) ∘ LLAfromECEF(wgs84)
    @test ECEFfromUTM(47, false, wgs84) == ECEFfromLLA(wgs84) ∘ LLAfromUTM(47, false, wgs84)

    @test UTMZfromECEF(wgs84) == UTMZfromLLA(wgs84) ∘ LLAfromECEF(wgs84)
    @test ECEFfromUTMZ(wgs84) == ECEFfromLLA(wgs84) ∘ LLAfromUTMZ(wgs84)

    @test ENUfromUTMZ(ecef, wgs84) == ENUfromLLA(ecef, wgs84) ∘ LLAfromUTMZ(wgs84)
    @test UTMZfromENU(ecef, wgs84) == UTMZfromLLA(wgs84) ∘ LLAfromENU(ecef, wgs84)

    @test ENUfromUTM(ecef, 47, false, wgs84) == ENUfromLLA(ecef, wgs84) ∘ LLAfromUTM(47, false, wgs84)
    @test UTMfromENU(ecef, 47, false, wgs84) == UTMfromLLA(47, false, wgs84) ∘ LLAfromENU(ecef, wgs84)
end


text:
# This file defines different sets of functions that are used for comparisons of some of our
# opacity calculations. Each set of functions is organized into it's own submodule.
#
# Tests are explicitly NOT defined in this file so that functions defined in this file can be used
# in scripts outside of this framework (to create plots in order to assess disagreement)

# ================================================================================================

# The first set of functions are defined to facillitate comparisons against the values shown in
# several panels of Figure 8.5 from Gray's 2005 edition of "The Observation and Analysis of Stellar
# Photospheres".

# For these panels, data was extracted using the WebPlotDigitizer tool. We've recorded values with
# greater precision than the values should actually have. The values are primarily meant to be used
# to make sure that our implementation give results to the correct order of magniutde and
# qualitatively has the expected wavelength dependences.

module Gray_opac_compare
using Korg, HDF5
using Korg: Interval, closed_interval, contained
using Interpolations: LinearInterpolation, Throw

# Load Gray05 data for a given panel. This returns a tuple holding two dictionaries:
# 1. The first holds the panel properties (i.e. the temperature and partial electron pressure that
#    was assumed for the calculation)
# 2. The second holds a 2-tuple for each source of opacity source. The first tuple entry specifies
#    wavelengths (in Ångstroms) and the second entry holds the opacity contribution at those
#    wavelengths.
function load_panel_data(panel_name, h5fname = joinpath(@__DIR__, "data/gray05_fig8.5.h5"))
    attribute_dict = Dict{String,Float64}()
    data_dict = Dict{String,Tuple{Array{Float64,1},Array{Float64,1}}}()

    h5open(h5fname, "r") do file
        # at the highest level, the file is divided into one group per panel
        panel_group = file[panel_name]

        # load the attributes from the panel group. This nominally holds the temperature and
        # partial electron pressure for which the values are computed.
        for key in keys(attributes(panel_group))
            attribute_dict[key] = read(attributes(panel_group)[key])
        end

        # In a given panel_group there is one subgroup per opacity source.
        for opacity_subgroup_name in keys(panel_group)
            opacity_subgroup = panel_group[opacity_subgroup_name]
            if !(length(keys(opacity_subgroup)) == 2 &&
                 haskey(opacity_subgroup, "lambda") &&
                 haskey(opacity_subgroup, "data"))
                msg = string("Each opacity subgroup is expected to just have a \"lambda\" ",
                             "dataset and a \"data\" dataset")
                throw(ErrorException(msg))
            end
            tmp = (read(opacity_subgroup["lambda"]), read(opacity_subgroup["data"]))
            if length(tmp[1]) != length(tmp[2])
                throw(ErrorException("Datasets within a given subgroup should have a fixed length"))
            end
            data_dict[opacity_subgroup_name] = tmp
        end
    end
    (attribute_dict, data_dict)
end

# To facillitate comparisons against opacity from H₂⁺ and He⁻, we need to have reasonable estimates
# of the ratio of free electrons to protons. In the future, we might want to hardcode the
# calculated values so that changes to the atomic data doesn't break all of these tests.

# from table D.2. This interpolates the log of the partition function.
# I updated cases were the table records 0.301 with log10(2.0)
function _gray05_H_I_partition_func(T)
    log_2 = log10(2.0)
    θ_vals = [  0.2,   0.4,   0.6,   0.8,   1.0,   1.2,   1.4,   1.6,   1.8,   2.0]
    table  = [0.368, 0.303, log_2, log_2, log_2, log_2, log_2, log_2, log_2, log_2]
    interpolator = LinearInterpolation(θ_vals, table, extrapolation_bc=Throw())
    10^interpolator(5040.0 / T)
end

# compute the ratio of the nₑ to n_H (the number density of all H I and H II)
# this function could be reused in the future for testing solutions to coupled Saha equations
function free_electrons_per_Hydrogen_particle(nₑ, T, abundances = Korg.solar_abundances,
                                              ionization_energies = Korg.ionization_energies,
                                              partition_funcs = Korg.partition_funcs)
    out = 0.0
    for element in 1:Korg.Natoms
        wII, wIII = Korg.saha_ion_weights(T, nₑ, element, Korg.ionization_energies,
                                          Korg.partition_funcs)

        nₑ_per_ndens_species = wII/(1 + wII + wIII) + 2wIII/(1 + wII + wIII)
        abundance = 10.0^(abundances[element]-12.0)
        out += abundance * nₑ_per_ndens_species
    end
    out
end

# At a given Temperature and number density of electrons, generate semi-realistic values for:
#    - number densities of H I
#    - number density of all hydgrogen particles (H I and H II)
#    - the mass density
#
# This is mostly used for comparisons against panels of figure 8.5 from Gray (2005). For most
# curves, we can simply assert that there is a fixed ionization fraction because we divide out the
# dependence on all of these terms during the calculation
function _semi_realisitic_dens(ne::F, fion::F= 0.02, HydrogenMassFrac::F= 0.76) where F <: Real
    # In the interest of semi-realistic numbers, Let's assume all free electrons come from Hydrogen
    nH_I = ne/fion
    nH = ne + nH_I
    ρ = nH * 1.67e-24/HydrogenMassFrac
    (nH_I, nH, ρ)
end

# Gray05 uses some confusing terminology and notation.
# - he defines the "atomic absorption coefficient" as the area per absorber and represent it with
#   the variable α. For reference, Rybicki & Lightman call this same quantity the cross-section and
#   represent it with the variable σ_ν.
# - he defines the "continuous absorption coefficient per neutral hydrogen atom" and denotes it
#   with κ. Rybicki & Lightman don't explicitly define this quantity, but it's compatible with
#   their nomenclature. R&L would express this quantity as α_ν/nH_I (in their notation). This is
#   NOT a cross-section in the general case (e.g. for continuum absorption of He, Gray still
#   defines the "continuous absorption coefficient per neutral hydrogen atom" as κ = α_ν/nH_I).
# - he uses κ_ν to denote the mass opacity (just like Rybicki & Lightman)
#
# Now we actually define the functions that compute the opacities in the form comparable with
# Gray05. These functons return the "continuous absorption coefficient per neutral hydrogen atom"
# divided by Pₑ, the partial electron pressure.
#
# For the sake of clarity, we label Rybicki & Lightman's α_ν = κ_ν*ρ as the "linear absorption
# coefficient" within these functions (note: throughout the rest of the codebase, we interchangably
# use the simpler term, "absorption coefficient", to refer to this same quantity)

# Combined H I bound-free and free-free opacity
function HI_coefficient(λ, T, Pₑ, H_I_ion_energy = 13.598)
    # λ should have units of Ångstroms
    # Pₑ is the partial pressure of the electrons in dyne/cm²
    nₑ =  Pₑ/(Korg.kboltz_cgs * T)
    ν = (Korg.c_cgs*1e8)/λ

    bf_coef = begin
        H_I_partition_val = 2.0 # implicitly in the implementation provided by Gray (2005)
        nH_I = nₑ * 100.0 # this is totally arbitrary
        bf_linear_absorption_coef = Korg.ContinuumAbsorption.H_I_bf([ν], T, nH_I/H_I_partition_val,
                                                                    H_I_ion_energy)[1]
        bf_linear_absorption_coef / (Pₑ * nH_I)
    end

    ff_coef = begin
        nH_total = nₑ * 100.0 # this is totally arbitrary
        ρ = nH_total * 1.67e-24/0.76 # this is totally arbitrary

        #Assume U_I(T) = 2.0 and U_II(T) = 1.0, as in Gray (2005)
        Us = Dict([Korg.species"H_I"=>(T -> 2.0), Korg.species"H_II"=>(T -> 1.0)])
        χs = [(Korg.RydbergH_eV, -1.0, -1.0)]
        wII, wIII = Korg.saha_ion_weights(T, nₑ, 1, χs, Us)

        nH_I = nH_total / (1 + wII)
        nH_II = nH_total * wII/(1 + wII)

        # compute the linear absorption coefficient  = dτ/ds = opacity*ρ
        ff_linear_absorption_coef = Korg.ContinuumAbsorption.H_I_ff([ν], T, nH_II, nₑ)[1]
        ff_linear_absorption_coef / (Pₑ * nH_I)
    end

    bf_coef + ff_coef
end

function Hminus_bf_coefficient(λ, T, Pₑ, ion_energy_H⁻ = 0.7552)
    # λ should have units of Ångstroms
    # Pₑ is the partial pressure of the electrons in dyne/cm²
    ne =  Pₑ/(Korg.kboltz_cgs * T)

    # the values of ndens, ρ and nH_I shouldn't actually matter since we will
    # just divide out all dependence on these variables.
    nH_I, nH, ρ = _semi_realisitic_dens(ne)

    partition_func = 2.0 # may want to include the temperature dependence of the partition function
    ν = (Korg.c_cgs*1e8)/λ

    linear_absorb_coef = Korg.ContinuumAbsorption.Hminus_bf([ν], T, nH_I/partition_func, ne,
                                                            ion_energy_H⁻)[1]
    linear_absorb_coef / (Pₑ * nH_I)
end

function Hminus_ff_coefficient(λ, T, Pₑ)
    # λ should have units of Ångstroms
    # Pₑ is the partial pressure of the electrons in dyne/cm²
    ne =  Pₑ/(Korg.kboltz_cgs * T)

    # the values of ndens, ρ and nH_I shouldn't actually matter since we will
    # just divide out all dependence on these variables.
    nH_I, nH, ρ = _semi_realisitic_dens(ne)

    partition_func = 2.0 # may want to include the temperature dependence of the partition function
    ν = (Korg.c_cgs*1e8)/λ

    linear_absorb_coef = Korg.ContinuumAbsorption.Hminus_ff([ν], T, nH_I/partition_func, ne)[1]
    linear_absorb_coef / (Pₑ * nH_I)
end

# computes the combine H₂⁺ free-free and bound-free absorption in units of cm^2 per H atom (not a
# typo)
# There seems to be an error in our function to compute H₂⁺
function H2plus_coefficient(λ, T, Pₑ)
    # λ should have units of Ångstroms
    # Pₑ is the partial pressure of the electrons in dyne/cm²
    ne =  Pₑ/(Korg.kboltz_cgs * T)

    # in the future, we may want to use the abundances, ionization_energies and partition functions
    # used in Gray (2005) since we are comparing to his plots
    ne_div_nH =
        free_electrons_per_Hydrogen_particle(ne, T)
    nH = ne/ne_div_nH

    wII, wIII = Korg.saha_ion_weights(T, ne, 0x01, Korg.ionization_energies,
                                         Korg.partition_funcs)
    nH_I = nH / (1 + wII)
    nH_II = nH * wII / (1 + wII)

    # set the partition function to 2.0 in nH_I_div_partition so that n(H I, n=1) = n(H I) for this
    # calculation, which is an assumption that Gray (2005) implicitly uses
    nH_I_div_partition = nH_I/2.0
    ρ = 1.0 # arbitrary value because we divide it out after

    ν = (Korg.c_cgs*1e8)/λ
    linear_absorb_coef = Korg.ContinuumAbsorption.H2plus_bf_and_ff([ν], T, nH_I, nH_II)[1]
    linear_absorb_coef / (Pₑ * nH_I)
end

# compute He⁻ free-free absorption in units of cm^2 per H atom (not a typo)
function Heminus_ff_coefficient(λ, T, Pₑ)
    # λ should have units of Ångstroms
    # Pₑ is the partial pressure of the electrons in dyne/cm²
    ne =  Pₑ/(Korg.kboltz_cgs * T)

    # the values of ndens, ρ and nH_I shouldn't actually matter since we will
    # just divide out all dependence on these variables.

    # In the interest of semi-realistic numbers, Let's assume all free
    # electrons come from Hydrogen
    nH_I, nH, ρ = _semi_realisitic_dens(ne)

    # we need to compute the number density of He I. Since this function exists to facillitate
    # comparisons against plots from Gray (2005), we should employ the abundances from
    # their book. They record that there is 8.51e-2 He particles per H particle in table 16.3
    nHe = nH*8.51e-2
    # Gray accounts for only neutral and singly ionized helium in his calculation.
    # We use the ionization energy tabulated in table D.1 (χ(He I) = 24.587).
    # Table D.2 of Gray (2005) records the partion function of He II to be
    # as 10^0.301, the closest value at the table's precision to 2.0
    UI, UII = (1.0, 10^-0.301)
    # Because this neglects He III, we don't use Korg.saha_ion_weights.
    wII = (2.0/ne * (UII/UI) * Korg.translational_U(Korg.electron_mass_cgs, T)
           * exp(-24.587/(Korg.kboltz_eV*T)))
    nHe_I = nHe * 1/(1+wII)

    ν = (Korg.c_cgs*1e8)/λ
    linear_absorb_coef = Korg.ContinuumAbsorption.Heminus_ff([ν], T, nHe_I/UI, ne)[1]
    linear_absorb_coef / (Pₑ * nH_I)
end

# There appears to be some errors in the H₂⁺ opacities, skipping them for now
const Gray05_opacity_form_funcs =
    Dict("H"          => (HI_coefficient,
                          Interval(0, Inf), "H I bf and ff"),
         "Hminus_bf"  => (Hminus_bf_coefficient,
                          closed_interval(2250.0, 15000.0),  "H⁻ bound-free"),
         "Hminus_ff"  => (Hminus_ff_coefficient,
                          closed_interval(2604.0, 113918.0), "H⁻ free-free"),
         "Heminus_ff" => (Heminus_ff_coefficient,
                          closed_interval(5063.0, 151878.0), "He⁻ free-free"),
         "H2plus"     => (H2plus_coefficient,
                          closed_interval(3847.0, 25000.0),   "H₂⁺ ff and bf"),
         )

# the absolute tolerances for values the opacity contributions from
#    - H⁻ bound-free
#    - H⁻ free-free
#    - He⁻ free-free
# There is a much larger discrepancy between the summed H I bf and ff opacity at large λ. But
# unlike for H₂⁺, this could simply be a consequence of the different approximations that are used
const Gray05_atols = Dict("a" => 0.05, "b" => 0.02, "c" => 0.04)

# this function is defined to make tests easier
function Gray05_comparison_vals(panel, opacity_func_name)
    panel_prop, panel_dict = load_panel_data(panel)
    orig_λ_vals, ref_data = panel_dict[opacity_func_name]
    temperature = panel_prop["temperature"]
    Pₑ = 10^panel_prop["log10(electron_pressure)"]

    func, bounds = Gray05_opacity_form_funcs[opacity_func_name][1:2]

    w = contained.(orig_λ_vals, Ref(bounds))
    calculated_vals = func.(orig_λ_vals[w], temperature, Pₑ)
    (calculated_vals*1e26, ref_data[w])
end

end # end the definition of the Gray_opac_compare module

# ================================================================================================

# Below, we define a set of functions that are used to compare our implementation of bound-free
# opacity calculations using data from the Opacity Project against calculations of the same
# opacities using the hydrogenic aproximation. This nominally means that we're comparing
# opacities from H I and He II
#
# In the future, it might also be nice to include comparisons for multi-electron data against
# experimental data.
module OP_compare
using Korg

function calc_hydrogenic_bf_absorption_coef(λ_vals,  T, ndens_species, species_name;
                                            use_OP_data = false)
    @assert species_name in ["H_I", "He_II"]

    if use_OP_data
        cross_sec_file = if species_name == "H_I"
            joinpath(@__DIR__, "data/TOPbase_cross_section_H_I.txt")
        else
            joinpath(@__DIR__, "data/TOPbase_cross_section_He_II.txt")
        end
        Korg.ContinuumAbsorption.absorption_coef_bf_TOPBase(λ_vals.*1e-8, [T], [ndens_species],
                                                            Korg.Species(species_name);
                                                            extrapolation_bc = 0.0,
                                                            cross_sec_file = cross_sec_file)[:, 1]
    else
        ν_vals = (Korg.c_cgs*1e8)./λ_vals # Hz
        ndens_div_partition = ndens_species/Korg.partition_funcs[Korg.Species(species_name)](T)
        κ_dflt_approach = if species_name == "H_I"
            H_I_ion_energy = 13.598
            Korg.ContinuumAbsorption.H_I_bf(ν_vals, T, ndens_div_partition, H_I_ion_energy)
        else
            He_II_ion_energy = 54.418
            Korg.ContinuumAbsorption.He_II_bf(ν_vals, T, ndens_div_partition, He_II_ion_energy)
        end
    end
end

# the following functions could all be refactored so that they are more concise... (they could
# honestly just be constants)

function _dflt_λ_vals(species_name)
    first, last = (species_name == "H_I") ? (80,80000) : (25, 40000)
    map((λ) -> Float32(λ), first:1.0:last)
end

function _λ_comp_intervals(species_name)
    # because the energy state data used by the Opacity Project is somewhat inaccurate (note the
    # disagreement gets worse at higher energy levels), we can only make meaningful comparisons
    # over specific wavelength intervals.

    if species_name == "H_I"
        [80.0       911.0;
         1000.0     3645.0;
         4000.0     8203.0;
         9300.0     14579.0;
         18000.0    22789.0]
    else
        [25.0       10.0^2.35;
         10.0^2.37  10.0^2.95;
         10.0^2.99  10.0^3.3;
         10.0^3.37  10.0^3.55;
         10.0^3.65  10.0^3.75]
    end
end

_hydrogenic_rtol(species_name) = (species_name == "H_I") ? 0.11 : 0.06
end # end the definition of OP_compare module


text:
"""
$(TYPEDEF)

Pose2 is a SE(2) mechanization of two Euclidean translations and one Circular rotation, used for general 2D SLAM.
"""
struct Pose2 <: IncrementalInference.InferenceVariable
  dims::Int
  labels::Vector{String}
  manifolds::Tuple{Symbol,Symbol,Symbol}
  Pose2(;labels::Vector{<:AbstractString}=String[]) = new(3, labels, (:Euclid, :Euclid, :Circular))
end


text:
function serialize(vis::CoreVisualizer, queue::CommandQueue)
    utime = time_ns()
    delete_cmds = Dict{String, Any}[]
    setgeometry_cmds = Dict{String, Any}[]
    settransform_cmds = Dict{String, Any}[]
    for path in queue.delete
        push!(delete_cmds, Dict("path" => path))
    end
    for path in queue.setgeometry
        visdata = vis.tree[path].data
        push!(setgeometry_cmds, serialize(path, visdata.geometries))
    end
    for path in queue.settransform
        visdata = vis.tree[path].data
        tform = serialize(visdata.transform)
        push!(settransform_cmds,
              Dict{String, Any}("path" => path,
                                "transform" => tform
              )
        )
    end
    data = Dict{String, Any}(
        "utime" => utime,
        "delete" => delete_cmds,
        "setgeometry" => setgeometry_cmds,
        "settransform" => settransform_cmds
    )
end

function serialize(path::AbstractVector, geomdatas::Vector{GeometryData})
    params = serialize.(geomdatas)
    if length(params) == 1
        Dict("path" => serialize(path),
             "geometry" => params[1])
    else
        Dict("path" => serialize(path),
             "geometries" => params)
    end
end

function serialize(geomdata::GeometryData)
    params = serialize(geomdata.geometry)
    params["color"] = serialize(geomdata.color)
    transform = compose(geomdata.transform,
                        intrinsic_transform(geomdata.geometry))
    if transform != IdentityTransformation()
        params["transform"] = serialize(transform)
    end
    params
end

intrinsic_transform(g) = IdentityTransformation()
intrinsic_transform(geomdata::GeometryData) = intrinsic_transform(geomdata.geometry)
intrinsic_transform(g::HyperRectangle) = Translation(center(g)...)
intrinsic_transform(g::HyperSphere) = Translation(center(g)...)
intrinsic_transform(g::HyperEllipsoid) = Translation(center(g)...)
intrinsic_transform(g::HyperCube) = Translation(center(g)...)
function intrinsic_transform(g::Cylinder{3})
    # DrakeVisualizer wants a cylinder to lie along the z axis
    R = rotation_between(SVector(0, 0, 1), g.extremity)
    Translation(center(g)) ∘ LinearMap(R)
end

serialize(color::Colorant) = (red(color),
                              green(color),
                              blue(color),
                              alpha(color))
serialize(p::Path) = string.(p)
serialize(v::Vector) = v
serialize(v::Vec) = convert(Vector, v)
serialize(v::Point) = convert(Vector, v)
serialize(v::StaticArray) = convert(Vector, v)
serialize(face::Face{N, T}) where {N, T} =
  raw.(convert(Face{N, GeometryTypes.OffsetInteger{-1, Int}}, face))
serialize(g::HyperRectangle) = Dict("type" => "box", "lengths" => serialize(widths(g)))
serialize(g::HyperSphere) = Dict("type" => "sphere", "radius" => radius(g))
serialize(g::HyperEllipsoid) = Dict("type" => "ellipsoid", "radii" => serialize(radii(g)))
serialize(g::Cylinder{3}) = Dict("type" => "cylinder",
                                      "length" => norm(g.extremity),
                                      "radius" => radius(g))
serialize(g::HyperCube) = Dict("type" => "box", "lengths" => widths(g))
serialize(g::GeometryPrimitive) = serialize(GLNormalMesh(g))
serialize(g::Triad) = Dict{String, Any}("type" => "triad",
                                        "scale" => g.scale,
                                        "tube" => g.tube)

function serialize(g::AbstractMesh)
    Dict("type" => "mesh_data",
         "vertices" => serialize.(vertices(g)),
         "faces" => serialize.(faces(g)))
end

function serialize(f::MeshFile)
    Dict("type" => "mesh_file",
         "filename" => f.filename,
         "scale" => SVector(1.0, 1.0, 1.0))
end

function serialize(g::PointCloud)
    params = Dict("type" => "pointcloud",
                  "points" => serialize.(g.points),
                  "channels" => Dict{String, Any}())
    for (channel, values) in g.channels
        params["channels"][string(channel)] = serialize.(values)
    end
    params
end

function serialize(g::PolyLine)
    params = Dict("type" => "line",
                  "points" => serialize.(g.points),
                  "radius" => g.radius,
                  "closed" => g.closed
        )
    if g.start_head !== nothing
        params["start_head"] = true
        params["head_radius"] = g.start_head.radius
        params["head_length"] = g.start_head.length
    end
    if g.end_head !== nothing
        params["end_head"] = true
        params["head_radius"] = g.end_head.radius
        params["head_length"] = g.end_head.length
    end
    params
end

function serialize(tform::Transformation)
    Dict{String, Vector{Float64}}("translation" => translation(tform),
                      "quaternion" => quaternion(tform))
end

quaternion(::IdentityTransformation) = SVector(1., 0, 0, 0)
quaternion(tform::AbstractAffineMap) = quaternion(transform_deriv(tform, SVector(0., 0, 0)))
quaternion(matrix::UniformScaling) = quaternion(IdentityTransformation())
quaternion(matrix::AbstractMatrix) = quaternion(Quat(matrix))
quaternion(quat::Quat) = SVector(quat.w, quat.x, quat.y, quat.z)

translation(tform::Transformation) = tform(SVector(0., 0, 0))


text:
using fci
using Test

@testset "fci.jl" begin
    include("test_compute_ss_terms_full.jl")
    include("test_precompute_spin_diag_terms.jl")
    include("test_excit_config.jl")
    include("test_get_all_configs.jl")
    include("test_get_index.jl")
    include("test_get_sigma.jl")
    include("test_make_lookup_table.jl")
    include("test_make_xy.jl")
end

max_stars_repo_name: nbraunsc/fci

max_issues_repo_name: nbraunsc/fci

max_forks_repo_name: nbraunsc/fci


text:
@testset "periodic" begin
    @testset "_periodic" begin
        @testset "$f" for f in (sin, cos)
            @testset "Typical 24-hour day" begin
                day_typical = ZonedDateTime(2015, 6, 1, tz"America/Winnipeg")
                hours_typical = day_typical .+ collect(Hour.(0:2))
                result = _periodic.(f, hours_typical, Day(1))
                expected = Dict(sin => [0., 0.2588, 0.5], cos => [1.0, 0.9659, 0.866])
                @test result ≈ expected[f] atol=1e-4
            end

            @testset "Spring daylight-saving time change (23 hours)" begin
                day_spring_dst = ZonedDateTime(2015, 3, 8, tz"America/Winnipeg")
                hours_spring_dst = day_spring_dst .+ collect(Hour.(0:2))
                result = _periodic.(f, hours_spring_dst, Day(1))
                expected = Dict(sin => [0., 0.2698, 0.5196], cos => [1., 0.9629, 0.8544])
                @test result ≈ expected[f] atol=1e-4
            end

            @testset "Fall daylight-saving time change (25 hours)" begin
                day_fall_dst = ZonedDateTime(2015, 11, 1, tz"America/Winnipeg")
                hours_fall_dst = day_fall_dst .+ collect(Hour.(0:2))
                result = _periodic.(f, hours_fall_dst, Day(1))
                expected = Dict(sin => [0., 0.2487, 0.4818], cos => [1., 0.9686, 0.8763])
                @test result ≈ expected[f] atol=1e-4
            end
        end

        @testset "phase shift" begin
            @test _periodic(sin, DateTime(2018), Week(1)) == 0
            @test _periodic(sin, DateTime(2018), Week(1), Day(7)) == 0
            @test _periodic(sin, DateTime(2018), Week(1), Day(-7)) == 0

            @test _periodic(sin, DateTime(2019), Week(1)) ≈ 0.7818 atol=1e-4
            @test _periodic(sin, DateTime(2019), Week(1), Day(7)) ≈ 0.7818 atol=1e-4
            @test _periodic(sin, DateTime(2019), Week(1), Day(5)) ≈ 0.4339 atol=1e-4

            # Using the phase shift we can make a Thursday the beginning of the week
            # (i.e. such that sin(t) == 0). The UNIX epoch 1970-01-01 fell on Thursday.
            # NOTE: these calculations can change depending on the period type given.
            @test _periodic(sin, DateTime(1970), Week(1), Day(3)) == 0
            @test _periodic(sin, DateTime(1970), Second(Week(1)), Day(5)) == 0

            # Phase shifts larger than the period should not change the result
            @test ==(
                _periodic(sin, DateTime(1970), Day(1), Day(0)),
                _periodic(sin, DateTime(1970), Day(1), Day(5)),
            )
        end
    end

    @testset "Real input" begin
        @testset "Constructors" begin
            @testset "Default" begin
                p = Periodic(sin, 5, 2)
                @test p isa Transform
                @test p.period == 5
                @test p.phase_shift == 2
            end

            @testset "No phase_shift" begin
                p = Periodic(sin, 5)
                @test p isa Transform
                @test p.period == 5
                @test p.phase_shift == 0
            end

            @testset "Non-positive period" for period in (0, -1)
                @test_throws DomainError Periodic(sin, period, 1)
                @test_throws DomainError Periodic(sin, period)
            end

            @testset "Different field subtypes" begin
                p = Periodic(sin, 4.0, 2)
                x = [2.0, 3.0, 5.0]

                @test FeatureTransforms.apply(x, p) ≈ [0.0, 1.0, -1.0] atol=1e-14
                @test p(x) ≈ [0.0, 1.0, -1.0] atol=1e-14
            end
        end
    end

    @testset "TimeType input" begin
        @testset "Constructors" begin
            @testset "Default" begin
                p = Periodic(sin, Day(5), Day(2))
                @test p isa Transform
                @test p.period == Day(5)
                @test p.phase_shift == Day(2)
                @test cardinality(p) == OneToOne()
            end

            @testset "No phase_shift" begin
                p = Periodic(sin, Day(5))
                @test p isa Transform
                @test p.period == Day(5)
                @test p.phase_shift == Day(0)
            end

            @testset "Non-positive period" for period in (Day(0), Day(-1))
                @test_throws DomainError Periodic(sin, period, Day(1))
                @test_throws DomainError Periodic(sin, period)
            end

            @testset "Different field subtypes" begin
                p = Periodic(sin, Week(1), Day(5))

                x = ZonedDateTime(2020, 1, 1, tz"EST") .+ (Day(0):Day(1):Day(5))
                # Use _periodic to get expected outputs because we test it elsewhere
                expected = _periodic.(sin, x, Week(1), Day(5))

                @test FeatureTransforms.apply(x, p) ≈ expected atol=1e-14
                @test p(x) ≈ expected atol=1e-14
            end
        end

        @testset "$f" for f in (sin, cos)
            p = Periodic(f, Day(5), Day(2))

            @testset "Basic" begin
                x = ZonedDateTime(2020, 1, 1, tz"EST") .+ (Day(0):Day(1):Day(5))
                # Use _periodic to get expected outputs because we test it elsewhere
                expected = _periodic.(f, x, Day(5), Day(2))

                @test FeatureTransforms.apply(x, p) ≈ expected atol=1e-14
                @test p(x) ≈ expected atol=1e-14
            end
        end
    end

    @testset "Type mismatch" begin
        @testset "Real data, Period transform" begin
            p = Periodic(sin, Day(5), Day(2))
            x = 0.:11.
            @test_throws MethodError FeatureTransforms.apply(x, p)
        end

        @testset "TimeType data, Real transform" begin
            p = Periodic(sin, 5, 2)
            x = ZonedDateTime(2020, 1, 1, tz"EST") .+ (Day(0):Day(1):Day(5))
            @test_throws MethodError FeatureTransforms.apply(x, p)
        end

        @test_throws ArgumentError Periodic(sin, 1, Day(3))
        @test_throws ArgumentError Periodic(sin, Day(1), 3)
    end
end


text:
module TurbulenceConvection

import ClimaCore
import LinearAlgebra
import DocStringExtensions
import StaticArrays
import StatsBase
import NCDatasets
import JSON
import Dierckx
import RootSolvers
import Statistics
import LambertW
import Thermodynamics
import Distributions
import FastGaussQuadrature
import CLIMAParameters
import OrdinaryDiffEq
import CloudMicrophysics
import SurfaceFluxes
import UnPack

const ODE = OrdinaryDiffEq
const CC = ClimaCore
const CCO = CC.Operators
const NC = NCDatasets
const SA = StaticArrays

const TD = Thermodynamics
const SF = SurfaceFluxes
const UF = SF.UniversalFunctions
const RS = RootSolvers

const CM = CloudMicrophysics
const CM0 = CloudMicrophysics.Microphysics_0M
const CM1 = CloudMicrophysics.Microphysics_1M
const liq_type = CM1.LiquidType()
const ice_type = CM1.IceType()
const rain_type = CM1.RainType()
const snow_type = CM1.SnowType()

const CP = CLIMAParameters
const CPP = CP.Planet
const CPSGS = CP.SubgridScale
const APS = CP.AbstractEarthParameterSet

const CPMP = CP.Atmos.Microphysics
const CPEDMF = CP.Atmos.EDMF

# For dispatching to inherited class
struct BaseCase end

up_sum(vals::AbstractArray) = reshape(sum(vals; dims = 1), size(vals, 2))

function parse_namelist(namelist, keys...; default = nothing, valid_options = nothing)
    param = namelist
    for k in keys
        if haskey(param, k)
            param = param[k]
            if valid_options ≠ nothing && !(param isa Dict)
                @assert param in valid_options
            end
        else
            if default == nothing
                error("No default value given for parameter ($(join(keys, ", "))).")
            else
                @info "Using default value, $default, for parameter ($(join(keys, ", ")))."
                return default
            end
        end
    end
    return param
end

include("ClimaParams.jl")
import .ClimaParams
const ICP = ClimaParams # internal clima parameters

#=
    debug_state(state, code_location::String)

A simple function for debugging the entire state,
specifically for when quantities that should remain
positive-definite become negative.

=#
function debug_state(state, code_location::String)
    prog_gm = center_prog_grid_mean(state)
    aux_gm = center_aux_grid_mean(state)
    prog_gm_f = face_prog_grid_mean(state)
    aux_gm_f = face_aux_grid_mean(state)

    prog_up = center_prog_updrafts(state)
    aux_up = center_aux_updrafts(state)
    prog_up_f = face_prog_updrafts(state)
    aux_up_f = face_aux_updrafts(state)

    prog_en = center_prog_environment(state)
    aux_en = center_aux_environment(state)
    aux_en_f = face_aux_environment(state)

    ######
    ###### Positive-definite variables
    ######

    vars_positive = [
        vec(prog_gm.θ_liq_ice),
        vec(prog_gm_f.w),
        vec(prog_up[1].ρarea),
        vec(prog_up[1].ρaθ_liq_ice),
        vec(prog_up_f[1].ρaw),
        vec(aux_en.area),
        vec(aux_en.θ_liq_ice),
    ]
    vars = vars_positive
    vars_conds = map(v -> any(v .< 0), vars)

    if any(vars_conds)
        @show code_location
        for (i, vc, v) in zip(1:length(vars), vars_conds, vars)
            vc || continue
            @show i, v
        end
        @show vars_conds
        error("Negative state for positive-definite field(s)")
    end

    ######
    ###### All listed variables
    ######
    vars = vars_positive
    vars_conds = map(v -> any(isnan.(v)) || any(isinf.(v)), vars)

    if any(vars_conds)
        @show code_location
        for (i, vc, v) in zip(1:length(vars), vars_conds, vars)
            vc || continue
            @show i, v
        end
        @show vars_conds
        error("Nan/Inf state for field(s)")
    end
end

include("parameters.jl")

include("Grid.jl")
include("dycore_api.jl")
include("NetCDFIO.jl")
include("diagnostics.jl")
include("Fields.jl")
include("ReferenceState.jl")
include("types.jl")
include("name_aliases.jl")
include("Operators.jl")

include("microphysics_coupling.jl")
include("turbulence_functions.jl")
include("utility_functions.jl")
include("TimeStepping.jl")
include("variables.jl")
include("EDMF_Precipitation.jl")
include("EDMF_Environment.jl")
include("EDMF_Updrafts.jl")
include("stochastic_closures.jl")
include("update_aux.jl")
include("Turbulence_PrognosticTKE.jl")
include("Forcing.jl")
include("Radiation.jl")
include("Surface.jl")
include("io.jl")
include("thermodynamics.jl")
include("surface_functions.jl")
include("closures/perturbation_pressure.jl")
include("closures/entr_detr.jl")
include("closures/nondimensional_exchange_functions.jl")
include("closures/mixing_length.jl")
include("closures/buoyancy_gradients.jl")

end


text:
#======
A. HAC
=======#

covariance(X, k::HAC, returntype::Type{T1} = CovarianceMatrix, factortype::Type{T2} = Cholesky; kwargs...) where {T1, T2}= covariance(X, k, HACCache(X, k), returntype, factortype; kwargs...)

function covariance(X::T,
                    k::HAC{F},
                    cache::HACCache = HACCache(X, k),
                    returntype::Type{T1} = CovarianceMatrix,
                    factortype::Type{T2} = Cholesky;
                    demean::Bool = true,
                    scale::Int = size(X,1)) where {T<:AbstractMatrix,
                                                   F,
                                                   T1<:Union{CovarianceMatrix, Matrix},
                                                   T2<:Union{Nothing, Factorization}}
    check_cache_consistenty(k, cache)
    demean!(cache, X, Val{demean})
    prewhiten!(cache)
    _covariance!(cache, k)
    finalize(T2, T1, k, cache, scale)
end

#=--
Implementation
--=#
function _covariance!(cache::HACCache, k::HAC{Optimal{T}}) where {T<:OptimalBandwidth}
    n, p = size(cache.XX)
    setupkernelweights!(k, p, eltype(cache.XX))
    optimal_bw!(cache, k, T())
    __covariance!(cache, k)
end

_covariance!(cache::HACCache, k::HAC{T}) where {T<:Fixed} = __covariance!(cache, k)

function __covariance!(cache::HACCache, k::HAC)
    n, p = size(cache.XX)
    fill!(cache.V, zero(eltype(cache.XX)))
    bw = first(k.bw)
    mul!(cache.V, cache.XX', cache.XX)
    triu!(cache.V)
    idxs = getcovindeces(k, n)
    @inbounds for j in idxs
        k_j = CovarianceMatrices.kernel(k, j/bw)
        LinearAlgebra.axpy!(k_j, CovarianceMatrices.Γ!(cache, j), cache.V)
    end
    LinearAlgebra.copytri!(cache.V, 'U')
    swhiten!(cache)
    rmul!(cache.V, 1/(n+isprewhiten(k)))
    nothing
end

#=======
B. HC
=======#
function covariance(X::T, k::K, returntype::Type{T1} = CovarianceMatrix, factortype::Type{T2} = Cholesky; demean::Bool = true, scale::Int = size(X,1)) where {T<:AbstractMatrix, K<:HC, T1<:Union{CovarianceMatrix, Matrix}, T2<:Union{Nothing, Factorization}}
    V = X'X
    finalize(V, T2, T1, k, scale)
end

function covariance(X::T, k::K, cache::HCCache, factortype::Type{T1} = Cholesky, returntype::Type{T2} = CovarianceMatrix; demean::Bool = true, scale::Int = size(X, 1)) where {T<:AbstractMatrix, K<:HC, T1<:Union{Nothing, Factorization}, T2<:Union{CovarianceMatrix, Matrix}}
    check_cache_consistenty(k, cache)
    _covariance!(cache, k)
    finalize(T2, T1, k, cache, scale)
end

#======
CRHC
=======#
function covariance(X::T, k::K, returntype::Type{T1} = CovarianceMatrix, factortype::Type{T2} = Cholesky; demean::Bool = true, scale::Int = size(X,1), sorted::Bool = false) where {T<:AbstractMatrix, K<:CRHC, T1<:Union{CovarianceMatrix, Matrix}, T2<:Union{Nothing, Factorization}}
    cache = CRHCCache(X, k.cl)
    covariance(X, k, cache, factortype, returntype, demean = demean, scale = scale, sorted = sorted)
end

function covariance(X::T, k::K, cache::CRHCCache, returntype::Type{T1} = CovarianceMatrix, factortype::Type{T2} = Cholesky; demean::Bool = true, scale::Int = size(X, 1), sorted::Bool = false ) where {T<:AbstractMatrix, K<:CRHC, T1<:Union{CovarianceMatrix, Matrix}, T2<:Union{Nothing, Factorization}}
    #check_cache_consistenty(k, cache)
    _covariance!(cache, X, k, sorted)
    finalize(cache.M, T2, T1, k, scale)
end

#=---
Implementation
---=#

function _covariance!(cache, X, k::CRHC, sorted::Bool)
    CovarianceMatrices.installsortedX!(cache, X, k, Val{sorted})
    bstarts = (searchsorted(cache.clus, j[2]) for j in enumerate(unique(cache.clus)))
    CovarianceMatrices.clusterize!(cache, bstarts)
    rmul!(cache.M, dof_adjustment(cache, k, bstarts))

end

function installsortedX!(cache, X, k, ::Type{Val{true}})
    copyto!(cache.q, X)
    copyto!(cache.clus, k.cl)
    nothing
end

function installsortedX!(cache, X, k, ::Type{Val{false}})
    n, p = size(cache.X)
    sortperm!(cache.clusidx, k.cl)
    cidx = cache.clusidx
    c  = k.cl
    cc = cache.clus
    XX = cache.q
    @inbounds for j in 1:p, i in eachindex(cache.clusidx)
        XX[i,j] = X[cidx[i], j]
    end
    @inbounds for i in eachindex(cache.clusidx)
        cc[i] = c[cidx[i]]
    end
    nothing
end

#======
Finalizers
======#
finalize(V, ::Type{T}, ::Type{Matrix}, k, scale::Real) where T = rmul!(V, 1/scale)
finalize(::Type{T}, ::Type{Matrix}, k, cache, scale::Real) where T = rmul!(cache.V, 1/scale)

function finalize(V, ::Type{SVD}, ::Type{CovarianceMatrix}, k, scale::Real)
    rmul!(V, 1/scale)
    CovarianceMatrix(svd(V), k, V)
end

function finalize(V, ::Type{Cholesky}, ::Type{CovarianceMatrix},  k, scale::Real)
    rmul!(V, 1/scale)
    CovarianceMatrix(cholesky(Hermitian(V)), k, V)
end

function finalize(::Type{SVD}, ::Type{CovarianceMatrix}, k, cache::AbstractCache, scale::Real)
    rmul!(cache.V, 1/scale)
    CovarianceMatrix(svd(cache.V), k, copy(cache.V))
end

function finalize(::Type{Cholesky}, ::Type{CovarianceMatrix},  k, cache::HACCache, scale::Real)
    rmul!(cache.V, 1/scale)
    CovarianceMatrix(cholesky(Hermitian(cache.V)), k, copy(cache.V))
end


text:
# using Printf

# for d in 0:max_diff_resolvable[]+100
#     tests_failed = 0
#     mean_diff = 0
#     max_diff = 0
#     mean_di = 0f0
#     max_di = 0f0
#     for lf1ui in 0x0000:Int16(2^15-d-1)
#         lf2ui = Int16(lf1ui + d)

#         lf1 = reinterpret(LogFixPoint16,lf1ui)
#         lf2 = reinterpret(LogFixPoint16,lf2ui)

#         r = Float32(lf1)+Float32(lf2)
#         result = LogFixPoint16(r)
#         result2 = lf1+lf2

#         ui1 = reinterpret(UInt16,result)
#         ui2 = reinterpret(UInt16,result2)

#         # filter out nans
#         failed = (ui1 == 0x8000 || ui2 == 0x8000) ? false : ui1 != ui2
#         tests_failed += failed

#         if failed
#             mean_diff += Int(ui1)-Int(ui2)
#             max_diff = max(abs(Int(ui1)-Int(ui2)),max_diff)
#             di = Float32.([prevfloat(result2),result2,nextfloat(result2)]) .- r
#             di = sort(abs.(di))[1:2]
#             dir = di[2]/sum(di)
#             mean_di += dir
#             max_di = max(max_di,dir)
#             # if dir > 0.7
#             #     println([lf1,lf2])
#             # end
#         end
#     end
#     if tests_failed > 0
#         ds = @sprintf("%4d",d)
#         ps = @sprintf("%5d",2^15-d-tests_failed)
#         fs = @sprintf("%5d",tests_failed)
#         ms = @sprintf("%2d",Int(mean_diff/tests_failed))
#         dis = @sprintf("%.6f",mean_di/tests_failed)
#         mdis = @sprintf("%.6f",max_di)
#         println("d=$ds, Pass: $ps, Fail: $fs, Diff: $ms, Mean dist: $dis, Max dist: $mdis")
#     end
# end

# # d=   1, Pass: 22623, Fail: 10144, Max diff: 1
# # d= 267, Pass: 32244, Fail:   257, Max diff: 1
# # d= 431, Pass: 29175, Fail:  3162, Max diff: 1
# # d= 512, Pass: 30577, Fail:  1679, Max diff: 1
# # d= 780, Pass: 31756, Fail:   232, Max diff: 1
# # d=1407, Pass: 30983, Fail:   378, Max diff: 1
# # d=1898, Pass: 16955, Fail: 13915, Max diff: 1
# # d=2364, Pass: 27703, Fail:  2701, Max diff: 1
# # d=2390, Pass: 23380, Fail:  6998, Max diff: 1
# # d=2535, Pass: 18473, Fail: 11760, Max diff: 1
# # d=3211, Pass: 25052, Fail:  4505, Max diff: 1
# # d=3493, Pass: 29253, Fail:    22, Max diff: 1
# # d=4201, Pass: 28545, Fail:    22, Max diff: 1
# # d=4578, Pass: 27461, Fail:   729, Max diff: 1
# # d=4579, Pass: 28188, Fail:     1, Max diff: 1
# # d=5389, Pass: 26733, Fail:   646, Max diff: 1
# # d=5390, Pass: 21643, Fail:  5735, Max diff: 1
# # d=5391, Pass: 23039, Fail:  4338, Max diff: 1
# # d=5392, Pass: 27367, Fail:     9, Max diff: 1

# for d in 0:max_diff_resolvable[]+100
#     tests_failed = 0
#     mean_diff = 0
#     max_diff = 0
#     mean_di = 0f0
#     max_di = 0f0
#     for lf1ui in 0x0000:Int16(2^15-d-1)
#         lf2ui = Int16(lf1ui + d)

#         lf1 = reinterpret(LogFixPoint16,lf1ui)
#         lf2 = reinterpret(LogFixPoint16,lf2ui)

#         r = Float32(lf2)-Float32(lf1)
#         result = LogFixPoint16(r)
#         result2 = lf2-lf1

#         ui1 = reinterpret(UInt16,result)
#         ui2 = reinterpret(UInt16,result2)

#         # filter out nans
#         failed = (ui1 == 0x8000 || ui2 == 0x8000) ? false : ui1 != ui2
#         tests_failed += failed

#         if failed
#             mean_diff += Int(ui1)-Int(ui2)
#             max_diff = max(abs(Int(ui1)-Int(ui2)),max_diff)
#             di = Float32.([prevfloat(result2),result2,nextfloat(result2)]) .- r
#             di = sort(abs.(di))[1:2]
#             dir = di[2]/sum(di)
#             mean_di += dir
#             max_di = max(max_di,dir)
#             # if dir > 0.7
#             #     println([lf1,lf2])
#             # end
#         end
#     end
#     if tests_failed > 0
#         ds = @sprintf("%4d",d)
#         ps = @sprintf("%5d",2^15-d-tests_failed)
#         fs = @sprintf("%5d",tests_failed)
#         ms = @sprintf("%2d",max_diff)
#         dis = @sprintf("%.6f",mean_di/tests_failed)
#         mdis = @sprintf("%.6f",max_di)
#         println("d=$ds, Pass: $ps, Fail: $fs, Diff: $ms, Mean dist: $dis, Max dist: $mdis")
#     end
# end

# er = fill(0f0,11760)
# for d in [2535]
#     i = 1
#     for lf1ui in 0x0000:Int16(2^15-d-1)
#         lf2ui = Int16(lf1ui + d)

#         lf1 = reinterpret(LogFixPoint16,lf1ui)
#         lf2 = reinterpret(LogFixPoint16,lf2ui)

#         result = Float32(lf1)+Float32(lf2)
#         if reinterpret(UInt16,LogFixPoint16(result)) != reinterpret(UInt16,lf1+lf2)
#             di = Float32.([lf1+lf2,nextfloat(lf1+lf2)]) .- result
#             er[i] = abs(di[1])/sum(abs.(di))
#             i += 1
#         end
#     end
# end

# for d in 0:10000
#     tests_failed = 0
#     max_diff = 0
#         for lf1ui in 0x0000:Int16(2^15-d-1)
#             lf2ui = Int16(lf1ui + d)

#             lf1 = reinterpret(LogFixPoint16,lf1ui)
#             lf2 = reinterpret(LogFixPoint16,lf2ui)

#             result = LogFixPoint16(Float32(lf2)-Float32(lf1))
#             tests_failed += reinterpret(UInt16,result) != reinterpret(UInt16,lf2-lf1)
#             max_diff = max(max_diff,abs(Int(reinterpret(UInt16,result))-Int(reinterpret(UInt16,lf2-lf1))))
#         end
#         if tests_failed > 0
#             ds = @sprintf("%4d",d)
#             ps = @sprintf("%5d",2^15-d-tests_failed)
#             fs = @sprintf("%5d",tests_failed)
#             println("d=$ds, Pass: $ps, Fail: $fs, Max diff: $max_diff")
#         end
#     end
# end

# # @testset "addition all" begin
# R = fill(Int32(0),1000,6)
# i = 1
# for _ in 1:10000000
#     lf1ui,lf2ui = rand(UInt16,2)

#     # abs to only test the addTable
#     lf1 = abs(reinterpret(LogFixPoint16,lf1ui))
#     lf2 = abs(reinterpret(LogFixPoint16,lf2ui))

#     result = LogFixPoint16(Float32(lf1)+Float32(lf2))
#     if reinterpret(UInt16,result) != reinterpret(UInt16,lf1+lf2)
#         R[i,1:end-1] = vcat([diff_val(lf1,lf2)],reinterpret.(UInt16,[lf1,lf2,result,lf1+lf2]))
#         R[i,end] = Int(reinterpret(UInt16,result)) - Int(reinterpret(UInt16,lf1+lf2))
#         global i += 1
#         global i = min(i,1000)
#     end
# end
# # end


text:
is_ch_type(::Val{:Enum8})  = true
is_ch_type(::Val{:Enum16})  = true
result_type(::Val{:Enum8}, args...)  = CategoricalVector{String}
result_type(::Val{:Enum16}, args...)  = CategoricalVector{String}

const ENUM_RE_ARG = r"""

           '((?:(?:[^'])|(?:\\'))+)'
           \s*=\s*
           (-?\d+)
           \s*$
       """x

function make_enum_map(::Type{BaseT}, args...) where {BaseT}
    map = Dict{String, BaseT}()
    for arg in args
        m = match(ENUM_RE_ARG, arg)
        isnothing(m) && error("Wrong enum argument $arg")
        map[m.captures[1]] = parse(BaseT, m.captures[2])
    end
    return map
end

function read_enum_data(sock::ClickHouseSock, num_rows::VarUInt,
                                        ::Type{BaseT}, args...) where {BaseT}
    map = make_enum_map(BaseT, args...)
    levels = collect(keys(map))
    enum_to_level = Dict(zip(values(map), 1:length(map)))
    data = chread(sock, Vector{BaseT}, num_rows)
    result = CategoricalVector{String, BaseT}(undef, 0, levels = levels)
    data .= getindex.(Ref(enum_to_level), data)
    result.refs = data
    return result
end

function write_enum_data(sock::ClickHouseSock, data::AbstractVector{String},
                                        ::Type{BaseT}, args...) where {BaseT}
    map = make_enum_map(BaseT, args...)
    d = Vector{BaseT}(undef, length(data))
    try
        d .= getindex.(Ref(map), data)
    catch exc
        if exc isa KeyError
            error("Value is not a valid enum variant: $(exc.key)")
        end
        rethrow()
    end
    chwrite(sock, d)
end

function write_enum_data(sock::ClickHouseSock, data::CategoricalVector{String},
                                        ::Type{BaseT}, args...) where {BaseT}
    map = make_enum_map(BaseT, args...)
    try

        level_to_enum = getindex.(Ref(map), levels(data))
        d = convert(Vector{BaseT}, getindex.(Ref(level_to_enum), data.refs))
        chwrite(sock, d)

    catch exc
        if exc isa KeyError
            error("Value is not a valid enum variant: $(exc.key)")
        end
        rethrow()
    end

end

read_col_data(sock::ClickHouseSock, num_rows::VarUInt, ::Val{:Enum8}, args...) =
                     read_enum_data(sock, num_rows, Int8, args...)
read_col_data(sock::ClickHouseSock, num_rows::VarUInt, ::Val{:Enum16}, args...) =
                    read_enum_data(sock, num_rows, Int16, args...)

write_col_data(sock::ClickHouseSock,
        data::AbstractVector{String},
         ::Val{:Enum8},
         args...) = write_enum_data(sock, data, Int8, args...)

write_col_data(sock::ClickHouseSock,
                data::AbstractVector{String},
                ::Val{:Enum16},
                 args...) = write_enum_data(sock, data, Int16, args...)

write_col_data(sock::ClickHouseSock,
             data::CategoricalVector{String},
             ::Val{:Enum8},
             args...) = write_enum_data(sock, data, Int8, args...)

write_col_data(sock::ClickHouseSock,
            data::CategoricalVector{String},
            ::Val{:Enum16}, args...) =
            write_enum_data(sock, data, Int16, args...)


text:
input = open("../inputs/day05.txt") do f
    read(f, String)
end

function range(s, low, high)
    lo = low
    hi = high
    for c in s
        if c == 'F' || c == 'L'
            hi = lo + (hi - lo) ÷ 2
        elseif c == 'B' || c == 'R'
            lo = lo + (hi - lo + 1) ÷ 2
        end
    end
    @assert (lo == hi)
    lo
end

function seat_id(s)
    row = range(s[1:7], 0, 127)
    col = range(s[8:10], 0, 7)
    row * 8 + col
end

all_seats = [ seat_id(ln) for ln in split(input, "\n")]
sort!(all_seats)

println("Part 1: ", max(all_seats...))

my_seat_id = first(
    [ all_seats[i] - 1
    for i in 2:length(all_seats)
    if all_seats[i - 1] + 1 != all_seats[i]])

println("Part 2: ", my_seat_id)


text:
# ======================================================================
# Description
# ========================================================================
# Calculates MOID between two keplerian orbits in AU
#
# Created: 6/16/20
# Original Author : Anivid Pedros Faura, anivid.pedrosfaura@colorado.edu
# Translated from .m to .jl by: Luke Bury, luke.bury@colorado.edu

# ========================================================================
# Libraries, Functions, and Necessary Paths
# ========================================================================
# # Libraries
# using DifferentialEquations # For using ODEProblem
# using LinearAlgebra # Allows for norm(), cross(), dot(), etc
# using Plots
using Printf # allows use of println(@sprintf())

# ========================================================================
# Dependent functions
# ========================================================================
function RefFrame(A,B)
    ### Reference frame rotation
    # 3-1-3 Tranformation (orbital elements)
    theta1 = A["Omega"]
    theta2 = A["i"]
    theta3 = A["argp"]
    RA = zeros(3,3) #from reference to A frame
    RA[1,1] = cos(theta3)*cos(theta1)-sin(theta3)*sin(theta1)*cos(theta2)
    RA[1,2] = sin(theta1)*cos(theta3)+cos(theta1)*cos(theta2)*sin(theta3)
    RA[1,3] = sin(theta2)*sin(theta3)
    RA[2,1] = -cos(theta1)*sin(theta3) - sin(theta1)*cos(theta2)*cos(theta3)
    RA[2,2] = -sin(theta3)*sin(theta1) + cos(theta3)*cos(theta2)*cos(theta1)
    RA[2,3] = sin(theta2)*cos(theta3)
    RA[3,1] = sin(theta2)*sin(theta1)
    RA[3,2] = -sin(theta2)*cos(theta1)
    RA[3,3] = cos(theta2)

    ### Preparing the orbit
    x = [cos(B["Omega"])*cos(B["argp"]) - sin(B["argp"])*cos(B["i"])*sin(B["Omega"])
         sin(B["Omega"])*cos(B["argp"]) + sin(B["argp"])*cos(B["i"])*cos(B["Omega"])
         sin(B["i"])*sin(B["argp"])]

    y = [-cos(B["Omega"])*sin(B["argp"]) - cos(B["argp"])*cos(B["i"])*sin(B["Omega"])
         -sin(B["Omega"])*sin(B["argp"]) + cos(B["argp"])*cos(B["i"])*cos(B["Omega"])
         sin(B["i"])*cos(B["argp"])];

    z = [sin(B["i"])*sin(B["Omega"])
        -sin(B["i"])*cos(B["Omega"])
        cos(B["i"])]

    xn = RA * x
    yn = RA * y
    zn = RA * z

    incliB = atan(sqrt(zn[1]^2 + zn[2]^2),zn[3])
    omegaB = -atan(zn[1],-zn[2])
    argpB = -atan(xn[3],yn[3])

    return incliB, omegaB, argpB
end

function WaterProcedure(incliB, omegaB, argpB, N, vtrueB, vL, vdis)
    if N<2
        vtrueB = zeros(4,1)
        vL = zeros(4,1)
        vdis = zeros(4,1)
        N = 4
        for j = 1:4
            vtrueB[j] = (0.25+0.5*j)*pi #evenly distributed 0.75 - 1.25 - 1.75 - 2.25
            a = (cos(omegaB)*cos(argpB+vtrueB[j]) - sin(omegaB)*sin(argpB+vtrueB[j])*cos(incliB))
            b = (sin(omegaB)*cos(argpB+vtrueB[j]) + cos(omegaB)*sin(argpB+vtrueB[j])*cos(incliB))
            vL[j] = atan(b,a)
            vdis[j] = 1e6; #set to something large
        end
    end

    return vtrueB, vL, vdis, N
end

# ========================================================================
# Calculating MOID
# ========================================================================
# -------------------------------------------------
# Setting up initial conditions
# -------------------------------------------------
# --------------------------
# Orbital parameters of body A
# --------------------------
A = Dict("sma" => 2.7691652, "e" => 0.0760091, "argp" => 73.59764*π/180, "Omega" => 80.30553*π/180, "i" => 10.59407*π/180)

# --------------------------
# Orbital parameters of body B
# --------------------------
B = Dict("sma" => 2.3655722, "e" => 0.127581, "argp" => 87.42605*π/180, "Omega" => 307.46872*π/180, "i" => 2.09575*π/180)

# -------------------------------------------------
# Calculating MOID
# -------------------------------------------------
# --------------------------
# New orbital elements: A frame
# --------------------------
# incliA = 0, omegaA = 0, argpA = 0
incliB, omegaB, argpB = RefFrame(A,B)

# --------------------------
# Scanning orbits
# --------------------------
# Scanning one full revolution of meridional plane to find local minima
# First guess for MOID is set to be big ~1e6AU

# Angle sweeping steps
cstep = 0.12 #rad - based on Wisnioski paper ~ 0.12rad
stepini = 0.07 #rad - for initial tuning step
stepfin = 1e-5 #rad - for final step of first tuning
stepmin = 1e-14 #rad - threshold step for secondtuning

# Initial Guess
trueB = -2 * cstep
moid = 1e6
dist_o = 1e6
vdis = ones(4,1)*1e6

for i in 1:2
    # First triplet
    rB = B["sma"] * (1-B["e"].^2) / (1+B["e"]*cos(trueB)) #compute the radius for B
    xB = rB * (cos(omegaB)*cos(argpB+trueB) - sin(omegaB)*sin(argpB+trueB)*cos(incliB))
    yB = rB * (sin(omegaB)*cos(argpB+trueB) + cos(omegaB)*sin(argpB+trueB)*cos(incliB))
    zB = rB * sin(argpB+trueB)*sin(incliB)

    rhoB = sqrt(xB^2+yB^2)
    L = atan(yB,xB)

    rA = A["sma"] * (1-A["e"].^2) / (1+A["e"]*cos(L)) #compute the radius for A
    rA2 = A["sma"] * (1-A["e"].^2) / (1-A["e"]*cos(L))

    if abs(rhoB-rA)>abs(rhoB+rA)
        rA = rA2
        L = L - pi
        diff = rhoB + rA2
    else
        diff = rhoB - rA
    end

    global D = zB^2+diff^2 # square of the distance
    # storing
    if i == 1
        global dist_oo = D
    else
        global dist_o = D
        global trueB_o = trueB
        global L_o = L
    end

    global trueB = trueB + cstep
end # for i in 1:2

# Full revolution
N = 0 # number of minima
dmin = D

vtrueB = fill(1,0).*0.0
vL = fill(1,0).*0.0
vdis = fill(1,0).*0.0

while trueB < (2*pi + cstep)
    rB = B["sma"] * (1-B["e"].^2) / (1+B["e"]*cos(trueB)) #compute the radius for B
    xB = rB * (cos(omegaB)*cos(argpB+trueB) - sin(omegaB)*sin(argpB+trueB)*cos(incliB))
    yB = rB * (sin(omegaB)*cos(argpB+trueB) + cos(omegaB)*sin(argpB+trueB)*cos(incliB))
    zB = rB * sin(argpB+trueB)*sin(incliB)

    rhoB = sqrt(xB^2+yB^2)
    L = atan(yB,xB)

    rA = A["sma"] * (1-A["e"].^2) / (1+A["e"]*cos(L)) #compute the radius for A
    rA2 = A["sma"] * (1-A["e"].^2) / (1-A["e"]*cos(L))

    if abs(rhoB-rA)>abs(rhoB+rA)
        rA = rA2
        L = L - pi
        diff = rhoB + rA2
    else
        diff = rhoB - rA
    end

    D = zB^2+diff^2 # square of the distance

    if (dist_o<=D) && (dist_o<=dist_oo)
        # global N = N + 1
        # vtrueB = trueB_o
        # vL = L_o
        # vdis = dist_o

        global vtrueB = append!(vtrueB,trueB_o)
        global vL = append!(vL, L_o)
        global vdis = append!(vdis, dist_o)
    end
    if dmin>D
        global dmin = D
    end

    global dist_oo = dist_o
    global trueB_o = trueB
    global L_o = L
    global dist_o = D
    global trueB = trueB + cstep
end

# -------------------------------------------------
# Water Procedure
# -------------------------------------------------
vtrueB, vL, vdis, N = WaterProcedure(incliB, omegaB, argpB, N, vtrueB, vL, vdis)

# -------------------------------------------------
# PARALLEL TUNING
# -------------------------------------------------
# Move objects separately along their orbits
# Smallest possible distance no longer meridional distance
rBt = fill(1,3).*NaN
rAt = fill(1,3).*NaN
xBt = fill(1,3).*NaN
yBt = fill(1,3).*NaN
zBt = fill(1,3).*NaN
xAt = fill(1,3).*NaN
yAt = fill(1,3).*NaN
k = 1
while k < N+2
    if k<=N
        global moid = vdis[k]
        global trueB_m = vtrueB[k]
        global L_m = vL[k]
        global myStep = stepini
        global threshold = stepfin #maybe here problem
    else
        if N == 2
            # in case of two minima are very close to each other(<1E-4 a.u.)
            # go to "water procedure"
            if (abs(vdis[1]-vdis[2])<1e-4)
                global N = 1
                global vtrueB, vL, vdis, N = WaterProcedure(incliB, omegaB, argpB, N)
                global k = 1
            else
                if (vdis[1]<moid)
                    global moid = vdis[1]
                    global trueB_m = vtrueB[1]
                    global L_m = vL[1]
                end
            end
        else
            # final tuning
            for i in 1:(N-1)
                if vdis[i]<moid
                    global moid = vdis[i]
                    global trueB_m = vtrueB[i]
                    global L_m = vL[i]
                end
            end
        end
        global myStep  = 2*stepini #inital state
        global threshold = stepmin #terminal state
    end

    rBt[2] = B["sma"] * (1-B["e"]^2) / (1+B["e"]*cos(trueB_m)) #compute the radius for B
    xBt[2] = rBt[2] * (cos(omegaB)*cos(argpB+trueB_m) - sin(omegaB)*sin(argpB+trueB_m)*cos(incliB))
    yBt[2] = rBt[2] * (sin(omegaB)*cos(argpB+trueB_m) + cos(omegaB)*sin(argpB+trueB_m)*cos(incliB))
    zBt[2] = rBt[2] * sin(argpB+trueB_m)*sin(incliB)

    rAt[2] = A["sma"] * (1-A["e"].^2) / (1+A["e"]*cos(L_m)) #compute the radius for A
    xAt[2] = rAt[2]*cos(L_m)
    yAt[2] = rAt[2]*sin(L_m)

    aleft = true
    aright = true
    bleft = true
    bright = true

    while (myStep>=threshold)
        lpoints = 0
        k1min   = 1
        k1max   = 3
        i1min   = 1
        i1max   = 3
        calc1   = false
        calc2   = false
        calc3   = false
        calc4   = false

        if bleft
            rBt[1] = B["sma"] * (1-B["e"]^2) / (1+B["e"]*cos(trueB_m-myStep)) #compute the radius for B
            xBt[1] = rBt[1] * (cos(omegaB)*cos(argpB+trueB_m-myStep) - sin(omegaB)*sin(argpB+trueB_m-myStep)*cos(incliB))
            yBt[1] = rBt[1] * (sin(omegaB)*cos(argpB+trueB_m-myStep) + cos(omegaB)*sin(argpB+trueB_m-myStep)*cos(incliB))
            zBt[1] = rBt[1] * sin(argpB+trueB_m-myStep)*sin(incliB)
            lpoints = lpoints + 1;
        end
        if bright
            rBt[3] = B["sma"] * (1-B["e"]^2) / (1+B["e"]*cos(trueB_m+myStep)) #compute the radius for B
            xBt[3] = rBt[3] * (cos(omegaB)*cos(argpB+trueB_m+myStep) - sin(omegaB)*sin(argpB+trueB_m+myStep)*cos(incliB))
            yBt[3] = rBt[3] * (sin(omegaB)*cos(argpB+trueB_m+myStep) + cos(omegaB)*sin(argpB+trueB_m+myStep)*cos(incliB))
            zBt[3] = rBt[3] * sin(argpB+trueB_m+myStep)*sin(incliB)
            lpoints = lpoints + 1;
        end
        if aleft
            rAt[1] = A["sma"] * (1-A["e"]^2) / (1+A["e"]*cos(L_m-myStep)) #compute the radius for A
            xAt[1] = rAt[1]*cos(L_m-myStep)
            yAt[1] = rAt[1]*sin(L_m-myStep)
            lpoints = lpoints + 1
        end
        if aright
            rAt[3] = A["sma"] * (1-A["e"]^2) / (1+A["e"]*cos(L_m+myStep)) #compute the radius for A
            xAt[3] = rAt[3]*cos(L_m+myStep)
            yAt[3] = rAt[3]*sin(L_m+myStep)
            lpoints = lpoints + 1
        end

        k1_t = 2
        i1_t = 2

        if lpoints == 1
            if aleft
                i1max = 1
            end
            if aright
                i1min = 3
            end
            if bright
                k1min = 3
            end
            if bleft
                k1max = 1
            end
        end

        if lpoints == 2
            if (aleft && bright)
                calc1 = true
            end
            if (aleft && bleft)
                calc2 = true
            end
            if (aright && bright)
                calc3 = true
            end
            if (aright && bleft)
                calc4 = true
            end
        end

        for k1 in k1min:k1max
            for i1 in i1min:i1max
                compute = true
                if lpoints == 2
                    if i1 != 1
                        if (k1 != 3 && calc1)||(k1 != 1 && calc2)
                            compute = false
                        end
                    end
                    if i1 != 3
                        if (k1 != 3 && calc3)||(k1 != 1 && calc4)
                            compute = false
                        end
                    end
                end
                if (k1 == 2 && i1 == 2)
                    compute = false
                end
                if (compute)
                    Dx = xBt[k1]-xAt[i1]
                    Dy = yBt[k1]-yAt[i1]
                    Dz = zBt[k1]
                    dist = Dx*Dx + Dy*Dy + Dz*Dz
                    if (dist<moid)
                        moid = dist
                        k1_t = k1
                        i1_t = i1
                    end
                end
            end
        end

        if (k1_t!=2) || (i1_t!=2)
            aleft=false
            aright=false
            bleft=false
            bright=false

            if i1_t!=2
                if i1_t == 1
                    aleft = true
                    L_m = L_m - myStep
                    rAt[3] = rAt[2]
                    xAt[3] = xAt[2]
                    yAt[3] = yAt[2]
                    rAt[2] = rAt[1]
                    xAt[2] = xAt[1]
                    yAt[2] = yAt[1]
                else
                    aright = true
                    L_m = L_m + myStep
                    rAt[1] = rAt[2]
                    xAt[1] = xAt[2]
                    yAt[1] = yAt[2]
                    rAt[2] = rAt[3]
                    xAt[2] = xAt[3]
                    yAt[2] = yAt[3]
                end
            end
            if k1_t!=2
                if k1_t == 1
                   bleft = true
                    trueB_m = trueB_m - myStep
                    rBt[3] = rBt[2]
                    xBt[3] = xBt[2]
                    yBt[3] = yBt[2]
                    zBt[3] = zBt[2]
                    rBt[2] = rBt[1]
                    xBt[2] = xBt[1]
                    yBt[2] = yBt[1]
                    zBt[2] = zBt[1]
                else
                    bright = true
                    trueB_m = trueB_m + myStep
                    rBt[1] = rBt[2]
                    xBt[1] = xBt[2]
                    yBt[1] = yBt[2]
                    zBt[1] = zBt[2]
                    rBt[2] = rBt[3]
                    xBt[2] = xBt[3]
                    yBt[2] = yBt[3]
                    zBt[2] = zBt[3]
                end
            end
        else
            aleft = true
            aright = true
            bleft = true
            bright = true
            myStep = myStep*0.15
        end
    end

    if k <= N
        vdis[k] = moid
        vtrueB[k] = trueB_m
        vL[k] = L_m
    end
    global k = k+1
end

moid = sqrt(moid)

println(@sprintf("\nThe MOID is: %1.16f\n", moid))

# ========================================================================
# Useful Formatting Structures
# ========================================================================
# -------------------------------------------------
#
# -------------------------------------------------

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

max_stars_repo_name: CCAR-ORCCA/MOID

max_issues_repo_name: CCAR-ORCCA/MOID

max_forks_repo_name: CCAR-ORCCA/MOID


text:
function DD.show_after(io::IO, mime::MIME"text/plain", A::AbstractRaster)

    if missingval(A) !== nothing
        printstyled(io, "with missingval: "; color=:light_black)
        print(io, string(missingval(A)), "\n")
    end
    if parent(A) isa DiskArrays.AbstractDiskArray
        if parent(A) isa FileArray
            printstyled(io, "\nfrom file:\n"; color=:light_black)
            print(io, filename(parent(A)))
        end
    else
        DD.print_array(io, mime, parent(A))
    end
end

function DD.show_after(io, mime, stack::AbstractRasterStack)
    if parent(stack) isa FileStack
        printstyled(io, "\nfrom file:\n"; color=:light_black)
        println(io, filename(stack))
    end
end

# Stack types can be enourmous. Just use nameof(T)
function Base.summary(io::IO, ser::AbstractRasterSeries{T,N}) where {T,N}
    if N == 0
        print(io, "0-dimensional ")
    elseif N == 1
        print(io, size(ser, 1), "-element ")
    else
        print(io, join(size(ser), "×"), " ")
    end
    print(io, string(nameof(typeof(ser)), "{$(nameof(T)),$N}"))
end

function Base.show(io::IO, mime::MIME"text/plain", A::AbstractRasterSeries{T,N}) where {T,N}
    lines = 0
    summary(io, A)
    DD.print_name(io, name(A))
    lines += Dimensions.print_dims(io, mime, dims(A))
    !(isempty(dims(A)) || isempty(refdims(A))) && println(io)
    lines += Dimensions.print_refdims(io, mime, refdims(A))
    println(io)
    ds = displaysize(io)
    ioctx = IOContext(io, :displaysize => (ds[1] - lines, ds[2]))
end

function Base.show(io::IO, mime::MIME"text/plain", lookup::AbstractProjected)
    LA.show_compact(io, mime, lookup)
    LA.print_index(io, mime, parent(lookup))
    print(io, " ")
    LA.print_order(io, lookup)
    print(io, " ")
    LA.print_span(io, lookup)
    print(io, " ")
    LA.print_sampling(io, lookup)
    if !(crs(lookup) isa Nothing)
        print(io, " crs: ", nameof(typeof(crs(lookup))))
    end
    if !(mappedcrs(lookup) isa Nothing)
        print(io, " mappedcrs: ", nameof(typeof(mappedcrs(lookup))))
    end
end


text:
@testset "4.median-of-two-sorted-arrays.jl" begin
    @test find_median_sorted_arrays([1, 3], [2]) == 2.0
    @test find_median_sorted_arrays([1, 2], [3, 4]) == 2.5
end


text:
function dlyap(A,Q)

	# Solving Discrete LYAPunov equation using Silvester equation
	#
	# Sylvester equation A*X + X*B + C = 0
	# if B = -(A')^(-1), C = Q*(A')^(-1)
	# => A*X - X*(A')^(-1) + Q*(A')^(-1) = 0
	# => A*X*A' - X + Q = 0  -> Discrete Lyapunov equation

	Atinv = Array{Float64,2}(inv(A'))
	B = -Atinv
	C = Q*Atinv
	X = sylvester(A, B, C)

	return X

end