import numpy as np

R_TERRE = 6371.0  # Rayon terrestre en km
ALTITUDE_SAT = 700.0  # Altitude des satellites en km
RAYON_ORBITE = R_TERRE + ALTITUDE_SAT
NOMBRE_SATELLITES = 3  # Taille de la constellation


# FONCTIONS GÉOMÉTRIQUES
def conversion_geographique_vers_cartesien(latitude_deg, longitude_deg, rayon=R_TERRE):
    """
    Convertit des coordonnées géographiques en coordonnées cartésiennes.
    """
    latitude_rad = np.radians(latitude_deg)
    longitude_rad = np.radians(longitude_deg)
    
    x = rayon * np.cos(latitude_rad) * np.cos(longitude_rad)
    y = rayon * np.cos(latitude_rad) * np.sin(longitude_rad)
    z = rayon * np.sin(latitude_rad)
    
    return np.array([x, y, z])

# Points terrestres à surveiller (coordonnées en degrés)
POINTS_TERRESTRES = np.array([
    conversion_geographique_vers_cartesien(0, 0),   
    conversion_geographique_vers_cartesien(20, 30), 
    conversion_geographique_vers_cartesien(-40, 60),
    conversion_geographique_vers_cartesien(20, -120), 
    conversion_geographique_vers_cartesien(55, -80)    
])


# FONCTIONS DE VISIBILITÉ
def est_visible_depuis_satellite(position_satellite, point_terrestre):
    """
    Détermine si un point terrestre est visible depuis un satellite.
    """
    # Vecteur normal à la surface terrestre au point considéré
    normale_locale = point_terrestre / np.linalg.norm(point_terrestre)
    
    # Vecteur allant du point terrestre vers le satellite
    vecteur_satellite_point = position_satellite - point_terrestre
    
    # Cosinus de l'angle entre la normale et le vecteur satellite→point
    cosinus_angle = np.dot(normale_locale, vecteur_satellite_point) / np.linalg.norm(vecteur_satellite_point)
    
    # Le point est visible si l'angle < 90° (cosinus > 0)
    return cosinus_angle > 0

def calculer_couverture(positions_satellites):
    """
    Calcule le nombre de points terrestres couverts par au moins un satellite.
    """
    nombre_points_couverts = 0
    
    for point in POINTS_TERRESTRES:
        # Vérifie si au moins un satellite voit ce point
        if any(est_visible_depuis_satellite(sat, point) for sat in positions_satellites):
            nombre_points_couverts += 1
            
    return nombre_points_couverts


# FONCTIONS DE CONVERSION SOLUTION ↔ POSITIONS
def conversion_solution_vers_positions(solution):
    """
    Convertit un vecteur solution en positions réelles de satellites.
    """
    positions = []
    
    for i in range(NOMBRE_SATELLITES):
        longitude = solution[i]
        latitude = solution[i + NOMBRE_SATELLITES]
        
        # Conversion coordo sphériques to cartésiennes
        x = RAYON_ORBITE * np.cos(latitude) * np.cos(longitude)
        y = RAYON_ORBITE * np.cos(latitude) * np.sin(longitude)
        z = RAYON_ORBITE * np.sin(latitude)
        
        positions.append(np.array([x, y, z]))
    
    return positions

def evaluer_solution(solution):
    """
    Fonction d'évaluation d'une solution (à maximiser).
    """
    positions = conversion_solution_vers_positions(solution)
    return calculer_couverture(positions)


# ALGORITHME 


def recherche_locale_tortue(solution_initiale, iterations_max=2500, pas_initial=0.4):
    """
    Implémente la recherche locale inspirée du comportement des tortues.
    """
    # Initialisation
    meilleure_solution = solution_initiale.copy()
    meilleure_valeur = evaluer_solution(meilleure_solution)
    pas_courant = pas_initial
    compteur_stagnation = 0
    
    print(f"  Debut recherche locale: couverture initiale = {meilleure_valeur}")
    
    # Boucle d'optimisation principale
    for iteration in range(iterations_max):
        # Génération d'un voisin par perturbation gaussienne
        voisin = meilleure_solution + np.random.normal(scale=pas_courant, size=meilleure_solution.shape)
        
        # Application des contraintes
        voisin[:NOMBRE_SATELLITES] %= (2 * np.pi)  # Longitudes périodiques
        voisin[NOMBRE_SATELLITES:] = np.clip(voisin[NOMBRE_SATELLITES:], -np.pi/2, np.pi/2)  # Latitudes bornées
        
        # Évaluation du voisin
        valeur_voisin = evaluer_solution(voisin)
        
        # Mise à jour si amélioration (MAXIMISATION)
        if valeur_voisin > meilleure_valeur:
            meilleure_solution = voisin.copy()
            meilleure_valeur = valeur_voisin
            compteur_stagnation = 0
        else:
            compteur_stagnation += 1
        
        # Mécanisme TIOA : ajustement adaptatif du pas
        if compteur_stagnation > 15:
            pas_courant *= 0.6
            compteur_stagnation = 0
            
            if pas_courant < 0.03:
                pas_courant = pas_initial
    
    print(f"  Fin recherche locale: couverture finale = {meilleure_valeur}")
    return meilleure_solution, meilleure_valeur


# ALGORITHME TIOA PRINCIPAL - STRATÉGIE MULTI-REDÉMARRAGE


def algorithme_optimisation_tortue(redemarrages=20):
    """
    Implémente l'algorithme TIOA complet avec stratégie multi-redémarrage.
    """
    # Initialisation des meilleures solutions
    meilleure_solution_globale = None
    meilleure_valeur_globale = -1
    
    print("ALGORITHME TIOA - OPTIMISATION DE CONSTELLATION SATELLITE")
    print(f"Objectif: Couvrir {len(POINTS_TERRESTRES)} points avec {NOMBRE_SATELLITES} satellites")
    print(f"Strategie: {redemarrages} redemarrages avec recherche locale adaptative")
    print()
    
    # Boucle de multi-redémarrage
    for restart in range(redemarrages):
        print(f"Redemarrage {restart + 1:02d}/{redemarrages}:")
        
        # Génération d'une solution initiale aléatoire
        solution_initiale = np.concatenate([
            np.random.uniform(0, 2 * np.pi, NOMBRE_SATELLITES),
            np.random.uniform(-0.8, 0.8, NOMBRE_SATELLITES)
        ])
        
        # Recherche locale à partir de cette initialisation
        meilleure_solution_locale, meilleure_valeur_locale = recherche_locale_tortue(solution_initiale)
        
        # Mise à jour de la meilleure solution globale
        if meilleure_valeur_locale > meilleure_valeur_globale:
            meilleure_solution_globale = meilleure_solution_locale.copy()
            meilleure_valeur_globale = meilleure_valeur_locale
            print(f"  >>> NOUVELLE SOLUTION GLOBALE: couverture = {meilleure_valeur_globale}")
        else:
            print(f"  Solution locale: couverture = {meilleure_valeur_locale}")
        
        print()
        
        # Critère d'arrêt anticipé si couverture parfaite
        if meilleure_valeur_globale == len(POINTS_TERRESTRES):
            print("*** COUVERTURE PARFAITE ATTEINTE - Arret anticipe ***")
            break
    
    return meilleure_solution_globale, meilleure_valeur_globale


# FONCTIONS D'AFFICHAGE ET D'ANALYSE


def analyser_visibilite(positions_satellites):
    """
    Analyse détaillée de la visibilité pour chaque point terrestre.
    """
    print("ANALYSE DE VISIBILITE:")
    
    for i, point in enumerate(POINTS_TERRESTRES):
        satellites_visibles = []
        for j, sat in enumerate(positions_satellites):
            if est_visible_depuis_satellite(sat, point):
                satellites_visibles.append(j + 1)
        
        if satellites_visibles:
            print(f"Point {i + 1}: visible par le(s) satellite(s) {satellites_visibles}")
        else:
            print(f"Point {i + 1}: NON COUVERT")

def afficher_resultats(solution_optimale, valeur_optimale):
    """
    Affiche les résultats détaillés de l'optimisation.
    """
    print("RESULTATS FINAUX DE L'OPTIMISATION")

    
    # Conversion en positions réelles
    positions_finales = conversion_solution_vers_positions(solution_optimale)
    
    # Affichage des positions satellites
    print("\nPOSITIONS OPTIMALES DES SATELLITES:")
    
    for i in range(NOMBRE_SATELLITES):
        longitude_rad = solution_optimale[i]
        latitude_rad = solution_optimale[i + NOMBRE_SATELLITES]
        longitude_deg = np.degrees(longitude_rad)
        latitude_deg = np.degrees(latitude_rad)
        x, y, z = positions_finales[i]
        
        print(f"Satellite {i + 1}:")
        print(f"  Longitude: {longitude_rad:.3f} rad ({longitude_deg:.1f} deg)")
        print(f"  Latitude:  {latitude_rad:.3f} rad ({latitude_deg:.1f} deg)")
        print(f"  Position:  ({x:.1f}, {y:.1f}, {z:.1f}) km")
        print()
    
    # Résumé de couverture
    print(f"COUVERTURE FINALE: {valeur_optimale} / {len(POINTS_TERRESTRES)} points")
    
    # Analyse détaillée
    analyser_visibilite(positions_finales)


if __name__ == "__main__":
    """
    Point d'entrée principal du programme.
    """
    # Lancement de l'optimisation TIOA
    solution_optimale, valeur_optimale = algorithme_optimisation_tortue(redemarrages=20)
    
    # Affichage des résultats
    afficher_resultats(solution_optimale, valeur_optimale)
    
    print("OPTIMISATION TERMINEE")