# Trajectoire d'une particule chargée dans des champs E & M, à l'aide d'équations vectorielles
from pylab import *
from scipy import integrate

m = 25.0              # Masse et
q = 5.0               # Charge de la particule
E = array([0, 0, .1]) # champ électrique, composantes Ex,Ey & Ez
B = array([0, 0,  5]) # champ magnétique

def solver(X, t0):             # X est un tableau à 6 éléments, t0 est nécessaire pour le solveur
v = array([X[3], X[4], X[5]])  # on fait le vecteur vitesse
a = q * (E + cross(v,B)) / m   # F = q v x B ; a = F/m
return [X[3], X[4], X[5], a[0], a[1], a[2]]

tmax = 50  # calcul sur 50 secondes au plus
x0 = [0, 0, 0, 0, 1, 0]      # position et vitesse à t = 0
t = arange(0, tmax, 0.01)    # tableau de coordonnées temporelles
pv = integrate.odeint(solver, x0, t)  # intégration pour la position et la vitesse

savetxt('xyz.txt',pv)                 # on enregistre 6 colonnes dans un fichier :  x,y,x, Vx, Vy, Vz

from mpl_toolkits.mplot3d import Axes3D
ax = Axes3D(figure())
ax.plot(pv[:,0], pv[:,1], pv[:,2])    # graphique 3d plot de x, y et z
ax.set_zlabel('Z axis')
show()