#JCourtin 06/2016
#D'après le BUP de Thierry ALHAHEL vol.99 Janvier 2005

#### LE VOILIER SOLAIRE 

from matplotlib import pyplot as plt
from math import *

#DEFINITION DES CONSTANTES:
c=3.0e8#m/s 
G=6.67e-11 #SI
MT=5.98e24#kg
RT=6.400e6#m
   
E=1.4e3#kW/m2
S=2.1022e6#m2
m=2000#kg
p=E/c  #pression radiation max
Fmax = 2*p*S
am=Fmax/m   # accélération de poussée max ~ 1 mg.

H=1 


################################ Partie à modifier #############################

##Fonctions dynamiques et cinématiques
#Calcul des forces
def Frad(r, phi, alpha, H):
    #???? force radiale
    return #force

def Ftan(r, phi, alpha, H):
    #???? force tangentielle
    return #force

#Calcul des accélérations
def dudt(r, phi, u, s, alpha, H):
    #???? accélération
    return #acc
    
def dsdt(r, phi, u, s, alpha, H):
    #???? accélération
    return #acc
    
#obtention de l'angle optimisé
def getAlpha(u,s):
    #ax = #???; ay = s#???
    """if (ax>=0):
        return asin(ay)        #Commenter cette partie
    else:
        return pi-asin(ay)""""
    return #effacer cette ligne
    
    
#energie mécanique
def getEnergie(r, phi, u, s):
    #E=#????
    return #E

##Conditions Initiales

R0   = 42250e3#m          #pourquoi cette valeur ?
phi0 = 3*pi/2 #rad

s0 = 3072.5 #m/s          #idem
u0 = 0.                   #idem

alpha0=getAlpha(u0, s0)
E0=getEnergie(R0, phi0, u0, s0)


## Intégration numérique avec Euler Explicite
t=[0.]; r=[R0]; phi=[phi0]; s=[s0]; u=[u0]; alpha=[alpha0]   #trajectoire

e=[E0]#énergie

Dt=1.#seconde
NH=150 #nb heure de simu
NInt=int(3600/Dt)  #Nb de calcul entre deux point tracés

N=int(NH*3600/Dt)

for n in range(N-1): #construit n+1 sur la base de n

    if (abs(r[n]*sin(phi[n]))<RT and cos(phi[n])>0.):  #expliquer
        H=0.
    else:
        H=1.
    """
    t   +=
    r   += 
    phi += 
    u   +=        #A compléter avec Euler
    s   += 
    alpha+=
    e   += 
    """
    
    
## FONCTIONS POUR TRACER LES GRAPHES ###########################################


##tracé de la trajectoire
def plotTrajectoire(r):
    index=[NInt*k for k in range(len(r)//NInt)] #sert à quoi ???
    #x=[compréhension de liste]
    #y=[compréhension de liste]
    
    plt.close()
    #plt.plot(x,y)
    #plt.plot(x,y,'o')
    #plt.show()
    
## Etude énergétique
def plotEnergie(t,e):
    index=[NInt*k for k in range(len(r)//NInt)]
    """
    t=[compréhension de liste]#temps en heure
    e=[compréhension de liste]                        #A compléter
    plt.close()
    plt.plot(t,e)
    plt.plot([0,200],[0,0])
    plt.show()
    """
    
def plotVlib(t,r,u,s):
    index=[NInt*k for k in range(len(r)//NInt)]
    """vLib=[compréhension de liste]
    t=[compréhension de liste]#temps en heure          #A compléter
    plt.close()
    plt.plot(t,vLib)
    plt.plot([0,200],[1,1])
    plt.show() 
    """
#### MAIN ################################# ne pas modifier ####################
if (__name__=="__main__"):
    
    
    #test1 : Trajectoire
    plotTrajectoire(r)
    
    #test2 : Energie
    #plotEnergie(t,e)
    
    #test 3 : vitesse de libération
    #plotVlib(t,r,u,s)