# -*- coding: utf-8 -*-
# Résolution numérique d'un problème à un corps soumis à un
# potentiel newtonien (en 1/r^2)
# d'équa. diff. vec{r}'' = -k/m vec{r}/r^3
# le vecteur de données considéré sera le vecteur (x, y , vx, vy)
# dont la dérivée est (vx, vy, -k/m*x/(x^2+y^2)^3/2, -k/m*y/(x^2+y^2)^3/2)
from numpy import *
from pylab import *
from PyQt5.QtGui import *
from PyQt5.QtWidgets import *
from PyQt5.QtCore import *
from point import Point
import os, tempfile, time
from matplotlib_widget import MyMplCanvas
import flottant as flottant
from rectangle_sensible import Rs
from debug import Debug
def rk4(derivs, y0, t):
"""
C'est le code de rk4 pris dans le module matplotplib.
Liste des paramètres d'entrée
derivs :
une fonction qui accepte en entrée un 6-uplet position,vitesse
et le papramètre temps, et qui renvoie en sortie un 6-uplet de dérivées.
y0 :
un 6-uplet représentant la position et la vitesse initiales
t :
une liste de dates régulièrement espacées pour lesquelles on veut
construire les points et vitessesde la trajectoire
Résultat de la fonction :
une liste de 6-uplets représentant les positions et vitesses aux
instants de la liste des dates données.
"""
Float=0.0
try: Ny = len(y0)
except TypeError:
yout = zeros( (len(t),), float)
else:
yout = zeros( (len(t), Ny), float)
yout[0] = y0
i = 0
for i in arange(len(t)-1):
thist = t[i]
dt = t[i+1] - thist
dt2 = dt/2.0
y0 = yout[i]
k1 = asarray(derivs(y0, thist))
k2 = asarray(derivs(y0 + dt2*k1, thist+dt2))
k3 = asarray(derivs(y0 + dt2*k2, thist+dt2))
k4 = asarray(derivs(y0 + dt*k3, thist+dt))
yout[i+1] = y0 + dt/6.0*(k1 + 2*k2 + 2*k3 + k4)
return yout
def euler(derivs, y0, t):
"""
C'est le code de la méthode d'Euler, qui est d'ordre 1 et très simple.
Liste des paramètres d'entrée
derivs :
une fonction qui accepte en entrée un 6-uplet position,vitesse
et le papramètre temps, et qui renvoie en sortie un 6-uplet de dérivées.
y0 :
un 6-uplet représentant la position et la vitesse initiales
t :
une liste de dates régulièrement espacées pour lesquelles on veut
construire les points et vitessesde la trajectoire
Résultat de la fonction :
une liste de 6-uplets représentant les positions et vitesses aux
instants de la liste des dates données.
"""
try: Ny = len(y0)
except TypeError:
yout = zeros( (len(t),), Float)
else:
yout = zeros( (len(t), Ny), Float)
yout[0] = y0
i = 0
for i in arange(len(t)-1):
thist = t[i]
dt = t[i+1] - thist
y0 = yout[i]
k1 = asarray(derivs(y0, thist))
yout[i+1] = y0 + dt*k1
return yout
def nouvChampGrav(x0,m):
"""Cette fonction renvoie une fonction anonyme représentant un champ
gravitationnel créé par un objet de masse m immobile aux coordonnées
spécifiées par x0.
le profil de la fonction résultat est :
(vecteur position -> vecteur accélération)
"""
return lambda x: (-6.67e-11*m*(x[0]-x0[0])/((x[0]-x0[0])**2+(x[1]-x0[1])**2+(x[2]-x0[2])**2)**1.5,
-6.67e-11*m*(x[1]-x0[1])/((x[0]-x0[0])**2+(x[1]-x0[1])**2+(x[2]-x0[2])**2)**1.5,
-6.67e-11*m*(x[2]-x0[2])/((x[0]-x0[0])**2+(x[1]-x0[1])**2+(x[2]-x0[2])**2)**(1.5))
class trajectoire(QObject):
def __init__(self,dt,t,y0, champ):
"""Les paramètres sont :
dt : intervalle de temps
t : durée totale de la simulation
y0 : vecteur à 6 composantes (position, vitesse initiales)
champ : une fonction donnant l'accélération à partir de la position
"""
QObject.__init__(self)
self.dt = dt
self.t = arange(0,t,self.dt)
self.y0 = y0
self.champ=champ
self.calcul()
def dessine(self):
x=[]
y=[]
for point in self.pv :
x.append(point[0])
y.append(point[1])
self.x = x
self.y = y
plot(self.x,self.y) # le dessin est en projection dans le plan x,y
axis('equal')
show()
def derivs(self,x,t):
"""calcul de la dérivée du vecteur à 6 composantes position,vitesse
"""
g = self.champ(x[0:3]) # le champ de gravité
return (x[3], # vitesse x
x[4], # vitesse y
x[5], # vitesse z
g[0], # acceleration x
g[1], # acceleration y
g[2]) # acceleration z
def calcul(self):
"""lance le calcul de la trajectoire à l'aide de l'algorithme de
Runge-Kutta qui est d'ordre 4 et rapide en même temps.
le résultat est dans self.pv, qui est une liste contenant
les 6-uplets position, vitesse.
"""
self.pv = rk4(self.derivs, self.y0, self.t)
def calculeNorme(nuplet):
s=0
for x in nuplet:
s+=x**2
return s**0.5
def projection(n):
if n >=0:
return lambda nuplet:nuplet[n]
else:
return 0.0
class Trajectoire(Rs):
def __init__(self, parent, rayon_astre, mainWin, debugger=Debug(0)):
Rs.__init__(self, parent, debugger=debugger)
self.parent = parent
self.mainWin=mainWin
self.milieuX=self.width()/2
self.milieuY=self.height()/2
self.rep= mainWin.rep
#self.setFrameShape(QFrame.Box)
self.setEchelle(6.4e6/25)
self.points={}
self.planetes={}
self.vitesse=[]
self.date=[]
self.dt=1
self.widget_vit_norm=None
self.widget_vitx=None
self.widget_vity=None
self.boum=-1.0
self.traj=None
def setEchelle(self,val, mode="mppx"):
""" Régle l'échelle, selon le mode choisi.
mode=mppx : échelle en mètre par pixel
mode=max : l'échelle sera ajustée pour que le point (0,val) soit dans
la fenêtre de trajectoire (à 95% du maximum)
"""
if mode=="mppx":
self.echelle=val
self.debug(9,"Échelle %s px/m (mode direct)" %self.echelle)
elif mode=="max":
self.echelle=val/self.milieuY/0.95
self.debug(9,"%s m pour %s px" %(val,self.milieuY))
self.debug(9,"Échelle %s px/m (mode max)" %self.echelle)
else:
self.debug(0,"Le mode %s ne convient pas pour setEchelle()" %mode)
def lance(self):
self.calcul_parametre()
if self.mainWin.ui.checkBox_efface.isChecked() :
self.efface()
self.debug(10,"efface la trajectoire depuis un lancé")
a=time.time()
self.traj= trajectoire(self.dt, self.t, self.pos+self.vit, self.gAstre)
self.debug(8, "%s, %s, %s, %s" %(self.dt, self.t, self.pos+self.vit, self.gAstre))
self.debug(5,"calcul en %s secondes" %(time.time()-a))
self.dessine_trajectoire()
#dessine les vitesses
self.grapheV()
def grapheV(self):
if self.widget_vit_norm != None:
self.widget_vit_norm.hide()
if self.widget_vitx != None:
self.widget_vitx.hide()
if self.widget_vity != None:
self.widget_vity.hide()
self.widget_vit_norme = MyMplCanvas(self.mainWin.ui.label_vit_norm,self.vitesse, calculeNorme, self.date, width=3, height=5, dpi=30, cliquable=True, titre="Norme de la vitesse")
self.widget_vitx = MyMplCanvas(self.mainWin.ui.label_vitx,self.vitesse, projection(0), self.date, width=3, height=5, dpi=30, cliquable=True, titre="Abscisse de la vitesse")
self.widget_vity = MyMplCanvas(self.mainWin.ui.label_vity,self.vitesse, projection(1), self.date, width=3, height=5, dpi=30, cliquable=True, titre="Ordonnée de la vitesse")
self.widget_vit_norme.show()
self.widget_vitx.show()
self.widget_vity.show()
def calcul_parametre(self):
self.masse_astre = self.mainWin.getMasseAstre()
self.distance_astre = self.mainWin.getDistanceAstre()
self.vitesse_x,self.vitesse_y = self.mainWin.getVitesse()
self.G=6.67259e-11
#calcul des énergies massiques (cinétique, potentielle, mécanique)
Ec_massique=0.5*(self.vitesse_x**2+self.vitesse_y**2)
Ep_massique=-self.masse_astre*self.G/self.distance_astre
Em_massique=Ec_massique+Ep_massique
if Em_massique >= 0:
self.debug(1,"L'énergie mécanique est excessive (%s J/kg), la trajectoire ne se fermera pas" %Em_massique)
t=self.tr("Ce n'est pas un satellite")
q=self.tr("L'énergie mécanique initiale est positive, l'objet lancé échappera à l'attraction de l'astre. Voulez-vous tracer une partie de la trajectoire ?")
ret=QMessageBox.question (self, t, q)
if ret==QMessageBox.Yes:
self.t=30*(self.distance_astre**3/self.masse_astre/self.G)**0.5
self.dt=self.t/100
else:
raise(ValueError)
else:
# calcul du grand axe a : Em = -k/2a pour une trajectoire elliptique
# donc a = -k/2Em
a = - self.masse_astre*self.G/2/Em_massique
# calcul de la période, en utilisant la troisième loi de Kepler
# T²=4pi²/MG*a³
self.t = 2*pi*(a**3/self.masse_astre/self.G)**0.5
self.dt = self.t/1000
self.debug(5,"masse astre %s" %self.masse_astre)
self.gAstre = nouvChampGrav((0,0,0),self.masse_astre)
self.pos = (0.0, self.distance_astre, 0)
self.vit = (self.vitesse_x,self.vitesse_y,0.0)
def dessine_trajectoire(self):
self.debug(9,"dessine la trajectoire")
self.dessine(self.traj.pv)
# on teste si toute la trajectoire tien bien là.
r=self.maxDistance(self.traj.pv)
if not self.surementVisible(r) and not self.boum > 0:
t=self.tr("Changement d'échelle")
q=self.tr("Un dépassement a été détecté. Voulez-vous changer d'échelle ?")
ret=QMessageBox.question (self, t, q)
if ret==QMessageBox.Yes:
self.setEchelle(r,"max")
self.efface()
self.dessine(self.traj.pv)
def maxDistance(self,points_vitesses):
"""renvoie la distance max entre le centre et le satellite en
projection sur le plan xy.
"""
r=0
for pv in points_vitesses:
if fabs(pv[0])>r: r=fabs(pv[0])
if fabs(pv[1])>r: r=fabs(pv[1])
return r
def surementVisible(self,r):
"""vrai si un cercle de rayon r est visible à coup sûr"""
return r/self.echelle < self.milieuY
def efface(self):
for k in self.points.keys():
objet=self.points[k]
objet.hide()
objet.clear()
self.points = {}
self.vitesse = []
def trace_point(self, x, y, couleur, type="petit", tau=None):
if tau:
cle=(x,y,tau)
else:
cle=(x,y)
if cle in self.points.keys():
# efface des points préexistants de même emplacement
objet=self.points[cle]
objet.hide()
objet.clear()
self.points[cle]=Point(self, (x,y), couleur, "", self.mainWin,type_de_point=type)
self.points[cle].show()
def dessine(self,points_vitesses):
self.boum=-1.0
OKtau=False
intervalle=self.mainWin.ui.intervale.text()
self.tau = int(flottant.traduit(intervalle))
if self.tau > 60: OKtau=True
self.vitesse = []
self.date=[]
tau_entier = 0
for i in range(0,len(points_vitesses),10):
# on ne trace qu'un point sur 10, soit 100 points
# sur les 1000 calculés
pv=points_vitesses[i]
pix_x=int(pv[0]/self.echelle+self.milieuX)
pix_y=int(-pv[1]/self.echelle+self.milieuY)
# pv est un hexuplet : 3 coordonnées de position, 3 de vitesse
self.vitesse.append((pv[3],pv[4]))
self.date.append(self.dt*i)
if ((pix_x-self.milieuX)**2+(pix_y-self.milieuY)**2)**(0.5) < int(self.mainWin.getRayonAstre()/self.echelle) :
self.trace_point(pix_x, pix_y, "red", type="boum")
self.boum=self.dt*i
break
else:
self.trace_point(pix_x, pix_y, "red")
if OKtau:
try :
#dessine un point tous les "tau" secondes si défini.
tau = i*self.dt/self.tau
if int(tau) > int(tau_entier) :
self.trace_point(pix_x, pix_y, "blue", type="gros", tau=tau)
tau_entier=int(tau)
except AttributeError:
pass
self.repaint()
self.update()
def paintEvent(self, event):
painter = QPainter()
painter.begin(self)
painter.setBrush(Qt.CrossPattern)
painter.setPen(Qt.green)
rayon=int(self.mainWin.getRayonAstre()/self.echelle)
painter.drawLine(self.milieuX, 0, self.milieuX, self.milieuY*2)
painter.drawLine(0,self.milieuY , 2*self.milieuX, self.milieuY)
img=self.astreImg
sourcerect=QRect(0,0,512,512)
self.milieuX=self.size().width()/2
self.milieuY=self.size().height()/2
targetrect=QRect(self.milieuX-rayon,self.milieuY-rayon,2*rayon,2*rayon)
#TODO : il faudrait redessiner le rectangle en fonction de l'échelle
#TODO : la taille n'est pas bien calculée quand on redimensionne !!!
painter.drawEllipse(targetrect)
painter.drawImage(targetrect,img,sourcerect)
painter.end()
def dir(self,choix):
return self.rep.chemin(choix)
def debug(self,level,msg):
self.mainWin.debug(level,msg)
def getPlanete(self,nom):
if nom not in self.planetes.keys():
handle, imageFile = tempfile.mkstemp(".png")
os.close(handle)
cmd="xplanet -latitude 90 -num_times 1 -glare 10 -body %s -radius 50 -searchdir %s -transpng %s -rotate -70" %(nom,self.dir("textures"),imageFile)
# la rotation de -70° permet au méridiens français de se trouver
# au centre de l'image. Cette valeur est empirique et dépend
# probablement de l'implémentation de Xplanet. Il n'est pas
# évident de jouer avec l'option -north qui seule permet de
# contrôler totalement l'orientation de l'image.
self.debug(8,"Création de l'image de planète par \"%s\"" %cmd)
os.system(cmd)
self.planetes[nom]=QImage(imageFile)
os.system("rm -f %s" %imageFile)
return self.planetes[nom]
def choisi_astre(self,nom):
self.astreImg= self.getPlanete(nom)
class Satellite(QLabel):
def __init__(self, parent, coord_sat, image):
QLabel.__init__(self, parent)
self.parent = parent
self.coord_sat=coord_sat
self.image = image
def paintEvent(self, event):
painter = QPainter()
painter.setBrush(Qt.CrossPattern)
painter.setPen(Qt.green)
painter.begin(self)
x, y = self.coord_sat.x(), self.coord_sat.y()
self.debug(5,"%s,%s" %(x,y))
painter.drawEllipse(x,y,30,16)
painter.drawPixmap(self.coord_sat, self.image)
painter.end()