# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/constants.py
# AU_M = 149597870700             # per IAU 2012 Resolution B2
# 
# DEG2RAD = 0.017453292519943296
# 
# Path: skyfield/earthlib.py
# def reverse_terra(xyz_au, gast, iterations=3):
#     """Deprecated conversion from GCRS -> lat,lon,t; neglects polar motion."""
#     x, y, z = xyz_au
#     R = sqrt(x*x + y*y)
# 
#     lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi
#     lat = arctan2(z, R)
# 
#     a = ERAD / AU_M
#     f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
#     e2 = 2.0*f - f*f
#     i = 0
#     C = 1.0
#     while i < iterations:
#         i += 1
#         C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
#         lat = arctan2(z + a * C * e2 * sin(lat), R)
#     elevation_m = ((R / cos(lat)) - a * C) * AU_M
#     return lat, lon, elevation_m


from numpy import einsum
from skyfield.api import Topos, load
from skyfield.constants import AU_M, DEG2RAD
from skyfield.earthlib import reverse_terra
"""How accurate is Skyfield's subpoint computation routine?

Let's call it with a varying number of iterations, and see how
accurately it can turn a Topos-generated position back into a latitude,
longitude, and elevation.

"""

def main():

# Path: skyfield/api.py
# N = E = +1.0
# S = W = -1.0
# 
# Path: skyfield/constants.py
# AU_M = 149597870700             # per IAU 2012 Resolution B2
# 
# DEG2RAD = 0.017453292519943296
# 
# Path: skyfield/earthlib.py
# def reverse_terra(xyz_au, gast, iterations=3):
#     """Deprecated conversion from GCRS -> lat,lon,t; neglects polar motion."""
#     x, y, z = xyz_au
#     R = sqrt(x*x + y*y)
# 
#     lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi
#     lat = arctan2(z, R)
# 
#     a = ERAD / AU_M
#     f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
#     e2 = 2.0*f - f*f
#     i = 0
#     C = 1.0
#     while i < iterations:
#         i += 1
#         C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
#         lat = arctan2(z + a * C * e2 * sin(lat), R)
#     elevation_m = ((R / cos(lat)) - a * C) * AU_M
#     return lat, lon, elevation_m


# Path: design/subpoint_accuracy.py
from numpy import einsum
from skyfield.api import Topos, load
from skyfield.constants import AU_M, DEG2RAD
from skyfield.earthlib import reverse_terra

"""How accurate is Skyfield's subpoint computation routine?

Let's call it with a varying number of iterations, and see how
accurately it can turn a Topos-generated position back into a latitude,
longitude, and elevation.

"""

def main():

    ts = load.timescale()