text:
export nutation_angles, mean_obliq, frame_tie, nutation, precession, cel_pole, e_tilt,
       ira_equinox, cio_location, cio_basis, era, tdb2tt, sidereal_time, wobble, spin,
       ter2cel, OnSurface, refract, equ2hor, CatEntry, Object, SkyPos, Observer, app_star,
       InSpace, place

"""
    nutation_angles(t)

Computes the values for nutation in longitude and nutation in obliquity for a given TDB julian date.

# Aruments
- `t::Real`: TDB time in Julian centuries since J2000.0

# Optional Arguments
- `accuracy::Symbol=:full`: Sets the accuracy level of `:full` or `:reduced`

# Returns
`(dpsi,deps)` where
- `dpsi`: Nutation in longitude in arcseconds
- `deps`: Nutation in obliquity in arcseconds
"""
function nutation_angles(t::Real; accuracy::Symbol=:full)
    @assert accuracy ∈ Set([:full, :reduced])
    t1 = t * 36525.0
    # High accuracy mode uses IAU 200A
    if accuracy == :full
        dpsi, deps = iau2000a(T0, t1)
    elseif accuracy == :reduced
        # Low accuracy mode uses the specially truncated version of IAU 2000A, called NU200K
        dpsi, deps = nu2000k(T0, t1)
    end
    # Convert output to arcseconds
    dpsi /= ASEC2RAD
    deps /= ASEC2RAD
    # Return
    return (dpsi, deps)
end

"""
    mean_obliq(jd_tdb)

Computes the mean obliquity of the ecliptic
"""
function mean_obliq(jd_tdb::Real)
    # Julian centuries since J2000.0 epoch
    t = (jd_tdb - T0) / 36525.0
    ε = ((((-0.0000000434 * t - 0.000000576) * t + 0.00200340) * t - 0.0001831) * t -
         46.836769) * t + 84381.406
    return ε
end

"""
    frame_tie(pos;direction)

Transforms a position vector from the dynamical reference system to ICRS, or vice versa.

# Argumenets
- `pos::AbstractVector`: Position vector, equitorial rectangular coordinates
- `direction::Symbol`: Transformation direction, `:dynamic2icrs` or `:icrs2dynamic`
"""
function frame_tie(pos::AbstractVector, direction::Symbol)
    @assert direction ∈ Set([:dynamic2icrs, :icrs2dynamic])
    if direction == :dynamic2icrs
        return frame_tie_rot * pos
    elseif direction == :icrs2dynamic
        return frame_tie_rot' * pos
    end
end

"""
    nutation(jd_tdb,pos)

Nutates equitorial rectangular coordinates from mean equator and equinox of epoch to true equator or
does the transformation in reverse given the flag `direction`.

# Arguments
- `jd_tdb::Real`: TDB Julian Date
- `pos::AbstractVector`: Position vector, geocentric equatorial rectangular coordinates, referred to mean equator and equinox of epoch

# Optional Arguments
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
- `direction::Symbol=:mean2true`: Perform the transformation in either the `:mean2true` or `:true2mean` direction
"""
function nutation(jd_tdb::Real, pos::AbstractVector; accuracy::Symbol=:full,
                  direction::Symbol=:mean2true)
    @assert direction ∈ Set([:mean2true, :true2mean])
    # Call e_tilt to get the obliquity and nutation angles
    oblm, oblt, _, psi, _ = e_tilt(jd_tdb; accuracy=accuracy)
    # Construct rotation matrix
    sobm, cobm = sincosd(oblm)
    sobt, cobt = sincosd(oblt)
    spsi, cpsi = sincos(psi * ASEC2RAD)
    rot = @SMatrix [cpsi -spsi*cobm -spsi*sobm
                    spsi*cobt cpsi * cobm * cobt+sobm * sobt cpsi * sobm * cobt-cobm * sobt
                    spsi*sobt cpsi * cobm * sobt-sobm * cobt cpsi * sobm * sobt+cobm * cobt]
    # Transform Position
    if direction == :mean2true
        return rot * pos
    elseif direction == :true2mean
        return rot' * pos
    end
end

"""
    precession(jd_tdb1,pos,jd_tdb2)

Precesses equatorial rectangular coordinates from one epoch to another. One of the two epochs must be J200.0.

# Arguments
- `jd_tdb1::Real`: TDB Julian Date of first epoch
- `jd_tdb2::Real`: TDB Julian Date of second epoch
- `pos::AbstractVector`: Position vector, geocentric equatorial rectangular coordinates, referred to mean dynamical equator and equinox of first epoch
"""
@memoize function precession(jd_tdb1::Real, pos::AbstractVector, jd_tdb2::Real)
    @assert (jd_tdb1 == T0) || (jd_tdb2 == T0) "One of jd_tdb1, jd_tdb2 must be T0"
    eps0 = 84381.406
    # t is time in TDB centuries between the two epochs
    t = (jd_tdb2 - jd_tdb1) / 36525.0
    if jd_tdb2 == T0
        t *= -1
    end

    psia = ((((-0.0000000951 * t + 0.000132851) * t - 0.00114045) * t - 1.0790069) * t +
            5038.481507) * t
    omegaa = ((((+0.0000003337 * t - 0.000000467) * t - 0.00772503) * t + 0.0512623) * t -
              0.025754) * t + eps0
    chia = ((((-0.0000000560 * t + 0.000170663) * t - 0.00121197) * t - 2.3814292) * t +
            10.556403) * t

    eps0 = eps0 * ASEC2RAD
    psia = psia * ASEC2RAD
    omegaa = omegaa * ASEC2RAD
    chia = chia * ASEC2RAD

    sa, ca = sincos(eps0)
    sb, cb = sincos(-psia)
    sc, cc = sincos(-omegaa)
    sd, cd = sincos(chia)

    # Construct precession rotation matrix
    rot = @SMatrix [cd * cb-sb * sd * cc cd * sb * ca + sd * cc * cb * ca-sa * sd * sc cd*sb*sa+sd*cc*cb*sa+ca*sd*sc
                    -sd * cb-sb * cd * cc -sd * sb * ca + cd * cc * cb * ca-sa * cd * sc -sd*sb*sa+cd*cc*cb*sa+ca*cd*sc
                    sb*sc -sc * cb * ca-sa * cc -sc * cb * sa+cc * ca]

    # Perform rotation from epoch to J2000.0
    if jd_tdb2 == T0
        return rot' * pos
    else
        return rot * pos
    end
end

"""
    cel_pole(tjd)

Computes the celestial pole offsets for high-precision applications. Each set of
offsets is a correction to the modeled position of the pole for a specific date, derived
from observations and published by the IERS.

This function differs from the C version, where this returns the corrections instead of mutating global state.

# Arguments
- `tjd::Real`: TDB or TT Julian date for pole offsets
- `type::Symbol=:angular`: Type of pole offset. `:angular` for corrections to angular
    coordinates of modeled pole referred to mean ecliptic of date, that is,
    delta-delta-psi and delta-delta-epsilon. `:positional` for corrections to components
    of modeled pole unit vector referred to GCRS axes, that is, dx and dy.
- `dople1`: Value of celestial pole offset in first coordinate.
- `dpole2`: Value of celestial pole offset is second coordinate.
"""
function cel_pole(tjd::Real, type::Symbol, dpole1::Real, dpole2::Real)
    @assert type ∈ Set([:angular, :positional])
    if type == :angular
        psi_cor = dpole1 * 1e-3
        eps_cor = dpole2 * 1e-3
    elseif type == :positional
        # Components of modeled pole unit vector referred to GCRS axes
        dx = dpole1
        dy = dpole2
        t = (tjd - T0) / 36525.0
        # Compute sine_e of mean obliquity of date
        mean_ob = mean_obliq(tjd)
        sin_e = sin(mean_ob * ASEC2RAD)
        # Trivial model of pole trajectory is GCRS allows computation of dz
        x = (2004.190 * t) * ASEC2RAD
        dz = -(x + 0.5 * x * x * x) * dx
        # Form pole offset vector in GCRS
        dp1 = [dx, dy, dz] .* 1e-3 .* ASEC2RAD
        # Precess pole offset vector to mean equator and equinox of date
        dp2 = frame_tie(dp1, :icrs2dynamic)
        dp3 = precession(T0, dp2, tjd)
        # Compute delta-delta-psi and delta-delta-epsilon in arcseconds
        psi_cor = (dp3[1] / sin_e) / ASEC2RAD
        eps_cor = dp3[2] / ASEC3RAD
    end
    return psi_cor, eps_cor
end

"""
    e_tilt(jd_tdb)

Computes quantities related to the orientation of the Earth's rotation axis at Julian date `jd_tdb`

# Arguments
- `jd_tdb::Real`: TDB Julian Date

# Optional Arguments
- `accuracy::Symbol=:full`: Sets the accuracy level of `:full` or `:reduced`

# Returns
`(mobl,tobl,ee,dpsi,deps)` where
- `mobl`: Mean obliquity of the ecliptic in degrees
- `tobl`: True obliquity of the ecliptic in degrees
- `ee`: Equation of the equinoxes in seconds of time
- `dpsi`: Nutation in longitude on arcseconds
- `deps`: Nutation in obliquity in arcseconds
"""
@memoize function e_tilt(jd_tdb::T; accuracy::Symbol=:full) where {T<:Real}
    # FIXME how do we include the PSI_COR and EPS_COR

    # Compute time in Julian centuries from epoch J2000.0
    t = (jd_tdb - T0) / 36525.0
    # Compute the nutation angles
    dp, de = nutation_angles(t; accuracy=accuracy)
    c_terms = ee_ct(jd_tdb, zero(T); accuracy=accuracy) / ASEC2RAD
    # Apply observed celestial pole offsets FIXME what do we do here?
    d_psi = dp # + PSI_COR
    d_eps = de # + EPS_COR
    # Compute mean obliquity of the ecliptic in arcseconds
    mean_ob = mean_obliq(jd_tdb)
    # Compute mean and true obliquity
    true_ob = mean_ob + d_eps
    # Convert obliquity values to degrees
    mean_ob /= 3600.0
    true_ob /= 3600.0
    # Compute equation of the equinoxes in seconds
    eq_eq = d_psi * cosd(mean_ob) + c_terms
    eq_eq /= 15.0
    # Return computed result
    return (mean_ob, true_ob, eq_eq, d_psi, d_eps)
end

"""
    ira_equinox(jd_tdb)

Computes the intermediate right ascension of the equinox at the input Julian date using an
analytical expression for the accumulated precession in right ascension. For the true equinox,
the result is the equation of the origins.

# Arguments
- `jd_tdb::Real`: TDB Julian Date

# Optional Arguments
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
- `equinox::Symbol=:mean`: Either `:mean` or `:true`
"""
@memoize function ira_equinox(jd_tdb::Real; accuracy::Symbol=:full,
                              equinox::Union{Symbol,Bool}=:mean)
    @assert equinox ∈ Set([:mean, :true])
    # Compute time in Julian centuries
    t = (jd_tdb - T0) / 36525.0
    # For the true equinox, obtain the equation of the equinoxes in seconds,
    # which includes the 'complementary terms'
    if equinox == :true
        _, _, eq_eq, _, _ = e_tilt(jd_tdb; accuracy=accuracy)
    elseif equinox == :mean
        eq_eq = 0.0
    end
    # Compute precession in RA in arcseconds taken from the reference
    prec_ra = 0.014506 +
              ((((-0.0000000368 * t - 0.000029956) * t - 0.00000044) * t + 1.3915817) * t +
               4612.156534) * t
    ra_eq = -(prec_ra / 15.0 + eq_eq) / 3600.0
    return ra_eq
end

"""
   cio_location(jd_tdb)

Computes the location of the celestial intermediate origin (CIO) for a given Julian date, as a
right ascension with respect to the true equinox of date.

# Arguments
- `jd_tdb::Real`: TDB Julian Date

# Optional Arguments
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
"""
@memoize function cio_location(jd_tdb::Real; accuracy::Symbol=:full)
    # NOTE: We are going to calculate the values instead of using the external dep for the time being
    # This results in ref_sys always being 2, true equator and equinox instead of interpolated GCRS
    ra_cio = -ira_equinox(jd_tdb; accuracy=accuracy, equinox=:true)
    return ra_cio
end

"""
    cio_basis(jd_tdb,ra_cio)

Computes the orthonomal basis vectors, with respect to the GCRS of the celestial intermediate system defined
by the celestial intermediate pole and origin.

# Arguments
- `jd_tdb::Real`: TDB Julian Date
- `ra_cio::Real`: Right ascension of the CIO at epoch (hours)

# Optional Arguments
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
"""
@memoize function cio_basis(jd_tdb::Real, ra_cio::Real; accuracy::Symbol=:full)
    # Compute unit vector z towards celestial pole
    z0 = [0.0, 0.0, 1.0]
    w1 = nutation(jd_tdb, z0; accuracy=accuracy, direction=:true2mean)
    w2 = precession(jd_tdb, w1, T0)
    z = frame_tie(w2, :dynamic2icrs)
    # Compute unit vectors x and y
    # First construct unit vector towards CIO in equator and equinox of date system
    w0 = @SVector [cosd(ra_cio * 15), sind(ra_cio * 15), 0]
    # Rotate the vector into the GCRS to form unit vector x
    w1 = nutation(jd_tdb, w0; accuracy=accuracy, direction=:true2mean)
    w2 = precession(jd_tdb, w1, T0)
    x = frame_tie(w2, :dynamic2icrs)
    # Compute unit vector y orthogonal to x and z
    y = z × x
    return x, y, z
end

"""
    era(jd_high)

This function returns the value of the earth rotation angle (theta in degrees) for a given UT1 Julian date.

# Arguments
- `jd_high::Real`: High-order part of UT1 Julian date
- `jd_low::Real=0.0`: Low-order part of UT1 Julian date
"""
function era(jd_high::Real, jd_low::Real=0.0)
    thet1 = 0.7790572732640 + 0.00273781191135448 * (jd_high - T0)
    thet2 = 0.00273781191135448 * jd_low
    thet3 = rem(jd_high, 1.0) + rem(jd_low, 1.0)
    theta = rem(thet1 + thet2 + thet3, 1.0) * 360.0
    if (theta < 0.0)
        theta += 360.0
    end
    return theta
end

"""
    tdb2tt(tdb_jd)

Computes the Terrestrial Time (TT) or Terrestrial Dynamic Time (TDT) Julian date
corresponding to a Barrycentric Dynamical Time (TDB) Julian Date)

# Arguments
- `tdb_jd::Real`: TDB Julian date
"""
function tdb2tt(tdb_jd::Real)
    t = (tdb_jd - T0) / 36525.0
    secdiff = 0.001657 * sin(628.3076 * t + 6.2401) +
              0.000022 * sin(575.3385 * t + 4.2970) +
              0.000014 * sin(1256.6152 * t + 6.1969) +
              0.000005 * sin(606.9777 * t + 4.0212) +
              0.000005 * sin(52.9691 * t + 0.4444) +
              0.000002 * sin(21.3299 * t + 5.5431) +
              0.000010 * t * sin(628.3076 * t + 4.2490)
    tt_jd = tdb_jd - secdiff / 86500.0
    return (tt_jd, secdiff)
end

"""
    sidereal_time(jd_ut1,0,delta_t)

Computes the Greenwich sidereal time, either mean or apparent, at Julian date `jd_high` + `jd_low` in hours.
See Chapter 5 in the NOVASC manual for information between equinox and CIO-based methods.

# Arguments
- `jd_high::Real`: high-order part of UT1 Julian date
- `jd_low::Real=0`: low-order part of UT1 Julian date
- `delta_t::Real=0`: Difference TT-UT1 at `jd_high` + `jd_low` in seconds

# Optional arguments
- `gst_type::Symbol=:mean`: Return results as mean (`:mean`) or apparent (`:apparent`) time
- `method::Symbol=:CIO`: Computation method, `:CIO`-based or `:equinox`-based
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
"""
@memoize function sidereal_time(jd_high::Real, jd_low::Real=0.0, delta_t::Real=0.0;
                                gst_type::Symbol=:mean, method::Symbol=:CIO,
                                accuracy::Symbol=:full)
    @assert method ∈ Set([:equinox, :CIO])
    @assert gst_type ∈ Set([:mean, :apparent])

    # Time argument for precession and nutation components of sidereal time is TDB. First approximation
    # is TDB = TT, then refine
    jd_ut = jd_high + jd_low
    jd_tt = jd_ut + (delta_t / 86400.0)
    jd_tdb = jd_tt
    _, secdiff = tdb2tt(jd_tdb)
    jd_tdb = jd_tt + (secdiff / 86400.0)
    t = (jd_tdb - T0) / 36525.0

    # Compute the Earth Rotation Angle. Time argument is UT1.
    theta = era(jd_high, jd_low)

    # Compute the equation of the equinoxes if needed, depending upon the input values of `gst_type` and `method`.
    # If not needed, set to zero
    if ((gst_type == :mean) && (method == :CIO)) ||
       ((gst_type == :apparent) && (method == :equinox))
        _, _, ee, _, _ = e_tilt(jd_tdb; accuracy=accuracy)
        eqeq = ee * 15
    else
        eqeq = 0
    end

    if method == :CIO
        # Obtain basis vectors, in the GCRS of the celestial intermediate system
        ra_cio = cio_location(jd_tdb; accuracy=accuracy)
        x, y, z = cio_basis(jd_tdb, ra_cio; accuracy=accuracy)
        # Compute the direction of the true equinox in the GCRS
        w1 = nutation(jd_tdb, [1, 0, 0]; direction=:true2mean, accuracy=accuracy)
        w2 = precession(jd_tdb, w1, T0)
        eq = frame_tie(w2, :dynamic2icrs)
        # Compute the hour angle of the equinox w.r.t. the TIO meridian
        ha_eq = theta - atand(eq ⋅ y, eq ⋅ x)
        # For mean, subtract the equation of the equinoxes
        ha_eq -= eqeq / 240
        ha_eq = rem(ha_eq, 360) / 15
        if ha_eq < 0
            ha_eq += 24
        end
        gst = ha_eq
    elseif method == :equinox
        # Precession-in-RA terms in mean sidereal time with coefficients in arcseconds
        st = eqeq +
             0.014506 +
             ((((-0.0000000368 * t - 0.000029956) * t - 0.00000044) * t + 1.3915817) * t +
              4612.156534) * t
        gst = rem((st / 3600.0 + theta), 360.0) / 15.0
        if gst < 0.0
            gst += 24.0
        end
    end
    return gst
end

"""
    wobble(tjd,xp,yp,pos)

Corrects a vector in the ITRS (rotating Earth-fixed system)
for polar motion, and also corrects the longitude origin
(by a tiny amount) to the Terrestrial Intermediate Origin
(TIO).  The ITRS vector is thereby transformed to the terrestrial
intermediate system, based on the true (rotational) equator and
TIO.  Because the true equator is the plane orthogonal to the
direction of the Celestial Intermediate Pole (CIP), the components
of the output vector are referred to z and x axes toward the CIP
and TIO, respectively.

# Arguments
- `tjd::Real`: TT or UT1 Julian Date
- `xp::Real`: Conventionally-defined X coordinate of the CIP in arcseconds
- `yp::Real`: Conventionally-defined Y coordinate of the CIP in arcseconds
- `pos::AbstractVector`: Position vector, geocentric equatorial rectangular coordinates

# Optional Arguments
- `direction::Symbol=:itrs2terr`: Either `:itrs2terr` or `:terr2itrs`
"""
function wobble(tjd::Real, xp::Real, yp::Real, pos::AbstractVector;
                direction::Symbol=:itrs2terr)
    @assert direction ∈ Set([:itrs2terr, :terr2itrs])
    xpole = xp * ASEC2RAD
    ypole = yp * ASEC2RAD
    # Compute approximate longitude of the TIO
    t = (tjd - T0) / 36525.0
    sprime = -47.0e-6 * t
    tiolon = -sprime * ASEC2RAD
    # Compute elements of rotation matrix
    sinx, cosx = sincos(xpole)
    siny, cosy = sincos(ypole)
    sinl, cosl = sincos(tiolon)

    xx = cosx * cosl
    yx = sinx * siny * cosl + cosy * sinl
    zx = -sinx * cosy * cosl + siny * sinl
    xy = -cosx * sinl
    yy = -sinx * siny * sinl + cosy * cosl
    zy = sinx * cosy * sinl + siny * cosl
    xz = sinx
    yz = -cosx * siny
    zz = cosx * cosy

    rot_mat = @SMatrix [xx yx zx; xy yy zy; xz yz zz]

    # Perform rotation
    if direction == :itrs2terr
        return rot_mat * pos
    elseif direction == :terr2itrs
        return rot_mat' * pos
    end
end

"""
    spin(angle,pos)

Transfroms a vector from one coordiate system to another with the same origin and axes rotated about the z-axis.

# Arguments
- `angle::Real`: Angle of rotation, positive counterclockwise when viewed from +z, in degrees
- `pos::AbstractVector`: Position Vector
"""
function spin(angle::Real, pos::AbstractVector)
    sinang, cosang = sincosd(angle)
    xx = cosang
    yx = sinang
    zx = 0.0
    xy = -sinang
    yy = cosang
    zy = 0.0
    xz = 0.0
    yz = 0.0
    zz = 1.0
    rot_mat = @SMatrix [xx yx zx; xy yy zy; xz yz zz]
    return rot_mat * pos
end

"""
    ter2cel(jd_ut_high,jd_ut_low,delta_t,vec)

Rotates a vector from the terrestrial to the celestial system. Specifically, it transforms
a vector in ITRS to GCRS by applying rotations for polar motion, Earth rotation, nutation, precession, and
the dynamical-to-GCRS frame tie.

# Arguments
- `jd_ut_high::Real`: High order part of UT1 Julian date
- `jd_ut_low::Real`: Low order part of UT1 Julian date
- `delta_t::Real`: Value of TT-UT1 at the input date
- `vec::AbstractVector`: Position vector, geocentric equatorial rectangular coordinates, reference axes set by `option`

# Optional Arguments
- `method::Symbol=:CIO`: Computation method, `:CIO`-based or `:equinox`-based
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
- `option::Symbol=:GRCS`: Sets position vector to be referred to `:GCRS` axes or `:equinox` of date and equator
- `xp::Real=0`: x coordinate of the celestial intermediate pole in arcseconds
- `yp::Real=0`: y coordinate of the celestial intermediate pole in arcseconds
"""
function ter2cel(jd_ut_high::Real, jd_ut_low::Real, delta_t::Real, vec::AbstractVector;
                 method::Symbol=:CIO, accuracy::Symbol=:full, option::Symbol=:GCRS,
                 xp::Real=0, yp::Real=0)
    @assert method ∈ Set([:CIO, :equinox])
    @assert option ∈ Set([:GCRS, :equinox])
    # Compute the TT Julian date
    jd_ut1 = jd_ut_high + jd_ut_low
    jd_tt = jd_ut1 + (delta_t / 86400.0)
    # Compute the TDB Julian date corresponding to the input UT1 Julian date
    jd_tdb = jd_tt
    _, secdiff = tdb2tt(jd_tdb)
    jd_tdb = jd_tt + secdiff / 86400.0
    # Apply polar motion
    if iszero(xp) && iszero(yp)
        v1 = vec
    else
        v1 = wobble(jd_tdb, xp, yp, vec; direction=:itrs2terr)
    end
    if method == :CIO
        # Obtain basis vectors, in the GCRS of the celestial intermediate system
        r_cio = cio_location(jd_tdb; accuracy=accuracy)
        x, y, z = cio_basis(jd_tdb, r_cio; accuracy=accuracy)
        # Compute and apply earth rotation angle, `theta` transforming the vector to the celestial intermediate system
        θ = era(jd_ut_high, jd_ut_low)
        v2 = spin(-θ, v1)
        # Transform the vector from the celestial intermediate system to the GCRS
        vec2 = hcat(x, y, z) * v2
    elseif method == :equinox
        gast = sidereal_time(jd_ut_high, jd_ut_low, delta_t; gst_type=:apparent,
                             method=:equinox, accuracy=accuracy)
        v2 = spin(-gast * 15, v1)
        if option == :equinox
            vec2 = v2
        elseif option == :GCRS
            # Apply precession, nutation, frame tie
            v3 = nutation(jd_tdb, v2; accuracy=accuracy, direction=:true2mean)
            v4 = precession(jd_tdb, v3, T0)
            vec2 = frame_tie(v4, :dynamic2icrs)
        end
    end
    return vec2
end

"""
    OneSurface(lat,lon,h,temp,p)

Structure holding the observer's location on the surface of the Earth.
Atmospheric parameters are optional and are only used by `equ2hor` in the
calculation of refraction.

# Arguments
- `latitude::Real`: ITRS latitude in degrees; north positive
- `longitude::Real`: ITRS longitude in degrees; east positive
- `height::Real`: Height of observer (meters)
- `temperature::Real=0`: Temperature (degrees Celsius)
- `pressure::Real=0`: Atmospheric pressure (millibars)
"""
struct OnSurface{T<:Real}
    latitude::T
    longitude::T
    height::T
    temperature::T
    pressure::T
end

OnSurface(lat::T, lon::T, h::T) where {T<:Real} = OnSurface(lat, lon, h, zero(T), zero(T))

"""
    refract(location,zd_obs)

Computes the atmospheric refraction in zenith distance.
This is approximate for optical wavelengths.

# Arguments
- `location::OnSurface`: The observation location
- `zd_obs::Real`: Observed zenith distance, in degrees

# Optional Arguments
- `ref_option::Symbol=:standard`: Either `:standard` to use standard atmospheric conditions, or `:location` to use the conditions included in `location`.
"""
function refract(location::OnSurface, zd_obs::Real; ref_option::Symbol=:standard)
    @assert ref_option ∈ Set([:standard, :location])
    # We only care about zenith distances between 0.1 and 91 degrees
    if 0.1 <= zd_obs <= 91
        if ref_option == :location
            p = location.pressure
            t = location.temperature
        elseif ref_option == :standard
            p = 1010.0 * exp(-location.height / 9.1e3)
            t = 10.0
        end
        h = 90.0 - zd_obs
        r = 0.016667 / tand(h + 7.31 / (h + 4.4))
        refr = r * (0.28 * p / (t + 273.0))
    else
        refr = 0.0
    end
    return refr
end

"""
    equ2hor(jd_ut1,delta_t,ra,dec)

This function transforms topocentric right ascension and
declination to zenith distance and azimuth.  It uses a method
that properly accounts for polar motion, which is significant at
the sub-arcsecond level.  This function can also adjust
coordinates for atmospheric refraction.

`xp` and `yp` are only needed if sub-arcsecond accuracy is needed.

# Arguments
- `jd_ut1::Real`: UT1 Julian Date
- `delta_t::Real`: Difference TT-UT1 ast `jd_ut1` in seconds
- `ra::Real`: Topocentric right ascension in hours, referred to true equator
- `dec::Real`: Topocentric declination in degrees, referred to true equatior
- `location::OnSurface`: Structure containing observer's location details

# Optional arguments
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
- `xp::Real=0`: x coordinate of the celestial intermediate pole in arcseconds
- `yp::Real=0`: y coordinate of the celestial intermediate pole in arcseconds
- `ref_option::Symbol=:none`: Whether to include refraction in the calculations.
Either `:none`, `:standard` to use standard atmospheric conditions, or `:location` to use the conditions
included in `location`.

# Returns
`(zd,az,rar,decr)` where
- `zd`: Topocentric zenith distance in degrees
- `az`: Topocentric azimuth (east of north) in degrees
- `rar`: Topocentric RA of object of interest in hours
- `decr`: Topocentric declination of object of interest in degrees
"""
function equ2hor(jd_ut1::Real, delta_t::Real, ra::Real, dec::Real, location::OnSurface;
                 accuracy::Symbol=:full, xp::Real=0.0, yp::Real=0.0,
                 ref_option::Symbol=:none)
    @assert ref_option ∈ Set([:none, :standard, :location])
    # Initialize things
    refr = 0.0
    rar = ra
    decr = dec
    # Trig it up
    sinlat, coslat = sincosd(location.latitude)
    sinlon, coslon = sincosd(location.longitude)
    sindc, cosdc = sincosd(dec)
    sinra, cosra = sincosd(ra * 15)
    # Set up orthonomarl basis vectors in local Earth-fixed system
    uze = @SVector [coslat * coslon, coslat * sinlon, sinlat]
    # Define vector towards local north in Earth-fixed system
    une = @SVector [-sinlat * coslon, -sinlat * sinlon, coslat]
    # Define vector towards local west in Earth-fixed system
    uwe = @SVector [sinlon, -coslon, 0.0]
    # Obtain vectors in celestial system and rotate earth-fixed orthonormal basis to celestial system
    uz = ter2cel(jd_ut1, 0.0, delta_t, uze; accuracy=accuracy, xp=xp, yp=yp,
                 option=:equinox, method=:equinox)
    un = ter2cel(jd_ut1, 0.0, delta_t, une; accuracy=accuracy, xp=xp, yp=yp,
                 option=:equinox, method=:equinox)
    uw = ter2cel(jd_ut1, 0.0, delta_t, uwe; accuracy=accuracy, xp=xp, yp=yp,
                 option=:equinox, method=:equinox)
    # Define unit vector `p` towards object in celestial system w.r.t equator and equinox of date
    p = [cosdc * cosra, cosdc * sinra, sindc]
    # Compute coordinates of object w.r.t orthonomral basis
    # Compuute projected components of `p` onto rotated Earth-fixed basis vectors
    pz = p[1] * uz[1] + p[2] * uz[2] + p[3] * uz[3]
    pn = p[1] * un[1] + p[2] * un[2] + p[3] * un[3]
    pw = p[1] * uw[1] + p[2] * uw[2] + p[3] * uw[3]
    # Compute azimuth and zenith distance
    proj = sqrt(pn^2 + pw^2)
    if proj > 0
        az = -atand(pw, pn)
    end
    if az < 0
        az += 360
    end
    if az >= 360
        az -= 360
    end
    zd = atand(proj, pz)
    # Apply atmospheric refraction if requested
    if ref_option != :none
        # Get refraction in zenith distance
        # This is iterative because refraction algorithms are always a function of observed zenith
        zd0 = zd
        while true
            zd1 = zd
            refr = refract(location, zd; ref_option=ref_option)
            zd = zd0 - refr
            abs(zd - zd1) > 3.0e-5 || break
        end
        # Apply refraction to celestial coordinates of object
        if (refr > 0) && (zd > 3.0e-4)
            # Shift position vector of object in celestial system to account for refraction
            sinzd, coszd = sincosd(zd)
            sinzd0, coszd0 = sincosd(zd0)
            # Compute refracted position vector
            pr = @. ((p - coszd0 * uz) / sinzd0) * sinzd + uz * coszd
            # Compute refracted ra and dec
            proj = sqrt(pr[1]^2 + pr[2]^2)
            if proj > 0.0
                rar = atand(pr[2], pr[1]) / 15.0
            end
            if rar < 0.0
                rar += 24.0
            end
            if rar >= 24.0
                rar -= 24.0
            end
            decr = atand(pr[3], proj)
        end
    end
    return zd, az, rar, decr
end

struct CatEntry
    starname::String
    catalog::String
    starnumber::Int
    ra::Real
    dec::Real
    promora::Real
    promodec::Real
    parallax::Real
    radialvelocity::Real
end

function show(io::IO, entry::CatEntry)
    println(io, "Name:                  ", entry.starname)
    println(io, "Catalog:               ", entry.catalog)
    println(io, "Number:                ", entry.starnumber)
    println(io, "Right Ascension:       ", entry.ra, " hr")
    println(io, "Declination:           ", entry.dec, " deg")
    println(io, "Proper motion in RA:   ", entry.promora, " mas/yr")
    println(io, "Proper motion in dec:  ", entry.promodec, " mas/yr")
    println(io, "Parallax:              ", entry.parallax, " mas")
    return println(io, "Radial Velocity:       ", entry.radialvelocity, " km/s")
end

"""
A celestial object of interest

# Arguments
- `type::Symbol`: `:major` planet, `:minor` planet, or `:extrasolar`
- `number::Int`: Object number. For `type=:major`, Earth would be `3`, etc.
For `:extrasolar` objects, set to `0` as the object is fully specified in `star`.
- `name::String`: Name of object
- `star::Union{Nothing,CatEntry}`: Astrometric data for any extrasolar celestial object

Convenience constructors supplied for extrasolar sources where you only need to supply the `CatEntry`.
"""
struct Object
    type::Symbol
    number::Int
    name::String
    star::Union{Nothing,CatEntry}
    function Object(type::Symbol, number::Int, name::String,
                    star::Union{Nothing,CatEntry}=nothing)
        @assert type ∈ Set([:major, :minor, :extrasolar])
        @assert !((type == :extrasolar) && isnothing(star)) "Extrasolar sources must have star data"
        @assert !((type == :extrasolar) && !iszero(number)) "Extrasolar number should be zero"
        @assert !(!(type == :extrasolar) && iszero(number)) "Non-Extrasolar number should be nonzero"
        return new(type, number, name, star)
    end
end

Object(type::Symbol, star::CatEntry) = Object(type, 0, star.starname, star)

function show(io::IO, obj::Object)
    println(io, "Type:          ",
            Dict(:major => "Major Planet", :minor => "Minor Planet",
                 :extrasolar => "Extrasolar")[obj.type])
    if !iszero(obj.number)
        println(io, "Number:        ", obj.number)
    end
    println(io, "Name:          ", obj.name)
    if !isnothing(obj.star)
        println(io, "Catalog Entry: ", obj.star.catalog, " ", obj.star.starnumber)
    end
end

"""
The sky position of a celestial object.

# Arguments
- `r̂::AbstractVector`: Unit vector towards object
- `ra::Real`: Apparent right ascension in hours
- `dec::Real`: Apparent declination in degrees
- `dis::Real`: True distance to solar systm body in AU
- `rv::Real`: Radial velocity in km/s
"""
struct SkyPos
    r̂::AbstractVector
    ra::Real
    dec::Real
    dis::Real
    rv::Real
end

"""
The observer's location and velocity on a near-Earth spacecraft

# Arguments
- `sc_pos::AbstractVector`: Geocentric position vector in km
- `sc_vel::AbstractVector`: Geocentric velocity vector (ẋ,ẏ,ż) in km/s
"""
struct InSpace
    sc_pos::AbstractVector
    sc_vel::AbstractVector
end

"""
The location of an observer.
Note, this is quite different from the C version as we will
use the union type to store the location instead of two potentially
nil pointers. As such the names of the arguments differ and is incompatible with C.

# Arguments
- `type::Symbol`: Location is either at `:geocenter`, on the `:surface` of Earth, or in `:space`.
In space implies on a near-earth spacecraft.
- `loc::Union{OnSurface,InSpace}`: Location of observer
"""
struct Observer
    type::Symbol
    loc::Union{OnSurface,InSpace,Nothing}
    function Observer(type::Symbol, loc::Union{OnSurface,InSpace,Nothing}=nothing)
        @assert type ∈ Set([:geocenter, :surface, :space])
        if type == :geocenter
            @assert isnothing(loc) "loc must be empty for geocentric observers"
        elseif type == :surface
            @assert typeof(loc) == OnSurface "Observer on Earth needs location data"
        elseif type == :space
            @assert typeof(loc) == InSpace "Observer in space needs location data"
        end
        return new(type, loc)
    end
end

function show(io::IO, obs::Observer)
    if obs.type == :geocenter
        println(io, "A geocentric observer")
    elseif obs.type == :surface
        println(io, "An observer on earth")
    elseif obs.type == :space
        prinlnt(io, "An observer in a near-Earth spacecraft")
    end
end

"""
    ephemeris(jd,cel_obj)

Retrieves the position and velocity of a solar system body from a fundamental ephemeris.

# Arguments
- `jd::AbstractVector`: TDB Julian date split into two parts, where the sum is the TDB Julian date
- `cel_obj::Object`: The celestial object of interest

# Optional arguments
- `origin::Symbol`: Either the solar system `:barycenter` or the center of mass of the `:sun`
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
"""
function ephemersis(jd::AbstractVector, cel_obj::Object; origin::Symbol=:barycenter,
                    accuracy::Symbol=:full)
    @assert origin ∈ Set([:barycenter, :sun])
    @assert accuracy ∈ Set([:full, :reduced])
    # Acess the ephemeris of the correct type
    if cel_obj.type == :major
        # Get the position and velocity of a major planet, Pluto, Sun, or Moon.
        # When high accuracy is specified, use function 'solarsystem_hp' rather
        # than 'solarsystem'
        if accuracy == :full
            pos, vel = solarsystem_hp(jd, cel_obj.number; origin=origin)
        elseif accuracy == :reduced
            pos, vel = solarsystem(sum(jd), cel_obj.number; origin=origin)
        end
    elseif cel_obj.type == :minor
        # Get the position and velocity of a minor planet
    elseif cel_obj.type == :extrasolar
        @error "Ephemersis only valid for major and minor planets"
    end
end

"""
     place(jd_tt,cel_obj,location,ΔT)

Computes the apparent direction of a star or solar system body at a specified time an in a specified coordinate system.

# Arguments
- `jd_tt::Real`: TT Julian date for place
- `cel_obj::Object`: The celestial object of interest
- `location::Observer`: The location of the observer
- `delta_t::Real`: Difference TT-UT1 at `jd_tt` in seconds

# Optional arguments
- `coord_sys::Symbol=:GCRS` Specifies the coordinate system of the output position.
Either `:GCRS`, `:true` equator and equinox, true equator and `:CIO` of date, or `:astrometric` coordinates.
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy
"""
function place(jd_tt::Real, cel_object::Object, location::Observer, delta_t::Real;
               coord_sys::Symbol=:GCRS, accuracy::Symbol=:full)
    #FIXME
    @assert coord_sys ∈ Set([:GCRS, :true, :CIO, :astrometric])
    # Create Eath and Sun `objects`
    earth = Object(:major, 3, "Earth")
    sun = Object(:major, 10, "Sun")
    if cel_object.type == :major && cel_object.numer == 3 && location.type != :space
        @error "Earth can only be an observed object when `location` is a near-Earth satellite" location
    end
    # Get position and velocity of Earth (geocenter) and Sun
    # Compute `jd_tdb`, the TDB Julian date corresponding to `jd_tt`
    jd_tdb = jd_tt
    x, secdif = tdb2tt(jd_tdb)
    jd_tdb = jd_tt + secdif / 86400.0
    # Get position and velocity of Earth w.r.t. barycenter of solar system in ICRS
    return jd = [jd_tdb, 0.0]
    #ephemeris(jd,earth,)

end

"""
    app_star(jd_tt,star)

Computes the apparent place of a star at date `jd_tt`, givne its catalog
mean place, proper motion, parallax, and radial velocity.

# Arguments
- `jd_tt::Real`: TT Julian date for apparent place
- `star::CatEntry`: Catalog entry structure containing catalog data for the object in the ICRS

# Optional arguments
- `accuracy::Symbol=:full`: Either `:full` or `:reduced` accuracy

# Returns
`(ra,dec)` where
- `ra`: Apparent right ascension in hours
- `dec`: Apparent declination in degrees
"""
function app_star(jd_tt::Real, star::CatEntry; accuracy::Symbol=:full)
    #FIXME NOTEST
    # Set up a structure of type `Object` containing the star data
    # Located outside solar system
    cel_obj = Object(:extrasolar, star)
    # Compute the apparent place with a call to function `place`
    sky_pos = place(jd_tt, cel_obj, Observer(:geocentric), 0.0; coord_sys=:true,
                    accuracy=accuracy)
    # Return ra and dec
    return (sky_pos.ra, sky_pos.dec)
end


text:
"""
    axes_types(::Type{T}, dim)

Returns the axis type along dimension `dim`.
"""
axes_types(x, dim) = axes_types(typeof(x), dim)
@inline axes_types(::Type{T}, dim) where {T} = axes_types(T, to_dims(T, dim))
@inline function axes_types(::Type{T}, dim::StaticInt{D}) where {T,D}
    if D > ndims(T)
        return OptionallyStaticUnitRange{One,One}
    else
        return _get_tuple(axes_types(T), dim)
    end
end
@inline function axes_types(::Type{T}, dim::Int) where {T}
    if dim > ndims(T)
        return OptionallyStaticUnitRange{One,One}
    else
        return axes_types(T).parameters[dim]
    end
end

"""
    axes_types(::Type{T}) -> Type

Returns the type of the axes for `T`
"""
axes_types(x) = axes_types(typeof(x))
axes_types(::Type{T}) where {T<:Array} = Tuple{Vararg{OneTo{Int},ndims(T)}}
function axes_types(::Type{T}) where {T}
    if parent_type(T) <: T
        return Tuple{Vararg{OptionallyStaticUnitRange{One,Int},ndims(T)}}
    else
        return axes_types(parent_type(T))
    end
end
axes_types(::Type{LinearIndices{N,R}}) where {N,R} = R
axes_types(::Type{CartesianIndices{N,R}}) where {N,R} = R
function axes_types(::Type{T}) where {T<:Union{Adjoint,Transpose}}
    P = parent_type(T)
    return Tuple{axes_types(P, static(2)), axes_types(P, static(1))}
end
function axes_types(::Type{T}) where {T<:PermutedDimsArray}
    return eachop_tuple(_get_tuple, to_parent_dims(T), axes_types(parent_type(T)))
end
function axes_types(::Type{T}) where {T<:AbstractRange}
    if known_length(T) === nothing
        return Tuple{OptionallyStaticUnitRange{One,Int}}
    else
        return Tuple{OptionallyStaticUnitRange{One,StaticInt{known_length(T)}}}
    end
end
function axes_types(::Type{T}) where {N,T<:Base.ReshapedArray{<:Any,N}}
    return Tuple{Vararg{OptionallyStaticUnitRange{One,Int},N}}
end

@inline function axes_types(::Type{T}) where {N,P,I,T<:SubArray{<:Any,N,P,I}}
    return eachop_tuple(_sub_axis_type, to_parent_dims(T), T)
end
@inline function _sub_axis_type(::Type{A}, dim::StaticInt) where {T,N,P,I,A<:SubArray{T,N,P,I}}
    return OptionallyStaticUnitRange{
        _int_or_static_int(known_first(axes_types(P, dim))),
        _int_or_static_int(known_length(_get_tuple(I, dim)))
    }
end

function axes_types(::Type{R}) where {T,N,S,A,R<:ReinterpretArray{T,N,S,A}}
    if _is_reshaped(R)
        if sizeof(S) === sizeof(T)
            return axes_types(A)
        elseif sizeof(S) > sizeof(T)
            return eachop_tuple(_reshaped_axis_type, to_parent_dims(R), R)
        else
            return eachop_tuple(axes_types, to_parent_dims(R), A)
        end
    else
        return eachop_tuple(_non_reshaped_axis_type, to_parent_dims(R), R)
    end
end

function _reshaped_axis_type(::Type{R}, dim::StaticInt) where {T,N,S,A,R<:ReinterpretArray{T,N,S,A}}
    return axes_types(parent_type(R), dim)
end
function _reshaped_axis_type(::Type{R}, dim::Zero) where {T,N,S,A,R<:ReinterpretArray{T,N,S,A}}
    return OptionallyStaticUnitRange{One,StaticInt{div(sizeof(S), sizeof(T))}}
end

function _non_reshaped_axis_type(::Type{R}, dim::StaticInt) where {T,N,S,A,R<:ReinterpretArray{T,N,S,A}}
    return axes_types(parent_type(R), dim)
end
function _non_reshaped_axis_type(::Type{R}, dim::One) where {T,N,S,A,R<:ReinterpretArray{T,N,S,A}}
    paxis = axes_types(A, dim)
    len = known_length(paxis)
    if len === nothing
        raxis = OptionallyStaticUnitRange{One,Int}
    else
        raxis = OptionallyStaticUnitRange{One,StaticInt{div(len * sizeof(S), sizeof(T))}}
    end
    return similar_type(paxis, Int, raxis)
end

#=
    similar_type(orignal_type, new_data_type)
=#
similar_type(::Type{OneTo{Int}}, ::Type{Int}, ::Type{OneTo{Int}}) = OneTo{Int}
similar_type(::Type{OneTo{Int}}, ::Type{Int}, ::Type{OptionallyStaticUnitRange{One,Int}}) = OneTo{Int}
similar_type(::Type{OneTo{Int}}, ::Type{Int}, ::Type{OptionallyStaticUnitRange{One,StaticInt{N}}}) where {N} = OptionallyStaticUnitRange{One,StaticInt{N}}

similar_type(::Type{OptionallyStaticUnitRange{One,Int}}, ::Type{Int}, ::Type{OneTo{Int}}) = OptionallyStaticUnitRange{One,Int}
similar_type(::Type{OptionallyStaticUnitRange{One,Int}}, ::Type{Int}, ::Type{OptionallyStaticUnitRange{One,Int}}) = OptionallyStaticUnitRange{One,Int}
similar_type(::Type{OptionallyStaticUnitRange{One,Int}}, ::Type{Int}, ::Type{OptionallyStaticUnitRange{One,StaticInt{N}}}) where {N} = OptionallyStaticUnitRange{One,StaticInt{N}}

similar_type(::Type{OptionallyStaticUnitRange{One,StaticInt{N}}}, ::Type{Int}, ::Type{OneTo{Int}}) where {N} = OptionallyStaticUnitRange{One,Int}
similar_type(::Type{OptionallyStaticUnitRange{One,StaticInt{N}}}, ::Type{Int}, ::Type{OptionallyStaticUnitRange{One,Int}}) where {N} = OptionallyStaticUnitRange{One,Int}
similar_type(::Type{OptionallyStaticUnitRange{One,StaticInt{N1}}}, ::Type{Int}, ::Type{OptionallyStaticUnitRange{One,StaticInt{N2}}}) where {N1,N2} = OptionallyStaticUnitRange{One,StaticInt{N2}}

"""
    axes(A, d)

Return a valid range that maps to each index along dimension `d` of `A`.
"""
axes(a, dim) = axes(a, to_dims(a, dim))
axes(a, dims::Tuple{Vararg{Any,K}}) where {K} = (axes(a, first(dims)), axes(a, tail(dims))...)
axes(a, dims::Tuple{T}) where {T} = (axes(a, first(dims)), )
axes(a, ::Tuple{}) = ()
function axes(a::A, dim::Integer) where {A}
    if parent_type(A) <: A
        return Base.axes(a, Int(dim))
    else
        return axes(parent(a), to_parent_dims(A, dim))
    end
end
function axes(A::CartesianIndices{N}, dim::Integer) where {N}
    if dim > N
        return static(1):static(1)
    else
        return getfield(axes(A), Int(dim))
    end
end
function axes(A::LinearIndices{N}, dim::Integer) where {N}
    if dim > N
        return static(1):static(1)
    else
        return getfield(axes(A), Int(dim))
    end
end

axes(A::SubArray, dim::Integer) = Base.axes(A, Int(dim))  # TODO implement ArrayInterface version
axes(A::ReinterpretArray, dim::Integer) = Base.axes(A, Int(dim))  # TODO implement ArrayInterface version
axes(A::Base.ReshapedArray, dim::Integer) = Base.axes(A, Int(dim))  # TODO implement ArrayInterface version

"""
    axes(A)

Return a tuple of ranges where each range maps to each element along a dimension of `A`.
"""
@inline function axes(a::A) where {A}
    if parent_type(A) <: A
        return Base.axes(a)
    else
        return axes(parent(a))
    end
end
axes(A::PermutedDimsArray) = permute(axes(parent(A)), to_parent_dims(A))
axes(A::Union{Transpose,Adjoint}) = _axes(A, parent(A))
_axes(A::Union{Transpose,Adjoint}, p::AbstractVector) = (One():One(), axes(p, One()))
_axes(A::Union{Transpose,Adjoint}, p::AbstractMatrix) = (axes(p, StaticInt(2)), axes(p, One()))
axes(A::SubArray) = Base.axes(A)  # TODO implement ArrayInterface version
axes(A::ReinterpretArray) = Base.axes(A)  # TODO implement ArrayInterface version
axes(A::Base.ReshapedArray) = Base.axes(A)  # TODO implement ArrayInterface version
axes(A::CartesianIndices) = A.indices
axes(A::LinearIndices) = A.indices

"""
    LazyAxis{N}(parent::AbstractArray)

A lazy representation of `axes(parent, N)`.
"""
struct LazyAxis{N,P} <: AbstractUnitRange{Int}
    parent::P

    LazyAxis{N}(parent::P) where {N,P} = new{N::Int,P}(parent)
    @inline function LazyAxis{:}(parent::P) where {P}
        if ndims(P) === 1
            return new{1,P}(parent)
        else
            return new{:,P}(parent)
        end
    end
end

@inline Base.parent(x::LazyAxis{N,P}) where {N,P} = axes(getfield(x, :parent), static(N))
@inline function Base.parent(x::LazyAxis{:,P}) where {P}
    return eachindex(IndexLinear(), getfield(x, :parent))
end

@inline parent_type(::Type{LazyAxis{N,P}}) where {N,P} = axes_types(P, static(N))
# TODO this approach to parent_type(::Type{LazyAxis{:}}) is a bit hacky. Something like
# LabelledArrays has a linear set of symbolic keys, which could be propagated through
# `to_indices` for key based indexing. However, there currently isn't a good way of handling
# that when the linear indices aren't linearly accessible through a child array (e.g, adjoint)
# For now we just make sure the linear elements are accurate.
parent_type(::Type{LazyAxis{:,P}}) where {P<:Array} = OneTo{Int}
@inline function parent_type(::Type{LazyAxis{:,P}}) where {P}
    if known_length(P) === nothing
        return OptionallyStaticUnitRange{StaticInt{1},Int}
    else
        return OptionallyStaticUnitRange{StaticInt{1},StaticInt{known_length(P)}}
    end
end

Base.keys(x::LazyAxis) = keys(parent(x))

Base.IndexStyle(::Type{T}) where {T<:LazyAxis} = IndexStyle(parent_type(T))

can_change_size(::Type{LazyAxis{N,P}}) where {N,P} = can_change_size(P)

known_first(::Type{T}) where {T<:LazyAxis} = known_first(parent_type(T))

known_length(::Type{LazyAxis{N,P}}) where {N,P} = known_size(P, N)
known_length(::Type{LazyAxis{:,P}}) where {P} = known_length(P)

@inline function known_last(::Type{T}) where {T<:LazyAxis}
    return _lazy_axis_known_last(known_first(T), known_length(T))
end
_lazy_axis_known_last(start::Int, length::Int) = (length + start) - 1
_lazy_axis_known_last(::Any, ::Any) = nothing

@inline function Base.first(x::LazyAxis{N})::Int where {N}
    if known_first(x) === nothing
        return offsets(getfield(x, :parent), static(N))
    else
        return known_first(x)
    end
end
@inline function Base.first(x::LazyAxis{:})::Int
    if known_first(x) === nothing
        return firstindex(getfield(x, :parent))
    else
        return known_first(x)
    end
end

@inline function Base.length(x::LazyAxis{N})::Int where {N}
    if known_length(x) === nothing
        return size(getfield(x, :parent), static(N))
    else
        return known_length(x)
    end
end
@inline function Base.length(x::LazyAxis{:})::Int
    if known_length(x) === nothing
        return lastindex(getfield(x, :parent))
    else
        return known_length(x)
    end
end

@inline function Base.last(x::LazyAxis)::Int
    if known_last(x) === nothing
        if known_first(x) === 1
            return length(x)
        else
            return (static_length(x) + static_first(x)) - 1
        end
    else
        return known_last(x)
    end
end

Base.to_shape(x::LazyAxis) = length(x)

@inline function Base.checkindex(::Type{Bool}, x::LazyAxis, i::Integer)
    if known_first(x) === nothing || known_last(x) === nothing
        return checkindex(Bool, parent(x), i)
    else  # everything is static so we don't have to retrieve the axis
        return (!(known_first(x) > i) || !(known_last(x) < i))
    end
end

@propagate_inbounds function Base.getindex(x::LazyAxis, i::Integer)
    @boundscheck checkindex(Bool, x, i) || throw(BoundsError(x, i))
    return Int(i)
end
@propagate_inbounds Base.getindex(x::LazyAxis, i::StepRange{T}) where {T<:Integer} = parent(x)[i]
@propagate_inbounds Base.getindex(x::LazyAxis, i::AbstractUnitRange{<:Integer}) = parent(x)[i]

Base.show(io::IO, x::LazyAxis{N}) where {N} = print(io, "LazyAxis{$N}($(parent(x))))")

"""
    lazy_axes(x)

Produces a tuple of axes where each axis is constructed lazily. If an axis of `x` is already
constructed or it is simply retrieved.
"""
@generated function lazy_axes(x::X) where {X}
    Expr(:block,
         Expr(:meta, :inline),
         Expr(:tuple, [:(LazyAxis{$dim}(x)) for dim in 1:ndims(X)]...)
    )
end
lazy_axes(x::LinearIndices) = axes(x)
lazy_axes(x::CartesianIndices) = axes(x)
@inline lazy_axes(x::MatAdjTrans) = reverse(lazy_axes(parent(x)))
@inline lazy_axes(x::VecAdjTrans) = (LazyAxis{1}(x), first(lazy_axes(parent(x))))
@inline lazy_axes(x::PermutedDimsArray) = permute(lazy_axes(parent(x)), to_parent_dims(A))

ext: jl


text:
#=

Implement a job scheduler which takes in a function f and an integer n, and calls f after n milliseconds.

=#

function job_scheduler(f::Function, n::Int)
    sleep(n/1000)
    f()
end

function job()
    return "Job Done!"
end

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
@reexport module _NumberTheory

import Base: ==
import Nemo
import Primes

using Base.Iterators: cycle
using Brainstorm: flatten, takewhile
using Brainstorm._Math: checked_add, checked_sub
using AutoHashEquals: @auto_hash_equals
using DataStructures: SortedDict
using IterTools: chain, drop, imap, take
using Lists: List, ListNode
using Nemo: eulerphi, factor, fmpz, isprime, is_prime, sigma, sqrtmod
using Pipe: @pipe
using Primes: ismersenneprime, isrieselprime, primes, primesmask
using Reexport: @reexport

include("primes.jl")
include("factor.jl")
include("mod.jl")
include("decimal.jl")
include("fibonacci.jl")
include("hailstone.jl")
include("triangle.jl")
include("pentagon.jl")
include("hexagon.jl")
include("confrac.jl")
include("diophantine.jl")

end


text:
# input population is a simplified version of
# PopulationHawkes. Does not connect to anything

# I need to abstract the type, because there is a variant to test weights
abstract type AbstractPopulationInput <: AbstractPopulation end

struct PopulationInput{PS<:PopulationState} <: AbstractPopulationInput
  state::PS
  spike_proposals::Vector{Float64} # memory allocation
  function PopulationInput(ps::PS) where PS<:PopulationState
    new{PS}(ps,fill(0.0,nneurons(ps)))
  end
end

abstract type SpikeGenerator end
struct InputUnit{SG<:SpikeGenerator} <: UnitType
  spike_generator::SG
end
@inline function compute_next_spike(t_now::Real,pop::AbstractPopulationInput,ineu::Integer)
  return compute_next_spike(t_now,pop.state.unittype.spike_generator,ineu)
end

# simplified constructor for population object
# (because it has no inputs or incoming weights, so it is pretty much a
# wrapper around its population state)
function PopulationInput(sg::SpikeGenerator,n::Integer)
  return PopulationInput(PopulationState(InputUnit(sg),n))
end

############
# A population with dummy weights, to test plasticity functions.
# N.B. those weights have no influence whatsoever on the dynamics!
struct PopulationInputTestWeights{N,PS<:PopulationState,
    TC<:NTuple{N,Connection},
    TP<:NTuple{N,PopulationState}} <: AbstractPopulationInput
  state::PS
  connections::TC
  pre_states::TP
  spike_proposals::Vector{Float64} # memory allocation
end
# and constructors
function PopulationInputTestWeights(state::PopulationState,
    (conn_pre::Tuple{C,PS} where {C<:Connection,PS<:PopulationState})...)
  connections = Tuple(getindex.(conn_pre,1))
  pre_states = Tuple(getindex.(conn_pre,2))
  spike_proposals = fill(-Inf,nneurons(state))
  return PopulationInputTestWeights(state,connections,pre_states,spike_proposals)
end
# one population only!
function PopulationInputTestWeights(state::PopulationState,conn::Connection)
  return PopulationInputTestWeights(state,(conn,state))
end

##############
# Different types of inputs :

##############
# Poisson with given rate
struct SGPoisson <: SpikeGenerator
  rates::Vector{Float64}
end
@inline function compute_next_spike(t_now::Float64,sg::SGPoisson,idxneu::Integer)
  Δt = -log(rand())/sg.rates[idxneu]
  return t_now+Δt
end

##############
# Spiketrains set externally
struct SGTrains <: SpikeGenerator
  trains::Vector{Vector{Float64}}
end
@inline function compute_next_spike(t_now::Float64,sg::SGTrains,idxneu::Integer)
  train = sg.trains[idxneu]
  t_now_plus = t_now+eps(1000.0) # add a small epsilon to move to next element
  idx = searchsortedfirst(train,t_now_plus)
  if idx > length(train)
    return Inf
  else
    return train[idx]
  end
end

##############
# Poisson varying according to a function
# for exact spike generation, a non-decreasing upper limit is needed, too
struct SGPoissonFunction <: SpikeGenerator
  ratefunction::Function # (t::Float64,idx::Int64) -> rate::Float64
  ratefunction_upper::Function # (t::Float64,idx::Int64) -> rate::Float64
end

# Thinning algorith, e.g.  Laub,Taimre,Pollet 2015
# pretty much the same as in `spike_generation_hawkes` ...
function _rand_by_thinning(t_start::Real,idx_neu::Integer,
    get_rate::Function,get_rate_upper::Function;
    Tmax=50.0,nowarning::Bool=false)
  t = t_start
  while (t-t_start)<Tmax # Tmax is upper limit, if rate too low
    (rup::Float64) = get_rate_upper(t,idx_neu)
    Δt = -log(rand())/rup # rand(Exponential())/rup
    t = t+Δt
    u = rand()*rup # rand(Uniform(0.0,rup))
    (_r::Float64) = get_rate(t,idx_neu)
    if u <= _r
      return t
    end
  end
  # if we are above Tmax, just return upper limit
  if !nowarning
    @warn "Upper limit reached, input firing below $(inv(Tmax)) Hz"
  end
  return Tmax + t_start
end
function compute_next_spike(t_now::Float64,sg::SGPoissonFunction,idx_neu::Integer)
  return _rand_by_thinning(t_now,idx_neu,sg.ratefunction,sg.ratefunction_upper)
end

# compatibility with ExpKernel stuff
# let's just do nothing, because inputs never depend on past states
function burn_spike!(::Real,ps::PopulationState,::Integer)
  typeassert(ps.unittype,InputUnit)
  return nothing
end


text:
using TheFix, Test

import TheFix.cleanse

A = ones(2, 3)
x = 1.0

TheFix.@safeword fix true
TheFix.@safeword FIX false

@testset "Fixed Errors" begin
    @testset "UndefVarError" begin
        @test cleanse(:xx; safe = false) == :x
        @test cleanse(:(exo(sim(1))); safe = false) == :(exp(sin(1)))
        @test cleanse(:(factorrial(factoreal(2))); safe = false) == :(factorial(factorial(2)))
    end
    @testset "DomainError" begin
        @test cleanse(:(sqrt(-1.0))) == :(sqrt(Complex(-1.0)))
        @test cleanse(:(sqrt.((-1.0, 2.0)))) == :(sqrt.(Complex.((-1.0, 2.0))))
        @test cleanse(:(2^(1-3))) == :(2 ^ float(1 - 3))
        @test cleanse(:(2 .^ (-2:2))) == :(2 .^ float.(-2:2))
        @test cleanse(:(2 .^ (-2, -1, 0, 1, 2))) == :(2 .^ float.((-2, -1, 0, 1, 2)))
        @test cleanse(:(log2(-2.0))) == :(log2(Complex(-2.0)))
        @test cleanse(:([2 1; 1 0]^(1 - 3))) == :([2 1; 1 0] ^ float(1 - 3))
        @test cleanse(:(exponent(0.0))) == :(exponent(0.0))
    end
    @testset "OverflowError" begin
        @test cleanse(:(factorial(21))) == :(factorial(big(21)))
        @test cleanse(:(factorrial(factoreal(4))); safe = false) == :(factorial(big(factorial(4))))
        @test cleanse(:(factorial.((21, 22)))) == :(factorial.(big.((21, 22))))
        @test cleanse(:(binomial(67, 30))) == :(binomial(widen(67), widen(30)))
        @test cleanse(:(gcd(typemin(Int), typemin(Int)))) == :(gcd(widen(typemin(Int)), widen(typemin(Int))))
        @test cleanse(:(gcd([typemin(Int), 1]))) == :(gcd(widen.([typemin(Int), 1])))
        @test cleanse(:(gcd.((typemin(Int), 1), (typemin(Int), 2)))) == :(gcd.(widen.((typemin(Int), 1)), widen.((typemin(Int), 2))))
        @test cleanse(:(gcd(typemin(Int), Int(0)))) == :(gcd(widen(typemin(Int)), widen(Int(0))))
        s = x -> 1102938470918723 + x*(102394812390847 + x*12349812309487)
        @test cleanse(:($s(1323457//20345))) == :($s(widen(1323457 // 20345)))
        @test cleanse(:($s.([1323457//20345, 1323458//20345]))) == :($s.(widen.([1323457 // 20345, 1323458 // 20345])))
        @test cleanse(:(-Rational{UInt}(1, 2))) == :(-(Rational{UInt}(1, 2)))
    end
    @testset "DivideError" begin
        @test cleanse(:(2 ÷ 0)) == :(2 / 0)
        @test cleanse(:(div(2, 0, RoundUp))) == :(2 / 0)
        @test cleanse(:(2 % 0)) == :(2 - 0 * (2 / 0))
        @test cleanse(:(mod(2, 0))) == :(2 - 0 * (2 / 0))
        @test cleanse(:(rem(2, 0, RoundDown))) == :(2 - 0 * round(2 / 0, RoundDown))
        @test cleanse(:(fld(2, 0))) == :(2 / 0)
        @test cleanse(:(divrem(2, 0))) == :((2 / 0, 2 - 0 * (2 / 0)))
        @test cleanse(:(divrem(2, 0, RoundNearest))) == :((2 / 0, 2 - 0 * round(2 / 0, RoundNearest)))
    end
    @testset "BoundsError" begin
        @test cleanse(:(A[7])) == :(A[clamp(7, extrema(eachindex(A))...)])
        @test cleanse(:(B[1:7]); safe = false) == :(A[intersect(1:7, eachindex(A))])
        @test cleanse(:(D[3, 3]); safe = false) == :(A[clamp(3, extrema(axes(A, 1))...), 3])
        @test cleanse(:(E[3, 4]); safe = false) == :(A[clamp(3, extrema(axes(A, 1))...), clamp(4, extrema(axes(A, 2))...)])
        @test cleanse(:(F[1:5, 1:5]); safe = false) == :(A[intersect(1:5, axes(A, 1)), intersect(1:5, axes(A, 2))])
    end
    @testset "Code blocks" begin
        @test cleanse(:(for k in 1:2 printline(k) end); safe = false) == :(for k in 1:2 println(k) end)
        @test cleanse(:(for k in 1:2
            printline(k)
        end); safe = false) == :(for k in 1:2
            println(k)
        end)
        @test cleanse(:(begin x = logarithm(4); y = sine(x) end); safe = false) == :(begin x = log(4); y = sin(x) end)
        @test cleanse(:(begin
            i = 1
            while i < 10
                global t = expp1(i)
                printline("t = $t")
                global i += 1
            end
        end); safe = false) == :(begin
            i = 1
            while i < 10
                global t = expm1(i)
                println("t = $t")
                global i += 1
            end
        end)
    end
end


text:
using LayerDicts
using Test

# Write your own tests here.
@testset "LayerDicts" begin

@testset "Constructors" begin
    # Parameters matching
    @test LayerDict(Dict{Int, Int}(), Dict{Int, Int}()) isa LayerDict{Int, Int}
    @test LayerDict((Dict{Int, Int}(), Dict{Int, Int}())) isa LayerDict{Int, Int}
    @test LayerDict([Dict{Int, Int}(), Dict{Int, Int}()]) isa LayerDict{Int, Int}

    # One parameter matching
    @test LayerDict(Dict{Int, Int}(), Dict{String, Int}()) isa LayerDict{Any, Int}
    @test LayerDict((Dict{Int, Int}(), Dict{String, Int}())) isa LayerDict{Any, Int}
    @test LayerDict([Dict{Int, Int}(), Dict{String, Int}()]) isa LayerDict{Any, Int}

    # One parameter promoting
    @test LayerDict(Dict{Int, Int}(), Dict{Float64, Int}()) isa LayerDict{Real, Int}
    @test LayerDict((Dict{Int, Int}(), Dict{Float64, Int}())) isa LayerDict{Real, Int}
    @test LayerDict([Dict{Int, Int}(), Dict{Float64, Int}()]) isa LayerDict{Real, Int}

    # Promoting to Any
    @test LayerDict(Dict{Any, Int}(), Dict{Float64, Any}(), Dict{Float64, Int}()) isa LayerDict{Any, Any}
    @test LayerDict((Dict{Any, Int}(), Dict{Float64, Any}(), Dict{Float64, Int}())) isa LayerDict{Any, Any}
    @test LayerDict([Dict{Any, Int}(), Dict{Float64, Any}(), Dict{Float64, Int}()]) isa LayerDict{Any, Any}

    # No dicts
    @test LayerDict() isa LayerDict{Any, Any}
    @test LayerDict(()) isa LayerDict{Any, Any}
    @test LayerDict(AbstractDict[]) isa LayerDict{Any, Any}
end

dict1 = Dict{Symbol, Int}(:foo => 1, :bar => 1)
dict2 = Dict{Symbol, Int}()
dict3 = Dict{Symbol, Int}(:bar => 3, :baz => 3)

ld = LayerDict{Symbol, Int}([dict1, dict2, dict3])

@testset "Getters" begin
    @testset "getindex" begin
        @test ld[:foo] == 1
        @test ld[:bar] == 1
        @test ld[:baz] == 3
        @test_throws KeyError ld[:quuz]
    end

    @testset "get (default value)" begin
        default = 4
        @test get(ld, :foo, default) == 1
        @test get(ld, :bar, default) == 1
        @test get(ld, :baz, default) == 3
        @test get(ld, :quuz, default) == 4
    end

    @testset "get (default function)" begin
        default = () -> 4
        @test get(default, ld, :foo) == 1
        @test get(default, ld, :bar) == 1
        @test get(default, ld, :baz) == 3
        @test get(default, ld, :quuz) == 4
    end

    @testset "haskey" begin
        @test haskey(ld, :foo)
        @test haskey(ld, :bar)
        @test haskey(ld, :baz)
        @test !haskey(ld, :quuz)
    end

    @testset "in" begin
        @test (:foo => 1) in ld
        @test (:bar => 1) in ld
        @test (:baz => 3) in ld
        @test !((:foo => 2) in ld)
        @test !((:bar => 3) in ld)
        @test !((:quuz => 4) in ld)
    end
end

@testset "Iterators" begin
    @testset "Pairs" begin
        pair_vec = collect(ld)
        @test length(pair_vec) == 3
        pair_set = Set(pair_vec)
        @test pair_set == Set([:foo => 1, :bar => 1, :baz => 3])
    end

    @testset "Keys" begin
        key_vec = collect(keys(ld))
        @test length(key_vec) == 3
        key_set = Set(key_vec)
        @test key_set == Set([:foo, :bar, :baz])
    end

    @testset "Values" begin
        val_vec = sort!(collect(values(ld)))
        @test val_vec == [1, 1, 3]
    end
end

end


text:
# mapping formula from a point (xi, eta) in reference space [-1,1]^2 to a point (x,y)
# in physical coordinate space for a quadrilateral element with straight sides
#     Alg. 95 from the blue book of Kopriva
function straight_side_quad_map(xi, eta, corner_points)

  x = (0.25 * (  corner_points[1,1] * (1.0 - xi) * (1.0 - eta)
               + corner_points[2,1] * (1.0 + xi) * (1.0 - eta)
               + corner_points[3,1] * (1.0 + xi) * (1.0 + eta)
               + corner_points[4,1] * (1.0 - xi) * (1.0 + eta)) )

  y = (0.25 * (  corner_points[1,2] * (1.0 - xi) * (1.0 - eta)
               + corner_points[2,2] * (1.0 + xi) * (1.0 - eta)
               + corner_points[3,2] * (1.0 + xi) * (1.0 + eta)
               + corner_points[4,2] * (1.0 - xi) * (1.0 + eta)) )

  return x, y
end

# Compute the metric terms for the straight sided quadrilateral mapping
#     Alg. 100 from the blue book of Kopriva
function straight_side_quad_map_metrics(xi, eta, corner_points)

  X_xi  = ( 0.25 * (  (1.0 - eta) * (corner_points[2,1] - corner_points[1,1])
                    + (1.0 + eta) * (corner_points[3,1] - corner_points[4,1])) )

  X_eta = ( 0.25 * (  (1.0 - xi) * (corner_points[4,1] - corner_points[1,1])
                    + (1.0 + xi) * (corner_points[3,1] - corner_points[2,1])) )

  Y_xi  = ( 0.25 * (  (1.0 - eta) * (corner_points[2,2] - corner_points[1,2])
                    + (1.0 + eta) * (corner_points[3,2] - corner_points[4,2])) )

  Y_eta = ( 0.25 * (  (1.0 - xi) * (corner_points[4,2] - corner_points[1,2])
                    + (1.0 + xi) * (corner_points[3,2] - corner_points[2,2])) )

  return X_xi, X_eta, Y_xi, Y_eta
end

# construct the (x,y) node coordinates in the volume of a straight sided element
function calc_node_coordinates!(node_coordinates, element, nodes, corners)

  for j in eachindex(nodes), i in eachindex(nodes)
    node_coordinates[:, i ,j ,element] .= straight_side_quad_map(nodes[i], nodes[j], corners)
  end

  return node_coordinates
end

# construct the metric terms for a straight sided element
function calc_metric_terms!(X_xi, X_eta, Y_xi, Y_eta, element, nodes, corners)

  for j in eachindex(nodes), i in eachindex(nodes)
    (X_xi[i, j, element],
     X_eta[i, j, element],
     Y_xi[i, j, element],
     Y_eta[i, j, element]) = straight_side_quad_map_metrics(nodes[i], nodes[j], corners)
  end

  return X_xi, X_eta, Y_xi, Y_eta
end

# construct the normals, their normailzation scalings for a straight sided element
function calc_normals_and_scaling!(normals, scaling, element, nodes, corners)

  # normals and boundary information for the left (local side 4) and right (local side 2)
  for j in eachindex(nodes)
    # side 2
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(1.0, nodes[j], corners)
    Jtemp = X_xi * Y_eta - X_eta * Y_xi
    scaling[j, 2, element]    = sqrt(Y_eta^2 + X_eta^2)
    normals[1, j, 2, element] = sign(Jtemp) * ( Y_eta / scaling[j, 2, element])
    normals[2, j, 2, element] = sign(Jtemp) * (-X_eta / scaling[j, 2, element])

    # side 4
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(-1.0, nodes[j], corners)
    Jtemp =  X_xi * Y_eta - X_eta * Y_xi
    scaling[j, 4, element]    =  sqrt(Y_eta^2 + X_eta^2)
    normals[1, j, 4, element] = -sign(Jtemp) * ( Y_eta / scaling[j, 4, element])
    normals[2, j, 4, element] = -sign(Jtemp) * (-X_eta / scaling[j, 4, element])
  end

  # normals and boundary information for the top (local side 3) and bottom (local side 1)
  for i in eachindex(nodes)
    # side 1
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], -1.0, corners)
    Jtemp =  X_xi * Y_eta - X_eta * Y_xi
    scaling[i, 1, element]    =  sqrt(Y_xi^2 + X_xi^2)
    normals[1, i, 1, element] = -sign(Jtemp) * (-Y_xi / scaling[i, 1, element])
    normals[2, i, 1, element] = -sign(Jtemp) * ( X_xi / scaling[i, 1, element])

    # side 3
    X_xi, X_eta, Y_xi, Y_eta = straight_side_quad_map_metrics(nodes[i], 1.0, corners)
    Jtemp = X_xi * Y_eta - X_eta * Y_xi
    scaling[i, 3, element]    = sqrt(Y_xi^2 + X_xi^2)
    normals[1, i, 3, element] = sign(Jtemp) * (-Y_xi / scaling[i, 3, element])
    normals[2, i, 3, element] = sign(Jtemp) * ( X_xi / scaling[i, 3, element])
  end

  return normals, scaling
end


text:
#Right now, we have some emacs keybindings for moving forward and back, (ctrl-f, ctrl-b, ctrl-a, ctrl-e)
# as well as ctrl-k, ctrl-j to get previous/next commands from command history

function after_search(b::Buffer)
    setmode(b, b.mode[2])
    after(b)
end

function search(b::Buffer)
    #This can be done faster for (very) large files,
    #also, no support for regex search
    txt = [unirange(b.text[y(b)], x(b)+1:length(b.text[y(b)]));
           b.text[y(b)+1:end]]

    inds = Base.search.(txt, b.state[:search])
    yoffset = findfirst(x->length(x)!=0, inds)
    if yoffset == 1
        b.cursor.pos .= [y(b)+yoffset-1, x(b)+inds[yoffset][1]]
        after(self)
    elseif yoffset>1
        b.cursor.pos .= [y(b)+yoffset-1, inds[yoffset][1]]
        after(self)
    else
        b.state[:log] = "NOT FOUND"
        escape(self)
    end
    b.state[:search_history] *= "\n"*b.state[:search]
    b.state[:search_n] = "$(parse(b.state[:search_n])+1)"
    b.state[:search_ind] = b.state[:search_n]
end

search_fn(c) = exec_fn(c, :search, search)
search_mode = Mode("search", Action(search_fn, x->true))


text:
#=
Grammar:
<exp>    ::= <term><exp'>
<exp'>   ::= +<term><exp'> | - <term><exp'> | ϵ
<term>   ::= <factor><term'>
<term'>  ::= *<factor><term'> | / <factor><term'> | ϵ
<factor> ::= (<exp>) |<ident>|<num>
<ident>  ::= A|B|...|Z
<num>    ::  Real, let's assume all Float64 numbers
=#

#=
   Recursive descent parser.
   We do LL(1). That is left-to-right leftmost derivation
=#

#=
 Let's construct one method for each non terminal.
 Let's assume that we have the method scan that updates our global token.
=#

global identLevel = 0
global current_token = ""

function print_and_ident(str::String)
  ident = ""
  for i in 1:identLevel
    ident *= "\t"
  end
  println("$ident $str")
end

function scan()
  res = read(stdin, Char)
  global current_token = string(res)
  return string(res)
end

function main()
  scan()
  res = 0
  res = exp(res)
  println("Result of res is: $res")
end

function exp(res)
  res = term(res)
  res = expPrime(res)
end

function term(res)
  res = factor(res)
  termPrime(res)
end

function expPrime(res)
  if current_token == "+"
    scan()
    res = res + term(res)
    expPrime(res)
  elseif current_token == "-"
    scan()
    res = res - term(res)
    expPrime(res)
  else #= ϵ =#
    return res
  end
end

function termPrime(res)
  if current_token == "*"
    scan()
    res = res * factor(res)
    termPrime(res)
  elseif current_token == "/"
    scan()
    res = res / factor(res)
    termPrime(res)
  else #= ϵ =#
    return res
  end
end

function factor(res)
  if current_token == "("
    scan()
    res = exp(res)
    if current_token != ")"
      throw("Right parenthesis expected")
    end
    scan()
    return res
  elseif nothing != match(r"[A-Z]+", current_token)
    res = parse(Float64, current_token)
    scan()
    return res
  elseif nothing != match(r"[0-9]+", current_token)
    res = parse(Float64, current_token)
    scan()
    return res
  else #= Faulty input =#
    throw("Unexpected token: $(current_token)")
  end
end

main()


text:
# This file contains the line dislocation from Paul Segall's book

export disp_antiplane_seg2, stress_antiplane_seg2
export disp_inplane_seg2, stress_inplane_seg2

function disp_antiplane_seg2(x1::T, x2::T, y1::T, d1::T, d2::T, s::T=one(T)) where T
    @assert abs(d1) > abs(d2) "d₁ denotes a deeper dislocation than d₂"
    -s / 2π * (
        atan((x1 - y1) / (x2 + d1)) - atan((x1 - y1) / (x2 - d1)) -
        atan((x1 - y1) / (x2 + d2)) + atan((x1 - y1) / (x2 - d2))
    )
end

function stress_antiplane_seg2(x1::T, x2::T, y1::T, d1::T, d2::T, G::T, s::T=one(T)) where T
    s13 = -s / 2π * G * (
        (d1 + x2)/((d1 + x2)^2 + (x1 - y1)^2) - (-(d1 - x2))/((d1 - x2)^2 + (x1 - y1)^2) -
        (d2 + x2)/((d2 + x2)^2 + (x1 - y1)^2) + (-(d2 - x2))/((d2 - x2)^2 + (x1 - y1)^2)
    )
    s23 = -s / 2π * G * (
        (-(x1 - y1))/((d1 + x2)^2 + (x1 - y1)^2) - (-(x1 - y1))/((d1 - x2)^2 + (x1 - y1)^2) -
        (-(x1 - y1))/((d2 + x2)^2 + (x1 - y1)^2) + (-(x1 - y1))/((d2 - x2)^2 + (x1 - y1)^2)
    )
    return s13, s23
end

function _disp_half_infinite_line_inplane(x1::T, x2::T, ξ1::T, ξ2::T, s1::T, s2::T, ν::T) where T
    θ1 = atan((x1 - ξ1) / (x2 - ξ2))
    θ2 = atan((x1 - ξ1) / (x2 + ξ2))
    r₁ = hypot(x1 - ξ1, x2 - ξ2)
    r₂ = hypot(x1 - ξ1, x2 + ξ2)
    r₁² = r₁ ^ 2
    r₂² = r₂ ^ 2
    r₁⁴ = r₁² ^ 2
    r₂⁴ = r₂² ^ 2

    u1 = (
        -s1/π/(1-ν) * ( (1-ν)/2*(θ2-θ1) + (x1-ξ1)*(x2-ξ2)/4r₁² - (x1-ξ1)*(x2+(3-4ν)*ξ2)/4r₂² + ξ2*x2*(x1-ξ1)*(x2+ξ2)/r₂⁴ ) +
        s2/π/(1-ν) * ( (1-2ν)/4*log(r₂/r₁) - (x2-ξ2)^2/4r₁² + (x2^2+ξ2^2-4(1-ν)*ξ2*(x2+ξ2))/4r₂² + x2*ξ2*(x2+ξ2)^2/r₂⁴ )
    )
    u2 = (
        -s1/π/(1-ν) * ( (1-2ν)/4*log(r₂/r₁) + (x2-ξ2)^2/4r₁² - ((x2+ξ2)^2-2*ξ2^2-2(1-2ν)*ξ2*(x2+ξ2))/4r₂² + x2*ξ2*(x2+ξ2)^2/r₂⁴ ) +
        s2/π/(1-ν) * ( (1-ν)/2*(θ1-θ2) + (x1-ξ1)*(x2-ξ2)/4r₁² - (x1-ξ1)*(x2+(3-4ν)*ξ2)/4r₂² - ξ2*x2*(x1-ξ1)*(x2+ξ2)/r₂⁴ )
    )
    return u1, u2
end

function disp_inplane_seg2(x1::T, x2::T, ξ1::T, ξ2::T, ζ1::T, ζ2::T, s1::T, s2::T, ν::T) where T
    @assert abs(ξ2) > abs(ζ2) "ξ₂ denotes a deeper dislocation than ζ₂"
    u1ξ, u2ξ = _disp_half_infinite_line_inplane(x1, x2, ξ1, ξ2, -s1, -s2, ν)
    u1ζ, u2ζ = _disp_half_infinite_line_inplane(x1, x2, ζ1, ζ2, s1, s2, ν)
    return u1ξ+u1ζ, u2ξ+u2ζ
end

function _stress_half_infinite_line_inplane(x1::T, x2::T, ξ1::T, ξ2::T, s1::T, s2::T, μ::T, ν::T) where T
    r₁ = hypot(x1 - ξ1, x2 - ξ2)
    r₂ = hypot(x1 - ξ1, x2 + ξ2)
    r₁² = r₁ ^ 2
    r₂² = r₂ ^ 2
    r₁⁴ = r₁² ^ 2
    r₂⁴ = r₂² ^ 2
    r₁⁶ = r₁² * r₁⁴
    r₂⁶ = r₂² * r₂⁴

    σ11 = (
        μ*s2/2π/(1-ν) * ( (x1-ξ1)*((x2-ξ2)^2-(x1-ξ1)^2)/r₁⁴ - (x1-ξ1)*((x2+ξ2)^2-(x1-ξ1)^2)/r₂⁴ + 4ξ2*(x1-ξ1)/r₂⁶*((2ξ2-x2)*(x2+ξ2)^2+(3x2+2ξ2)*(x1-ξ1)^2) ) +
        μ*s1/2π/(1-ν) * ( (x2-ξ2)*((x2-ξ2)^2+3(x1-ξ1)^2)/r₁⁴ - (x2+ξ2)*((x2+ξ2)^2+3(x1-ξ1)^2)/r₂⁴ + 2ξ2/r₂⁶*(6x2*(x2+ξ2)*(x1-ξ1)^2-(x2-ξ2)*(x2+ξ2)^3-(x1-ξ1)^4) )
    )

    σ22 = (
        -μ*s2/2π/(1-ν) * ( (x1-ξ1)*(3(x2-ξ2)^2+(x1-ξ1)^2)/r₁⁴ - (x1-ξ1)*(3(x2+ξ2)^2+(x1-ξ1)^2)/r₂⁴ - 4ξ2*x2*(x1-ξ1)/r₂⁶*(3(x2+ξ2)^2-(x1-ξ1)^2) ) +
        μ*s1/2π/(1-ν) * ( (x2-ξ2)*((x2-ξ2)^2-(x1-ξ1)^2)/r₁⁴ - (x2+ξ2)*((x2+ξ2)^2-(x1-ξ1)^2)/r₂⁴ - 2ξ2/r₂⁶*(6x2*(x2+ξ2)*(x1-ξ1)^2-(3x2+ξ2)*(x2+ξ2)^3+(x1-ξ1)^4)  )
    )

    σ12 = (
        μ*s2/2π/(1-ν) * ( (x2-ξ2)*((x2-ξ2)^2-(x1-ξ1)^2)/r₁⁴ - (x2+ξ2)*((x2+ξ2)^2-(x1-ξ1)^2)/r₂⁴ + 2ξ2/r₂⁶*(6x2*(x2+ξ2)*(x1-ξ1)^2-(x1-ξ1)^4+(ξ2-x2)*(x2+ξ2)^3) ) +
        μ*s1/2π/(1-ν) * ( (x1-ξ1)*((x2-ξ2)^2-(x1-ξ1)^2)/r₁⁴ - (x1-ξ1)*((x2+ξ2)^2-(x1-ξ1)^2)/r₂⁴ + 4ξ2*x2*(x1-ξ1)/r₂⁶*(3(x2+ξ2)^2-(x1-ξ1)^2) )
    )

    return σ11, σ22, σ12
end

function stress_inplane_seg2(x1::T, x2::T, ξ1::T, ξ2::T, ζ1::T, ζ2::T, s1::T, s2::T, G::T, ν::T) where T
    @assert abs(ξ2) > abs(ζ2) "ξ₂ denotes a deeper dislocation than ζ₂"
    σ11ξ, σ22ξ, σ12ξ = _stress_half_infinite_line_inplane(x1, x2, ξ1, ξ2, -s1, -s2, G, ν)
    σ11ζ, σ22ζ, σ12ζ = _stress_half_infinite_line_inplane(x1, x2, ζ1, ζ2, s1, s2, G, ν)
    return σ11ξ+σ11ζ, σ22ξ+σ22ζ, σ12ξ+σ12ζ
end


text:
# A bright soliton in 1D.

push!(LOAD_PATH, "../src")

using Fluxion

t = UnknownDimension("t")
x = BasicDimension("x", RegularGrid(0, 20, 100))
psi = UnknownField("psi", [x, t])

eq = equal(diff(psi, t), 0.5im * diff(psi, x, x) + (1im * abs(psi)^2 - 0.5im) * psi)

function xdensity(psi, t)
    abs(psi)^2
end

k = fourier_space(x)
function kdensity(psi, t)
    psi_k = fourier(psi, x => k)
    abs(psi_k)^2
end

psi0 = 1 / cosh(10 - x)

results = integrate(
    eq, psi0, 0,
    samplers=Dict(
        :xdensity => (xdensity, linspace(0, 5, 101)),
        :kdensity => (kdensity, linspace(0, 5, 101))
        ))

plot_results(results)


text:
# We have a similar QP set-up as before

# Simply include this file

# QP Set-Up
# By definition, the QP has the following Properties
# minimize_x (1/2) xT Q x + cT x
# subject to Ax ≤ B
#
# where (all are Real)
# Q is an nxn Symmetric Matrix (that is positive definite for convex)
# x is an nx1 vector
# c is an nx1 vector
# A is an mxn Matrix
# b is an mx1 vector

using LinearAlgebra, SparseArrays

include("../objective/QP_Linear_objectives.jl")
include("../constraints/constraintManager.jl")
include("../other_utils/utils.jl")

# --------------------------
# Lagrangian
# φ(y) = f(x) + (ρ/2) ||c(y)||_2^2 + λ c(y)
#      = f(x) + (ρ/2) c(y)'c(y)    + λ c(y)
# --------------------------
# Where y = [x; s; t]
# --------------------------
# lambdaInit = zeros(size(fObj(x0)))
# rhoInit = 1

@doc raw"""
    augLag(obj::objectiveFunc, cM::constraintManager, rho::Float64)

This is a mutable struct that stores an SOCP Augmented Lagrangian.

For a single SOCP constraint, the lagrangian is
```math
φ(y) = f(y) + (ρ/2) ||c(y)||_2^2 + λ c(y)
```
```math
φ(y) = f(y) + (ρ/2) c(y)'c(y)    + λ c(y)
```

Where `y` is represented a primal vector struct.

Generally, we initialize `ρ` to 1.

See also: [`evalAL`](@ref), [`evalGradAL`](@ref), [`evalHessAl`](@ref)
"""
mutable struct augLag
    obj::objectiveFunc
    cM::constraintManager
    rho::Float64
end

function getPrimals(al::augLag, primalDualVec)
    if typeof(al.cM) == constraintManager_Dynamics
        # println("Adjusting Primal Vec")
        affLambdaSize = size(al.cM.affineLambdaList, 1)
        primalSize = size(primalDualVec, 1) - affLambdaSize
        primals = primalDualVec[1:primalSize]
    else
        # There are no duals
        primals = primalDualVec
    end

    return primals
end

"""
    evalAL(alQP::augLag, y)

Evaluates the Augmented Lagrangian with the primals `y`. If linear dynamics are
present, it will special case the dynamics.

returns a real number
"""
function evalAL(alQP::augLag, y)
    # φ(y) = f(x) + (ρ/2) c(y)'c(y)    + λ c(y)
    # φ(y) = [(1/2) xT Q x + cT x] + (ρ/2) (c(y))'(c(y)) + λ (c(y))

    primals = getPrimals(alQP, y)

    fCurr = fObjQP(alQP.obj, primals)[1]
    cCurr = evalConstraints(alQP.cM, primals, alQP.rho)

    # println("f(x) = $fCurr")
    # println(" with size -> $(size(fCurr))")
    # println("c(x) = $cCurr")
    # println(" with size -> $(size(cCurr))")

    if typeof(alQP.cM) == constraintManager_Dynamics
        # φ(y) + λ d(y)
        duals = alQP.cM.affineLambdaList
        cAff = evalAffineEq(alQP.cM, primals)
        # println("Duals (λ) = $(size(duals))")
        # println("Ax - b = $(size(cAff))")
        cAffVal = duals'cAff
        # println("λ'(Ax-b) = $(size(cAffVal)) = $cAffVal")
        # println("-> $(cAffVal[1])")
        cCurr += cAffVal[1]
    end

    return fCurr + cCurr
end

"""
    evalGradAL(alQP::augLag, y)

Evaluates the gradient of the augmented lagrangian with the primals `y`. The
format is adjusted depending on the type of constraint manager.

For a base constraint manager (all constraints in the augmented lagrangian),
the gradient is just ∇φ = ∇f + Σ ∇ci.

For a constraint manager that separates the equality constraints, the gradient
also includes the derivatives with respect to the duals of the equality
constraints. So ∇φ = [∇f + Σ ∇ci; Ax - b]

returns a vector of size `y`
"""
function evalGradAL(alQP::augLag, y)
    # ∇φ(y) = ∇f(x) + J(c(y))'(ρ c(y) + λ)
    # y = [x; s; t]
    # ∇φ(y) is (n+m+1)x1

    primals = getPrimals(alQP, y)

    gradfCurr = dfdxQP(alQP.obj, primals)
    gradCCurr = evalGradConstraints(alQP.cM, primals, alQP.rho)

    gradPhiPrimal = gradfCurr + gradCCurr

    if typeof(alQP.cM) == constraintManager_Base
        return gradPhiPrimal
    elseif typeof(alQP.cM) == constraintManager_Dynamics
        duals = alQP.cM.affineLambdaList
        Aaff = affineBlock(alQP.cM)
        return [gradPhiPrimal + Aaff'duals; evalAffineEq(alQP.cM, primals)]
    end
end

"""
    evalHessAl(alQP::augLag, y)

Evaluates the Hessian of the Augmented Lagrangian of an SOCP with the
primals `y`. The format is adjusted depending on the type of constraint manager.

For a base constraint manager (all constraints in the augmented lagrangian),
the hessian is just Hφ = Hf + Σ Hci.

For a constraint manager that separates the equality constraints, the gradient
also includes the derivatives with respect to the duals of the equality
constraints. So Hφ = [(Hf + Σ Hci)     A'; A    0]

returns a Symmetric Matrix of size `y`×`y` or `y + b`×`y + b`
"""
function evalHessAl(alQP::augLag, y)

    primals = getPrimals(alQP, y)

    hessf = hessQP(alQP.obj)
    hessC = evalHessConstraints(alQP.cM, primals, alQP.rho)

    hessPhiPrimals = hessf + hessC

    if typeof(alQP.cM) == constraintManager_Base
        return hessPhiPrimals
    elseif typeof(alQP.cM) == constraintManager_Dynamics
        Aaff = affineBlock(alQP.cM)
        blockSize = size(Aaff, 1)
        return [hessPhiPrimals Aaff'; Aaff spzeros(blockSize, blockSize)]
    end
end

"""
    calcNormGradResiduals(alQP::augLag, yList)

Calculates the 2-norm of the residuals of an augmented lagrangian at each of
the inputs `y`. The 2-norm in safe in that it is bounded below (e.g. by
10^(-20)) and can be plotted on a semilog or log-log plot.

See also: [`safeNorm`](@ref)
"""
function calcNormGradResiduals(alQP::augLag, yList)
    #=
    Calculate the residuals where the AL is the merit function.
    Returns the norm of the gradient of the AL at each point in xArr
    =#

    primals = getPrimals(alQP, yList[1])
    pSize = size(primals, 1)
    resArr = [evalGradAL(alQP, y[1:pSize]) for y in yList]

    return safeNorm(resArr)
end

"""
    calcALArr(alQP::augLag, yList::Array{SOCP_primals, 1})

Evaluates the augmented lagrangian at each input `y`.
"""
function calcALArr(alQP::augLag, yList)
    #=
    Calculate the value of the Augmented Lagrangian at each point
    =#
    return [evalAL(alQP, y) for y in yList]
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
import PWS
using Test

system = PWS.gene_expression_system(dtimes=0:0.2:1.0)
conf = PWS.generate_configuration(system)
ensemble = PWS.MarginalEnsemble(system)

# Test the shooting moves
jump_problem = ensemble.jump_problem
original_traj = conf.s_traj
new_traj = copy(conf.s_traj)
@test original_traj == new_traj
PWS.shoot_forward!(new_traj, original_traj, jump_problem, 0.5)
@test original_traj != new_traj
@test issorted(new_traj.t) && allunique(new_traj.t)
@test length(new_traj.t) == length(new_traj.u) == length(new_traj.i)

new_traj = copy(conf.s_traj)
@test original_traj == new_traj
PWS.shoot_backward!(new_traj, original_traj, jump_problem, 0.5)
@test original_traj != new_traj
@test issorted(new_traj.t) && allunique(new_traj.t)
@test length(new_traj.t) == length(new_traj.u) == length(new_traj.i)

mchain = PWS.chain(ensemble; θ=1.0)
@test PWS.energy(conf, mchain) < 0

# Test Thermodynamic Integration

alg = PWS.TIEstimate(100, 16, 1000)
result = PWS.simulate(alg, conf, ensemble)
@test all(PWS.log_marginal(result) .>= 0)


text:
function getLowestCount(chars)
    dict = Dict()
    for char in chars
       dict[char] = get(dict, char, 0) + 1
    end
    return sort(collect(dict), by= tuple -> last(tuple))[1]
end

f = open("../day6.txt")
    lines = readlines(f)
    chars = map(s -> split(strip(s), ""), lines)
    reshaped = hcat(chars...)
    columns = map(x -> getLowestCount(reshaped[x,:]), 1:length(reshaped[:,1]))
    answer = *(map(c -> c[1], columns)...)
    println(answer)
close(f)


text:
using Revise
using Revise354

rlogger = Revise.debug_logger()

entr([], [Revise354]) do
    for log in rlogger.logs
        println(log)
    end
    Revise354.test()
end