#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sat Jun  4 16:22:00 2022

@author: macbookraphael
"""
import numpy as np
import matplotlib.pyplot as plt
import csv
from math import sqrt, pi
from matplotlib.patches import Polygon

def surface(xA, yA, xB, yB):
    """ calcule la surface d'un triangle OAB (O est l'origine du repère) 
    xA, yA, xB, yB sont les coordonnées des sommets A et B """
    a = sqrt(xA**2 + yA**2)
    b = sqrt(xB**2 + yB**2)
    c = sqrt((xB - xA)**2 + (yB - yA)**2)
    p = (a + b + c) / 2 
    s = (1/4) * sqrt(p * (p-a) * (p-b) * (p-c)) # formule de Héron
    return s

def jpl_horizons(donnees):
    """ donnees est le fichier de données éphémérides du satellite 
    Utilise une source de données issue de JPL HORIZONS """
    with open(donnees) as source:
        donnees = list(csv.DictReader(source, delimiter=","))
    dates = [donnee["Calendar Date (TDB)"][6:-14] for donnee in donnees]
    dates_Julien = [float(donnee["JDTDB"].replace(',','.')) for donnee in donnees]
    dt = dates_Julien[1] - dates_Julien[0]
    t = np.array([dt*i for i in range(len(dates))])
    x = np.array([float(donnee["X"].replace(',','.')) for donnee in donnees])
    y = np.array([float(donnee["Y"].replace(',','.')) for donnee in donnees])
    return dates, t, x, y

def kepler_2(donnees, nom, pas):
    """ Reprend la démarche vue avec Halimède pour un satellites dont l'éphéméride est dans le fichier donnees """
    
    dates, t, x, y = jpl_horizons(donnees)

    fig, ax = plt.subplots()

    plt.scatter(0, 0)
    plt.plot(x, y, 'r.')
    plt.axis("equal")
    plt.title(f"Position de {nom} (en km) dans le référentiel de Neptune\n")

    for i in range(0, len(x), pas):
        xA, yA, xB, yB = x[i], y[i], x[i+1], y[i+1]
        triangle = Polygon([(0,0),(xA, yA),(xB, yB)], facecolor='0.8', edgecolor='0.5')
        ax.add_patch(triangle)

    plt.show()

    for i in range(0, len(x), pas):
        xA, yA, xB, yB = x[i], y[i], x[i+1], y[i+1]
        dateA, dateB = dates[i], dates[i+1]
        a, b = t[i], t[i+1]
        s = surface(xA, yA, xB, yB)
        print(f"Entre {dateA} et {dateA} (durée {b-a} j): aire balayée = {s:#.2g} km^2")

donnees = "Lune.csv"
nom = "nereid"
pas = 100
kepler_2(donnees, nom, pas)