# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/constants.py
# AU_KM = 149597870.700
# 
# AU_M = 149597870700             # per IAU 2012 Resolution B2
# 
# Path: skyfield/sgp4lib.py
# def TEME_to_ITRF(jd_ut1, rTEME, vTEME, xp=0.0, yp=0.0, fraction_ut1=0.0):
#     """Deprecated: use the TEME and ITRS frame objects instead."""
#     theta, theta_dot = theta_GMST1982(jd_ut1, fraction_ut1)
#     angular_velocity = multiply.outer(_zero_zero_minus_one, theta_dot)
# 
#     R = rot_z(-theta)
# 
#     if len(rTEME.shape) == 1:
#         rPEF = (R).dot(rTEME)
#         vPEF = (R).dot(vTEME) + _cross(angular_velocity, rPEF)
#     else:
#         rPEF = mxv(R, rTEME)
#         vPEF = mxv(R, vTEME) + _cross(angular_velocity, rPEF)
# 
#     if xp == 0.0 and yp == 0.0:
#         rITRF = rPEF
#         vITRF = vPEF
#     else:
#         W = (rot_x(yp)).dot(rot_y(xp))
#         rITRF = (W).dot(rPEF)
#         vITRF = (W).dot(vPEF)
#     return rITRF, vITRF
# 
# Path: skyfield/timelib.py
# def julian_date(year, month=1, day=1, hour=0, minute=0, second=0.0):
#     """Given a proleptic Gregorian calendar date and time, build a Julian date.
# 
#     The difference between a “Julian day” and a “Julian date” is that
#     the “day” is the integer part, while the “date” includes a fraction
#     indicating the time.
# 
#     """
#     return julian_day(year, month, day) - 0.5 + (
#         second + minute * 60.0 + hour * 3600.0) / DAY_S
# 
# Path: skyfield/constants.py
# DEG2RAD = 0.017453292519943296


from numpy import array
from skyfield import api
from skyfield.api import EarthSatellite, load
from skyfield.constants import AU_KM, AU_M
from skyfield.sgp4lib import TEME_to_ITRF
from skyfield.timelib import julian_date
from ..constants import DEG2RAD

line1 = '1 25544U 98067A   18184.80969102  .00001614  00000-0  31745-4 0  9993'
line2 = '2 25544  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106'

# Here are numbers from HORIZONS, which I copied into the test below:
#
#Ephemeris / WWW_USER Wed Jul  4 19:16:45 2018 Pasadena, USA      / Horizons
#...
#2458303.500000000 = A.D. 2018-Jul-04 00:00:00.0000 TDB
# X = 2.633404251158200E-05 Y = 1.015087620439817E-05 Z = 3.544778677556393E-05
# VX=-1.751248694205384E-03 VY= 4.065407557020968E-03 VZ= 1.363540232307603E-04
#2458304.500000000 = A.D. 2018-Jul-05 00:00:00.0000 TDB
# X =-2.136440257814821E-05 Y =-2.084170814514480E-05 Z =-3.415494123796893E-05
# VX= 2.143876266215405E-03 VY=-3.752167957502106E-03 VZ= 9.484159290242074E-04

# TODO: try with array of dates

def test_iss_against_horizons():
    ts = api.load.timescale()
    s = EarthSatellite(line1, line2)

    hp = array([
        [2.633404251158200E-5, 1.015087620439817E-5, 3.544778677556393E-5],
        [-2.136440257814821E-5, -2.084170814514480E-5, -3.415494123796893E-5],
    ]).T
    hv = array([
        [-1.751248694205384E-3, 4.065407557020968E-3, 1.363540232307603E-4],
        [2.143876266215405E-3, -3.752167957502106E-3, 9.484159290242074E-4],
    ]).T


# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/constants.py
# AU_KM = 149597870.700
# 
# AU_M = 149597870700             # per IAU 2012 Resolution B2
# 
# Path: skyfield/sgp4lib.py
# def TEME_to_ITRF(jd_ut1, rTEME, vTEME, xp=0.0, yp=0.0, fraction_ut1=0.0):
#     """Deprecated: use the TEME and ITRS frame objects instead."""
#     theta, theta_dot = theta_GMST1982(jd_ut1, fraction_ut1)
#     angular_velocity = multiply.outer(_zero_zero_minus_one, theta_dot)
# 
#     R = rot_z(-theta)
# 
#     if len(rTEME.shape) == 1:
#         rPEF = (R).dot(rTEME)
#         vPEF = (R).dot(vTEME) + _cross(angular_velocity, rPEF)
#     else:
#         rPEF = mxv(R, rTEME)
#         vPEF = mxv(R, vTEME) + _cross(angular_velocity, rPEF)
# 
#     if xp == 0.0 and yp == 0.0:
#         rITRF = rPEF
#         vITRF = vPEF
#     else:
#         W = (rot_x(yp)).dot(rot_y(xp))
#         rITRF = (W).dot(rPEF)
#         vITRF = (W).dot(vPEF)
#     return rITRF, vITRF
# 
# Path: skyfield/timelib.py
# def julian_date(year, month=1, day=1, hour=0, minute=0, second=0.0):
#     """Given a proleptic Gregorian calendar date and time, build a Julian date.
# 
#     The difference between a “Julian day” and a “Julian date” is that
#     the “day” is the integer part, while the “date” includes a fraction
#     indicating the time.
# 
#     """
#     return julian_day(year, month, day) - 0.5 + (
#         second + minute * 60.0 + hour * 3600.0) / DAY_S
# 
# Path: skyfield/constants.py
# DEG2RAD = 0.017453292519943296


# Path: skyfield/tests/test_earth_satellites.py
from numpy import array
from skyfield import api
from skyfield.api import EarthSatellite, load
from skyfield.constants import AU_KM, AU_M
from skyfield.sgp4lib import TEME_to_ITRF
from skyfield.timelib import julian_date
from ..constants import DEG2RAD


line1 = '1 25544U 98067A   18184.80969102  .00001614  00000-0  31745-4 0  9993'
line2 = '2 25544  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106'

# Here are numbers from HORIZONS, which I copied into the test below:
#
#Ephemeris / WWW_USER Wed Jul  4 19:16:45 2018 Pasadena, USA      / Horizons
#...
#2458303.500000000 = A.D. 2018-Jul-04 00:00:00.0000 TDB
# X = 2.633404251158200E-05 Y = 1.015087620439817E-05 Z = 3.544778677556393E-05
# VX=-1.751248694205384E-03 VY= 4.065407557020968E-03 VZ= 1.363540232307603E-04
#2458304.500000000 = A.D. 2018-Jul-05 00:00:00.0000 TDB
# X =-2.136440257814821E-05 Y =-2.084170814514480E-05 Z =-3.415494123796893E-05
# VX= 2.143876266215405E-03 VY=-3.752167957502106E-03 VZ= 9.484159290242074E-04

# TODO: try with array of dates

def test_iss_against_horizons():
    ts = api.load.timescale()
    s = EarthSatellite(line1, line2)

    hp = array([
        [2.633404251158200E-5, 1.015087620439817E-5, 3.544778677556393E-5],
        [-2.136440257814821E-5, -2.084170814514480E-5, -3.415494123796893E-5],
    ]).T
    hv = array([
        [-1.751248694205384E-3, 4.065407557020968E-3, 1.363540232307603E-4],
        [2.143876266215405E-3, -3.752167957502106E-3, 9.484159290242074E-4],
    ]).T


    two_meters = 2.0 / AU_M