text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to ECI to ECI transformations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Get the current EOP Data.
#
# TODO: The EOP data obtained from IERS website does not match the values in the
# examples in [1]. However, this should be enough, because 1) the individual
# functions at the low level are tested using the same values of [1], and 2) the
# difference is smaller than 30 cm.

eop_iau1980  = read_iers_eop("./eop_IAU1980.txt",  :IAU1980)
eop_iau2000a = read_iers_eop("./eop_IAU2000A.txt", :IAU2000A)

# File: ./src/transformations/eci_to_eci.jl
# =========================================

################################################################################
#                                  IAU-76/FK5
################################################################################

# Functions: rECItoECI
# --------------------

# GCRF <=> J2000
# ==============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#   v_gcrf =   -4.7432201570 i +    0.7905364970 j +    5.5337557270 k [km/s]
#
# one gets:
#
#   r_j2000 = 5102.50960000   i + 6123.01152000   j + 6378.13630000   k [km]
#   v_j2000 =   -4.7432196000 i +    0.7905366000 j +    5.5337561900 k [km/s]
#
################################################################################

@testset "Function rECItoECI GCRF <=> J2000" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## GCRF => J2000
    ## =============

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]
    v_gcrf = [-4.7432201570; 0.7905364970; 5.5337557270]

    ## DCM
    ## ---

    D_J2000_GCRF = rECItoECI(GCRF(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_J2000_GCRF*r_gcrf
    v_j2000 = D_J2000_GCRF*v_gcrf

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7

    ## Quaternion
    ## ----------

    q_J2000_GCRF = rECItoECI(Quaternion, GCRF(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = vect(conj(q_J2000_GCRF)*r_gcrf*q_J2000_GCRF)
    v_j2000 = vect(conj(q_J2000_GCRF)*v_gcrf*q_J2000_GCRF)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7

    ## J2000 => GCRF
    ## =============

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_GCRF_J2000 = rECItoECI(J2000(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_GCRF_J2000*r_j2000
    v_gcrf = D_GCRF_J2000*v_j2000

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7

    ## Quaternion
    ## ----------

    q_GCRF_J2000 = rECItoECI(Quaternion, J2000(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(conj(q_GCRF_J2000)*r_j2000*q_GCRF_J2000)
    v_gcrf = vect(conj(q_GCRF_J2000)*v_j2000*q_GCRF_J2000)

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7
end

# GCRF <=> MOD
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#   v_gcrf =   -4.7432201570 i +    0.7905364970 j +    5.5337557270 k [km/s]
#
# one gets:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#   v_mod =   -4.7462630520 i +    0.7860140450 j +    5.5317905620 k [km/s]
#
################################################################################

@testset "Function rECItoECI GCRF <=> MOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## GCRF => MOD
    ## ===========

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]
    v_gcrf = [-4.7432201570; 0.7905364970; 5.5337557270]

    ## DCM
    ## ---

    D_MOD_GCRF = rECItoECI(GCRF(), MOD(), JD_UTC, eop_iau1980)

    r_mod = D_MOD_GCRF*r_gcrf
    v_mod = D_MOD_GCRF*v_gcrf

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    ## Quaternion
    ## ----------

    q_MOD_GCRF = rECItoECI(Quaternion, GCRF(), MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(conj(q_MOD_GCRF)*r_gcrf*q_MOD_GCRF)
    v_mod = vect(conj(q_MOD_GCRF)*v_gcrf*q_MOD_GCRF)

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    ## MOD => GCRF
    ## ===========

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]
    v_mod = [-4.7462630520; 0.7860140450; 5.5317905620]

    ## DCM
    ## ---

    D_GCRF_MOD = rECItoECI(MOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_GCRF_MOD*r_mod
    v_gcrf = D_GCRF_MOD*v_mod

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7

    ## Quaternion
    ## ----------

    q_GCRF_MOD = rECItoECI(Quaternion, MOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(conj(q_GCRF_MOD)*r_mod*q_GCRF_MOD)
    v_gcrf = vect(conj(q_GCRF_MOD)*v_mod*q_GCRF_MOD)

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7
end

# GCRF <=> TOD
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#   v_gcrf =   -4.7432201570 i +    0.7905364970 j +    5.5337557270 k [km/s]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#   v_tod =   -4.7460883850 i +    0.7860783240 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECItoECI GCRF <=> TOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## GCRF => TOD
    ## ===========

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]
    v_gcrf = [-4.7432201570; 0.7905364970; 5.5337557270]

    ## DCM
    ## ---

    D_TOD_GCRF = rECItoECI(GCRF(), TOD(), JD_UTC, eop_iau1980)

    r_tod = D_TOD_GCRF*r_gcrf
    v_tod = D_TOD_GCRF*v_gcrf

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TOD_GCRF = rECItoECI(Quaternion, GCRF(), TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(conj(q_TOD_GCRF)*r_gcrf*q_TOD_GCRF)
    v_tod = vect(conj(q_TOD_GCRF)*v_gcrf*q_TOD_GCRF)

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    ## TOD => GCRF
    ## ===========

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]
    v_tod = [-4.7460883850; 0.7860783240; 5.5319312880]

    ## DCM
    ## ---

    D_GCRF_TOD = rECItoECI(TOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_GCRF_TOD*r_tod
    v_gcrf = D_GCRF_TOD*v_tod

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7

    ## Quaternion
    ## ----------

    q_GCRF_TOD = rECItoECI(Quaternion, TOD(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(conj(q_GCRF_TOD)*r_tod*q_GCRF_TOD)
    v_gcrf = vect(conj(q_GCRF_TOD)*v_tod*q_GCRF_TOD)

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7
end

# GCRF <=> TEME
# =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_gcrf = 5102.50895790   i + 6123.01140070   j + 6378.13692820   k [km]
#   v_gcrf =   -4.7432201570 i +    0.7905364970 j +    5.5337557270 k [km/s]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECItoECI GCRF <=> TEME" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## GCRF => TEME
    ## ===========

    r_gcrf = [5102.50895790; 6123.01140070; 6378.13692820]
    v_gcrf = [-4.7432201570; 0.7905364970; 5.5337557270]

    ## DCM
    ## ---

    D_TEME_GCRF = rECItoECI(GCRF(), TEME(), JD_UTC, eop_iau1980)

    r_teme = D_TEME_GCRF*r_gcrf
    v_teme = D_TEME_GCRF*v_gcrf

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TEME_GCRF = rECItoECI(Quaternion, GCRF(), TEME(), JD_UTC, eop_iau1980)

    r_teme = vect(conj(q_TEME_GCRF)*r_gcrf*q_TEME_GCRF)
    v_teme = vect(conj(q_TEME_GCRF)*v_gcrf*q_TEME_GCRF)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## TEME => GCRF
    ## ============

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_GCRF_TEME = rECItoECI(TEME(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = D_GCRF_TEME*r_teme
    v_gcrf = D_GCRF_TEME*v_teme

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7

    ## Quaternion
    ## ----------

    q_GCRF_TEME = rECItoECI(Quaternion, TEME(), GCRF(), JD_UTC, eop_iau1980)

    r_gcrf = vect(conj(q_GCRF_TEME)*r_teme*q_GCRF_TEME)
    v_gcrf = vect(conj(q_GCRF_TEME)*v_teme*q_GCRF_TEME)

    @test r_gcrf[1] ≈ +5102.50895790 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01140070 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13692820 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201570  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364970  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557270  atol=1e-7
end

# J2000 <=> MOD
# =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC     = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#   v_j2000 =  -4.7432196000 i + 0.7905366000    j + 5.5337561900    k [km/s]
#
# one gets:
#
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#   v_mod =   -4.7462630520 i +    0.7860140450 j +    5.5317905620 k [km/s]
#
# If not EOP correction is used, then we get:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#   v_mod =   -4.7462624950 i +    0.7860141490 j +    5.5317910250 k [km/s]
#
################################################################################

@testset "Function rECItoECI J2000 <=> MOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## J2000 => MOD
    ## ============

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_MOD_J2000 = rECItoECI(J2000(), MOD(), JD_UTC, eop_iau1980)

    r_mod = D_MOD_J2000*r_j2000
    v_mod = D_MOD_J2000*v_j2000

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    ## Quaternion
    ## ----------

    q_MOD_J2000 = rECItoECI(Quaternion, J2000(), MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(conj(q_MOD_J2000)*r_j2000*q_MOD_J2000)
    v_mod = vect(conj(q_MOD_J2000)*v_j2000*q_MOD_J2000)

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    ## MOD => J2000
    ## ============

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]
    v_mod = [-4.7462630520; 0.7860140450; 5.5317905620]

    ## DCM
    ## ---

    D_J2000_MOD = rECItoECI(MOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_J2000_MOD*r_mod
    v_j2000 = D_J2000_MOD*v_mod

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7

    ## Quaternion
    ## ----------

    q_J2000_MOD = rECItoECI(Quaternion, MOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = vect(conj(q_J2000_MOD)*r_mod*q_J2000_MOD)
    v_j2000 = vect(conj(q_J2000_MOD)*v_mod*q_J2000_MOD)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7

    # No EOP Corrections
    # ==========================================================================

    ## J2000 => MOD
    ## ============

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_MOD_J2000 = rECItoECI(J2000(), MOD(), JD_UTC)

    r_mod = D_MOD_J2000*r_j2000
    v_mod = D_MOD_J2000*v_j2000

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## Quaternion
    ## ----------

    q_MOD_J2000 = rECItoECI(Quaternion, J2000(), MOD(), JD_UTC)

    r_mod = vect(conj(q_MOD_J2000)*r_j2000*q_MOD_J2000)
    v_mod = vect(conj(q_MOD_J2000)*v_j2000*q_MOD_J2000)

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## MOD => J2000
    ## ============

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]
    v_mod = [-4.7462624950; 0.7860141490; 5.5317910250]

    ## DCM
    ## ---

    D_J2000_MOD = rECItoECI(MOD(), J2000(), JD_UTC)

    r_j2000 = D_J2000_MOD*r_mod
    v_j2000 = D_J2000_MOD*v_mod

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9

    ## Quaternion
    ## ----------

    q_J2000_MOD = rECItoECI(Quaternion, MOD(), J2000(), JD_UTC)

    r_j2000 = vect(conj(q_J2000_MOD)*r_mod*q_J2000_MOD)
    v_j2000 = vect(conj(q_J2000_MOD)*v_mod*q_J2000_MOD)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9
end

# J2000 <=> TOD
# =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#   v_j2000 =  -4.7432196000 i + 0.7905366000    j + 5.5337561900    k [km/s]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#   v_tod =   -4.7460883850 i +    0.7860783240 j +    5.5319312880 k [km/s]
#
# If not EOP correction is used, then we get:
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#   v_tod =   -4.7460885670 i +    0.7860772220 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECItoECI J2000 <=> TOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## J2000 => TOD
    ## ============

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_TOD_J2000 = rECItoECI(J2000(), TOD(), JD_UTC, eop_iau1980)

    r_tod = D_TOD_J2000*r_j2000
    v_tod = D_TOD_J2000*v_j2000

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TOD_J2000 = rECItoECI(Quaternion, J2000(), TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(conj(q_TOD_J2000)*r_j2000*q_TOD_J2000)
    v_tod = vect(conj(q_TOD_J2000)*v_j2000*q_TOD_J2000)

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    ## TOD => J2000
    ## ============

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]
    v_tod = [-4.7460883850; 0.7860783240; 5.5319312880]

    ## DCM
    ## ---

    D_J2000_TOD = rECItoECI(TOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_J2000_TOD*r_tod
    v_j2000 = D_J2000_TOD*v_tod

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7

    ## Quaternion
    ## ----------

    q_J2000_TOD = rECItoECI(Quaternion, TOD(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = vect(conj(q_J2000_TOD)*r_tod*q_J2000_TOD)
    v_j2000 = vect(conj(q_J2000_TOD)*v_tod*q_J2000_TOD)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7

    # No EOP corrections
    # ==========================================================================

    ## J2000 => TOD
    ## ============

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_TOD_J2000 = rECItoECI(J2000(), TOD(), JD_UTC)

    r_tod = D_TOD_J2000*r_j2000
    v_tod = D_TOD_J2000*v_j2000

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TOD_J2000 = rECItoECI(Quaternion, J2000(), TOD(), JD_UTC)

    r_tod = vect(conj(q_TOD_J2000)*r_j2000*q_TOD_J2000)
    v_tod = vect(conj(q_TOD_J2000)*v_j2000*q_TOD_J2000)

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## TOD => J2000
    ## ============

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]
    v_tod = [-4.7460885670; 0.7860772220; 5.5319312880]

    ## DCM
    ## ---

    D_J2000_TOD = rECItoECI(TOD(), J2000(), JD_UTC)

    r_j2000 = D_J2000_TOD*r_tod
    v_j2000 = D_J2000_TOD*v_tod

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9

    ## Quaternion
    ## ----------

    q_J2000_TOD = rECItoECI(Quaternion, TOD(), J2000(), JD_UTC)

    r_j2000 = vect(conj(q_J2000_TOD)*r_tod*q_J2000_TOD)
    v_j2000 = vect(conj(q_J2000_TOD)*v_tod*q_J2000_TOD)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-7
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-7
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-7

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-9
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-9
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-9
end

# J2000 <=> TEME
# ==============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_j2000 = 5102.50960000  i + 6123.01152000   j + 6378.13630000   k [km]
#   v_j2000 =  -4.7432196000 i + 0.7905366000    j + 5.5337561900    k [km/s]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECItoECI J2000 <=> TEME" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## J2000 => TEME
    ## =============

    r_j2000 = [5102.50960000; 6123.01152000; 6378.13630000]
    v_j2000 = [-4.7432196000; 0.7905366000; 5.5337561900]

    ## DCM
    ## ---

    D_TEME_J2000 = rECItoECI(J2000(), TEME(), JD_UTC, eop_iau1980)

    r_teme = D_TEME_J2000*r_j2000
    v_teme = D_TEME_J2000*v_j2000

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TEME_J2000 = rECItoECI(Quaternion, J2000(), TEME(), JD_UTC)

    r_teme = vect(conj(q_TEME_J2000)*r_j2000*q_TEME_J2000)
    v_teme = vect(conj(q_TEME_J2000)*v_j2000*q_TEME_J2000)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## TEME => J2000
    ## =============

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_J2000_TEME = rECItoECI(TEME(), J2000(), JD_UTC, eop_iau1980)

    r_j2000 = D_J2000_TEME*r_teme
    v_j2000 = D_J2000_TEME*v_teme

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7

    ## Quaternion
    ## ----------

    q_J2000_TEME = rECItoECI(Quaternion, TEME(), J2000(), JD_UTC)

    r_j2000 = vect(conj(q_J2000_TEME)*r_teme*q_J2000_TEME)
    v_j2000 = vect(conj(q_J2000_TEME)*v_teme*q_J2000_TEME)

    @test r_j2000[1] ≈ +5102.50960000 atol=1e-4
    @test r_j2000[2] ≈ +6123.01152000 atol=1e-4
    @test r_j2000[3] ≈ +6378.13630000 atol=1e-4

    @test v_j2000[1] ≈ -4.7432196000  atol=1e-7
    @test v_j2000[2] ≈ +0.7905366000  atol=1e-7
    @test v_j2000[3] ≈ +5.5337561900  atol=1e-7
end

# MOD <=> TOD
# ===========

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#   v_mod =   -4.7462630520 i +    0.7860140450 j +    5.5317905620 k [km/s]
#
# one gets:
#
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#   v_tod =   -4.7460883850 i +    0.7860783240 j +    5.5319312880 k [km/s]
#
# If not EOP correction is used, then we get:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#   v_mod =   -4.7462624950 i +    0.7860141490 j +    5.5317910250 k [km/s]
#
#   r_tod = 5094.51478040   i + 6127.36646120   j + 6380.34453270   k [km]
#   v_tod =   -4.7460885670 i +    0.7860772220 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECItoECI MOD <=> TOD" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## MOD => TOD
    ## ==========

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]
    v_mod = [-4.7462630520; 0.7860140450; 5.5317905620]

    ## DCM
    ## ---

    D_TOD_MOD = rECItoECI(MOD(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = D_TOD_MOD*r_mod
    v_tod = D_TOD_MOD*v_mod

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TOD_MOD = rECItoECI(Quaternion, MOD(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(conj(q_TOD_MOD)*r_mod*q_TOD_MOD)
    v_tod = vect(conj(q_TOD_MOD)*v_mod*q_TOD_MOD)

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    ## TOD => MOD
    ## ==========

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]
    v_tod = [-4.7460883850; 0.7860783240; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TOD = rECItoECI(TOD(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = D_MOD_TOD*r_tod
    v_mod = D_MOD_TOD*v_tod

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    ## Quaternion
    ## ----------

    q_MOD_TOD = rECItoECI(Quaternion, TOD(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(conj(q_MOD_TOD)*r_tod*q_MOD_TOD)
    v_mod = vect(conj(q_MOD_TOD)*v_tod*q_MOD_TOD)

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    # No EOP corrections
    # ==========================================================================

    ## MOD => TOD
    ## ==========

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]
    v_mod = [-4.7462624950; 0.7860141490; 5.5317910250]

    ## DCM
    ## ---

    D_TOD_MOD = rECItoECI(MOD(), JD_UTC, TOD(), JD_UTC)

    r_tod = D_TOD_MOD*r_mod
    v_tod = D_TOD_MOD*v_mod

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TOD_MOD = rECItoECI(Quaternion, MOD(), JD_UTC, TOD(), JD_UTC)

    r_tod = vect(conj(q_TOD_MOD)*r_mod*q_TOD_MOD)
    v_tod = vect(conj(q_TOD_MOD)*v_mod*q_TOD_MOD)

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## TOD => MOD
    ## ==========

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]
    v_tod = [-4.7460885670; 0.7860772220; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TOD = rECItoECI(TOD(), JD_UTC, MOD(), JD_UTC)

    r_mod = D_MOD_TOD*r_tod
    v_mod = D_MOD_TOD*v_tod

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## Quaternion
    ## ----------

    q_MOD_TOD = rECItoECI(Quaternion, TOD(), JD_UTC, MOD(), JD_UTC)

    r_mod = vect(conj(q_MOD_TOD)*r_tod*q_MOD_TOD)
    v_mod = vect(conj(q_MOD_TOD)*v_tod*q_MOD_TOD)

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9
end

# MOD <=> TEME
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_mod = 5094.02837450   i + 6127.87081640   j + 6380.24851640   k [km]
#   v_mod =   -4.7462630520 i +    0.7860140450 j +    5.5317905620 k [km/s]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
# If not EOP correction is used, then we get the same result by using:
#
#   r_mod = 5094.02901670   i + 6127.87093630   j + 6380.24788850   k [km]
#   v_mod =   -4.7462624950 i +    0.7860141490 j +    5.5317910250 k [km/s]
#
################################################################################

@testset "Function rECItoECI MOD <=> TEME" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## MOD => TEME
    ## ===========

    r_mod = [5094.02837450; 6127.87081640; 6380.24851640]
    v_mod = [-4.7462630520; 0.7860140450; 5.5317905620]

    ## DCM
    ## ---

    D_TEME_MOD = rECItoECI(MOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = D_TEME_MOD*r_mod
    v_teme = D_TEME_MOD*v_mod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TEME_MOD = rECItoECI(Quaternion, MOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = vect(conj(q_TEME_MOD)*r_mod*q_TEME_MOD)
    v_teme = vect(conj(q_TEME_MOD)*v_mod*q_TEME_MOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## TEME => MOD
    ## ===========

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TEME = rECItoECI(TEME(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = D_MOD_TEME*r_teme
    v_mod = D_MOD_TEME*v_teme

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    ## Quaternion
    ## ----------

    q_MOD_TEME = rECItoECI(Quaternion, TEME(), JD_UTC, MOD(), JD_UTC, eop_iau1980)

    r_mod = vect(conj(q_MOD_TEME)*r_teme*q_MOD_TEME)
    v_mod = vect(conj(q_MOD_TEME)*v_teme*q_MOD_TEME)

    @test r_mod[1] ≈ +5094.02837450 atol=1e-4
    @test r_mod[2] ≈ +6127.87081640 atol=1e-4
    @test r_mod[3] ≈ +6380.24851640 atol=1e-4

    @test v_mod[1] ≈ -4.7462630520  atol=1e-7
    @test v_mod[2] ≈ +0.7860140450  atol=1e-7
    @test v_mod[3] ≈ +5.5317905620  atol=1e-7

    # No EOP corrections
    # ==========================================================================

    ## MOD => TEME
    ## ===========

    r_mod = [5094.02901670; 6127.87093630; 6380.24788850]
    v_mod = [-4.7462624950; 0.7860141490; 5.5317910250]

    ## DCM
    ## ---

    D_TEME_MOD = rECItoECI(MOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = D_TEME_MOD*r_mod
    v_teme = D_TEME_MOD*v_mod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TEME_MOD = rECItoECI(Quaternion, MOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = vect(conj(q_TEME_MOD)*r_mod*q_TEME_MOD)
    v_teme = vect(conj(q_TEME_MOD)*v_mod*q_TEME_MOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-7
    @test r_teme[2] ≈ +6127.64465950 atol=1e-7
    @test r_teme[3] ≈ +6380.34453270 atol=1e-7

    @test v_teme[1] ≈ -4.7461314870  atol=1e-9
    @test v_teme[2] ≈ +0.7858180410  atol=1e-9
    @test v_teme[3] ≈ +5.5319312880  atol=1e-9

    ## TEME => MOD
    ## ===========

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_MOD_TEME = rECItoECI(TEME(), JD_UTC, MOD(), JD_UTC)

    r_mod = D_MOD_TEME*r_teme
    v_mod = D_MOD_TEME*v_teme

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9

    ## Quaternion
    ## ----------

    q_MOD_TEME = rECItoECI(Quaternion, TEME(), JD_UTC, MOD(), JD_UTC)

    r_mod = vect(conj(q_MOD_TEME)*r_teme*q_MOD_TEME)
    v_mod = vect(conj(q_MOD_TEME)*v_teme*q_MOD_TEME)

    @test r_mod[1] ≈ +5094.02901670 atol=1e-7
    @test r_mod[2] ≈ +6127.87093630 atol=1e-7
    @test r_mod[3] ≈ +6380.24788850 atol=1e-7

    @test v_mod[1] ≈ -4.7462624950  atol=1e-9
    @test v_mod[2] ≈ +0.7860141490  atol=1e-9
    @test v_mod[3] ≈ +5.5317910250  atol=1e-9
end

# TOD <=> TEME
# ============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_tod = 5094.51620300   i + 6127.36527840   j + 6380.34453270   k [km]
#   v_tod =   -4.7460883850 i +    0.7860783240 j +    5.5319312880 k [km/s]
#
# one gets:
#
#   r_teme = 5094.18016210   i + 6127.64465950   j + 6380.34453270   k [km]
#   v_teme =   -4.7461314870 i +    0.7858180410 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function rECItoECI TOD <=> TEME" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## TOD => TEME
    ## ===========

    r_tod = [5094.51620300; 6127.36527840; 6380.34453270]
    v_tod = [-4.7460883850; 0.7860783240; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_TOD = rECItoECI(TOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = D_TEME_TOD*r_tod
    v_teme = D_TEME_TOD*v_tod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TEME_TOD = rECItoECI(Quaternion, TOD(), JD_UTC, TEME(), JD_UTC, eop_iau1980)

    r_teme = vect(conj(q_TEME_TOD)*r_tod*q_TEME_TOD)
    v_teme = vect(conj(q_TEME_TOD)*v_tod*q_TEME_TOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## TEME => TOD
    ## ===========

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TOD_TEME = rECItoECI(TEME(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = D_TOD_TEME*r_teme
    v_tod = D_TOD_TEME*v_teme

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TOD_TEME = rECItoECI(Quaternion, TEME(), JD_UTC, TOD(), JD_UTC, eop_iau1980)

    r_tod = vect(conj(q_TOD_TEME)*r_teme*q_TOD_TEME)
    v_tod = vect(conj(q_TOD_TEME)*v_teme*q_TOD_TEME)

    @test r_tod[1] ≈ +5094.51620300 atol=1e-4
    @test r_tod[2] ≈ +6127.36527840 atol=1e-4
    @test r_tod[3] ≈ +6380.34453270 atol=1e-4

    @test v_tod[1] ≈ -4.7460883850  atol=1e-7
    @test v_tod[2] ≈ +0.7860783240  atol=1e-7
    @test v_tod[3] ≈ +5.5319312880  atol=1e-7

    # No EOP corrections
    # ==========================================================================

    ## TOD => TEME
    ## ===========

    r_tod = [5094.51478040; 6127.36646120; 6380.34453270]
    v_tod = [-4.7460885670; 0.7860772220; 5.5319312880]

    ## DCM
    ## ---

    D_TEME_TOD = rECItoECI(TOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = D_TEME_TOD*r_tod
    v_teme = D_TEME_TOD*v_tod

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_TEME_TOD = rECItoECI(Quaternion, TOD(), JD_UTC, TEME(), JD_UTC)

    r_teme = vect(conj(q_TEME_TOD)*r_tod*q_TEME_TOD)
    v_teme = vect(conj(q_TEME_TOD)*v_tod*q_TEME_TOD)

    @test r_teme[1] ≈ +5094.18016210 atol=1e-4
    @test r_teme[2] ≈ +6127.64465950 atol=1e-4
    @test r_teme[3] ≈ +6380.34453270 atol=1e-4

    @test v_teme[1] ≈ -4.7461314870  atol=1e-7
    @test v_teme[2] ≈ +0.7858180410  atol=1e-7
    @test v_teme[3] ≈ +5.5319312880  atol=1e-7

    ## TEME => TOD
    ## ===========

    r_teme = [5094.18016210; 6127.64465950; 6380.34453270]
    v_teme = [-4.7461314870; 0.7858180410; 5.5319312880]

    ## DCM
    ## ---

    D_TOD_TEME = rECItoECI(TEME(), JD_UTC, TOD(), JD_UTC)

    r_tod = D_TOD_TEME*r_teme
    v_tod = D_TOD_TEME*v_teme

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9

    ## Quaternion
    ## ----------

    q_TOD_TEME = rECItoECI(Quaternion, TEME(), JD_UTC, TOD(), JD_UTC)

    r_tod = vect(conj(q_TOD_TEME)*r_teme*q_TOD_TEME)
    v_tod = vect(conj(q_TOD_TEME)*v_teme*q_TOD_TEME)

    @test r_tod[1] ≈ +5094.51478040 atol=1e-7
    @test r_tod[2] ≈ +6127.36646120 atol=1e-7
    @test r_tod[3] ≈ +6380.34453270 atol=1e-7

    @test v_tod[1] ≈ -4.7460885670  atol=1e-9
    @test v_tod[2] ≈ +0.7860772220  atol=1e-9
    @test v_tod[3] ≈ +5.5319312880  atol=1e-9
end

################################################################################
#                                IAU-2006/2010
################################################################################

## GCRF <=> CIRS
## =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
# NOTE: It seems that the results in the Example 3-14 is computed **without**
# the `dX` and `dY` corrections, whereas in the Table 3-6 they are computed
# **with** the corrections.
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_cirs = -5100.01840470   i + 6122.78636480   j + 6380.34453270   k [km]
#   v_cirs =    -4.7453803300 i -    0.7903414530 j +    5.5319312880 k [km/s]
#
# one gets the following (this is the result in Table 3-6):
#
#   r_gcrf = 5102.50895290  i + 6123.01139910 j + 6378.13693380 k [km]
#   v_gcrf =  -4.7432201610 i + 0.7905364950  j + 5.5337557240  k [km/s]
#
################################################################################

@testset "Function rECItoECI GCRF <=> CIRS" begin
    JD_UTC = DatetoJD(2004, 4, 6, 7, 51, 28.386009)

    ## GCRF => CIRS
    ## ============

    r_gcrf = [5102.50895290; 6123.01139910; 6378.13693380]
    v_gcrf = [-4.7432201610; 0.7905364950; 5.5337557240]

    ## DCM
    ## ---

    D_CIRS_GCRF = rECItoECI(GCRF(), CIRS(), JD_UTC, eop_iau2000a)

    r_cirs = D_CIRS_GCRF*r_gcrf
    v_cirs = D_CIRS_GCRF*v_gcrf

    @test r_cirs[1] ≈ +5100.01840470 atol=1e-4
    @test r_cirs[2] ≈ +6122.78636480 atol=1e-4
    @test r_cirs[3] ≈ +6380.34453270 atol=1e-4

    @test v_cirs[1] ≈ -4.7453803300  atol=1e-7
    @test v_cirs[2] ≈ +0.7903414530  atol=1e-7
    @test v_cirs[3] ≈ +5.5319312880  atol=1e-7

    ## Quaternion
    ## ----------

    q_CIRS_GCRF = rECItoECI(Quaternion, GCRF(), CIRS(), JD_UTC, eop_iau2000a)

    r_cirs = vect(conj(q_CIRS_GCRF)*r_gcrf*q_CIRS_GCRF)
    v_cirs = vect(conj(q_CIRS_GCRF)*v_gcrf*q_CIRS_GCRF)

    @test r_cirs[1] ≈ +5100.01840470 atol=1e-4
    @test r_cirs[2] ≈ +6122.78636480 atol=1e-4
    @test r_cirs[3] ≈ +6380.34453270 atol=1e-4

    @test v_cirs[1] ≈ -4.7453803300  atol=1e-7
    @test v_cirs[2] ≈ +0.7903414530  atol=1e-7
    @test v_cirs[3] ≈ +5.5319312880  atol=1e-7

    ## CIRS => GCRF
    ## ============

    r_cirs = [+5100.01840470; +6122.78636480; +6380.34453270]
    v_cirs = [-4.7453803300; +0.7903414530; +5.5319312880]

    ## DCM
    ## ---

    D_GCRF_CIRS = rECItoECI(CIRS(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = D_GCRF_CIRS*r_cirs
    v_gcrf = D_GCRF_CIRS*v_cirs

    @test r_gcrf[1] ≈ +5102.50895290 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201610  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364950  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557240  atol=1e-7

    ## Quaternion
    ## ----------

    q_GCRF_CIRS = rECItoECI(Quaternion, CIRS(), GCRF(), JD_UTC, eop_iau2000a)

    r_gcrf = vect(conj(q_GCRF_CIRS)*r_cirs*q_GCRF_CIRS)
    v_gcrf = vect(conj(q_GCRF_CIRS)*v_cirs*q_GCRF_CIRS)

    @test r_gcrf[1] ≈ +5102.50895290 atol=1e-4
    @test r_gcrf[2] ≈ +6123.01139910 atol=1e-4
    @test r_gcrf[3] ≈ +6378.13693380 atol=1e-4

    @test v_gcrf[1] ≈ -4.7432201610  atol=1e-7
    @test v_gcrf[2] ≈ +0.7905364950  atol=1e-7
    @test v_gcrf[3] ≈ +5.5337557240  atol=1e-7
end


text:
using Zygote
using Random
using DelimitedFiles
using LinearAlgebra
using GalacticOptim,Optim,NLopt#,BlackBoxOptim
using DifferentialEquations#,DiffEqFlux
#using ForwardDiff
using Plots
using LaTeXStrings
using GLOQ
#using ReverseDiff
#include("../src/GLOQ.jl")
pyplot()
# BlackBoxOptim somehow downgrade some packages and as a result breaks the auto-differentiation with Zygote
# we should avoid it.
# NLopt seems to be really fast. Their LBFGS with box constrains is faster than Optim's FMinbox(LBFGS)

N_states = 4;
freqs = [4.10817777; 3.88303940; 3.61937188]
omegas = 2.0*pi.*freqs
gamma1   = [2.59451982e-05; 4.54296537e-05; 0.0]
gamma2   = [2.58005604e-05; 9.00000000e-05; 0.0]

# driving frequency
omr_t1_01 = 2.0*pi*4.1081786
omr_t1_12 = 2.0*pi*3.8830355
omr_ramsey_01 = 2.0*pi*(4.1081786-5e-4)
omr_ramsey_12 = 2.0*pi*(3.8830355-2.0e-3)
omr_echo_01 = 2.0*pi*4.1081786
omr_echo_12 = 2.0*pi*3.8830355

# charge noise
char_noise = 2.0*pi.*[-1.017e-5,1.058e-4,0]

# width of the control signal
width = 21.25
TC_01 = 2.5*width

width = 17.00
TC_12 = 2.5*width

# in nano second
T_T1_01 = 200.0*1000.0
T_T1_12 = 100.0*1000.0
T_Ramsey_01 = 40.0*1000.0
T_Ramsey_12 = 20.0*1000.0
T_Echo_01 = 40.0*1000.0
T_Echo_12 = 20.0*1000.0

# data size used for characterization
n1use = 201
n2use = 201
n3use = 201
n4use = 281
n5use = 201
n6use = 201

# time points evaluated
dt_t1_01 = T_T1_01/200
dt_t1_12 = T_T1_12/200
dt_ramsey_01 = T_Ramsey_01/200
dt_ramsey_12 = T_Ramsey_12/400
dt_echo_01 = T_Echo_01/200
dt_echo_12 = T_Echo_12/200

t_t1_01 = zeros(Float64,n1use)
for i = 1:n1use
	t_t1_01[i] = (i-1)*dt_t1_01
end
T_evaluated_t1_01 = t_t1_01[end]

t_t1_12 = zeros(Float64,n2use)
for i = 1:n2use
	t_t1_12[i] = (i-1)*dt_t1_12
end
T_evaluated_t1_12 = t_t1_12[end]

t_ramsey_01 = zeros(Float64,n3use)
for i = 1:n3use
	t_ramsey_01[i] = (i-1)*dt_ramsey_01
end
T_evaluated_ramsey_01 = t_ramsey_01[end]

t_ramsey_12 = zeros(Float64,n4use)
for i = 1:n4use
	t_ramsey_12[i] = (i-1)*dt_ramsey_12
end
T_evaluated_ramsey_12 = t_ramsey_12[end]

t_echo_01 = zeros(Float64,n5use)
for i = 1:n5use
	t_echo_01[i] = (i-1)*dt_echo_01
end
T_evaluated_echo_01 = t_echo_01[end]

t_echo_12 = zeros(Float64,n6use)
for i = 1:n6use
	t_echo_12[i] = (i-1)*dt_echo_12
end
T_evaluated_echo_12 = t_echo_12[end]

data_t1_01 = readdlm("data/CD8_Alice_0-1_T1.txt",)
data_t1_01 = transpose(data_t1_01[:,1:n1use])
data_t1_12 = readdlm("data/CD8_Alice_1-2_T1.txt")
data_t1_12 = transpose(data_t1_12[:,1:n2use])
data_ramsey_01 = readdlm("data/Alice_ramsey_0-1_20200722.txt")
data_ramsey_01 = transpose(data_ramsey_01[:,1:n3use])
data_ramsey_12 = readdlm("data/Alice_ramsey_1-2_20200723.txt")
data_ramsey_12 = transpose(data_ramsey_12[:,1:n4use])
data_echo_01 = readdlm("data/CD8_Alice_0-1_echo.txt")
data_echo_01 = transpose(data_echo_01[:,1:n5use])
data_echo_12 = readdlm("data/CD8_Alice_1-2_echo.txt")
data_echo_12 = transpose(data_echo_12[:,1:n6use])

#println(size(data_t1_01))
#println(size(data_t1_12))
#println(size(data_ramsey_01))
#println(size(data_ramsey_12))
#println(size(data_echo_01))
#println(size(data_echo_12))

# initial conditions
u0_t1_01 = [1.0;0.0;0.0;0.0]
u0_t1_12 = [0.0;1.0;0.0;0.0]
u0_ramsey_01 = [1.0;0.0;0.0;0.0]
u0_ramsey_12 = [0.0;1.0;0.0;0.0]
u0_echo_01 = [1.0;0.0;0.0;0.0]
u0_echo_12 = [0.0;1.0;0.0;0.0]

v0_t1_01 = [0.0;0.0;0.0;0.0]
v0_t1_12 = [0.0;0.0;0.0;0.0]
v0_ramsey_01 = [0.0;0.0;0.0;0.0]
v0_ramsey_12 = [0.0;0.0;0.0;0.0]
v0_echo_01 = [0.0;0.0;0.0;0.0]
v0_echo_12 = [0.0;0.0;0.0;0.0]

########################################################
# Define loss function
########################################################
function loss(p)#,p_keywords)
	_rho_Ramsey_u,_rho_Ramsey_v = GLOQ.RamseyParityForwardSolve(u0_ramsey_12,v0_ramsey_12,
				     [p[1:2];omegas[3]],omr_ramsey_12,
					 [p[7:8];char_noise[3]],
					 [p[3:4];0.0],[p[5:6];0.0],
					 initial_state,
					 TC_12,t_ramsey_12,N_states)
	_population_Ramsey12 = GLOQ.get_population(_rho_Ramsey_u)
	_loss = sum(abs2,_population_Ramsey12[:,1:3]-data_ramsey_12).*dt_ramsey_12/1000
	return _loss
end

global _function_call = 0
global _iteration_number = 0
function loss_gala(p,p_keywords)
	# Ramsey experiments
	_rho_Ramsey_01_u,_rho_Ramsey_01_v = GLOQ.RamseyParityForwardSolve(u0_ramsey_01,v0_ramsey_01,
					 [2.0*pi.*p[1:2];omegas[3]],omr_ramsey_01,
					 [2.0*pi.*p[7:8];char_noise[3]],
					 [p[3:4];0.0],[p[5:6];0.0],
					 0,
					 TC_01,t_ramsey_01,N_states)
	_population_Ramsey_01 = GLOQ.get_population(_rho_Ramsey_01_u)

	_rho_Ramsey_12_u,_rho_Ramsey_12_v = GLOQ.RamseyParityForwardSolve(u0_ramsey_12,v0_ramsey_12,
				     [2.0*pi.*p[1:2];omegas[3]],omr_ramsey_12,
					 [2.0*pi.*p[7:8];char_noise[3]],
					 [p[3:4];0.0],[p[5:6];0.0],
					 1,
					 TC_12,t_ramsey_12,N_states)
	_population_Ramsey_12 = GLOQ.get_population(_rho_Ramsey_12_u)

	# Echo experiments
	_rho_Echo_01_u,_rho_Echo_01_v = GLOQ.EchoParityForwardSolve(u0_echo_01,v0_echo_01,
					 [2.0*pi.*p[1:2];omegas[3]],omr_echo_01,
					 [2.0*pi.*p[7:8];char_noise[3]],
					 [p[3:4];0.0],[p[5:6];0.0],
					 0,
					 TC_01,t_echo_01,N_states)
	_population_Echo_01 = GLOQ.get_population(_rho_Echo_01_u)

	_rho_Echo_12_u,_rho_Echo_12_v = GLOQ.EchoParityForwardSolve(u0_echo_12,v0_echo_12,
					 [2.0*pi.*p[1:2];omegas[3]],omr_echo_12,
					 [2.0*pi.*p[7:8];char_noise[3]],
					 [p[3:4];0.0],[p[5:6];0.0],
					 1,
					 TC_12,t_echo_12,N_states)
	_population_Echo_12 = GLOQ.get_population(_rho_Echo_12_u)

	# T1-decay experiments
	_rho_T1_01_u,_rho_T1_01_v = GLOQ.T1ParityForwardSolve(u0_t1_01,v0_t1_01,
					 [2.0*pi.*p[1:2];omegas[3]],omr_t1_01,
					 [2.0*pi.*p[7:8];char_noise[3]],
					 [p[3:4];0.0],[p[5:6];0.0],
					 0,
					 TC_01,t_t1_01,N_states)
	_population_T1_01 = GLOQ.get_population(_rho_T1_01_u)

	_rho_T1_12_u,_rho_T1_12_v = GLOQ.T1ParityForwardSolve(u0_t1_12,v0_t1_12,
					 [2.0*pi.*p[1:2];omegas[3]],omr_t1_12,
					 [2.0*pi.*p[7:8];char_noise[3]],
					 [p[3:4];0.0],[p[5:6];0.0],
					 1,
					 TC_12,t_t1_12,N_states)
	_population_T1_12 = GLOQ.get_population(_rho_T1_12_u)

	################################################################################
	_loss = sum(abs2,_population_Ramsey_01[:,1:3]-data_ramsey_01)*dt_ramsey_01/T_evaluated_ramsey_01+
			sum(abs2,_population_Ramsey_12[:,1:3]-data_ramsey_12)*dt_ramsey_12/T_evaluated_ramsey_12
			sum(abs2,_population_Echo_01[:,1:3]-data_echo_01).*dt_echo_01/T_evaluated_echo_01+
			sum(abs2,_population_Echo_12[:,1:3]-data_echo_12).*dt_echo_12/T_evaluated_echo_12#+
			#sum(abs2,_population_T1_01[:,1:3]-data_t1_01).*dt_t1_01/T_evaluated_t1_01+
			#sum(abs2,_population_T1_12[:,1:3]-data_t1_12).*dt_t1_12/T_evaluated_t1_12

	global _function_call
	_function_call += 1
	if(_function_call%5==0)
		println("Function call: ",_function_call," Loss = ",_loss)
	end
	return _loss
end

# initial guess
p_initial = [4.10817777;3.88303940; # transition frequencies
			 gamma1[1];gamma1[2]; # parameters for the decoupling
			 gamma2[1];gamma2[2]; # parameters for the dephasing
			 -1.017e-5;1.058e-4] # charge noise
#loss_gala(p_initial,[])
# upper and lower bound for optimizaiton
perturbation = [0.02; 0.02;  # freqs
		 		 0.9; 0.9;   # gam1
		 	 	 0.9; 0.9;	   # gma2
				-2.0; 1.0]   # charge noise

lower_bound = (1.0.-perturbation).*p_initial
upper_bound = (1.0.+perturbation).*p_initial
#@time Zygote.gradient(loss, p_initial)

########################################################
# Optimization
########################################################

optimization_function_az = OptimizationFunction(loss_gala, GalacticOptim.AutoZygote())
prob_az = GalacticOptim.OptimizationProblem(optimization_function_az, p_initial,
										    lb = lower_bound, ub = upper_bound)

_function_call = 0
_iteration_number = 0

@time sol = solve(prob_az , Opt(:LD_LBFGS,length(p_initial)),
				  maxiters=50,
				  ftol_rel=1e-3)

########################################################
# Visualization
########################################################
# Echo 0-1
rho_Echo01_u,rho_Echo01_v = GLOQ.EchoParityForwardSolve(u0_echo_01,v0_echo_01,
				 [2.0*pi.*sol.u[1:2];omegas[3]],omr_echo_01,
				 [2.0*pi.*sol.u[7:8];char_noise[3]],
				 [sol.u[3:4];0.0],[sol.u[5:6];0.0],
				 0,
				 TC_01,t_echo_01,N_states)
optim_echo_01 = GLOQ.get_population(rho_Echo01_u)
optim_echo_01 = optim_echo_01[:,1:3]

# Echo 1-2
rho_Echo12_u,rho_Echo12_v = GLOQ.EchoParityForwardSolve(u0_echo_12,v0_echo_12,
				 [2.0*pi.*sol.u[1:2];omegas[3]],omr_echo_12,
				 [2.0*pi.*sol.u[7:8];char_noise[3]],
				 [sol.u[3:4];0.0],[sol.u[5:6];0.0],
				 1,
				 TC_12,t_echo_12,N_states)
optim_echo_12 = GLOQ.get_population(rho_Echo12_u)
optim_echo_12 = optim_echo_12[:,1:3]

# Ramsey 0-1
rho_Ramsey01_u,rho_Ramsey01_v = GLOQ.RamseyParityForwardSolve(u0_ramsey_01,v0_ramsey_01,
				 [2.0*pi.*sol.u[1:2];omegas[3]],omr_ramsey_01,
				 [2.0*pi.*sol.u[7:8];char_noise[3]],
				 [sol.u[3:4];0.0],[sol.u[5:6];0.0],
				 0,
				 TC_01,t_ramsey_01,N_states)
optim_ramsey_01 = GLOQ.get_population(rho_Ramsey01_u)
optim_ramsey_01 = optim_ramsey_01[:,1:3]

# Ramsey 1-2
rho_Ramsey12_u,rho_Ramsey12_v = GLOQ.RamseyParityForwardSolve(u0_ramsey_12,v0_ramsey_12,
				 [2.0*pi.*sol.u[1:2];omegas[3]],omr_ramsey_12,
				 [2.0*pi.*sol.u[7:8];char_noise[3]],
				 [sol.u[3:4];0.0],[sol.u[5:6];0.0],
				 1,
				 TC_12,t_ramsey_12,N_states)
optim_ramsey_12 = GLOQ.get_population(rho_Ramsey12_u)
optim_ramsey_12 = optim_ramsey_12[:,1:3]

# T1 0-1
rho_T101_u,rho_T101_v = GLOQ.T1ParityForwardSolve(u0_t1_01,v0_t1_01,
				 [2.0*pi.*sol.u[1:2];omegas[3]],omr_t1_01,
				 [2.0*pi.*sol.u[7:8];char_noise[3]],
				 [sol.u[3:4];0.0],[sol.u[5:6];0.0],
				 0,
				 TC_01,t_t1_01,N_states)
optim_t1_01 = GLOQ.get_population(rho_T101_u)
optim_t1_01 = optim_t1_01[:,1:3]

# T1 1-2
rho_T112_u,rho_T112_v = GLOQ.T1ParityForwardSolve(u0_t1_12,v0_t1_12,
				 [2.0*pi.*sol.u[1:2];omegas[3]],omr_t1_12,
				 [2.0*pi.*sol.u[7:8];char_noise[3]],
				 [sol.u[3:4];0.0],[sol.u[5:6];0.0],
				 1,
				 TC_12,t_t1_12,N_states)
optim_t1_12 = GLOQ.get_population(rho_T112_u)
optim_t1_12 = optim_t1_12[:,1:3]

# Plotting
fnt = Plots.font("Helvetica",16)
lfnt = Plots.font("Helvetica",12)
Plots.default(titlefont=fnt,
			  guidefont=fnt,
			  tickfont=fnt,
			  legendfont=lfnt,
			  linewidth=2);

# Echo 0-1
fig_e01 = plot(t_echo_01./1000.0,data_echo_01,line=(:dash),
		   size=(1200,600),legend=:outerright,
		   title="Echo 0-1",
		   label=[L"$\rho_{00}$, data" L"$\rho_{11}$, data" L"$\rho_{22}$, data"])
plot!(fig_e01,t_echo_01./1000.0,optim_echo_01,legend=:outerright,
		  label=[L"$\rho_{00}$, sim" L"$\rho_{11}$, sim" L"$\rho_{22}$, sim"])
display(fig_e01)

fig_e01_error = plot(t_ramsey_01./1000.0,optim_echo_01-data_echo_01,
		 size=(1200,600),legend=:outerright,
		 title = "Echo 0-1",
		 label=[L"$\rho_{00}$" L"$\rho_{11}$" L"$\rho_{22}$"])
display(fig_e01_error)

# Echo 1-2
fig_e12 = plot(t_echo_12./1000.0,data_echo_12,line=(:dash),
		   size=(1200,600),legend=:outerright,
		   title="Echo 1-2",
		   label=[L"$\rho_{00}$, data" L"$\rho_{11}$, data" L"$\rho_{22}$, data"])
plot!(fig_e12,t_echo_12./1000.0,optim_echo_12,legend=:outerright,
		  label=[L"$\rho_{00}$, sim" L"$\rho_{11}$, sim" L"$\rho_{22}$, sim"])
display(fig_e12)

fig_e12_error = plot(t_echo_12./1000.0,optim_echo_12-data_echo_12,
		 size=(1200,600),legend=:outerright,
		 title = "Echo 1-2",
		 label=[L"$\rho_{00}$" L"$\rho_{11}$" L"$\rho_{22}$"])
display(fig_e12_error)

# Ramsey 0-1
fig_r01 = plot(t_ramsey_01./1000.0,data_ramsey_01,line=(:dash),
		   size=(1200,600),legend=:outerright,
		   title="Ramsey 0-1",
		   label=[L"$\rho_{00}$, data" L"$\rho_{11}$, data" L"$\rho_{22}$, data"])
plot!(fig_r01,t_ramsey_01./1000.0,optim_ramsey_01,legend=:outerright,
		  label=[L"$\rho_{00}$, sim" L"$\rho_{11}$, sim" L"$\rho_{22}$, sim"])
display(fig_r01)

fig_r01_error = plot(t_ramsey_01./1000.0,optim_ramsey_01-data_ramsey_01,
		 size=(1200,600),legend=:outerright,
		 title = "Ramsey 0-1",
		 label=[L"$\rho_{00}$" L"$\rho_{11}$" L"$\rho_{22}$"])
display(fig_r01_error)

# Ramsey 1-2
fig_r12 = plot(t_ramsey_12./1000.0,data_ramsey_12,line=(:dash),
		   size=(1200,600),legend=:outerright,
		   title="Ramsey 1-2",
		   label=[L"$\rho_{00}$, data" L"$\rho_{11}$, data" L"$\rho_{22}$, data"])
plot!(fig_r12,t_ramsey_12./1000.0,optim_ramsey_12,legend=:outerright,
		  label=[L"$\rho_{00}$, sim" L"$\rho_{11}$, sim" L"$\rho_{22}$, sim"])
display(fig_r12)

fig_r12_error = plot(t_ramsey_12./1000.0,optim_ramsey_12-data_ramsey_12,
				   size=(1200,600),legend=:outerright,
		   		   title="Ramsey 1-2",
		 	   	   label=[L"$\rho_{00}$" L"$\rho_{11}$" L"$\rho_{22}$"])
display(fig_r12_error)

# T1 0-1
fig_t01 = plot(t_t1_01./1000.0,data_t1_01,line=(:dash),
		   size=(1200,600),legend=:outerright,
		   title="T1 0-1",
		   label=[L"$\rho_{00}$, data" L"$\rho_{11}$, data" L"$\rho_{22}$, data"])
plot!(fig_t01,t_t1_01./1000.0,optim_t1_01,legend=:outerright,
		  label=[L"$\rho_{00}$, sim" L"$\rho_{11}$, sim" L"$\rho_{22}$, sim"])
display(fig_t01)

fig_t01_error = plot(t_t1_01./1000.0,optim_t1_01-data_t1_01,
		 size=(1200,600),legend=:outerright,
		 title = "T1 0-1",
		 label=[L"$\rho_{00}$" L"$\rho_{11}$" L"$\rho_{22}$"])
display(fig_t01_error)

# T1 1-2
fig_t12 = plot(t_t1_12./1000.0,data_t1_12,line=(:dash),
		   size=(1200,600),legend=:outerright,
		   title="T1 1-2",
		   label=[L"$\rho_{00}$, data" L"$\rho_{11}$, data" L"$\rho_{22}$, data"])
plot!(fig_t12,t_t1_12./1000.0,optim_t1_12,legend=:outerright,
		  label=[L"$\rho_{00}$, sim" L"$\rho_{11}$, sim" L"$\rho_{22}$, sim"])
display(fig_t12)

fig_t12_error = plot(t_t1_12./1000.0,optim_t1_12-data_t1_12,
		 size=(1200,600),legend=:outerright,
		 title = "T1 1-2",
		 label=[L"$\rho_{00}$" L"$\rho_{11}$" L"$\rho_{22}$"])
display(fig_t12_error)

# output
println("Optimized results:",sol.u)

ext: jl


text:
function redirect_output(fun::Function)
    buf = IOBuffer()
    fun(IOContext(buf, :color => true))
    String(take!(buf))
end

function test_output(fun::Function, reference="")
    res = redirect_output(fun)
    test_esc_str_equal(res, reference)
end

@testset "Solver trace" begin
    @test format(ColumnSeparator(), 1) == "│"

    test_output() do io
        s = SolverTrace(10, io=io)
    end

    test_output("$(crayon"underline")       $(crayon"reset")\n") do io
        s = SolverTrace(10, io=io)
        print_header(s)
    end

    test_output("""[$(crayon"fg:(255,0,0)") 1/10$(crayon"reset")]
[$(crayon"fg:(255,57,0)") 2/10$(crayon"reset")]
[$(crayon"fg:(255,113,0)") 3/10$(crayon"reset")]
[$(crayon"fg:(255,170,0)") 4/10$(crayon"reset")]
[$(crayon"fg:(255,227,0)") 5/10$(crayon"reset")]
[$(crayon"fg:(227,255,0)") 6/10$(crayon"reset")]
[$(crayon"fg:(170,255,0)") 7/10$(crayon"reset")]
[$(crayon"fg:(113,255,0)") 8/10$(crayon"reset")]
[$(crayon"fg:(57,255,0)") 9/10$(crayon"reset")]
[$(crayon"fg:(0,255,0)")10/10$(crayon"reset")]
""") do io
    N = 10
    s = SolverTrace(N, io=io, progress_meter=false)
    for i = 1:N
        SolverTraces.next!(s)
    end
end

    test_output(
        """[$(crayon"fg:(255,0,0)") 1/20$(crayon"reset")]
[$(crayon"fg:(255,27,0)") 2/20$(crayon"reset")]
[$(crayon"fg:(255,81,0)") 4/20$(crayon"reset")]
[$(crayon"fg:(255,134,0)") 6/20$(crayon"reset")]
[$(crayon"fg:(255,188,0)") 8/20$(crayon"reset")]
[$(crayon"fg:(255,242,0)")10/20$(crayon"reset")]
[$(crayon"fg:(215,255,0)")12/20$(crayon"reset")]
[$(crayon"fg:(161,255,0)")14/20$(crayon"reset")]
[$(crayon"fg:(107,255,0)")16/20$(crayon"reset")]
[$(crayon"fg:(54,255,0)")18/20$(crayon"reset")]
[$(crayon"fg:(0,255,0)")20/20$(crayon"reset")]
"""
    ) do io
        N = 20
        s = SolverTrace(N, io=io)
        for i = 1:N
            SolverTraces.next!(s)
        end
    end

    test_output(
        "$(crayon"underline")Hello   World$(crayon"reset")\n   1  │    1 \n   2  │    4 \n   3  │    9 \n   4  │   16 \n   5  │   25 \n   6  │   36 \n   7  │   49 \n   8  │   64 \n   9  │   81 \n  10  │  100 \n"
    ) do io
        N = 10
        s =  SolverTrace(N,
                         ScalarColumn(identity, "Hello", Int),
                         ColumnSeparator(),
                         ScalarColumn(i -> i^2, "World", Int),
                         io=io)
        print_header(s)
        for i = 1:N
            SolverTraces.next!(s)
        end
    end

    @testset "Callbacks" begin
        N = 20
        iterations = Int[]
        iterations2 = Int[]
        s = SolverTrace(N, callbacks = ((i -> push!(iterations, i)),))
        for i = 1:N
            SolverTraces.next!(s) do
                push!(iterations2, i)
            end
        end

        @test iterations == iterations2 == [1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
    end
end


text:
struct RangeIterator
    k::Int
    d::Int
    r::Int
end

"""
    RangeIterator(n::Int,k::Int)

Returns an iterator splitting the range `1:n` into `min(k,n)` parts of (almost) equal size.
"""
RangeIterator(n::Int, k::Int) = RangeIterator(min(n,k),divrem(n,k)...)
Base.length(it::RangeIterator) = it.k
endpos(it::RangeIterator, i::Int) = i*it.d+min(i,it.r)
Base.iterate(it::RangeIterator, i::Int=1) = i>it.k ? nothing : (endpos(it,i-1)+1:endpos(it,i), i+1)

function threaded_range_loop(fun::Function, n::Integer)
    @sync for r in RangeIterator(n, Threads.nthreads())
        Threads.@spawn for i in r
            fun(i)
        end
    end
end

threaded_range_loop(fun::Function, v::AbstractVector) =
    threaded_range_loop(i -> fun(v[i]), length(v))

threaded_enumerate_range_loop(fun::Function, v::AbstractVector) =
    threaded_range_loop(i -> fun(i,v[i]), length(v))


text:
mutable struct Die
    face::Vector{String}
end

function roll!(die, dir)
    face = die.face
    if dir == 'E'
        face[1], face[3], face[4], face[6] = face[4], face[1], face[6], face[3]
    elseif dir == 'W'
        face[1], face[3], face[4], face[6] = face[3], face[6], face[1], face[4]
    elseif dir == 'N'
        face[1], face[2], face[5], face[6] = face[2], face[6], face[1], face[5]
    else
        face[1], face[2], face[5], face[6] = face[5], face[1], face[6], face[2]
    end
end

findindex(a, v) = findfirst(x -> x == a, v)
Base.copy(d::Die) = Die(copy(d.face))

function issamedice(die1, die2)
    if die1.face == die2.face
        return true
    end

    for i in "EEENEEENEEENNEEESEEESEEE"
        roll!(die2, i)
        if die1.face == die2.face
            return true
        end
    end

    return false
end

die1 = Die(split(readline()))
die2 = Die(split(readline()))
if issamedice(die1, die2)
    println("Yes")
else
    println("No")
end


text:
using Random, NelderMead, Test

Random.seed!(0)

@testset "NelderMead tests" begin
  @testset "ND Rosenbrock" begin
    rosenbrock1d(x) =  (1.0 - x[1])^2
    function rosenbrocknd(x)
      ithterm(i) = (1 - x[i])^2 + 100 * (x[i+1] - x[i]^2)^2
      return mapreduce(ithterm, +, 1:length(x)-1)
    end

    function test_solution(solution, N)
      minimum, minimiser, returncode, numiters = solution

      @test minimum <= stopval
      @test isapprox(ones(N), minimiser, rtol=1.0e-6)
      @test returncode == :STOPVAL_REACHED
    end

    stopval = eps()

    objective(x, N) = N == 1 ? rosenbrock1d(x) : rosenbrocknd(x)

    @testset "Single optimise - initial position and simplex size" begin
      for N ∈ 1:10
        solution = NelderMead.optimise(x->objective(x, N), zeros(N), ones(N),
                                      stopval=stopval, maxiters=100_000)
        test_solution(solution, N)
      end
    end
    @testset "Single optimise - vertex locations" begin
      for N ∈ 1:10
        solution = NelderMead.optimise(x->objective(x, N), [rand(N) for i ∈ 1:N+1],
                                    stopval=stopval, maxiters=100_000)
        test_solution(solution, N)
      end
    end
    @testset "Grid optimise" begin
      for N ∈ 1:3
        gridsize = [rand(1:3) for i ∈ 1:N]
        solutions = NelderMead.optimise(x->objective(x, N), zeros(N), ones(N), gridsize,
                                    stopval=stopval, maxiters=100_000)
        @assert length(solutions) >= 1
        test_solution.(solutions, N)
      end
    end
    @testset "Errors are caught" begin
      N = 2
      @test_throws ArgumentError NelderMead.optimise(x->objective(x, N), [rand(1),
                                                  rand(2), rand(3)])
      @test_throws ArgumentError NelderMead.optimise(x->objective(x, N), rand(2),
                                                  rand(3))
      @test_throws ArgumentError NelderMead.optimise(x->objective(x, N), [0.0],
                                                  [0.0], [1])
      @test_throws ArgumentError NelderMead.optimise(x->objective(x, N), rand(1),
                                                  rand(2), [1, 2, 3])
      @test_throws ArgumentError NelderMead.optimise(x->objective(x, N), rand(2),
                                                  rand(2), [0, 0])
    end

    @testset "Stretched grid causes endless loop" begin
      N = 2
      function stretchedobjective(x)
        x[2] /= 10.0
        return objective(x, N)
      end
      solution = NelderMead.optimise(stretchedobjective, zeros(N), ones(N),
                                  maxiters=1000, xtol_rel=2eps())
      _, _, returncode, numiters = solution
      @test returncode == :ENDLESS_LOOP
    end

    @testset "Nan causing endless loop" begin
      N = 2
      function nanobjective(x)
        output = objective(x, N)
        x[2] > 1.1 && (output += NaN)
        return output
      end
      solution = NelderMead.optimise(nanobjective, zeros(N), ones(N) / 2,
                                  maxiters=1000, xtol_rel=2eps())
      _, _, returncode, numiters = solution
      @test returncode == :ENDLESS_LOOP
    end

  end
end


text:
# OpenCL.Program

type Program <: CLObject
    id::CL_program
    binary::Bool

    function Program(program_id::CL_program;
                     retain::Bool=false, binary::Bool=false)
        if retain
            @check api.clRetainProgram(program_id)
        end
        p = new(program_id, binary)
        finalizer(p, _finalize)
        return p
    end
end

function _finalize(p::Program)
    if p.id != C_NULL
        @check api.clReleaseProgram(p.id)
        p.id = C_NULL
    end
end

Base.show(io::IO, p::Program) = begin
    ptr_val = convert(UInt, Base.pointer(p))
    ptr_address = "0x$(hex(ptr_val, Sys.WORD_SIZE>>2))"
    print(io, "OpenCL.Program(@$ptr_address)")
end

Base.pointer(p::Program) = p.id

Base.getindex(p::Program, pinfo::Symbol) = info(p, pinfo)

function Program(ctx::Context; source=nothing, binaries=nothing)
    local program_id::CL_program
    if source !== nothing && binaries !== nothing
        throw(ArgumentError("Program be source or binary"))
    end
    if source !== nothing
        byte_source = [String(source)]
        err_code = Ref{CL_int}()
        program_id = api.clCreateProgramWithSource(ctx.id, 1, byte_source, C_NULL, err_code)
        if err_code[] != CL_SUCCESS
            throw(CLError(err_code[]))
        end
        return Program(program_id, binary=false)

    elseif binaries !== nothing
        ndevices = length(binaries)
        device_ids = Array(CL_device_id, ndevices)
        bin_lengths = Array(Csize_t, ndevices)
        binary_status = Array(CL_int, ndevices)
        binary_ptrs= Array(Ptr{UInt8}, ndevices)
        try
            for (i, (dev, bin)) in enumerate(binaries)
                device_ids[i] = dev.id
                bin_lengths[i] = length(bin)
                binary_ptrs[i] = Base.unsafe_convert(Ptr{UInt8}, pointer(bin))
            end
            err_code = Array(CL_int, 1)
            program_id = api.clCreateProgramWithBinary(ctx.id, ndevices, device_ids, bin_lengths,
                                                       binary_ptrs, binary_status, err_code)
            if err_code[1] != CL_SUCCESS
                throw(CLError(err_code[1]))
            end
            for status in binary_status
                if status != CL_SUCCESS
                    throw(CLError(status))
                end
            end
        catch err
            throw(err)
        end
        return Program(program_id, binary=true)
    end
end

#TODO: build callback...
function build!(p::Program; options="", raise=true)
    opts = String(options)
    ndevices = 0
    device_ids = C_NULL
    @check api.clBuildProgram(p.id, cl_uint(ndevices), device_ids, opts, C_NULL, C_NULL)
    if raise
        for (dev, status) in info(p, :build_status)
            if status == CL_BUILD_ERROR
                #TODO: throw(CLBuildError(self.logs[dev], self.logs)
                error("$p build error on device $dev")
            end
        end
    end
    return p
end

let
    num_devices(p::Program) = begin
        ret = Ref{CL_uint}()
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_NUM_DEVICES, sizeof(ret), ret, C_NULL)
        return ret[]
    end

    devices(p::Program) = begin
        ndevices = num_devices(p)
        device_ids = Array(CL_device_id, ndevices)
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_DEVICES,
                                    sizeof(CL_device_id) * ndevices, device_ids, C_NULL)
        return [Device(device_ids[i]) for i in 1:ndevices]
    end

    build_status(p::Program) = begin
        status_dict = Dict{Device, CL_build_status}()
        status = Array(CL_build_status, 1)
        err_code = Array(CL_int, 1)
        for d in devices(p)
            @check api.clGetProgramBuildInfo(p.id, d.id, CL_PROGRAM_BUILD_STATUS,
                                             sizeof(CL_build_status), status, C_NULL)
            status_dict[d] = status[1]
        end
        return status_dict
    end

    build_logs(p::Program) = begin
        logs = Dict{Device, String}()
        for d in devices(p)
            log_len = Ref{Csize_t}()
            @check api.clGetProgramBuildInfo(p.id, d.id, CL_PROGRAM_BUILD_LOG,
                                             0, C_NULL, log_len)
            if log_len[] == 0
                logs[d] = ""
                continue
            end
            log_bytestring = Array(CL_char, log_len[])
            @check api.clGetProgramBuildInfo(p.id, d.id, CL_PROGRAM_BUILD_LOG,
                                            log_len[], log_bytestring, C_NULL)
            logs[d] = CLString(log_bytestring)
        end
        return logs
    end

    binaries(p::Program) = begin
        binary_dict = Dict{Device, Array{UInt8}}()
        slen = Ref{Csize_t}()
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_BINARY_SIZES,
                                    0, C_NULL, slen)

        sizes = zeros(Csize_t, slen[])
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_BINARY_SIZES,
                                    slen[], sizes, C_NULL)
        bins = Array(Ptr{UInt8}, length(sizes))
        # keep a reference to the underlying binary arrays
        # as storing the pointer to the array hides the additional
        # reference from julia's garbage collector
        bin_arrays = Any[]
        for (i, s) in enumerate(sizes)
            if s > 0
                bin = Array(UInt8, s)
                bins[i] = pointer(bin)
                push!(bin_arrays, bin)
            else
                bins[i] = Base.unsafe_convert(Ptr{UInt8}, C_NULL)
            end
        end
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_BINARIES,
                                    length(sizes) * sizeof(Ptr{UInt8}),
                                    bins, C_NULL)
        bidx = 1
        for (i, d) in enumerate(devices(p))
            if sizes[i] > 0
                binary_dict[d] = bin_arrays[bidx]
                bidx += 1
            end
        end
        return binary_dict
    end

    source(p::Program) = begin
        p.binary && throw(CLError(-45))
        src_len = Ref{Csize_t}()
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_SOURCE, 0, C_NULL, src_len)
        src_len[] <= 1 && return nothing
        src = Array(Cchar, src_len[])
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_SOURCE, src_len[], src, C_NULL)
        return CLString(src)
    end

    context(p::Program) = begin
        ret = Ref{CL_context}()
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_CONTEXT,
                                    sizeof(CL_context), ret, C_NULL)
        return Context(ret[], retain = true)
    end

    reference_count(p::Program) = begin
        ret = Ref{CL_uint}()
        @check api.clGetProgramInfo(p.id, CL_PROGRAM_REFERENCE_COUNT,
                                    sizeof(CL_uint), ret, C_NULL)
        return ret[]
    end

    const info_map = Dict{Symbol, Function}(
        :reference_count => reference_count,
        :devices => devices,
        :context => context,
        :num_devices => num_devices,
        :source => source,
        :binaries => binaries,
        :build_log => build_logs,
        :build_status => build_status,
    )

    function info(p::Program, pinfo::Symbol)
        try
            func = info_map[pinfo]
            func(p)
        catch err
            if isa(err, KeyError)
                throw(ArgumentError("OpenCL.Program has no info for $pinfo"))
            else
                throw(err)
            end
        end
    end
end

#OpenCL 1.2
#TODO: create_program_with_built_in_kernels(ctx, devices, kernel_names)
#TODO: link_program(ctx, programs; options=[], devices=None)
#TODO: unload_platform_compiler(platform)


text:
include("../src/geom_module.jl")

module pointtest
using ..geom
using PyPlot
PyPlot.pygui(true)

function test_rot_mat()

    phi = 0.0
    rot_mat = geom.RotationMatrix(phi)

    phi = pi
    rot_mat = geom.RotationMatrix(phi)

    phi = pi/2
    rot_mat = geom.RotationMatrix(phi)

    println(rot_mat)

end

function test_point()

    p1 = geom.Point{Float64}(0, 1)
    p2 = geom.Point{Float64}(1, 0)
    p3 = geom.Point{Float64}(1, 1)
    p4 = geom.Point{Float64}(0, 0)

    println("at first index: $(p1[1])" )
    println("at second index: $(p1[2])" )

    println("at first index: $(p1[1] == p1.x)" )
    println("at first index: $(p1[2] == p1.y)" )

end

function test_point2()

    p1 = geom.Point{Float64}(2.0, 1)
    p2 = geom.Point{Float64}(1, -1.0)

    p_sum = p1 + p2
    p_neg = p1 - p2
    p_prod = p1 * 10

    println(" p sum = $(p_sum) ")
    println(" p sum = $(p_sum - geom.Point(p1.x + p2.x, p1.y + p2.y) ) ")
    println(" p prod = $(p_prod) ")
    println(" p_neg = $(p_neg) ")

    rot_mat = geom.RotationMatrix(pi/2)

    p_rot = geom.rotate(p1, rot_mat)
    println(" p_rot = $(p_rot) ")

end

function test_point3()

    p1 = geom.Point(10., 3.)

    p_shift = geom.shift(p1, 10., 5.)
    p_scale = geom.scale(p1, 3)

    println(" p_shift = $(p_shift) ")
    println(" p_scale = $(p_scale) ")

end

end # module

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Synchrony, Base.Test
datadir = joinpath(Pkg.dir("Synchrony"), "test", "data")

# Test Morlet wavelet bases
#
# Reference implementation is from Torrence, C., and G. P. Compo. “A
# Practical Guide to Wavelet Analysis.” Bulletin of the American
# Meteorological Society 79, no. 1 (1998): 61–78.
#
# Our bases differ from Torrence and Compo's by a factor corresponding
# to the length of the FFT, which eliminates the need to normalize
# the output of the inverse FFT.
d1 = [wavebases(MorletWavelet([0.1], 5), 1024) * 1024; zeros(511, 1)]
d2 = map(Float64, readdlm(joinpath(datadir, "morlet_bases_f_0.1_k0_5.txt")))
@test_approx_eq d1 d2

# Test Morlet wavelet by comparing output for a 0.1 Hz oscillation
# embedded in white noise
test_in = map(Float64, readdlm(joinpath(datadir, "wavelet_test_in.txt"))[:])
foi = map(Float64, readdlm(joinpath(datadir, "wavelet_test_foi.txt"))[:])
d1 = cwt(test_in, MorletWavelet(foi, 5))
d2 = complex(map(Float64, readdlm(joinpath(datadir, "wavelet_test_out_re.txt"), '\t')),
             map(Float64, readdlm(joinpath(datadir, "wavelet_test_out_im.txt"), '\t')))
coi_periods = map(Float64, readdlm(joinpath(datadir, "wavelet_test_coi.txt"))[:])
for j = 1:length(foi)
    period = 1/foi[j]
    for i = 1:length(coi_periods)
        @test isnan(real(d1[i, j])) == (period > coi_periods[i]) || abs(period - coi_periods[i]) <= eps(coi_periods[i])
    end
end
d2[isnan(real(d1))] = NaN
@test_approx_eq d1 d2

# Test tstd
for w in (MorletWavelet([0.0003]), MorseWavelet([0.0003], 8, 3))
    wb = [wavebases(w, 100000, 3); zeros(49999)]
    p = abs2(ifftshift(bfft(wb)))
    x = (0:length(p)-1)/3
    μ = sum(p.*x)
    σ = sqrt(sum(p.*abs2(x.-μ)))
    @test_approx_eq_eps σ tstd(w) σ*1e-5
end

# Test fstd
for w in (MorletWavelet([0.75]), MorseWavelet([0.75], 8, 3))
    x = frequencies(100000, 3)[1]
    p = abs2(wavebases(w, 100000, 3))*100000
    μ = sum(p.*x)
    σ = sqrt(sum(p.*abs2(x.-μ)))
    @test_approx_eq_eps σ fstd(w) σ*1e-5
end

# Test missing data handling in wavelet transform
y = ones(1000)
y[499] = 0
y[500] = 1/3
y[501] = 2/3
foi = 2.^(0:0.25:7)
w = MorletWavelet(foi, 5)
cois = ceil(Int, tstd(w)*2*1000)
z1 = cwt(y, w, 1000)
y[500] = NaN
z2 = cwt(y, w, 1000)
for k = 1:length(foi), i = 1:length(y)
    @test isnan(real(z2[i, k])) == (i <= cois[k] || i >= length(y) - cois[k] + 1 ||
                                    500 - cois[k] <= i <= 500 + cois[k])
end
@test_approx_eq z1[!isnan(real(z2))] z2[!isnan(real(z2))]
y[501] = NaN
z2 = cwt(y, w, 1000)
for k = 1:length(foi), i = 1:length(y)
    @test isnan(real(z2[i, k])) == (i <= cois[k] || i >= length(y) - cois[k] + 1 ||
                                    500 - cois[k] <= i <= 501 + cois[k])
end
@test_approx_eq z1[!isnan(real(z2))] z2[!isnan(real(z2))]


text:
nodes_OMIB = [
    Bus(
        1, #number
        "Bus 1", #Name
        "REF", #BusType (REF, PV, PQ)
        0, #Angle in radians
        1.06, #Voltage in pu
        (min = 0.94, max = 1.06), #Voltage limits in pu
        69,
        nothing,
        nothing,
    ), #Base voltage in kV
    Bus(2, "Bus 2", "PV", 0, 1.045, (min = 0.94, max = 1.06), 69, nothing, nothing),
]

branch_OMIB = [PSY.Line(
    "Line1", #name
    true, #available
    0.0, #active power flow initial condition (from-to)
    0.0, #reactive power flow initial condition (from-to)
    Arc(from = nodes_OMIB[1], to = nodes_OMIB[2]), #Connection between buses
    0.01, #resistance in pu
    0.05, #reactance in pu
    (from = 0.0, to = 0.0), #susceptance in pu
    18.046, #rate in MW
    1.04,
)]  #angle limits (-min and max)

@testset "Dynamic Machines" begin
    Basic = BaseMachine(
        0.0, #R
        0.2995, #Xd_p
        1.05, #eq_p
        615.0,
    )  #MVABase
    @test Basic isa PowerSystems.DynamicComponent

    oneDoneQ = OneDOneQMachine(
        0.0, #R
        0.8979, #Xd
        0.646, #Xq
        0.2995, #Xd_p
        0.04, #Xq_p
        7.4, #Td0_p
        0.033, #Tq0_p
        615.0,
    )   #MVABase
    @test oneDoneQ isa PowerSystems.DynamicComponent

    AndersonFouad = AndersonFouadMachine(
        0.0, #R
        0.8979, #Xd
        0.646, #Xq
        0.2995, #Xd_p
        0.646, #Xq_p
        0.23, #Xd_pp
        0.4, #Xq_pp
        7.4, #Td0_p
        0.01, #Tq0_p #Data not available in Milano: Used 0.01
        0.03, #Td0_pp
        0.033, #Tq0_pp
        615.0,
    ) #MVABase
    @test AndersonFouad isa PowerSystems.DynamicComponent

    KundurMachine = SimpleFullMachine(
        0.003, #R on Example 3.1 and 4.1 of Kundur
        0.0006, #R_f
        0.0284, #R_1d or RD in Machowski
        0.0062, #R_1q or RQ on Machowski
        1.81, #L_d
        1.76, #L_q
        1.66, #L_ad or k*M_f or k*M_D in Machowski
        1.61, #L_aq or k*M_Q in Machowski
        1.66, #L_f1d or L_fD in Machowski. Assumed to be equal to L_ad
        1.825, #L_ff
        0.1713, #L_1d or L_D in Machowski
        0.7525, #L_1q or L_Q in Machowski
        555.0,
    ) #MVABase
    @test KundurMachine isa PowerSystems.DynamicComponent

    KundurFullMachine = FullMachine(
        0.003, #R on Example 3.1 and 4.1 of Kundur
        0.0006, #R_f
        #0.003, #R_f
        0.0284, #R_1d or RD in Machowski
        0.0062, #R_1q or RQ on Machowski
        1.81, #L_d
        1.76, #L_q
        1.66, #L_ad or k*M_f or k*M_D in Machowski
        1.61, #L_aq or k*M_Q in Machowski
        1.66, #L_f1d or L_fD in Machowski. Assumed to be equal to L_ad
        1.825, #L_ff
        0.1713, #L_1d or L_D in Machowski
        0.7525, #L_1q or L_Q in Machowski
        555.0,
    ) #MVABase
    @test KundurFullMachine isa PowerSystems.DynamicComponent

    Mach2_benchmark = OneDOneQMachine(
        0.0, #R
        1.3125, #Xd
        1.2578, #Xq
        0.1813, #Xd_p
        0.25, #Xq_p
        5.89, #Td0_p
        0.6, #Tq0_p
        100.0,
    )   #MVABase
    @test Mach2_benchmark isa PowerSystems.DynamicComponent
end

################ Shaft Data #####################
@testset "Dynamic Shaft" begin
    BaseShaft = SingleMass(
        5.148, #H
        2.0,
    ) #D
    @test BaseShaft isa PowerSystems.DynamicComponent

    FiveShaft = FiveMassShaft(
        5.148,  #H
        0.3348, #H_hp
        0.7306, #H_ip
        0.8154, #H_lp
        0.0452, #H_ex,
        2.0,    #D
        0.5180, #D_hp
        0.2240, #D_ip
        0.2240, #D_lp
        0.1450, #D_ex
        0.0518, #D_12
        0.0224, #D_23
        0.0224, #D_34
        0.0145, #D_45
        33.07,  #K_hp
        28.59,  #K_ip
        44.68,  #K_lp
        21.984,
    ) #K_ex
    @test FiveShaft isa PowerSystems.DynamicComponent
end

################# PSS Data #####################
@testset "Dynamic PSS" begin
    no_pss = PSSFixed(0.0)
    @test no_pss isa PowerSystems.DynamicComponent
end
################ TG Data #####################
@testset "Dynamic Turbine Governor Constructors" begin
    fixed_tg = TGFixed(1.0) #eff
    @test fixed_tg isa PowerSystems.DynamicComponent

    typeI_tg = TGTypeI(
        0.02, #R
        0.1, #Ts
        0.45, #Tc
        0.0, #T3
        0.0, #T4
        50.0, #T5
        0.3, #P_min
        1.2,
    ) #P_max
    @test typeI_tg isa PowerSystems.DynamicComponent

    typeII_tg = TGTypeII(
        0.05, #R
        0.3, #T1
        0.1, #T2
        1.0, #τ_max
        0.1,
    ) #τ_min
    @test typeII_tg isa PowerSystems.DynamicComponent
end
################ AVR Data #####################
@testset "Dynamic AVR Constructors" begin
    proportional_avr = AVRSimple(5000.0) #Kv
    @test proportional_avr isa PowerSystems.DynamicComponent

    fixed_avr = AVRFixed(1.05) #Emf
    @test fixed_avr isa PowerSystems.DynamicComponent

    typeI_avr = AVRTypeI(
        200.0, #Ka
        1.0, #Ke
        0.0012, #Kf
        0.02, #Ta
        0.19, #Te
        1.0, #Tf
        0.001, #Tr
        9.9, #Vr_max
        0.0, #Vr_min
        0.0006, #Ae
        0.9,
    )
    @test typeI_avr isa PowerSystems.DynamicComponent

    gen2_avr_benchmark = AVRTypeII(
        20.0, #K0 - Gain
        0.2, #T1 - 1st pole
        0.063, #T2 - 1st zero
        0.35, #T3 - 2nd pole
        0.01, #T4 - 2nd zero
        0.314, #Te - Field current time constant
        0.001, #Tr - Measurement time constant
        5.0, #Vrmax
        -5.0, #Vrmin
        0.0039, #Ae - 1st ceiling coefficient
        1.555,
    ) #Be - 2nd ceiling coefficient
    @test gen2_avr_benchmark isa PowerSystems.DynamicComponent
end
######################### Generators ########################
@testset "Dynamic Generators" begin
    #Components for the test
    Basic = BaseMachine(
        0.0, #R
        0.2995, #Xd_p
        1.05, #eq_p
        615.0,
    )  #MVABase

    BaseShaft = SingleMass(
        5.148, #H
        2.0,
    ) #D

    fixed_avr = AVRFixed(1.05) #Emf

    proportional_avr = AVRSimple(5000.0) #Kv

    fixed_tg = TGFixed(1.0) #eff

    no_pss = PSSFixed(0.0)

    oneDoneQ = OneDOneQMachine(
        0.0, #R
        0.8979, #Xd
        0.646, #Xq
        0.2995, #Xd_p
        0.04, #Xq_p
        7.4, #Td0_p
        0.033, #Tq0_p
        615.0,
    )   #MVABase

    Gen1AVR = DynamicGenerator(
        1, #Number
        "TestGen",
        nodes_OMIB[2],#bus
        1.0, # ω_ref,
        1.05,
        0.4,
        0.0,
        Basic,
        BaseShaft,
        proportional_avr, #avr
        fixed_tg, #tg
        no_pss,
    )
    @test Gen1AVR isa PowerSystems.Component

    Gen1AVRnoAVR = DynamicGenerator(
        1, #Number
        "TestGen",
        nodes_OMIB[2],#bus
        1.0, # ω_ref,
        1.05,
        0.4,
        0.0,
        Basic,
        BaseShaft,
        fixed_avr, #avr
        fixed_tg, #tg
        no_pss,
    )
    @test Gen1AVRnoAVR isa PowerSystems.Component

    Gen2AVRnoAVR = DynamicGenerator(
        1, #Number
        "TestGen",
        nodes_OMIB[2],#bus
        1.0, # ω_ref,
        1.02,
        0.4,
        0.0,
        oneDoneQ,
        BaseShaft,
        fixed_avr, #avr
        fixed_tg, #tg
        no_pss,
    )
    @test Gen2AVRnoAVR isa PowerSystems.Component

    Gen2AVR = DynamicGenerator(
        1, #Number
        "TestGen",
        nodes_OMIB[2],#bus
        1.0, # ω_ref,
        1.02,
        0.4,
        0.0,
        oneDoneQ,
        BaseShaft,
        proportional_avr, #avr
        fixed_tg, #tg
        no_pss,
    )
    @test Gen2AVR isa PowerSystems.Component

    sys = PSY.System(100)
    for bus in nodes_OMIB
        PSY.add_component!(sys, bus)
    end
    for lines in branch_OMIB
        PSY.add_component!(sys, lines)
    end
    PSY.add_component!(sys, Gen1AVR)

    @test collect(PSY.get_components(DynamicGenerator, sys))[1] == Gen1AVR

end


text:
# Autogenerated wrapper script for Capstone_jll for x86_64-unknown-freebsd
export cstool, libcapstone

JLLWrappers.@generate_wrapper_header("Capstone")
JLLWrappers.@declare_executable_product(cstool)
JLLWrappers.@declare_library_product(libcapstone, "libcapstone.so.4")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        cstool,
        "bin/cstool",
    )

    JLLWrappers.@init_library_product(
        libcapstone,
        "lib/libcapstone.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()