earthTest.py
# -------------------------------------------------------------------------
#
# PYTHON for PHYSICS DUMMIES
# Problème 1
#
# Script de test
# Vincent Legat
#
# -------------------------------------------------------------------------
#
import numpy as np
# ============================================================
# FONCTIONS A MODIFIER [begin]
#
# -1- earthData() : paramètres physiques :-)
#
# G : constante universelle de la gravité en [UA^3/kg day^2]
# M : masse du soleil en [kg]
# m : masse de la Terre en [kg]
# v : vitesse de la Terre autour du soleil en [UA/day]
#
def earthData():
G = 1.0
M = 1.0
m = 1.0
v = 1.0
T = 2*np.pi
#
# Modifier les cinq valeurs précédentes :-)
#
return G,M,m,v,T
#
# -2- earthGravity(x,data) : champs gravitationnel
#
# g : vecteur numpy de taille deux contenant le champs 2D
#
def earthGravity(x,data):
G,M,m,v,T = data()
g = np.zeros(2)
#
# A complèter :-)
#
return g
#
# -3- earthEuler(n,g,dat) : intégration de la trajectoire
# par la méthode d'Euler explicite
# Il est possible de faire mieux : on le verra plus tard :-)
# Voir cours LEPL1104
#
# X : tableau numpy de taille (n+1)x2 contenant
# contenant les coordonnées discrètes de la trajectoire
#
def earthEuler(n,g,data):
G,M,m,v,T = data()
X = np.zeros((n+1,2)); X[0,:] = [1.0,0.0]
#
# A complèter :-)
#
return X
#
# FONCTIONS A MODIFIER [end]
# ============================================================
#
# -4- Le programme main :-)
# Calcul de constantes physiques via earthData
# Intégration de la trajectoire via earthGravity
# Et un joli plot !
#
def main():
G,M,m,v,T = earthData()
n = 10000;
print(" === Gravitational constant = %9.2e [UA^3/kg day^2]" % G)
print(" === Mass of the Sun = %9.2e [kg]" % M)
print(" === Mass of the Earth = %9.2e [kg]" % m)
print(" === Speed of the Earth = %9.2e [UA/day]" % v)
print(" === Period of the Earth = %4.1f [day]" % T)
print(" === Number of time steps = %d " % n)
X = earthEuler(n,earthGravity,earthData)
print(" === Final position of the Earth after one revolution = [%4.1f,%4.1f] [UA]" % (X[-1,0],X[-1,1]))
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['toolbar'] = 'None'
plt.figure("La Terre tourne autour du Soleil :-)")
plt.plot([0],[0],'ro',markerSize=20)
plt.plot(X[:,0],X[:,1],'-b')
plt.plot(X[-1,0],X[-1,1],'bo',markerSize=5)
plt.xlim(-1.5, 1.5)
plt.ylim(-1.5, 1.5)
plt.gca().set_aspect('equal',adjustable='box')
plt.axis("off")
plt.show()
main()