<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 La physique de l'ascenseur spatial
Nous allons d'abord charger les librairies Python de Mathematiques (numpy) et de représentation graphique (matplotlib). 
Sans numpy nous n'aurions pas de fonction racine (sqrt). 
Quant à matplotlib, il s'agit d'une librairie qui fonctionne un peu comme la calculatrice du Lyçée et nous permet de tracer des courbes.
Step1:  Définissons ensuite les constantes utilisées
Step2:  1. Movement orbital
Nous devons calculer la vitesse de rotation de la terre et la hauteur de l'orbite géostationnaire.
Nous en profiterons aussi pour calculer l'accélération à la surface de la terre.
1.1. Calcul de la vitesse de rotation de la terre
Soit Wt la vitesse angulaire de rotation de la terre, en radians par seconde
Step3:  1.2 Fonctions accélération
Nous allons définir des fonctions pour l'accélération gravitationnelle et l'accélération centrifuge.
Step4:  Exemple
Step5:  Exemple
Step6:  1.3 Calcul de la hauteur de l'orbite geostationnaire
Nous pouvons calculer la valeur exacte de l'orbite geostationnaire. Cette valeur doit
Step7:  Représentation graphique
Step8:  Nous observons que la Force totale s'annule pour le rayon de l'orbite géostationnaire. Cette courbe verte nous montre la force par unité de longueur et par unité de masse exercée sur chaque point du cable. La force est négative (attractive) au desssous de l'orbite géostationnaire et elle est positive au dessus.
2. Le cable
Notre but dans cette deuxième partie est de 
2.1 L'idée d'equilibre
Step9:  2.2 La modélisation du cable
Nous voulons comprendre pourquoi le cable est capable de se maintenir verticalement sans exercer de force sur la surface terrestre. 
Intuitivement, il semble que le cable sera en équilibre lorsque les deux parties exerceront la même force, l'une positive et l'autre négative, sur le point de rencontre, Rg. 
Le schéma suivant montre la courbe de l'accélération totale encadrée par des valeurs par défaut et par excès.
Step10:  Solution de l'intégrale
Step11:  Encadrer l'accélération
Step12:  Approximation de la Tension
Nous allons utiliser l'encadrement de l'accélération pour réaliser l'expérience de l'approximation de la Tension. 
Nous pouvons calculer les approximations par défaut et par excès de la Tension pour différentes valeurs du nombre de maillons.
Step13:  2.3 Calcul du contrepoids et de la longueur du cable
Step14:  Ci-dessous, en deux plots séparés
Step15:  3. Quels matériaux por le cable?
Nous pouvons maintenant calculer la résistance du cable pour différents matériaux pour savoir s'il existe un matériau capable de supporter la tension maximale. 
Nos allons multiplier la courbe de l'accélération de l'accélération par la densité pour obtenir une vraie Tension en Pascals.  
Valeurs de densité et de tension maximale por différents matériaux
Step16:  Calcul de résistance du cable
<ASSISTANT_TASK:>
Python Code:
import numpy as np
import matplotlib.pyplot as plt
#print(plt.style.available)
plt.style.use('presentation')
Explanation: La physique de l'ascenseur spatial
Nous allons d'abord charger les librairies Python de Mathematiques (numpy) et de représentation graphique (matplotlib). 
Sans numpy nous n'aurions pas de fonction racine (sqrt). 
Quant à matplotlib, il s'agit d'une librairie qui fonctionne un peu comme la calculatrice du Lyçée et nous permet de tracer des courbes.
End of explanation
G  = 6.67E-11    # Constante de gravitation universelle en m(3)*s(-2)*kg(-1)
Mt = 5.98E24     # Masse de la terre en kg  
Rt = 6378E3      # Rayon  de la terre en m
Explanation: Définissons ensuite les constantes utilisées:
End of explanation
Wt = (2*np.pi)/(24*60*60)
print("Vitesse de rotation de la terre = ",Wt,"Radians/seconde")
print("Vitesse de rotation de la terre = ",round (Wt,7),"Radians/seconde")
Explanation: 1. Movement orbital
Nous devons calculer la vitesse de rotation de la terre et la hauteur de l'orbite géostationnaire.
Nous en profiterons aussi pour calculer l'accélération à la surface de la terre.
1.1. Calcul de la vitesse de rotation de la terre
Soit Wt la vitesse angulaire de rotation de la terre, en radians par seconde:
End of explanation
def AccGravitationnelle (r):
    AccGravitationnelle = (G*Mt)/(r**2)
    return AccGravitationnelle
def AccCentrifuge (r):
    AccCentrifuge = r * (Wt**2)
    return AccCentrifuge
def AccTotale (r):
    AccTotale = AccGravitationnelle(r) - AccCentrifuge(r)
    return AccTotale
Explanation: 1.2 Fonctions accélération
Nous allons définir des fonctions pour l'accélération gravitationnelle et l'accélération centrifuge.
End of explanation
print("Accélération gravitationnelle à la surface de la terre:",AccGravitationnelle(Rt), "(m.s-2)")
print("Accélération gravitationnelle à la surface de la terre:",round(AccGravitationnelle(Rt),2), "(m.s-2)")
Explanation: Exemple: Calcul de l'accélération gravitationnelle à la surface de la terre:
End of explanation
print("Accélération centrifuge à la surface de la terre:",round(AccCentrifuge(Rt),9), "(m.s-2)")
Explanation: Exemple: Calcul de l'accélération centrifuge à la surface de la terre:
End of explanation
Rg =((G*Mt)/(Wt**2))**(1/3)
Hg = (Rg - Rt)
print("Rayon de l'orbite geostationnaire:", int (Rg/1000), 'km')
print("Hauteur de l'orbite geostationnaire:", int (Hg/1000), 'km')
Explanation: 1.3 Calcul de la hauteur de l'orbite geostationnaire
Nous pouvons calculer la valeur exacte de l'orbite geostationnaire. Cette valeur doit
End of explanation
r = np.arange(6E6, 160E6, 1E3)
r_km = r*1E-3
plt.figure(figsize=(16, 12), dpi= 100)
plt.plot(r_km,  AccGravitationnelle(r), color="cyan",   label='Accélération Gravitationnelle')
plt.plot(r_km,  AccCentrifuge(r),       color="orange", label='Accélération Centrifuge' )
plt.plot(r_km,  AccTotale(r),           color="green",  label='Accélération Totale' )
plt.axvline(x=Rt*1E-3, linestyle=":", color="brown")
plt.axvline(x=Rg*1E-3, linestyle=":", color="brown")
plt.xlabel("Rayon de l'orbite en km", fontsize=18)
plt.ylabel("Accélération en m.s-2", fontsize=18)
plt.axis([0, 160E3, -1, 11])
plt.grid()
plt.title("Accélérations sur le fil en fonction du rayon de l'orbite", fontsize=22)
plt.legend(loc='upper right', fontsize=16, fancybox='true')
plt.savefig('images/Accelerations.png')
plt.show()
Explanation: Représentation graphique:
Nous allons tracer les courbes des deux accélérations, gravitationelle et centrifuge.
Elles devraient se croiser dans un point d'abcisse r = 42241 km.
Le rayon r prendra des valeurs entre 6.000 km et 160.000 km
End of explanation
Rf = 150000E3
r = np.arange(6E6, 160E6, 1E3)
r_km = r*1E-3
plt.figure(figsize=(16, 12), dpi= 80)
plt.plot(r_km,  AccTotale(r), color="green", label='Accélération Totale' )
#Lignes verticales délimitent les trois rayons:
plt.axvline(x=Rt*1E-3, linestyle=":", color="brown")
plt.axvline(x=Rg*1E-3, linestyle=":", color="brown")
plt.axvline(x=Rf*1E-3, linestyle=":", color="brown")
plt.fill_between(r_km, AccTotale(r), 0, where=(r_km > Rt*1E-3)&(r_km < Rg*1E-3), color='orange', alpha=.4)
plt.fill_between(r_km, AccTotale(r), 0, where=(r_km > Rg*1E-3)&(r_km < Rf*1E-3), color="green", alpha=.4)
plt.xlabel("Rayon de l'orbite en km")
plt.ylabel("Accélération en m.s-2")
plt.axis([0, 160E3, -1, 11])
plt.grid()
plt.title("Accélération sur le fil en fonction du rayon de l'orbite")
plt.legend(loc='upper right')
plt.savefig('images/Equilibre.png')
plt.show()
Explanation: Nous observons que la Force totale s'annule pour le rayon de l'orbite géostationnaire. Cette courbe verte nous montre la force par unité de longueur et par unité de masse exercée sur chaque point du cable. La force est négative (attractive) au desssous de l'orbite géostationnaire et elle est positive au dessus.
2. Le cable
Notre but dans cette deuxième partie est de 
2.1 L'idée d'equilibre:
Nous voulons comprendre pourquoi le cable est capable de se maintenir verticalement sans exercer de force sur la surface terrestre. 
Intuitivement, il semble que le cable sera en équilibre lorsque les deux parties exerceront la même force. 
La valeur exacte de Rf est pour l'instant inconnue. La documentation indique une valeur de 150.000 km que nous utiliserons pour tracer les courbes. Nous vérifierons cette valeur plus tard.
End of explanation
NombreMaillons = 16
LongueurMaillon = 1E-3*(Rg-Rt)/NombreMaillons
print("Nombre de maillons:", NombreMaillons, "; Longueur d'un maillon:",round(LongueurMaillon,0), "km")
fig, ax = plt.subplots(figsize=(16, 8))
ax.set_xticks(np.arange(Rt*1E-3, Rg*1E-3, LongueurMaillon))
r = np.arange(Rt, Rg, 1E3)
r_km = r*1E-3
plt.plot(r_km, AccTotale(r), color="blue", label="Accélération")
r = np.arange(Rt, Rg, LongueurMaillon*1E3)
r_km = r*1E-3
plt.step(r_km, AccTotale(r), where='pre', color="orange", label="Approximation par défault")
plt.step(r_km, AccTotale(r), where='post', color="green", label="Approximation par excés")
plt.axvline(x=Rt*1E-3, linestyle=":", color="brown")
plt.axvline(x=Rg*1E-3, linestyle=":", color="brown")
ax.grid(markevery=LongueurMaillon)
plt.ylim(-1, 25)
plt.legend(loc='upper right')
plt.title("Expérience sur l'approximation de l'intégrale de l'accélération")
plt.savefig('images/Integrale.png')
plt.show()
Explanation: 2.2 La modélisation du cable
Nous voulons comprendre pourquoi le cable est capable de se maintenir verticalement sans exercer de force sur la surface terrestre. 
Intuitivement, il semble que le cable sera en équilibre lorsque les deux parties exerceront la même force, l'une positive et l'autre négative, sur le point de rencontre, Rg. 
Le schéma suivant montre la courbe de l'accélération totale encadrée par des valeurs par défaut et par excès.
End of explanation
def IntegraleAcc (r):
    IntegraleAcc = G*Mt*( -1/r - (r**2)/(2*Rg**3) + 1/Rt + (Rt**2)/(2*Rg**3))
    return IntegraleAcc
Explanation: Solution de l'intégrale:
Nous allons introduire une autre fonction que nous appelerons IntegraleAcc. 
Cette fonction est une primitive de la fonction AccTotale, sa dérivée en r est égale à AccTotale. La valeur de cette primitive pour r=Rt est nulle.
End of explanation
NombreMaillons = 16
print(Rt)
LongueurMaillon = (Rg-Rt)/NombreMaillons
print("Nombre de maillons:", NombreMaillons, ",Longueur d'un maillon:",round(LongueurMaillon*1E-3,0), "km")
# Maillons est la liste 
Maillons_pre  = np.linspace(Rt, Rg,NombreMaillons+1)
Maillons      = Maillons_pre + LongueurMaillon/2
Maillons_post = Maillons_pre + LongueurMaillon
Acc_pre  = AccTotale(Maillons_pre)
Acc      = AccTotale(Maillons)
Acc_post = AccTotale(Maillons_post)
# Obtenir les valeurs maximum et minimum de l'accélération
AccMax = np.maximum (Acc_pre,Acc_post)
AccMin = np.minimum (Acc_pre,Acc_post)
# Représentation graphique:
r = np.arange(Rt, Rg, 1E3)
r_km = r*1E-3
plt.figure(figsize=(16, 10), dpi= 80, facecolor='w', edgecolor='k')
plt.plot(Maillons, AccMax, color='red', label='Approximation par excès')
plt.plot(r, AccTotale(r), color='blue', label="Acceleration Totale")
plt.plot(Maillons, AccMin, color='green', label="Approximation par excès")
plt.plot(Maillons, AccMax,'ro')
plt.plot(Maillons, Acc,   'bo')
plt.plot(Maillons, AccMin,'go')
plt.xticks(Maillons)
plt.ylim(-1, 25)
plt.grid()
plt.axvline(x=Rt, linestyle=":", color="brown")
plt.axvline(x=Rg, linestyle=":", color="brown")
plt.legend()
plt.title("Encadrer l'accélération")
plt.savefig('images/AccMaxMin.png')
plt.show()
Explanation: Encadrer l'accélération
End of explanation
NombreMaillons = 16
print(Rt)
LongueurMaillon = (Rg-Rt)/NombreMaillons
print("Nombre de maillons:", NombreMaillons, ",Longueur d'un maillon:",round(LongueurMaillon*1E-3,0), "km")
# Maillons est la liste 
Maillons_pre  = np.linspace(Rt, Rg,NombreMaillons+1)
Maillons      = Maillons_pre + LongueurMaillon/2
Maillons_post = Maillons_pre + LongueurMaillon
Acc_pre  = AccTotale(Maillons_pre)
Acc      = AccTotale(Maillons)
Acc_post = AccTotale(Maillons_post)
# Obtenir les valeurs maximum et minimum de l'accélération
AccMax = np.maximum (Acc_pre,Acc_post)
AccMin = np.minimum (Acc_pre,Acc_post)
#Surface de chaque maillon.
Aire_max     = LongueurMaillon*AccMax
Aire_min     = LongueurMaillon*AccMin
Int_min      = [0]*(NombreMaillons+1)
Int_max      = [0]*(NombreMaillons+1)
#Calculer les approximations
for i in range(NombreMaillons+1):
    Int_min[i]  = sum(Aire_min[0:i] ) 
    Int_max[i]  = sum(Aire_max[0:i] )
# Représentation graphique:
r = np.arange(Rt, Rg, 1E3)
r_km = r*1E-3
plt.figure(figsize=(16, 10), dpi= 80, facecolor='w', edgecolor='k')
plt.plot(r, IntegraleAcc(r), label="Fonction primitive de l'accelération Totale")
plt.plot(Maillons, Int_min,   label='Intégrale, approximation min')
plt.plot(Maillons, Int_max,   label='Intégrale, approximation max')
plt.xticks(Maillons)
plt.axvline(x=Rt, linestyle=":", color="brown")
plt.axvline(x=Rg, linestyle=":", color="brown")
plt.grid()
plt.legend(loc='lower center')
plt.title("Expérience sur l'approximation de la l'intégrale de l'accélération")
plt.savefig('images/Tension-16.png')
plt.show()
Explanation: Approximation de la Tension
Nous allons utiliser l'encadrement de l'accélération pour réaliser l'expérience de l'approximation de la Tension. 
Nous pouvons calculer les approximations par défaut et par excès de la Tension pour différentes valeurs du nombre de maillons.
End of explanation
# Rf: rayon du bout final du cable, tension nulle
# Valeur théorique trouvée dans la documentation, "Physics of the Space Elevator"
Rf = (Rt/2)*( np.sqrt(1+8*(Rg/Rt)**3) - 1)
LongueurCable = Rf-Rt
print ("Rayon final: Rf=",int(Rf*1E-3),"km")
print ("Longueur totale du cable, LongueurCable =", int(LongueurCable*1E-3),"km")
Explanation: 2.3 Calcul du contrepoids et de la longueur du cable:
Nous allons calculer la valeur exacte de Rf avec la formule que nous avons trouvée dans l'article 'Physics of the space elevator'. Cette valeur doit correspondre sur le graphique ci-dessus au deuxième point pour lequel la Tension s'annule.
End of explanation
r = np.arange(6E6, 160E6, 1E3)
r_km = r*1E-3
plt.figure(figsize=(16, 12), dpi= 100, facecolor='w', edgecolor='k')
#Première figure: 2 rangées, 1 colonne, figure 1
plt.subplot(211)
plt.title("Accélération Totale et son intégrale en fonction du rayon")
plt.axvline(x=Rt*1E-3, linestyle=":", color="brown")
plt.text(Rt*1E-3, -3.2, 'Rt', fontsize ='16', color="brown")
plt.axvline(x=Rg*1E-3, linestyle=":", color="brown")
plt.text(Rg*1E-3, -3.2, 'Rg', fontsize ='16', color="brown")
plt.axvline(x=Rf*1E-3, linestyle=":", color="brown")
plt.text(Rf*1E-3, -3.2, 'Rf', fontsize ='16', color="brown")
plt.ylabel("Accélération Totale (m.s-2)")
plt.plot(r_km,  AccTotale(r), color="green",  linewidth='1.5', label="Accélération Totale")
plt.legend()
plt.grid()
#Deuxième figure: 2 rangées, 1 colonne, figure 2
plt.subplot(212)
plt.axvline(x=Rt*1E-3, linestyle=":", color="brown")
plt.axvline(x=Rg*1E-3, linestyle=":", color="brown")
plt.axvline(x=Rf*1E-3, linestyle=":", color="brown")
plt.plot(r_km, IntegraleAcc(r), label="Integrale de l'accélération")
plt.xlabel("Rayon de l'orbite en km")
plt.ylabel("Intégrale de l'accélération Totale (m2.s-2)")
plt.legend()
plt.grid()
plt.savefig('images/acctotale-integrale.png')
plt.show()
Explanation: Ci-dessous, en deux plots séparés:
End of explanation
densite_acier  = 7900  # kg/m3
densite_kevlar = 1440  # kg/m3
densite_cnt    = 1300  # kg/m3
Tmax_acier  = 5E9      # Pascals (N/m2)
Tmax_kevlar = 3.6E9    # Pascals (N/m2)
Tmax_cnt    = 130E9    # Pascals (N/m2)
Explanation: 3. Quels matériaux por le cable?
Nous pouvons maintenant calculer la résistance du cable pour différents matériaux pour savoir s'il existe un matériau capable de supporter la tension maximale. 
Nos allons multiplier la courbe de l'accélération de l'accélération par la densité pour obtenir une vraie Tension en Pascals.  
Valeurs de densité et de tension maximale por différents matériaux:
Nous avons extrait ces valeurs de la documentation 'Physics of the Space Elevator'
End of explanation
r = np.arange(Rt, Rg, 1E3)
r_km = r*1E-3
plt.figure(figsize=(16, 12), dpi= 80, facecolor='w', edgecolor='k')
plt.axvline(x=Rg*1E-3, linestyle=":", color="brown")
#Acier:
plt.plot(r_km, IntegraleAcc(r)*densite_acier, color="red", label="Tension Acier")
plt.axhline(y=Tmax_acier, linestyle=":", color="red", linewidth=3, label="Résistance Acier" )
#Kevlar:
plt.plot(r_km, IntegraleAcc(r)*densite_kevlar, color="orange", label="Tension Kevlar")
plt.axhline(y=Tmax_kevlar, linestyle=":", color="orange",linewidth=3, label="Résistance Kevlar")
#Carbon nanotubes:
plt.plot(r_km, IntegraleAcc(r)*densite_cnt, color="green", label="Tension CNT")
plt.axhline(y=Tmax_cnt, linestyle=":", color="green",linewidth=3, Label="Résistance CNT")
plt.xlabel("Rayon de l'orbite en km")
plt.ylabel("Tension en Pascals (Pa), (N.m-2)")
plt.grid()
plt.title("Tension et résistance du cable")
plt.legend(loc="upper left")
plt.savefig('images/materiaux.png')
plt.show()
print("Tension maximum du cable au point Rg:" )
print("Tension Maximum par unité de masse:", int(IntegraleAcc(Rg))*1E-9,"Giga Pascals/kg")
print("Tension maximum Rg, Acier:",  int(IntegraleAcc(Rg)*densite_acier*1E-9),  "Giga Pascals")
print("Tension maximum Rg, Kevlar:", int(IntegraleAcc(Rg)*densite_kevlar*1E-9), "Giga Pascals")
print("Tension maximum Rg, CNT:",    int(IntegraleAcc(Rg)*densite_cnt*1E-9),    "Giga Pascals")
Explanation: Calcul de résistance du cable:
Les calculs ci-dessous tracent les courbes de Tension entre Rt et Rg et les valeurs maximales de Tension (Tensile strength en Anglais) que peuvent supporter les trois matériaux:
End of explanation
<END_TASK>