import json
import math
import random
from datetime import datetime, timedelta, timezone

import numpy as np
from skyfield.api import load, wgs84

SIDEREAL_DAY_SECONDS = 86164.0905
EARTH_OMEGA = 2.0 * math.pi / SIDEREAL_DAY_SECONDS

MASTER_JSON = "master_ephemeris.json"
ts = load.timescale()
eph = load("de440s.bsp")
earth = eph["earth"]

with open(MASTER_JSON, "r") as f:
    master = json.load(f)

START_TIME = datetime.fromisoformat(master["start_time_utc"])
END_TIME = datetime.fromisoformat(master["end_time_utc"])

COEFF_RANGE = float(master["coeff_range"])
INT16_SCALE = 32767.0 / COEFF_RANGE

THETA0 = float(master["earth_theta0_rad"])

R_SLOW = np.asarray(
    master["earth_r_slow"],
    dtype=np.float64
)

BODY_NAMES = list(master["bodies"].keys())


def random_datetime(start, end):
    return start + timedelta(seconds=random.random() * (end - start).total_seconds())


def chebyshev_eval_int16(coefficients, u):
    """
    this functions should really line in its own file.
    """
    b1 = 0.0
    b2 = 0.0
    for i in range(len(coefficients) - 1, 0, -1):
        c = float(coefficients[i]) / INT16_SCALE
        b0 = (2.0 * u * b1 - b2 + c)
        b2 = b1
        b1 = b0
    c0 = float(coefficients[0]) / INT16_SCALE
    return u * b1 - b2 + c0


def master_icrf_vector(body_name, utc_datetime):
    """
    build geocentric ECI/ICRF direction
    """
    body = master["bodies"][body_name]
    segment_days = int(body["segment_days"])
    segment_seconds = segment_days * 86400.0
    elapsed_seconds = (utc_datetime - START_TIME).total_seconds()
    # if elapsed_seconds < 0:
    #     raise ValueError("time before ephemeris start")
    # if utc_datetime > END_TIME:
    #     raise ValueError("time after ephemeris end")
    segment_index = int(elapsed_seconds // segment_seconds)
    segments = body["segments"]
    if segment_index >= len(segments):
        segment_index = len(segments) - 1
    segment_start_seconds = (segment_index * segment_seconds)
    seconds_into_segment = (elapsed_seconds - segment_start_seconds)
    fraction = (seconds_into_segment / segment_seconds)
    u = (2.0 * fraction - 1.0)
    segment = segments[segment_index]

    x = chebyshev_eval_int16(segment["x"], u)
    y = chebyshev_eval_int16(segment["y"], u)
    z = chebyshev_eval_int16(segment["z"], u)

    vector = np.array([x, y, z], dtype=np.float64)
    vector /= np.linalg.norm(vector)
    return vector


def icrf_to_ecef(icrf, utc_datetime):
    dt_seconds = (utc_datetime - START_TIME).total_seconds()
    theta = (THETA0 + EARTH_OMEGA * dt_seconds)
    theta %= 2.0 * math.pi
    c = math.cos(theta)
    s = math.sin(theta)
    rz = np.array([
        [c, s, 0.0],
        [-s, c, 0.0],
        [0.0, 0.0, 1.0]], dtype=np.float64)
    ecef = (rz @ R_SLOW @ icrf)
    ecef /= np.linalg.norm(ecef)
    return ecef


def moon_topocentric_ecef(
        moon_icrf_unit,
        utc_datetime,
        latitude_deg,
        longitude_deg
):
    """
    I found the moon error was too large. After some trail and error and research,
    some of the error comes from topocentric parallax. That's when the moon is technicality
    overhead (same half side of the sphere as the observer) but the observer is still blocked
    by the curved earth. Think if you put a very big building at lat 0, lon 0
    and you are on the north pole. technicality the building top is radially above you,
    but the horizon is still in your way. The other planets do not suffer from this as they
    are so far away
    Use:
        average Moon distance = 384400 km
        Earth radius          = 6371 km
    Using these values makes the moon parallax issue drop to ~zero
    """

    AU_KM = 149597870.7
    MOON_DISTANCE_AU = 384400.0 / AU_KM
    EARTH_RADIUS_AU = 6371.0 / AU_KM
    moon_ecef_unit = icrf_to_ecef(moon_icrf_unit, utc_datetime)
    moon_ecef = (moon_ecef_unit * MOON_DISTANCE_AU)
    lat = math.radians(latitude_deg)
    lon = math.radians(longitude_deg)
    cos_lat = math.cos(lat)
    sin_lat = math.sin(lat)
    cos_lon = math.cos(lon)
    sin_lon = math.sin(lon)
    observer_ecef = np.array([
        EARTH_RADIUS_AU * cos_lat * cos_lon,
        EARTH_RADIUS_AU * cos_lat * sin_lon,
        EARTH_RADIUS_AU * sin_lat], dtype=np.float64)
    topocentric = (moon_ecef - observer_ecef)
    # need direction after the subtraction.
    topocentric /= np.linalg.norm(topocentric)
    return topocentric


def skyfield_icrf_vector(body_name, utc_datetime):
    body_info = master["bodies"][body_name]
    skyfield_name = body_info.get("skyfield_name", body_name)
    target = eph[skyfield_name]
    t = ts.from_datetime(utc_datetime)
    observation = earth.at(t).observe(target)
    xyz = np.asarray(observation.position.au, dtype=np.float64)
    xyz /= np.linalg.norm(xyz)
    return xyz


def ecef_to_altaz(ecef, latitude_deg, longitude_deg):
    """may only need to do a dot product with perfect marble sphere earth"""
    lat = math.radians(latitude_deg)
    lon = math.radians(longitude_deg)
    sin_lat = math.sin(lat)
    cos_lat = math.cos(lat)
    sin_lon = math.sin(lon)
    cos_lon = math.cos(lon)
    east = np.array([-sin_lon, cos_lon, 0.0])
    north = np.array([-sin_lat * cos_lon, -sin_lat * sin_lon, cos_lat])
    up = np.array([cos_lat * cos_lon, cos_lat * sin_lon, sin_lat])
    east_value = float(np.dot(ecef, east))
    north_value = float(np.dot(ecef, north))
    up_value = float(np.dot(ecef, up))
    up_value = np.clip(up_value, -1.0, 1.0)
    altitude = math.degrees(math.asin(up_value))
    azimuth = math.degrees(math.atan2(east_value, north_value))
    if azimuth < 0.0:
        azimuth += 360.0
    return altitude, azimuth


def master_altaz(body_name, utc_datetime, latitude_deg, longitude_deg):
    icrf = master_icrf_vector(body_name, utc_datetime)
    if body_name == "moon":
        ecef = moon_topocentric_ecef(icrf, utc_datetime, latitude_deg, longitude_deg)
    else:
        ecef = icrf_to_ecef(icrf, utc_datetime)
    return ecef_to_altaz(ecef, latitude_deg, longitude_deg)


def skyfield_altaz(body_name, utc_datetime, latitude_deg, longitude_deg):
    body_info = master["bodies"][body_name]
    skyfield_name = body_info.get("skyfield_name", body_name)
    target = eph[skyfield_name]
    observer = wgs84.latlon(latitude_deg, longitude_deg)
    t = ts.from_datetime(utc_datetime)
    astrometric = (earth + observer).at(t).observe(target)
    apparent = astrometric.apparent()
    alt, az, distance = apparent.altaz()
    return alt.degrees, az.degrees, distance.au


## test
# body = random.choice(BODY_NAMES)
body = "moon"
utc_datetime = random_datetime(START_TIME, END_TIME)
latitude = random.uniform(-89.0, 89.0)  ## if you live on the north pole, please let me know if you have buit this
longitude = random.uniform(-180.0, 180.0)

utc_datetime = datetime.now(timezone.utc)
latitude = 38.861852
longitude = -77.048736

approx_icrf = master_icrf_vector(body, utc_datetime)

truth_icrf = skyfield_icrf_vector(body, utc_datetime)

approx_alt, approx_az = master_altaz(body, utc_datetime, latitude, longitude)

truth_alt, truth_az, distance_au = skyfield_altaz(body, utc_datetime, latitude, longitude)


def angular_error_deg(a, b):
    a = np.asarray(a, dtype=np.float64)
    b = np.asarray(b, dtype=np.float64)

    a /= np.linalg.norm(a)
    b /= np.linalg.norm(b)

    dot = np.dot(a, b)
    dot = np.clip(dot, -1.0, 1.0)

    return math.degrees(math.acos(dot))


print(f"testing {body}")
icrf_error = angular_error_deg(approx_icrf, truth_icrf)
print(f"angular error {icrf_error}")
print(f"skyfield overhead {truth_alt > 0.0}, approx overhead {approx_alt > 0.0}")

# utc_datetime = random_datetime(START_TIME, END_TIME)
# latitude = random.uniform(-89.0, 89.0)
# longitude = random.uniform(-180.0, 180.0)


print(f"time:      {utc_datetime.isoformat()}")
print(f"latitude:  {latitude:.6f}")
print(f"longitude: {longitude:.6f}")
print()

print(
    f"{'body':<10} "
    f"{'approx alt':>11} "
    f"{'sky alt':>11} "
    f"{'approx OH':>10} "
    f"{'sky OH':>8} "
    f"{'match':>7} "
    f"{'ICRF err':>10}"
)
print("-" * 75)

for body in BODY_NAMES:
    approx_icrf = master_icrf_vector(body, utc_datetime)
    truth_icrf = skyfield_icrf_vector(body, utc_datetime)

    approx_alt, approx_az = master_altaz(body, utc_datetime, latitude, longitude)

    truth_alt, truth_az, distance_au = skyfield_altaz(body, utc_datetime, latitude, longitude)

    icrf_error = angular_error_deg(approx_icrf, truth_icrf)

    approx_overhead = approx_alt > 0.0
    truth_overhead = truth_alt > 0.0
    match = approx_overhead == truth_overhead

    print(
        f"{body:<10} "
        f"{approx_alt:>11.4f} "
        f"{truth_alt:>11.4f} "
        f"{str(approx_overhead):>10} "
        f"{str(truth_overhead):>8} "
        f"{str(match):>7} "
        f"{icrf_error:>9.5f}°"
    )

# print(master_icrf_vector("moon", utc_datetime))
# print(icrf_to_ecef(master_icrf_vector("moon", utc_datetime), utc_datetime))
# approx_alt, approx_az = master_altaz("moon", utc_datetime, latitude, longitude)
# print(approx_alt)