"""Καύσιμα με ρητά αποθέματα κατά τον Κανονισμό (ΕΕ) 965/2012 (CAT.OP.MPA.181, AMC1: βασικό σχήμα καυσίμου).

Μοντέλο καυσίμου ταξιδιού (Breguet με σταθερό «κόστος» ανόδου/προσέγγισης):
    trip(W0, r) = W0 · (1 − e^(−r/k)) + O
k (ν.μ.) και O (kg) βαθμονομούνται ανά τύπο, είτε στο διάγραμμα φορτίου-εμβέλειας του κατασκευαστή (αφού αφαιρεθούν
τα αποθέματα που υποθέτει το διάγραμμα) είτε σε δημοσιευμένα καύσιμα ταξιδιού για γνωστές αποστάσεις (ATR, De Havilland).
Ρυθμός καυσίμου σε πλεύση (kg/h) = W · V / k· σε αναμονή = συντελεστής αναμονής × ρυθμός πλεύσης.

Αποθέματα (policy), ανά πτήση με βάρος απογείωσης W0 και απόσταση r:
    απρόβλεπτα = max(cont_pct% του ταξιδιού, cont_min_hold λεπτά αναμονής)
    εναλλακτικό = trip(W_άφιξης, alt_nm)            (ματαιωμένη προσέγγιση, άνοδος, πλεύση, κάθοδος, προσέγγιση)
    τελικό απόθεμα = final_hold_min λεπτά αναμονής (με το βάρος άφιξης στο εναλλακτικό)
    + προαιρετικά flight_min λεπτά πτήσης με ρυθμό πλεύσης (π.χ. «45 min flight» της Embraer)
Το καύσιμο τροχοδρόμησης δεν μετρά στο βάρος απογείωσης (καίγεται πριν από την απογείωση).
"""
import math


def trip(w0, r, k, o):
    return w0 * (1 - math.exp(-r / k)) + o if r > 0 else o


def cruise_ff(w, v, k):
    """kg/h σε πλεύση (dW/dt = W·V/k)."""
    return w * v / k


def breakdown(w0, r, k, o, v, hold_factor, pol):
    """Καύσιμα (kg) για πτήση r ν.μ. με βάρος απογείωσης w0 και αποθέματα pol."""
    tr = trip(w0, r, k, o)
    w_dest = w0 - tr
    hold = lambda w: hold_factor * cruise_ff(w, v, k) / 60          # kg ανά λεπτό αναμονής
    cont = max(pol.get("cont_pct", 0) / 100 * tr, pol.get("cont_min_hold", 0) * hold(w_dest))
    alt = trip(w_dest, pol["alt_nm"], k, o) if pol.get("alt_nm") else 0.0
    w_alt = w_dest - alt
    final = pol.get("final_hold_min", 0) * hold(w_alt) + pol.get("flight_min", 0) * cruise_ff(w_alt, v, k) / 60
    return {"trip": tr, "contingency": cont, "alternate": alt, "final_reserve": final,
            "total": tr + cont + alt + final}


def max_range(w0, fuel, k, o, v, hold_factor, pol, r_max=8000.0):
    """Μέγιστη απόσταση (ν.μ.) ώστε ταξίδι + αποθέματα ≤ fuel. None αν δεν φτάνει ούτε για r = 0."""
    if breakdown(w0, 0, k, o, v, hold_factor, pol)["total"] > fuel:
        return None
    lo, hi = 0.0, r_max
    if breakdown(w0, hi, k, o, v, hold_factor, pol)["total"] <= fuel:
        return hi
    for _ in range(80):
        mid = (lo + hi) / 2
        if breakdown(w0, mid, k, o, v, hold_factor, pol)["total"] <= fuel:
            lo = mid
        else:
            hi = mid
    return lo


def _solve_o(w0, r, fuel, k, v, hold_factor, pol):
    """O τέτοιο ώστε το σύνολο για r να ισούται με fuel (το σύνολο αυξάνει με το O)."""
    lo, hi = -fuel, fuel
    for _ in range(80):
        mid = (lo + hi) / 2
        if breakdown(w0, r, k, mid, v, hold_factor, pol)["total"] < fuel:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2


def calibrate(points, v, hold_factor, pol):
    """points: [(w0, r_nm, fuel_kg), (w0, r_nm, fuel_kg)] με το ίδιο σύνολο αποθεμάτων pol.
    Επιστρέφει (k, O) με O ≥ 0· αν δύο σημεία δεν δίνουν O ≥ 0, O = 0 και k από το σημείο μεγαλύτερης απόστασης."""
    (wa, ra, fa), (wb, rb, fb) = sorted(points, key=lambda p: p[1])
    def resid(k):
        o = _solve_o(wa, ra, fa, k, v, hold_factor, pol)
        return breakdown(wb, rb, k, o, v, hold_factor, pol)["total"] - fb, o
    lo, hi = 500.0, 400000.0
    rlo, rhi = resid(lo)[0], resid(hi)[0]
    if rlo * rhi < 0:
        for _ in range(100):
            mid = math.sqrt(lo * hi)
            rm = resid(mid)[0]
            if (rm > 0) == (rlo > 0):
                lo, rlo = mid, rm
            else:
                hi = mid
        k = math.sqrt(lo * hi)
        o = resid(k)[1]
        if o >= 0:
            return k, o
    # εναλλακτικά: O = 0, k από το μακρύτερο σημείο
    lo, hi = 500.0, 400000.0
    for _ in range(100):
        mid = math.sqrt(lo * hi)
        if breakdown(wb, rb, mid, 0.0, v, hold_factor, pol)["total"] > fb:
            lo = mid
        else:
            hi = mid
    return math.sqrt(lo * hi), 0.0