text:
using Pkg
 Pkg.activate(".")

verbose = true
 if verbose   println("# Loading RvSpecML")    end
 using RvSpectML
include("shared/scripts.jl")
 if verbose   println("# Loading other packages")    end
 using DataFrames, Query, Statistics, Dates

include("read_neid_solar_data_20190918.jl")

order_list_timeseries = extract_orders(all_spectra,pipeline_plan)

line_list_df = prepare_line_list_pass1(linelist_for_ccf_filename, all_spectra, pipeline_plan,  v_center_to_avoid_tellurics=ccf_mid_velocity, Δv_to_avoid_tellurics = 21e3)

mask_scale_factors = [ 1.0,  2, 4, 6, 8, 10, 12, 14, 16 ]
 rms_rvs = zeros(length(mask_scale_factors))
 rms_binned_rvs = zeros(length(mask_scale_factors))
 println("Starting old Tophat CCFs...")
 for (i,mask_scale_factor) in enumerate(mask_scale_factors)
   println("# mask_scale_factor = ", mask_scale_factor)
   local (ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_scale_factor=mask_scale_factor, mask_type=:tophat,ccf_mid_velocity=ccf_mid_velocity, recalc=true, use_old = false)
   alg_fit_rv = RVFromCCF.MeasureRvFromCCFGaussian(frac_of_width_to_fit=3)
   local rvs_ccf = calc_rvs_from_ccf_total(ccfs, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false, alg_fit_rv=alg_fit_rv)
   rms_rvs[i] = std(rvs_ccf.-mean(rvs_ccf))
   rms_binned_rvs[i] = std(RvSpectML.bin_rvs_consecutive(rvs_ccf.-mean(rvs_ccf), 4))
 end

 using Plots
   plt1 = plot(mask_scale_factors, rms_rvs, color=1, label="New Tophat")
   plt2 = plot(mask_scale_factors, rms_binned_rvs, color=1, label=:none)
   plot(plt1,plt2,layout=(2,1) )
   xlabel!("Mask width scale parameter")
   ylabel!("RMS (m/s)")

#=
# Commented out since redundant now that old and new match
mask_scale_factors2 = [ 1.0,  2, 4, 6, 8, 10, 12 ]
 rms_rvs2 = zeros(length(mask_scale_factors2))
 rms_binned_rvs2 = zeros(length(mask_scale_factors2))
 println("Starting old Tophat CCFs...")
 for (i,mask_scale_factor) in enumerate(mask_scale_factors2)
   println("# mask_scale_factor = ", mask_scale_factor)
   (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:tophat, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true, use_old = true)
   rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
   rms_rvs2[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
   rms_binned_rvs2[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
 end
 scatter!(plt1,mask_scale_factors2, rms_rvs2, color=2, label="Old Tophat")
 scatter!(plt2,mask_scale_factors2, rms_binned_rvs2, color=2, label=:none)
 plot(plt1,plt2,layout=(2,1) )
 #xlabel!("Mask width scale parameter")
 #ylabel!("RMS (m/s)")
=#

mask_scale_factors3 = [ 1.0,  2, 4, 6, 8, 10, 12 ]
 rms_rvs3 = zeros(length(mask_scale_factors3))
 rms_binned_rvs3 = zeros(length(mask_scale_factors3))
  println("Starting Half Cos...")
  for (i,mask_scale_factor) in enumerate(mask_scale_factors3)
    println("# mask_scale_factor = ", mask_scale_factor)
    (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:halfcos, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true, use_old = false)
    rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
    rms_rvs3[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
    rms_binned_rvs3[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
  end
  plot!(plt1,mask_scale_factors3, rms_rvs3, color=3, label="Half Cos")
  plot!(plt2,mask_scale_factors3, rms_binned_rvs3, color=3, label=:none)
  plot(plt1,plt2,layout=(2,1) )

mask_scale_factors4 = [ 1.0,  2.0, 3.0, 4.0, 4.5, 5.0, 5.5, 6.0 ]
 rms_rvs4 = zeros(length(mask_scale_factors4))
   rms_binned_rvs4 = zeros(length(mask_scale_factors4))
   rms_rvs4_quad = zeros(length(mask_scale_factors4))
   rms_rvs4_cent = zeros(length(mask_scale_factors4))
   println("Starting Gaussian...")
    for (i,mask_scale_factor) in enumerate(mask_scale_factors4)
      println("# mask_scale_factor = ", mask_scale_factor)
      (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:gaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true, use_old = false)
      rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
      rms_rvs4[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
      rms_binned_rvs4[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
      rvs_ccf_quad = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false, alg_fit_rv= RVFromCCF.MeasureRvFromCCFQuadratic())
      #rvs_ccf_cent = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false, alg_fit_rv= RVFromCCF.MeasureRvFromCCFCentroid())
      rms_rvs4_quad[i] = std(rvs_ccf_quad.-mean(rvs_ccf_quad))
      #rms_rvs4_cent[i] = std(rvs_ccf_cent.-mean(rvs_ccf_cent))
    end

scatter!(plt1,mask_scale_factors4, rms_rvs4, color=4, label="Gaussian- Gaussian")
    scatter!(plt2,mask_scale_factors4, rms_binned_rvs4, color=4, label=:none)
    plot!(plt1,mask_scale_factors4, rms_rvs4_quad, color=5, label="Gaussian- Quadratic")
    #scatter!(plt1,mask_scale_factors4, rms_rvs4_cent, color=6, label="Gaussian- Centroid")
    plot(plt1,plt2,layout=(2,1) )

#=
# Commented out since this is slow, so people only run it when they really want to
mask_scale_factors5 = [ 4.0, 6, 8 ]
 rms_rvs5 = zeros(length(mask_scale_factors5))
   rms_binned_rvs5 = zeros(length(mask_scale_factors5))
   println("Starting Super Gaussian...")
    for (i,mask_scale_factor) in enumerate(mask_scale_factors5)
      println("# mask_scale_factor = ", mask_scale_factor)
      (ccfs_expr, v_grid_expr) = ccf_total(order_list_timeseries, line_list_df, pipeline_plan, mask_type=:supergaussian, mask_scale_factor=mask_scale_factor,  ccf_mid_velocity=ccf_mid_velocity, recalc=true, use_old = false)
      rvs_ccf_expr = calc_rvs_from_ccf_total(ccfs_expr, pipeline_plan, v_grid=v_grid_expr, times = order_list_timeseries.times, recalc=true, bin_consecutive=4, bin_nightly=false)
      rms_rvs5[i] = std(rvs_ccf_expr.-mean(rvs_ccf_expr))
      rms_binned_rvs5[i] = std(bin_rvs_consecutive(rvs_ccf_expr.-mean(rvs_ccf_expr), 4))
    end

plot!(plt1,mask_scale_factors5, rms_rvs5, color=5, label="Super-Gaussian")
    plot!(plt2,mask_scale_factors5, rms_binned_rvs5, color=5, label=:none)
    plot(plt1,plt2,layout=(2,1) )
=#

plt1 = plot(mask_scale_factors, rms_rvs, color=1, label="New Tophat")
 plt2 = plot(mask_scale_factors, rms_binned_rvs, color=1, label=:none)
 plot(plt1,plt2,layout=(2,1) )
 xlabel!("Mask width scale parameter")
 ylabel!("RMS (m/s)")
 #scatter!(plt1,mask_scale_factors2, rms_rvs2, color=2, label="Old Tophat")
 #scatter!(plt2,mask_scale_factors2, rms_binned_rvs2, color=2, label=:none)
 plot!(plt1,mask_scale_factors3, rms_rvs3, color=3, label="Half Cos")
 plot!(plt2,mask_scale_factors3, rms_binned_rvs3, color=3, label=:none)
 plot(plt1,plt2,layout=(2,1) )
 scatter!(plt1,mask_scale_factors4, rms_rvs4, color=4, label="Gaussian- Gaussian")
 scatter!(plt2,mask_scale_factors4, rms_binned_rvs4, color=4, label=:none)
 scatter!(plt1,mask_scale_factors4, rms_rvs4_quad, color=5, label="Gaussian- Quadratic")
 #scatter!(plt1,mask_scale_factors4, rms_rvs4_cent, color=6, label="Gaussian- Centroid")
 #plot!(plt1,mask_scale_factors5, rms_rvs5, color=5, label="Super-Gaussian")
 #plot!(plt2,mask_scale_factors5, rms_binned_rvs5, color=5, label=:none)

ylims!(plt1,0.8,1.0)
 ylims!(plt2,0,0.5)
 plot(plt1,plt2,layout=(2,1) )

savefig("neid_mask_width_comp.png")


text:
using Pkg

function importSim(simFile)
  simName = split(basename(simFile), ".")[1]
  if (!startswith(simFile, "/"))
    simFile = "$(pwd())/$simFile"
  end
  include(simFile)
  simName
end

function activatePkg()
  Pkg.activate(abspath(PROGRAM_FILE, "..", ".."))
end

ext: jl


text:
mutable struct GFProposeState
    choices::DynamicChoiceMap
    weight::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFProposeState(params::Dict{Symbol,Any})
    GFProposeState(choicemap(), 0., AddressVisitor(), params)
end

function traceat(state::GFProposeState, dist::Distribution{T},
              args, key) where {T}
    local retval::T

    # check that key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # sample return value
    retval = random(dist, args...)

    # update assignment
    set_value!(state.choices, key, retval)

    # update weight
    state.weight += logpdf(dist, retval, args...)

    retval
end

function traceat(state::GFProposeState, gen_fn::GenerativeFunction{T,U},
              args, key) where {T,U}
    local retval::T

    # check that key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # get subtrace
    (submap, weight, retval) = propose(gen_fn, args)

    # update assignment
    set_submap!(state.choices, key, submap)

    # update weight
    state.weight += weight

    retval
end

function splice(state::GFProposeState, gen_fn::DynamicDSLFunction, args::Tuple)
    prev_params = state.params
    state.params = gen_fn.params
    retval = exec(gen_fn, state, args)
    state.params = prev_params
    retval
end

function propose(gen_fn::DynamicDSLFunction, args::Tuple)
    state = GFProposeState(gen_fn.params)
    retval = exec(gen_fn, state, args)
    (state.choices, state.weight, retval)
end

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
function logsumexp{I<:Number}(x::AbstractArray{I})
    m = maximum(x)
    r = 0.0
    for i = 1:endof(x)
        r += exp(x[i] - m)
    end
    return log(r) + m
end


text:
using Pkg;
Pkg.Registry.add(RegistrySpec(; url="https://github.com/grlap/Registry"));
exit()

max_stars_repo_name: grlap/nghttpjl

max_issues_repo_name: grlap/nghttpjl


text:
abstract type FormatLine end

"""
A line in the vtree file format
"""
abstract type VtreeFormatLine <: FormatLine end

const VtreeFormatLines = AbstractVector{<:VtreeFormatLine}

struct VtreeCommentLine{T<:AbstractString} <: VtreeFormatLine
    comment::T
end

struct VtreeHeaderLine <: VtreeFormatLine end

struct VtreeInnerLine <: VtreeFormatLine
    node_id::UInt32
    left_id::UInt32
    right_id::UInt32
end

struct VtreeLeafLine <: VtreeFormatLine
    node_id::UInt32
    variable::Var
end

function compile_vtree_format_lines(lines::VtreeFormatLines,
                                    ::Type{V}=PlainVtree)::V where V<:Vtree
    compile_vtree_format_lines_m(lines, V)[1]
end

function compile_vtree_format_lines_m(lines::VtreeFormatLines,
                                      ::Type{V}=PlainVtree) where  V<:Vtree

    id2node = Dict{UInt32, V}()
    root = nothing

    compile_line(::Union{VtreeHeaderLine,VtreeCommentLine}) = () # do nothing

    compile_line(ln::VtreeLeafLine) = begin
        n = V(ln.variable)
        id2node[ln.node_id] = root = n
    end

    compile_line(ln::VtreeInnerLine) = begin
        left_node = id2node[ln.left_id]
        right_node = id2node[ln.right_id]
        n = V(left_node,right_node)
        id2node[ln.node_id] = root = n
    end

    foreach(compile_line, lines)
    root, id2node
end


text:
module VisualStringDistances

using LinearAlgebra
using UnbalancedOptimalTransport: UnbalancedOptimalTransport, DiscreteMeasure, KL,
                                  sinkhorn_divergence!

using DelimitedFiles
using StaticArrays

export printglyph, word_measure, visual_distance

include("compat.jl")
include("glyphs.jl")
include("glyphcoordinates.jl")
include("glue.jl")

"""
    visual_distance(::Type{T}, s::Union{Char,AbstractString},
                         t::Union{Char,AbstractString}; D=KL(one(T)), ϵ=T(0.1),
                         normalize=nothing) where {T}

Computes a measure of distance between the strings `s` and `t` in terms of their visual representation
as rendered by GNU Unifont and quantified by an unbalanced Sinkhorn divergence from UnbalancedOptimalTransport.jl.

* The keyword argument `D` chooses the `UnbalancedOptimalTransport.AbstractDivergence` used to penalize the creation
  or destruction of "mass" (black pixels). For `D = VisualStringDistances.KL(ρ)` for some number `ρ ≥ 0`,
  the distance is non-negative and zero if and only if the two visual representations of the strings
  are the same, as is generally desired.
* The keyword argument `ϵ` sets the "entropic regularization" in the Sinkhorn divergence; see the
  [documentation](https://ericphanson.github.io/UnbalancedOptimalTransport.jl/stable/optimal_transport/)
  there for more information. In short, smaller `ϵ` computes a quantity more directly related to the cost
  of moving mass, but takes longer to compute.
* The keyword argument `normalize` can be chosen to be a function which returns a normalizing constant
  given the maximum length of the two strings. The choice `normalize=identity` thus divides the result
  by the maximum length of the two strings. The choice `normalize=sqrt` has been found to give
  a good balance in some settings.

One may use [`printglyph`](@ref) to see the visual representation of the strings as rendered by GNU Unifont.

!!! note
    At the time of this writing, GNU Unifont is capable of rendering 57086 different unicode characters.
    However, it renders some unicode characters with the same graphical representation; specifically,
    689 distinct unicode characters have duplicate representations. Here's a set of six duplicates, for
    example:

    * 'Ꮋ': Unicode U+13BB (category Lu: Letter, uppercase)
    * 'Н': Unicode U+041D (category Lu: Letter, uppercase)
    * 'ꓧ': Unicode U+A4E7 (category Lo: Letter, other)
    * 'Ⲏ': Unicode U+2C8E (category Lu: Letter, uppercase)
    * 'Η': Unicode U+0397 (category Lu: Letter, uppercase)
    * 'H': ASCII/Unicode U+0048 (category Lu: Letter, uppercase)

    The visual distance between these, therefore, is returned as zero (up to numerical error).

## Example

```julia
julia> using VisualStringDistances

julia> printglyph("abc")
------------------------
------------------------
------------------------
---------#--------------
---------#--------------
---------#--------------
--####---#-###----####--
-#----#--##---#--#----#-
------#--#----#--#------
--#####--#----#--#------
-#----#--#----#--#------
-#----#--#----#--#------
-#---##--##---#--#----#-
--###-#--#-###----####--
------------------------
------------------------

julia> printglyph("def")
------------------------
------------------------
------------------------
------#-------------##--
------#------------#----
------#------------#----
--###-#---####-----#----
-#---##--#----#--#####--
-#----#--#----#----#----
-#----#--######----#----
-#----#--#---------#----
-#----#--#---------#----
-#---##--#----#----#----
--###-#---####-----#----
------------------------
------------------------

julia> visual_distance("abc", "def")
31.57060117541754

julia> visual_distance("abc", "abe")
4.979840716647487
```

"""
function visual_distance(::Type{T}, s::Union{Char,AbstractString},
                         t::Union{Char,AbstractString}; D=KL(one(T)), ϵ=T(0.1),
                         normalize=nothing) where {T}
    d = sinkhorn_divergence!(D, word_measure(T, s), word_measure(T, t), ϵ)
    if normalize !== nothing
        d = d / normalize(max(length(s), length(t)))
    end
    return d
end

# `Float64` default.
function visual_distance(s::Union{Char,AbstractString}, t::Union{Char,AbstractString};
                         D=KL(1.0), ϵ=0.1, normalize=nothing)
    visual_distance(Float64, s, t; D=D, ϵ=ϵ, normalize=normalize)
end

end


text:
# Exercise 1a
using LinearAlgebra

function trap_weights(x)
    return step(x) * [0.5; ones(length(x) - 2); 0.5]
end
x = range(0.0, 1.0, length = 10)
ω = trap_weights(x)
f(x) = x^2
dot(f.(x), ω)

struct UniformTrapezoidal
    count::Int
    Δ::Float64
end

Base.iterate(S::UniformTrapezoidal, state=1) = state > S.count ? nothing : (if state == 1 || state == count; Δ*0.5 else Δ*1.0 end, state+1)

function integration(f, x::AbstractArray; count = 10)
    Δ = step(x)
    int = similar(x)
    for (i,value) in enumerate(UniformTrapezoidal(count,Δ))
        int[i] = f(x[i]) * value
    end
    return sum(int)
end

count = 100
x = range(0.0, 1.0, length = count)
@show Δ = step(x)
f(x) = x^2

integral = integration(f,x,count = count)

struct Squares
    count::Int
end

Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)

for i in Squares(7)
    println(i)
end

struct SquaresVector <: AbstractArray{Int, 1}
    count::Int
end
Base.size(S::SquaresVector) = (S.count,)
Base.IndexStyle(::Type{<:SquaresVector}) = IndexLinear()
Base.getindex(S::SquaresVector, i::Int) = i*i

s = SquaresVector(4)
s[2]

# Exercise 1b

struct TrapezoidalVector <: AbstractArray{Int, 1}
    count::Int
end
Base.size(S::TrapezoidalVector) = (S.count,)
Base.IndexStyle(::Type{<:TrapezoidalVector}) = IndexLinear()
Base.getindex(S::TrapezoidalVector, i::Int) = 1
Base.setindex!(S::TrapezoidalVector, v, 2) = 2

V = TrapezoidalVector(4)

max_stars_repo_name: shanemcmiken/quantecon-notebooks-julia

max_issues_repo_name: shanemcmiken/quantecon-notebooks-julia

max_forks_repo_name: shanemcmiken/quantecon-notebooks-julia


text:
##################
## LLA <-> ECEF ##
##################

struct LLAfromECEF <: Transformation
    el::Ellipsoid

    a::Float64      # major axis
    f::Float64      # flattening
    e2::Float64     # Eccentricity squared
    e2m::Float64    # 1 - e2
    e2a::Float64    # |e2|
    e4a::Float64    # e2^2

    function LLAfromECEF(el,a,f,e2,e2m,e2a,e4a)
        if !(isfinite(a) && a > 0)
            error("Major radius is not positive")
        end
        if !(isfinite(f) && f < 1)
            error("Minor radius is not positive")
        end
        return new(el,a,f,e2,e2m,e2a,e4a)
    end
end
Base.show(io::IO, trans::LLAfromECEF) = print(io, "LLAfromECEF($(trans.el))")

"""
    LLAfromECEF(ellipsoid)

Construct a `Transformation` object to convert from ECEF coordinates
to LLA coordinates. Pre-caches ellipsoidal parameters for efficiency.
"""
function LLAfromECEF(el::Ellipsoid)
    a = el.a
    b = el.b
    f = el.f
    e2 = f*(2-f)
    e2m = (1-f)*(1-f)  #1 - e2
    e2a = abs(e2)
    e4a = e2*e2

    return LLAfromECEF(el, a, f, e2, e2m, e2a, e4a)
end

LLAfromECEF(datum::Datum) = LLAfromECEF(ellipsoid(datum))

function (trans::LLAfromECEF)(ecef::ECEF)
    # Ported to Julia by Andy Ferris, 2016 and re-released under MIT license.
    #/**
    # * \file Geocentric.cpp
    # * \brief Implementation for GeographicLib::Geocentric class
    # *
    # * Copyright (c) Charles Karney (2008-2015) <charles@karney.com> and licensed
    # * under the MIT/X11 License.  For more information, see
    # * http://geographiclib.sourceforge.net/
    # **********************************************************************/
    R = hypot(ecef.x, ecef.y)
    if R == 0
        slam = 0.0
        clam = 1.0
    else
        slam = ecef.y / R
        clam = ecef.x / R
    end
    h = hypot(R, ecef.z)    # Distance to center of earth

    if (trans.e4a == 0)
        # Treat the spherical case.  Dealing with underflow in the general case
        # with _e2 = 0 is difficult.  Origin maps to north pole, same as the
        # ellipsoidal case below
        if h == 0
            sphi = 1.0
            cphi = 0.0
        else
            sphi = ecef.z / h
            cphi = R / h
        end
        h -= trans.a
    else # Ellipsoidal
        # Treat prolate spheroids by swapping r and z here and by switching
        # the arguments to phi = atan(...) at the end.
        p = (R / trans.a) * (R / trans.a)
        q = trans.e2m * (ecef.z / trans.a) * (ecef.z / trans.a)
        r = (p + q - trans.e4a) / 6
        if (trans.f < 0)
            tmp = p
            p = q
            q = tmp
        end

        if ( !(trans.e4a * q == 0 && r <= 0) )

            # Avoid possible division by zero when r = 0 by multiplying
            # equations for s and t by r^3 and r, resp.
            S = trans.e4a * p * q / 4 # S = r^3 * s
            r2 = r * r
            r3 = r * r2
            disc = S * (2 * r3 + S)
            u = r
            if (disc >= 0)
                T3 = S + r3
                # Pick the sign on the sqrt to maximize abs(T3).  This minimizes
                # loss of precision due to cancellation.  The result is unchanged
                # because of the way the T is used in definition of u.
                T3 += (T3 < 0 ? -sqrt(disc) : sqrt(disc)) # T3 = (r * t)^3
                # N.B. cbrt always returns the real root.  cbrt(-8) = -2.
                T = cbrt(T3) # T = r * t
                # T can be zero; but then r2 / T -> 0.
                u += T + (T != 0 ? r2 / T : 0.0)
            else
                # T is complex, but the way u is defined the result is real.
                ang = atan(sqrt(-disc), -(S + r3))
                # There are three possible cube roots.  We choose the root which
                # avoids cancellation.  Note that disc < 0 implies that r < 0.
                u += 2 * r * cos(ang / 3)
            end

            v = sqrt(u*u + trans.e4a * q) # guaranteed positive
            # Avoid loss of accuracy when u < 0.  Underflow doesn't occur in
            # e4 * q / (v - u) because u ~ e^4 when q is small and u < 0.
            uv = (u < 0 ? trans.e4a * q / (v - u) : u + v) # u+v, guaranteed positive
            # Need to guard against w going negative due to roundoff in uv - q.
            w = max(0.0, trans.e2a * (uv - q) / (2 * v))
            # Rearrange expression for k to avoid loss of accuracy due to
            # subtraction.  Division by 0 not possible because uv > 0, w >= 0.
            k = uv / (sqrt(uv + w*w) + w)
            k1 = (trans.f >= 0 ? k : k - trans.e2)
            k2 = (trans.f >= 0 ? k + trans.e2 : k)
            d = k1 * R / k2
            H = hypot(ecef.z/k1, R/k2)
            sphi = (ecef.z/k1) / H
            cphi = (R/k2) / H
            h = (1 - trans.e2m/k1) * hypot(d, ecef.z)
        else  # e4 * q == 0 && r <= 0
            # This leads to k = 0 (oblate, equatorial plane) and k + e^2 = 0
            # (prolate, rotation axis) and the generation of 0/0 in the general
            # formulas for phi and h.  using the general formula and division by 0
            # in formula for h.  So handle this case by taking the limits:
            # f > 0: z -> 0, k      ->   e2 * sqrt(q)/sqrt(e4 - p)
            # f < 0: R -> 0, k + e2 -> - e2 * sqrt(q)/sqrt(e4 - p)
            zz = sqrt((trans.f >= 0 ? trans.e4a - p : p) / trans.e2m)
            xx = sqrt( trans.f <  0 ? trans.e4a - p : p)
            H = hypot(zz, xx)
            sphi = zz / H
            cphi = xx / H
            if (ecef.z < 0)
                sphi = -sphi # for tiny negative Z (not for prolate)
            end
            h = - trans.a * (trans.f >= 0 ? trans.e2m : 1.0) * H / trans.e2a
        end
    end
    lat = 180 * atan(sphi, cphi) / pi
    lon = 180 * atan(slam, clam) / pi

    return LLA(lat, lon, h)
end

"""
    ECEFfromLLA(ellipsoid_or_datum)

Construct a `Transformation` object to convert from LLA coordinates
to ECEF coordinates.
"""
struct ECEFfromLLA <: Transformation
    el::Ellipsoid
end
Base.show(io::IO, trans::ECEFfromLLA) = print(io, "ECEFfromLLA($(trans.el))")

ECEFfromLLA(datum::Datum) = ECEFfromLLA(ellipsoid(datum))

function (trans::ECEFfromLLA)(lla::LLA)
    ϕdeg, λdeg, h = lla.lat, lla.lon, lla.alt

    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)
    sinλ, cosλ = sind(λdeg), cosd(λdeg)

    N = trans.el.a / sqrt(1 - trans.el.e2 * sinϕ^2)  # Radius of curvature (meters)

    x = (N + h) * cosϕ * cosλ
    y = (N + h) * cosϕ * sinλ
    z = (N * (1 - trans.el.e2) + h) * sinϕ

    return ECEF(x, y, z)
end

Base.inv(trans::LLAfromECEF) = ECEFfromLLA(trans.el)
Base.inv(trans::ECEFfromLLA) = LLAfromECEF(trans.el)

##################
## ENU <-> NED ##
##################

"""
    ENUfromNED(ned::NED)

Construct a `Transformation` object to convert from local `NED` coordinates
to local `ENU` coordinates centered at the same origin. This is a simple
permutation of coordinates and sign change for the altitude.
"""
struct ENUfromNED <: Transformation end     # singleton type

Base.show(io::IO, ::ENUfromNED) = print(io, "ENUfromNED()")

function (::ENUfromNED)(ned::NED)
    ENU(ned.e, ned.n, -ned.d)
end

"""
    NEDfromENU(enu::ENU)

Construct a `Transformation` object to convert from local `ENU` coordinates
to local `NED` coordinates centered at the same origin. This is a simple
permutation of coordinates and sign change for the altitude.
"""
struct NEDfromENU <: Transformation end     # singleton type

Base.show(io::IO, ::NEDfromENU) = print(io, "NEDfromENU()")

function (::NEDfromENU)(enu::ENU)
    NED(enu.n, enu.e, -enu.u)
end

Base.inv(::ENUfromNED) = NEDfromENU()
Base.inv(::NEDfromENU) = ENUfromNED()

##################
## ECEF <-> ENU ##
##################

"""
    ENUfromECEF(origin, datum)
    ENUfromECEF(origin::UTM, zone, isnorth, datum)
    ENUfromECEF(origin::ECEF, lat, lon)

Construct a `Transformation` object to convert from global `ECEF` coordinates
to local `ENU` coordinates centered at the `origin`. This object pre-caches both the
ECEF coordinates and latitude and longitude of the origin for maximal efficiency.
"""
struct ENUfromECEF{T} <: Transformation
    origin::ECEF{T}
    lat::T
    lon::T
end
ENUfromECEF(origin::LLA, datum) = ENUfromECEF(ECEFfromLLA(datum)(origin), origin.lat, origin.lon)
function ENUfromECEF(origin::ECEF, datum)
    origin_lla = LLAfromECEF(datum)(origin)
    ENUfromECEF(origin, origin_lla.lat, origin_lla.lon)
end
Base.show(io::IO, trans::ENUfromECEF) = print(io, "ENUfromECEF($(trans.origin), lat=$(trans.lat)°, lon=$(trans.lon)°)")
Base.isapprox(t1::ENUfromECEF, t2::ENUfromECEF; kwargs...) = isapprox(t1.origin, t2.origin; kwargs...) && isapprox(t1.lat, t2.lat; kwargs...) && isapprox(t1.lon, t2.lon; kwargs...)

function (trans::ENUfromECEF)(ecef::ECEF)
    ϕdeg, λdeg = trans.lat, trans.lon

    ∂x = ecef.x - trans.origin.x
    ∂y = ecef.y - trans.origin.y
    ∂z = ecef.z - trans.origin.z

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # R = [     -sinλ       cosλ  0.0
    #      -cosλ*sinϕ -sinλ*sinϕ cosϕ
    #       cosλ*cosϕ  sinλ*cosϕ sinϕ]
    #
    # east, north, up = R * [∂x, ∂y, ∂z]
    east  = ∂x * -sinλ      + ∂y * cosλ       + ∂z * 0.0
    north = ∂x * -cosλ*sinϕ + ∂y * -sinλ*sinϕ + ∂z * cosϕ
    up    = ∂x * cosλ*cosϕ  + ∂y * sinλ*cosϕ  + ∂z * sinϕ

    return ENU(east, north, up)
end

"""
    ECEFfromENU(origin, datum)
    ECEFfromENU(origin::UTM, zone, isnorth, datum)
    ECEFfromENU(origin::ECEF, lat, lon)

Construct a `Transformation` object to convert from local `ENU` coordinates
centred at `origin` to global `ECEF` coodinates. This object pre-caches both the
ECEF coordinates and latitude and longitude of the origin for maximal efficiency.
"""
struct ECEFfromENU{T} <: Transformation
    origin::ECEF{T}
    lat::T
    lon::T
end
ECEFfromENU(origin::LLA, datum) = ECEFfromENU(ECEFfromLLA(datum)(origin), origin.lat, origin.lon)
function ECEFfromENU(origin::ECEF, datum)
    origin_lla = LLAfromECEF(datum)(origin)
    ECEFfromENU(origin, origin_lla.lat, origin_lla.lon)
end
Base.show(io::IO, trans::ECEFfromENU) = print(io, "ECEFfromENU($(trans.origin), lat=$(trans.lat)°, lon=$(trans.lon)°)")
Base.isapprox(t1::ECEFfromENU, t2::ECEFfromENU; kwargs...) = isapprox(t1.origin, t2.origin; kwargs...) && isapprox(t1.lat, t2.lat; kwargs...) && isapprox(t1.lon, t2.lon; kwargs...)

function (trans::ECEFfromENU)(enu::ENU)
    ϕdeg, λdeg = trans.lat, trans.lon

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # Rᵀ = [-sinλ -cosλ*sinϕ cosλ*cosϕ
    #        cosλ -sinλ*sinϕ sinλ*cosϕ
    #         0.0       cosϕ      sinϕ]
    # Δx, Δy, Δz = Rᵀ * [east, north, up]
    Δx = -sinλ * enu.e + -cosλ*sinϕ * enu.n + cosλ*cosϕ * enu.u
    Δy =  cosλ * enu.e + -sinλ*sinϕ * enu.n + sinλ*cosϕ * enu.u
    Δz =   0.0 * enu.e +       cosϕ * enu.n +      sinϕ * enu.u

    X = trans.origin.x + Δx
    Y = trans.origin.y + Δy
    Z = trans.origin.z + Δz

    return ECEF(X,Y,Z)
end

Base.inv(trans::ECEFfromENU) = ENUfromECEF(trans.origin, trans.lat, trans.lon)
Base.inv(trans::ENUfromECEF) = ECEFfromENU(trans.origin, trans.lat, trans.lon)

##################
## ECEF <-> NED ##
##################

"""
    NEDfromECEF(origin, datum)
    NEDfromECEF(origin::UTM, zone, isnorth, datum)
    NEDfromECEF(origin::ECEF, lat, lon)

Construct a composite transformation NEDfromENU() ∘ ENUfromECEF(origin, datum)
to convert from global `ECEF` coordinates to local `NED` coordinates centered at the `origin`.
This object pre-caches both the ECEF coordinates and latitude and longitude of the origin for maximal efficiency.
"""
NEDfromECEF(origin, datum) = NEDfromENU() ∘ ENUfromECEF(origin, datum)

"""
    ECEFfromNED(origin, datum)
    ECEFfromNED(origin::UTM, zone, isnorth, datum)
    ECEFfromNED(origin::ECEF, lat, lon)

Construct a composite transformation ECEFfromENU(origin,datum) ∘ ENUfromNED()
to convert from local `NED` coordinates centred at `origin` to global `ECEF` coodinates.
This object pre-caches both the ECEF coordinates and latitude and longitude of the origin for maximal efficiency.
"""
ECEFfromNED(origin, datum) = ECEFfromENU(origin,datum) ∘ ENUfromNED()

#################
## LLA <-> ENU ##
#################

"""
    ENUfromLLA(origin, datum)

Creates composite transformation `ENUfromECEF(origin, datum) ∘ ECEFfromLLA(datum)`.
"""
ENUfromLLA(origin, datum) = ENUfromECEF(origin, datum) ∘ ECEFfromLLA(datum)

"""
    LLAfromENU(origin, datum)

Creates composite transformation `LLAfromECEF(datum) ∘ ECEFfromENU(origin, datum)`.
"""
LLAfromENU(origin, datum) = LLAfromECEF(datum) ∘ ECEFfromENU(origin, datum)

#################
## LLA <-> NED ##
#################

"""
    NEDfromLLA(origin, datum)

Creates composite transformation `NEDfromECEF(origin, datum) ∘ ECEFfromLLA(datum)`.
"""
NEDfromLLA(origin, datum) = NEDfromECEF(origin, datum) ∘ ECEFfromLLA(datum)

"""
    LLAfromNED(origin, datum)

Creates composite transformation `LLAfromECEF(datum) ∘ ECEFfromNED(origin, datum)`.
"""
LLAfromNED(origin, datum) = LLAfromECEF(datum) ∘ ECEFfromNED(origin, datum)

#################
## LLA <-> UTM ##
#################

"""
    LLAfromUTM(zone, isnorth::Bool, datum)

Construct a `Transformation` object to convert from `UTM` coordinates in
the specified zone and hemisphere (`isnorth = true` or `false`) to global `LLA`
coordinates. Pre-caches ellipsoidal parameters for efficiency and performs
Charles Karney's accurate 6th-order series expansion algorithm.

(See also `LLAfromUTMZ`)
"""
struct LLAfromUTM{Datum,Order} <: Transformation
    zone::UInt8
    isnorth::Bool # hemisphere
    tm::TransverseMercator{Order}
    datum::Datum
end
LLAfromUTM(zone::UInt8, h, d) = LLAfromUTM(UInt8(zone), h, TransverseMercator(d), d)
LLAfromUTM(zone::Integer, h, d) = LLAfromUTM(UInt8(zone), h, d)
Base.show(io::IO, trans::LLAfromUTM) = print(io, "LLAfromUTM(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")), $(trans.datum))")

function (trans::LLAfromUTM)(utm::UTM)
    if trans.zone == 0
        # Do inverse steriographic projection
        k0 = 0.994
        x = (utm.x - 2e6)
        y = (utm.y - 2e6)
        (lat,lon,γ,k) = polarst_inv(trans.isnorth, k0, trans.tm, x, y)

        return LLA(lat, lon, utm.z)
    else
        lon_ref = Float64(utm_meridian(trans.zone))
        k0 = 0.9996 # Horizontal scaling factor
        x = (utm.x - 5e5) # Convention has 500km offset for easting
        y = (utm.y - (trans.isnorth ? 0.0 : 1e7)) # Northing offset for southern hemisphere
        #(lat,lon,k,γ) = transverse_mercator_inv(lat_ref, lon_ref, x, y, ellipsoid(trans.datum))
        (lat,lon,γ,k) = transverse_mercator_reverse(lon_ref, x, y, k0, trans.tm)
        return LLA(lat, lon, utm.z) # Note: scaling not applied to vertical dimension
    end
end

"""
    UTMfromLLA(zone, isnorth::Bool, datum)

Construct a `Transformation` object to convert from global `LLA` coordinates
to `UTM` coordinates in the specified zone and hemisphere (`isnorth = true` or `false`). Pre-caches
ellipsoidal parameters for efficiency and performs Charles Karney's accurate
6th-order series expansion algorithm.

(See also `UTMZfromLLA`)
"""
struct UTMfromLLA{Datum,Order} <: Transformation
    zone::UInt8
    isnorth::Bool # true = north, false = south
    tm::TransverseMercator{Order}
    datum::Datum
end
UTMfromLLA(zone::UInt8, h, d) = UTMfromLLA(zone, h, TransverseMercator(d), d)
UTMfromLLA(zone::Integer, h, d) = UTMfromLLA(UInt8(zone), h, d)
Base.show(io::IO, trans::UTMfromLLA) = print(io, "UTMfromLLA(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")), $(trans.datum))")

function (trans::UTMfromLLA)(lla::LLA)
    if trans.zone == 0
        # Do polar steriographic projection
        k0 = 0.994
        (x,y,γ,k) = polarst_fwd(trans.isnorth, k0, trans.tm, lla.lat, lla.lon)
        x = x + 2e6
        y = y + 2e6

        return UTM(x, y, lla.alt)
    else
        lon_ref = Float64(utm_meridian(trans.zone))
        k0 = 0.9996 # Horizontal scaling factor
        #(x,y,k,γ) = transverse_mercator(lat_ref, lon_ref, lla.lat, lla.lon, ellipsoid(trans.datum))
        (x,y,γ,k) = transverse_mercator_forward(lon_ref, lla.lat, lla.lon, k0, trans.tm)
        x = 5e5 + x # Convention has 500km offset for easting
        y = (trans.isnorth ? 0.0 : 1e7) + y # Northing offset for southern hemisphere
        # also, k = k * k0
        return UTM(x, y, lla.alt) # Note: scaling not applied to vertical dimension
    end
end

Base.inv(trans::LLAfromUTM) = UTMfromLLA(trans.zone, trans.isnorth, trans.tm, trans.datum)
Base.inv(trans::UTMfromLLA) = LLAfromUTM(trans.zone, trans.isnorth, trans.tm, trans.datum)

##################
## ECEF <-> UTM ##
##################

"""
    UTMfromECEF(zone, isnorth, datum)

Creates composite transformation `UTMfromLLA(zone, isnorth, datum) ∘ LLAfromECEF(datum)`.
"""
UTMfromECEF(zone, isnorth, datum) = UTMfromLLA(zone, isnorth, datum) ∘ LLAfromECEF(datum)

"""
    ECEFfromUTM(zone, isnorth, datum)

Creates composite transformation `ECEFfromLLA(datum) ∘ LLAfromUTM(zone, isnorth, datum)`.
"""
ECEFfromUTM(zone, isnorth, datum) = ECEFfromLLA(datum) ∘ LLAfromUTM(zone, isnorth, datum)

##################
## LLA <-> UTMZ ##
##################

"""
    LLAfromUTMZ(datum)

Construct a `Transformation` object to convert from global `UTMZ`
coordinates to global `LLA` coordinates. Pre-caches
ellipsoidal parameters for efficiency and performs Charles Karney's accurate
6th-order series expansion algorithm.

(See also `LLAfromUTM`)
"""
struct LLAfromUTMZ{Datum,Order} <: Transformation
    tm::TransverseMercator{Order}
    datum::Datum
end
LLAfromUTMZ(datum) = LLAfromUTMZ(TransverseMercator(datum), datum)
Base.show(io::IO, trans::LLAfromUTMZ) = print(io, "LLAfromUTMZ($(trans.datum))")

function (trans::LLAfromUTMZ)(utm::UTMZ)
    if utm.zone == 0
        # Do inverse steriographic projection
        k0 = 0.994
        x = (utm.x - 2e6)
        y = (utm.y - 2e6)
        (lat,lon,γ,k) = polarst_inv(utm.isnorth, k0, trans.tm, x, y)

        return LLA(lat, lon, utm.z)
    else
        lon_ref = Float64(utm_meridian(utm.zone))
        k0 = 0.9996 # Horizontal scaling factor
        x = (utm.x - 5e5) # Convention has 500km offset for easting
        y = (utm.y - (utm.isnorth ? 0.0 : 1e7)) # Northing offset for southern hemisphere
        #(lat,lon,k,γ) = transverse_mercator_inv(lat_ref, lon_ref, x, y, ellipsoid(trans.datum))
        (lat,lon,γ,k) = transverse_mercator_reverse(lon_ref, x, y, k0, trans.tm)
        return LLA(lat, lon, utm.z) # Note: scaling not applied to vertical dimension
    end
end

"""
    UTMZfromLLA(datum)

Construct a `Transformation` object to convert from global `LLA` coordinates
to global `UTMZ` coordinates. The zone and hemisphere is automatically calculated
following the standard definitions (including exceptions in Norway). Pre-caches
ellipsoidal parameters for efficiency and performs Charles Karney's accurate
6th-order series expansion algorithm.

(See also `UTMfromLLA`)
"""
struct UTMZfromLLA{Datum,Order} <: Transformation
    tm::TransverseMercator{Order}
    datum::Datum
end
UTMZfromLLA(datum) = UTMZfromLLA(TransverseMercator(datum), datum)
Base.show(io::IO, trans::UTMZfromLLA) = print(io, "UTMZfromLLA($(trans.datum))")

function (trans::UTMZfromLLA)(lla::LLA)
    (zone, isnorth) = utm_zone(lla)
    if zone == 0
        # Do polar steriographic projection
        k0 = 0.994
        (x,y,γ,k) = polarst_fwd(isnorth, k0, trans.tm, lla.lat, lla.lon)
        x = x + 2e6
        y = y + 2e6

        return UTMZ(x, y, lla.alt, 0, isnorth)
    else
        lon_ref = Float64(utm_meridian(zone))
        k0 = 0.9996 # Horizontal scaling factor
        #(x,y,k,γ) = transverse_mercator(lat_ref, lon_ref, lla.lat, lla.lon, ellipsoid(trans.datum))
        (x,y,γ,k) = transverse_mercator_forward(lon_ref, lla.lat, lla.lon, k0, trans.tm)
        x = 5e5 + x # Convention has 500km offset for easting
        y = (isnorth ? 0.0 : 1e7) + y # Northing offset for southern hemisphere
        # also, k = k * k0
        return UTMZ(x, y, lla.alt, zone, isnorth) # Note: scaling not applied to vertical dimension
    end
end

Base.inv(trans::LLAfromUTMZ) = UTMZfromLLA(trans.tm, trans.datum)
Base.inv(trans::UTMZfromLLA) = LLAfromUTMZ(trans.tm, trans.datum)

###################
## UTM <-> UTMZ ##
###################

"""
    UTMZfromUTM(zone, isnorth, datum)

Transformation to append the UTM/UPS zone and hemisphere to a raw `UTM` data point.
"""
struct UTMZfromUTM{Datum} <: Transformation
    zone::Int
    isnorth::Bool
    datum::Datum
end
UTMZfromUTM(zone::Integer, h, d::D) where {D} = UTMZfromUTM{D}(UInt8(zone), h, d)
Base.show(io::IO, trans::UTMZfromUTM) = print(io, "UTMZfromUTM(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")), $(trans.datum))")
#Base.show(io::IO, trans::UTMZfromUTM) = print(io, "UTMZfromUTM(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")))")

"""
    UTMfromUTMZ(zone, isnorth, datum)

Transformation to remove the zone and hemisphere from `UTMZ` data point, and
automatically convert the data to the specified zone if necessary.
"""
struct UTMfromUTMZ{Datum} <: Transformation
    zone::Int
    isnorth::Bool
    datum::Datum
end
UTMfromUTMZ(zone::Integer, h, d::D) where {D} = UTMfromUTMZ{D}(UInt8(zone), h, d)
Base.show(io::IO, trans::UTMfromUTMZ) = print(io, "UTMfromUTMZ(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")), $(trans.datum))")
#Base.show(io::IO, trans::UTMfromUTMZ) = print(io, "UTMfromUTMZ(zone=$(trans.zone == 0 ? "polar" : trans.zone) ($(trans.isnorth ? "north" : "south")))")

(trans::UTMZfromUTM)(utm::UTM) = UTMZ(utm.x, utm.y, utm.z, trans.zone, trans.isnorth)
function (trans::UTMfromUTMZ)(utm::UTMZ)
    if trans.zone == utm.zone && trans.isnorth == utm.isnorth
        UTM(utm.x, utm.y, utm.z)
    else
        # Should this be an error or an automatic transformation to the correct zone?
        #error("Incorrect UTM zone")
        UTMfromLLA(trans.zone, trans.isnorth, trans.datum)(LLAfromUTMZ(trans.datum)(utm))
    end
end

Base.inv(trans::UTMfromUTMZ) = UTMZfromUTM(trans.zone, trans.isnorth, trans.datum)
Base.inv(trans::UTMZfromUTM) = UTMfromUTMZ(trans.zone, trans.isnorth, trans.datum)

###################
## ECEF <-> UTMZ ##
###################

"""
    UTMZfromECEF(datum)

Creates composite transformation `UTMZfromLLA(datum) ∘ LLAfromECEF(datum)`.
"""
UTMZfromECEF(datum) = UTMZfromLLA(datum) ∘ LLAfromECEF(datum)

"""
    ECEFfromUTMZ(datum)

Creates composite transformation `ECEFfromLLA(datum) ∘ LLAfromUTMZ(datum)`.
"""
ECEFfromUTMZ(datum) = ECEFfromLLA(datum) ∘ LLAfromUTMZ(datum)

##################
## ENU <-> UTMZ ##
##################

ENUfromECEF(origin::UTMZ, datum) = ENUfromECEF(LLAfromUTMZ(datum)(origin), datum)
ECEFfromENU(origin::UTMZ, datum) = ECEFfromENU(LLAfromUTMZ(datum)(origin), datum)

"""
    ENUfromUTMZ(origin, datum)

Creates composite transformation `ENUfromLLA(origin, datum) ∘ LLAfromUTMZ(datum)`.
"""
ENUfromUTMZ(origin, datum) = ENUfromLLA(origin, datum) ∘ LLAfromUTMZ(datum)

"""
    UTMZfromENU(origin, datum)

Creates composite transformation `UTMZfromLLA(datum) ∘ LLAfromENU(origin, datum)`.
"""
UTMZfromENU(origin, datum) = UTMZfromLLA(datum) ∘ LLAfromENU(origin, datum)

#################
## ENU <-> UTM ##
#################
ENUfromECEF(origin::UTM, zone::Integer, isnorth::Bool, datum) = ENUfromECEF(LLAfromUTM(zone, isnorth, datum)(origin), datum)
ECEFfromENU(origin::UTM, zone::Integer, isnorth::Bool, datum) = ECEFfromENU(LLAfromUTM(zone, isnorth, datum)(origin), datum)

# Assume origin and utm point share the same zone and hemisphere
UTMfromENU(origin::UTM, zone::Integer, isnorth::Bool, datum) = UTMfromLLA(zone, isnorth, datum) ∘ LLAfromENU(UTMZ(origin, zone, isnorth), datum)
ENUfromUTM(origin::UTM, zone::Integer, isnorth::Bool, datum) = ENUfromLLA(UTMZ(origin, zone, isnorth), datum) ∘ LLAfromUTM(zone, isnorth, datum)

"""
    UTMfromENU(origin, zone, isnorth, datum)

Creates composite transformation `UTMfromLLA(zone, isnorth, datum) ∘ LLAfromENU(origin, datum)`.
If `origin` is a `UTM` point, then it is assumed it is in the given specified zone and hemisphere.
"""
UTMfromENU(origin, zone::Integer, isnorth::Bool, datum) = UTMfromLLA(zone, isnorth, datum) ∘ LLAfromENU(origin, datum)

"""
    ENUfromUTM(origin, zone, isnorth, datum)

Creates composite transformation `UTMfromLLA(zone, isnorth, datum) ∘ LLAfromENU(origin, datum)`.
If `origin` is a `UTM` point, then it is assumed it is in the given specified zone and hemisphere.
"""
ENUfromUTM(origin, zone::Integer, isnorth::Bool, datum) = ENUfromLLA(origin, datum) ∘ LLAfromUTM(zone, isnorth, datum)

##################
## NED <-> UTMZ ##
##################

NEDfromECEF(origin::UTMZ, datum) = NEDfromECEF(LLAfromUTMZ(datum)(origin), datum)
ECEFfromNED(origin::UTMZ, datum) = ECEFfromNED(LLAfromUTMZ(datum)(origin), datum)

"""
    NEDfromUTMZ(origin, datum)

Creates composite transformation `ENUfromLLA(origin, datum) ∘ LLAfromUTMZ(datum)`.
"""
NEDfromUTMZ(origin, datum) = NEDfromLLA(origin, datum) ∘ LLAfromUTMZ(datum)

"""
    UTMZfromNED(origin, datum)

Creates composite transformation `UTMZfromLLA(datum) ∘ LLAfromNED(origin, datum)`.
"""
UTMZfromNED(origin, datum) = UTMZfromLLA(datum) ∘ LLAfromNED(origin, datum)

#################
## NED <-> UTM ##
#################
NEDfromECEF(origin::UTM, zone::Integer, isnorth::Bool, datum) = NEDfromECEF(LLAfromUTM(zone, isnorth, datum)(origin), datum)
ECEFfromNED(origin::UTM, zone::Integer, isnorth::Bool, datum) = ECEFfromNED(LLAfromUTM(zone, isnorth, datum)(origin), datum)

# Assume origin and utm point share the same zone and hemisphere
UTMfromNED(origin::UTM, zone::Integer, isnorth::Bool, datum) = UTMfromLLA(zone, isnorth, datum) ∘ LLAfromNED(UTMZ(origin, zone, isnorth), datum)
NEDfromUTM(origin::UTM, zone::Integer, isnorth::Bool, datum) = NEDfromLLA(UTMZ(origin, zone, isnorth), datum) ∘ LLAfromUTM(zone, isnorth, datum)

"""
    UTMfromNED(origin, zone, isnorth, datum)

Creates composite transformation `UTMfromLLA(zone, isnorth, datum) ∘ LLAfromNED(origin, datum)`.
If `origin` is a `UTM` point, then it is assumed it is in the given specified zone and hemisphere.
"""
UTMfromNED(origin, zone::Integer, isnorth::Bool, datum) = UTMfromLLA(zone, isnorth, datum) ∘ LLAfromNED(origin, datum)

"""
    NEDfromUTM(origin, zone, isnorth, datum)

Creates composite transformation `UTMfromLLA(zone, isnorth, datum) ∘ LLAfromNED(origin, datum)`.
If `origin` is a `UTM` point, then it is assumed it is in the given specified zone and hemisphere.
"""
NEDfromUTM(origin, zone::Integer, isnorth::Bool, datum) = NEDfromLLA(origin, datum) ∘ LLAfromUTM(zone, isnorth, datum)


text:
using Test
using AbstractLogic

@test LogicalCombo().logical == Bool[]
@test LogicalCombo().domain == []
@test LogicalCombo().keys == Symbol[]
@test LogicalCombo().commands == String[]

@test size(expand(LogicalCombo(), a = 1:3, b= 1:4, c = 1:2), 1) == 3*4*2
@test size(expand(LogicalCombo(), ["a","b","c"], ["a1","a2","a3"]), 1) == 3^3
@test size(expand(LogicalCombo(),[:a => 1:3, :b => 1:4, :c => 1:5]), 1) == 3*4*5
@test size(expand(LogicalCombo(), a = 1:3, b = 1:4, c = 1:5), 1) == 3*4*5
@test size(expand(LogicalCombo(a = 1:3, b = 1:4), c = 1:5), 1) == 3*4*5
@test_throws "key :a already defined!" expand(LogicalCombo(a=1), a = 1:3)
@test size(expand(LogicalCombo()), 1) == 0

@test nfeasible(LogicalCombo()) == 0
@test LogicalCombo() |> collect |> size == (1,1)

@test string(LogicalCombo()) == string(LogicalCombo(Symbol[], [], Bool[]))

x = LogicalCombo([:x=>1:2, :y=>1:4])
@test  x |> nfeasible == 8

@test x[1:2] == [true, true]
@test x[1] == x[2]
@test x[5,"x"] == x[2,"y"]
@test x[2,"x"] != x[2,"y"]
@test x[:,:x] == x[:,"x"]
@test x[1:2,:] == x[[1,2],:]
@test x[:,:x,:] == x[:,:x,true]
@test x[1,:x,true] == 1
@test x[1:3] = true

@test x[:,:,true] != x[:,:,!true]
@test x[:,:] |> size == (8,2)
@test x[:,:, varnames=true] |> size == (9,2)

x[1:2] = [false, false]; @test x[1:2] == [false, false]
x[3] = false; @test x[3] == false

@test range(x) |> length == 2
@test range(x, "x") |> length == 1

@test length(AbstractLogic.sample(LogicalCombo())) == 0
@test length(AbstractLogic.sample(x)) == 2
@test length(AbstractLogic.sample(x,3)) == 6
@test length(AbstractLogic.sample(x,2,false)) == 4


text:
using Hyperspecialize
println(:(@concretize Real Float32))
println(macroexpand(Main, :(@concretize Real Float32)))
println(:(@widen Real Float32))
println(macroexpand(Main, :(@widen Real Float32)))
println(:(@concretization Real))
println(macroexpand(Main, :(@concretization Real)))
println(:(@concretize (foo, Real) Float32))
println(macroexpand(Main, :(@concretize (foo, Real) Float32)))
println(:(@widen (foo, Real) Float32))
println(macroexpand(Main, :(@widen (foo, Real) Float32)))
println(:(@concretization (foo, Real)))
println(macroexpand(Main, :(@concretization (foo, Real))))
println(:(@replicable bar(x::@hyperspecialize((foo, Real))) = x))
println(macroexpand(Main, :(@replicable bar(x::@hyperspecialize((foo, Real))) = x)))