1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
from pylab import *

m = 2.0              # Masse de la particule
q = 5.0              # Charge
t = 0.0
dt = 0.001
N = 5000             # nombre de pas pour l'intégration

E = array([0.0, 0.0, 0.0])
B = array([0.0, 0.0, 5.0])

pos = zeros([N,3])   # tableau de N positions
pos[0] = [5,0,0]     # position initiale, x, y, et z
vel = zeros([N,3])
vel[0] = [20,0,0]    # vitesse initiale. vx, vy et vz

for k in range(N-1):
F = q * (E + cross(vel[k],B))        # F = q * [E (v x B)]
vel[k+1] = vel[k] + F/m * dt         # a = F/m; dv = a.dt
pos[k+1] = pos[k] + vel[k] * dt      # dx = v.dt
t = t + dt

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