"""Shared constants and helpers for the seta_haystack study."""
import os
import ssl
import time
import urllib.request
import urllib.parse

import numpy as np

HERE = os.path.dirname(os.path.abspath(__file__))
STUDY = os.path.dirname(HERE)
DATA = os.path.join(STUDY, "data")
FIGS = os.path.join(STUDY, "figures")
RAW = os.environ.get("SETA_RAW", "fermi_raw/seta_haystack")
for _d in (DATA, FIGS, RAW):
    os.makedirs(_d, exist_ok=True)

AU_KM = 1.495978707e8
SEED = 20261012

# ----------------------------------------------------------------------------
# downloads
# ----------------------------------------------------------------------------

def fetch(url, dest, timeout=900, min_bytes=100, retries=3):
    """Download url to dest unless dest already exists (non-empty)."""
    if os.path.exists(dest) and os.path.getsize(dest) >= min_bytes:
        return dest
    last = None
    for k in range(retries):
        try:
            req = urllib.request.Request(url, headers={"User-Agent": "seta-haystack-study/1.0"})
            with urllib.request.urlopen(req, timeout=timeout) as r, open(dest + ".part", "wb") as f:
                while True:
                    b = r.read(1 << 20)
                    if not b:
                        break
                    f.write(b)
            os.replace(dest + ".part", dest)
            return dest
        except Exception as ex:  # noqa
            last = ex
            time.sleep(5 * (k + 1))
    raise RuntimeError(f"download failed for {url}: {last}")


def fetch_text(url, timeout=120, retries=3):
    last = None
    for k in range(retries):
        try:
            req = urllib.request.Request(url, headers={"User-Agent": "seta-haystack-study/1.0"})
            with urllib.request.urlopen(req, timeout=timeout) as r:
                return r.read().decode("utf-8", "replace")
        except Exception as ex:  # noqa
            last = ex
            time.sleep(3 * (k + 1))
    raise RuntimeError(f"fetch failed for {url}: {last}")


# ----------------------------------------------------------------------------
# photometry
# ----------------------------------------------------------------------------

def H_from_D(D_m, p):
    """Absolute magnitude of a sphere of diameter D_m (metres), geometric albedo p."""
    return 5.0 * np.log10(1329.0e3 / (np.asarray(D_m, float) * np.sqrt(p)))


def D_from_H(H, p):
    """Diameter in metres."""
    return 1329.0e3 / np.sqrt(p) * 10 ** (-np.asarray(H, float) / 5.0)


def phase_HG(alpha_rad, G=0.15):
    """HG phase function correction in magnitudes (>=0), Bowell et al. 1989."""
    t = np.tan(np.clip(alpha_rad, 0, np.pi * 0.999) / 2.0)
    phi1 = np.exp(-3.33 * t ** 0.63)
    phi2 = np.exp(-1.87 * t ** 1.22)
    return -2.5 * np.log10(np.clip((1 - G) * phi1 + G * phi2, 1e-12, None))


def trailing_loss(rate_deg_day, t_exp=None, seeing=None):
    """Loss of limiting magnitude (mag, positive = harder to detect) versus sky rate w,
    using the empirical fit for the Catalina Sky Survey 1.5-m telescope by Nesvorny et al.
    (2024, NEOMOD2, eqs 3-4): Vlim(w) = Vlim(0) + A w for w < w1, and
    Vlim(0) + A w1 - 2.5 log10[1 + C (w - w1)] above, with A = 0.052, C = 0.192,
    w1 = 3.6 deg/day.  (Their fit covers w < 10 deg/day; we extrapolate beyond.)
    t_exp and seeing are unused placeholders kept for interface stability."""
    w = np.asarray(rate_deg_day, float)
    A, C, w1 = 0.052, 0.192, 3.6
    gain = np.where(w < w1, A * w, A * w1 - 2.5 * np.log10(1.0 + C * np.clip(w - w1, 0, None)))
    return -gain


def trailing_loss_theory(rate_deg_day, t_exp=30.0, seeing=2.5):
    """Alternative: Veres & Chesley (2017) detection trailing loss,
    dm = 1.25 log10(1 + 0.42 x^2), x = rate*t_exp/(24*seeing)."""
    x = np.asarray(rate_deg_day) * t_exp / (24.0 * seeing)
    return 1.25 * np.log10(1.0 + 0.42 * x * x)


# ----------------------------------------------------------------------------
# two-body propagation (heliocentric, ecliptic frame, units AU / days)
# ----------------------------------------------------------------------------
K_GAUSS = 0.01720209895


def kepler_E(M, e, n_iter=14):
    """Solve Kepler's equation by Newton iteration (robust starter for high e)."""
    M = np.mod(M + np.pi, 2 * np.pi) - np.pi
    e = np.broadcast_to(e, M.shape)
    E = np.where(e < 0.8, M + e * np.sin(M), np.pi * np.where(M >= 0, 1.0, -1.0))
    for _ in range(n_iter):
        E = E - (E - e * np.sin(E) - M) / (1 - e * np.cos(E))
    return E


def helio_xyz(a, e, inc, Om, w, M):
    """Heliocentric ecliptic xyz (AU). Angles in rad. Arrays broadcast."""
    E = kepler_E(M, e)
    xp = a * (np.cos(E) - e)
    yp = a * np.sqrt(1 - e * e) * np.sin(E)
    cO, sO = np.cos(Om), np.sin(Om)
    cw, sw = np.cos(w), np.sin(w)
    ci, si = np.cos(inc), np.sin(inc)
    x = (cO * cw - sO * sw * ci) * xp + (-cO * sw - sO * cw * ci) * yp
    y = (sO * cw + cO * sw * ci) * xp + (-sO * sw + cO * cw * ci) * yp
    z = (sw * si) * xp + (cw * si) * yp
    return x, y, z


def observe(xyz_obj, xyz_earth):
    """Return r, delta, phase angle, solar elongation (rad) for object and Earth
    heliocentric positions (tuples of arrays)."""
    x, y, z = xyz_obj
    ex, ey, ez = xyz_earth
    dx, dy, dz = x - ex, y - ey, z - ez
    r = np.sqrt(x * x + y * y + z * z)
    d = np.sqrt(dx * dx + dy * dy + dz * dz)
    re = np.sqrt(ex * ex + ey * ey + ez * ez)
    cos_alpha = (x * dx + y * dy + z * dz) / (r * d)
    cos_elong = -(ex * dx + ey * dy + ez * dz) / (re * d)
    alpha = np.arccos(np.clip(cos_alpha, -1, 1))
    elong = np.arccos(np.clip(cos_elong, -1, 1))
    return r, d, alpha, elong, (dx / d, dy / d, dz / d)


def sky_rate_deg_day(unit, dt_days=1.0):
    """Angular rate (deg/day) from unit vectors sampled along the last axis."""
    ux, uy, uz = unit
    dot = ux[..., 1:] * ux[..., :-1] + uy[..., 1:] * uy[..., :-1] + uz[..., 1:] * uz[..., :-1]
    ang = np.degrees(np.arccos(np.clip(dot, -1, 1))) / dt_days
    ang = np.concatenate([ang[..., :1], ang], axis=-1)
    return ang