import os
from urllib.error import HTTPError

import numpy as np
import pytest

from astropy import units as u
from astropy.constants import c
from astropy.coordinates.builtin_frames import TETE
from astropy.coordinates.earth import EarthLocation
from astropy.coordinates.funcs import get_sun
from astropy.coordinates.representation import (
    CartesianRepresentation,
    UnitSphericalRepresentation,
)
from astropy.coordinates.sky_coordinate import SkyCoord
from astropy.coordinates.solar_system import (
    BODY_NAME_TO_KERNEL_SPEC,
    _get_apparent_body_position,
    get_body,
    get_body_barycentric,
    get_body_barycentric_posvel,
    get_moon,
    solar_system_ephemeris,
)
from astropy.tests.helper import assert_quantity_allclose
from astropy.time import Time
from astropy.units import allclose as quantity_allclose
from astropy.utils.compat.optional_deps import HAS_JPLEPHEM, HAS_SKYFIELD
from astropy.utils.data import download_file, get_pkg_data_filename
from astropy.utils.exceptions import AstropyDeprecationWarning

if HAS_SKYFIELD:
    from skyfield.api import Loader

KITT_PEAK = EarthLocation.from_geodetic(
    lon=-111.6 * u.deg, lat=31.963333333333342 * u.deg, height=2120 * u.m
)


@pytest.fixture(scope="module")
def skyfield_ephemeris(tmp_path_factory):
    load = Loader(tmp_path_factory.mktemp("skyfield_ephemeris"))
    try:
        planets = load("de421.bsp")
        ts = load.timescale()
    except OSError as e:
        if os.environ.get("CI", False) and "timed out" in str(e):
            pytest.xfail("Timed out in CI")
        else:
            raise
    yield planets, ts
    planets.close()


@pytest.fixture(scope="module")
def horizons_ephemeris():
    """
    Test positions generated by JPL Horizons accessed on
    2016-03-28, with refraction turned on.
    """
    geocentric_apparent_frame = TETE(obstime=Time("1980-03-25 00:00"))
    t = Time("2014-09-25T00:00", location=KITT_PEAK)
    kitt_peak_apparent_frame = TETE(obstime=t, location=t.location)
    return {
        "geocentric": {
            "mercury": SkyCoord(
                ra="22h41m47.78s",
                dec="-08d29m32.0s",
                distance=c * 6.323037 * u.min,
                frame=geocentric_apparent_frame,
            ),
            "moon": SkyCoord(
                ra="07h32m02.62s",
                dec="+18d34m05.0s",
                distance=c * 0.021921 * u.min,
                frame=geocentric_apparent_frame,
            ),
            "jupiter": SkyCoord(
                ra="10h17m12.82s",
                dec="+12d02m57.0s",
                distance=c * 37.694557 * u.min,
                frame=geocentric_apparent_frame,
            ),
            "sun": SkyCoord(
                ra="00h16m31.00s",
                dec="+01d47m16.9s",
                distance=c * 8.294858 * u.min,
                frame=geocentric_apparent_frame,
            ),
        },
        "kitt_peak": {
            "mercury": SkyCoord(
                ra="13h38m58.50s",
                dec="-13d34m42.6s",
                distance=c * 7.699020 * u.min,
                frame=kitt_peak_apparent_frame,
            ),
            "moon": SkyCoord(
                ra="12h33m12.85s",
                dec="-05d17m54.4s",
                distance=c * 0.022054 * u.min,
                frame=kitt_peak_apparent_frame,
            ),
            "jupiter": SkyCoord(
                ra="09h09m55.55s",
                dec="+16d51m57.8s",
                distance=c * 49.244937 * u.min,
                frame=kitt_peak_apparent_frame,
            ),
        },
    }


@pytest.mark.remote_data
@pytest.mark.skipif(not HAS_SKYFIELD, reason="requires skyfield")
@pytest.mark.parametrize("body", ["mercury", "jupiter barycenter", "moon"])
def test_positions_skyfield(body, skyfield_ephemeris):
    """
    Test positions against those generated by skyfield.
    """
    planets, ts = skyfield_ephemeris
    t = Time("1980-03-25 00:00")
    frame = TETE(obstime=t)

    skyfield_t = ts.from_astropy(t)
    skyfield_coords = planets["earth"].at(skyfield_t).observe(planets[body]).apparent()
    ra, dec, dist = skyfield_coords.radec(epoch="date")

    skyfield_coords = SkyCoord(
        ra.to(u.deg), dec.to(u.deg), distance=dist.to(u.km), frame=frame
    )
    # planet positions w.r.t true equator and equinox
    astropy_coords = get_body(
        "jupiter" if body == "jupiter barycenter" else body, time=t, ephemeris="de430"
    ).transform_to(frame)

    assert astropy_coords.separation(skyfield_coords) < 1 * u.arcsec
    assert astropy_coords.separation_3d(skyfield_coords) < 10 * u.km


@pytest.mark.parametrize(
    ("body", "sep_tol", "dist_tol", "location"),
    (
        ("mercury", 7.0 * u.arcsec, 1000 * u.km, "geocentric"),
        ("jupiter", 78.0 * u.arcsec, 76000 * u.km, "geocentric"),
        ("moon", 20.0 * u.arcsec, 80 * u.km, "geocentric"),
        ("sun", 5.0 * u.arcsec, 11.0 * u.km, "geocentric"),
        ("mercury", 7.0 * u.arcsec, 500 * u.km, "kitt_peak"),
        ("jupiter", 78.0 * u.arcsec, 82000 * u.km, "kitt_peak"),
    ),
)
def test_erfa_planet(body, sep_tol, dist_tol, location, horizons_ephemeris):
    """Test predictions using erfa/plan94.

    Accuracies are maximum deviations listed in erfa/plan94.c, for Jupiter and
    Mercury, and that quoted in Meeus "Astronomical Algorithms" (1998) for the Moon.
    """
    if location == "kitt_peak":
        # Add uncertainty in position of Earth
        dist_tol += 1300 * u.km

    horizons = horizons_ephemeris[location][body]
    astropy = get_body(body, horizons.frame.obstime, ephemeris="builtin").transform_to(
        horizons.frame
    )

    assert astropy.separation(horizons) < sep_tol
    assert_quantity_allclose(astropy.distance, horizons.distance, atol=dist_tol)


@pytest.mark.remote_data
@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
@pytest.mark.parametrize(
    "body,location",
    (
        ("mercury", "geocentric"),
        ("jupiter", "geocentric"),
        ("sun", "geocentric"),
        ("moon", "geocentric"),
        ("mercury", "kitt_peak"),
        ("jupiter", "kitt_peak"),
        ("moon", "kitt_peak"),
    ),
)
def test_de432s_planet(body, location, horizons_ephemeris):
    horizons = horizons_ephemeris[location][body]
    astropy = get_body(body, horizons.frame.obstime, ephemeris="de432s").transform_to(
        horizons.frame
    )

    assert astropy.separation(horizons) < 5 * u.arcsec
    assert_quantity_allclose(astropy.distance, horizons.distance, atol=20 * u.km)


@pytest.mark.remote_data
@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
@pytest.mark.parametrize("bodyname", ("mercury", "jupiter"))
def test_custom_kernel_spec_body(bodyname):
    """
    Checks that giving a kernel specifier instead of a body name works
    """
    t = Time("2014-09-25T00:00", location=KITT_PEAK)
    coord_by_name = get_body(bodyname, t, ephemeris="de432s")
    coord_by_kspec = get_body(BODY_NAME_TO_KERNEL_SPEC[bodyname], t, ephemeris="de432s")

    assert_quantity_allclose(coord_by_name.ra, coord_by_kspec.ra)
    assert_quantity_allclose(coord_by_name.dec, coord_by_kspec.dec)
    assert_quantity_allclose(coord_by_name.distance, coord_by_kspec.distance)


@pytest.mark.remote_data
@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
def test_horizons_consistency_with_precision():
    """
    A test to compare at high precision against output of JPL horizons.

    Tests ephemerides, and conversions from ICRS to GCRS to TETE. We are aiming for
    better than 2 milli-arcsecond precision.

    We use the Moon since it is nearby, and moves fast in the sky so we are
    testing for parallax, proper handling of light deflection and aberration.
    """
    moon_data = np.loadtxt(get_pkg_data_filename("data/jpl_moon.dat"))
    loc = EarthLocation.from_geodetic(
        -67.787260 * u.deg, -22.959748 * u.deg, 5186 * u.m
    )
    times = Time("2020-04-06 00:00") + np.arange(0, 24, 1) * u.hour

    apparent_frame = TETE(obstime=times, location=loc)
    with solar_system_ephemeris.set("de430"):
        astropy = get_body("moon", times, loc).transform_to(apparent_frame)
    # JPL Horizons has a known offset (frame bias) of 51.02 mas in RA.
    usrepr = UnitSphericalRepresentation(
        moon_data[:, 0] * u.deg + 51.02376467 * u.mas, moon_data[:, 1] * u.deg
    )
    horizons = apparent_frame.realize_frame(usrepr)
    assert_quantity_allclose(astropy.separation(horizons), 0 * u.mas, atol=1.5 * u.mas)


@pytest.mark.remote_data
@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
@pytest.mark.parametrize(
    "time",
    (Time("1960-01-12 00:00"), Time("1980-03-25 00:00"), Time("2010-10-13 00:00")),
)
def test_get_sun_consistency(time):
    """
    Test that the sun from JPL and the builtin get_sun match
    """
    sun_jpl_gcrs = get_body("sun", time, ephemeris="de432s")
    assert get_sun(time).separation(sun_jpl_gcrs) < 0.1 * u.arcsec


def test_get_body_nonscalar_regression():
    """
    Test that the builtin ephemeris works with non-scalar times.

    See Issue #5069.
    """
    times = Time(["2015-08-28 03:30", "2015-09-05 10:30"])
    # the following line will raise an Exception if the bug recurs.
    get_body("moon", times, ephemeris="builtin")


def test_barycentric_pos_posvel_same():
    # Check that the two routines give identical results.
    ep1 = get_body_barycentric("earth", Time("2016-03-20T12:30:00"))
    ep2, _ = get_body_barycentric_posvel("earth", Time("2016-03-20T12:30:00"))
    np.testing.assert_array_equal(ep1.xyz, ep2.xyz)


def test_earth_barycentric_velocity_rough():
    # Check that a time near the equinox gives roughly the right result.
    ep, ev = get_body_barycentric_posvel("earth", Time("2016-03-20T12:30:00"))
    assert_quantity_allclose(ep.xyz, [-1.0, 0.0, 0.0] * u.AU, atol=0.01 * u.AU)
    expected = (
        u.Quantity([0.0 * u.one, np.cos(23.5 * u.deg), np.sin(23.5 * u.deg)])
        * -30.0
        * u.km
        / u.s
    )
    assert_quantity_allclose(ev.xyz, expected, atol=1.0 * u.km / u.s)


def test_earth_barycentric_velocity_multi_d():
    # Might as well test it with a multidimensional array too.
    t = Time("2016-03-20T12:30:00") + np.arange(8.0).reshape(2, 2, 2) * u.yr / 2.0
    ep, ev = get_body_barycentric_posvel("earth", t)
    # note: assert_quantity_allclose doesn't like the shape mismatch.
    # this is a problem with np.testing.assert_allclose.
    assert quantity_allclose(
        ep.get_xyz(xyz_axis=-1),
        [[-1.0, 0.0, 0.0], [+1.0, 0.0, 0.0]] * u.AU,
        atol=0.06 * u.AU,
    )
    expected = u.Quantity([0.0 * u.one, np.cos(23.5 * u.deg), np.sin(23.5 * u.deg)]) * (
        [[-30.0], [30.0]] * u.km / u.s
    )
    assert quantity_allclose(ev.get_xyz(xyz_axis=-1), expected, atol=2.0 * u.km / u.s)


@pytest.mark.remote_data
@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
@pytest.mark.parametrize(
    ("body", "pos_tol", "vel_tol"),
    (
        pytest.param("mercury", 1000.0 * u.km, 1.0 * u.km / u.s, id="mercury"),
        pytest.param("jupiter", 100000.0 * u.km, 2.0 * u.km / u.s, id="jupiter"),
        pytest.param("earth", 10 * u.km, 10 * u.mm / u.s, id="earth"),
        pytest.param("moon", 18 * u.km, 50 * u.mm / u.s, id="moon"),
    ),
)
def test_barycentric_velocity_consistency(body, pos_tol, vel_tol):
    # Tolerances are about 1.5 times the rms listed for plan94 and epv00,
    # except for Mercury (which nominally is 334 km rms), and the Moon
    # (which nominally is 6 km rms).
    t = Time("2016-03-20T12:30:00")
    ep, ev = get_body_barycentric_posvel(body, t, ephemeris="builtin")
    dp, dv = get_body_barycentric_posvel(body, t, ephemeris="de432s")
    assert_quantity_allclose(ep.xyz, dp.xyz, atol=pos_tol)
    assert_quantity_allclose(ev.xyz, dv.xyz, atol=vel_tol)
    # Might as well test it with a multidimensional array too.
    t = Time("2016-03-20T12:30:00") + np.arange(8.0).reshape(2, 2, 2) * u.yr / 2.0
    ep, ev = get_body_barycentric_posvel(body, t, ephemeris="builtin")
    dp, dv = get_body_barycentric_posvel(body, t, ephemeris="de432s")
    assert_quantity_allclose(ep.xyz, dp.xyz, atol=pos_tol)
    assert_quantity_allclose(ev.xyz, dv.xyz, atol=vel_tol)


@pytest.mark.remote_data
@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
@pytest.mark.parametrize(
    "time",
    (Time("1960-01-12 00:00"), Time("1980-03-25 00:00"), Time("2010-10-13 00:00")),
)
def test_url_or_file_ephemeris(time):
    # URL for ephemeris de432s used for testing:
    url = "http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de432s.bsp"

    # Pass the ephemeris directly as a URL.
    coord_by_url = get_body("earth", time, ephemeris=url)

    # Translate the URL to the cached location on the filesystem.
    # Since we just used the url above, it should already have been downloaded.
    filepath = download_file(url, cache=True)

    # Get the coordinates using the file path directly:
    coord_by_filepath = get_body("earth", time, ephemeris=filepath)

    # Using the URL or filepath should give exactly the same results:
    np.testing.assert_array_equal(coord_by_url.ra, coord_by_filepath.ra)
    np.testing.assert_array_equal(coord_by_url.dec, coord_by_filepath.dec)
    np.testing.assert_array_equal(coord_by_url.distance, coord_by_filepath.distance)


@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
def test_ephemeris_non_existing_url(monkeypatch):
    def request_invalid_url(*args, **kwargs):
        raise HTTPError(code=404, msg="Not Found", fp=None, hdrs=None, url="")

    monkeypatch.setattr("urllib.request.OpenerDirector.open", request_invalid_url)
    with pytest.raises(HTTPError, match="^HTTP Error 404: Not Found$"):
        get_body(
            "earth",
            time=Time("1960-01-12 00:00"),
            ephemeris="https://www.astropy.org/path/to/nonexisting/file.bsp",
        )


@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
@pytest.mark.parametrize(
    "ephemeris,expected_error",
    [
        pytest.param(
            "de001",
            pytest.raises(HTTPError, match="^HTTP Error 404: Not Found$"),
            marks=pytest.mark.remote_data,
            id="non_existing_JPL_ephemeris_version",
        ),
        pytest.param(
            "not_an_ephemeris",
            pytest.raises(ValueError, match="^Malformed URL: 'not_an_ephemeris'$"),
            marks=pytest.mark.remote_data,
            id="invalid_string",
        ),
        pytest.param(
            "/path/to/nonexisting/file.bsp",
            pytest.raises(
                ValueError, match="^Malformed URL: '/path/to/nonexisting/file.bsp'$"
            ),
            id="missing_local_file",
        ),
    ],
)
def test_ephemeris_wrong_input(ephemeris, expected_error):
    with expected_error:
        get_body("earth", Time("1960-01-12 00:00"), ephemeris=ephemeris)


@pytest.mark.skipif(not HAS_JPLEPHEM, reason="requires jplephem")
def test_ephemeris_local_file_not_ephemeris():
    # NOTE: This test currently leaves the file open (ResourceWarning).
    # To fix this issue, an upstream fix is required in jplephem
    # package.
    with pytest.warns(ResourceWarning), pytest.raises(ValueError, match="^file starts"):
        get_body("earth", Time("1960-01-12 00:00"), ephemeris=__file__)


def test_get_body_accounts_for_location_on_Earth():
    """Regression test for #10271"""
    t = Time(58973.534052125986, format="mjd")
    # GCRS position of ALMA at this time
    obs_p = CartesianRepresentation(
        5724535.74068625, -1311071.58985697, -2492738.93017009, u.m
    )

    icrs_sun_from_alma = _get_apparent_body_position("sun", t, "builtin", obs_p)
    icrs_sun_from_geocentre = _get_apparent_body_position(
        "sun", t, "builtin", CartesianRepresentation(0, 0, 0, u.m)
    )

    difference = (icrs_sun_from_alma - icrs_sun_from_geocentre).norm()
    assert_quantity_allclose(difference, 0.13046941 * u.m, atol=1 * u.mm)


def test_get_moon_deprecation():
    time_now = Time.now()
    with pytest.warns(
        AstropyDeprecationWarning, match=r'Use get_body\("moon"\) instead\.$'
    ):
        moon = get_moon(time_now)
    assert moon == get_body("moon", time_now)