#!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as ode
from mpl_toolkits import mplot3d

q=1.e-9
m=1.e-9
qm=q/m
b=[0.,0.,0.1]
e=[0.,0.,0.]
betta=0.0
v0=10.

def func(y,t):
   yp=np.zeros_like(y)
   yp[0]=y[3]
   yp[1]=y[4]
   yp[2]=y[5]
   yp[3]=qm*(y[4]*b[2]-y[5]*b[1]+e[0])-betta*y[3]
   yp[4]=qm*(y[5]*b[0]-y[3]*b[2]+e[1])-betta*y[4]
   yp[5]=qm*(y[3]*b[1]-y[4]*b[0]+e[2])-betta*y[5]
   return yp

y0=[0.,0.,0.,0.,10.,20.]
time=np.arange(0,600,5)
res=ode.odeint(func,y0,time)
ax=qm*(res[:,4]*b[2]-res[:,5]*b[1]+e[0])-betta*res[:,3]
ay=qm*(res[:,5]*b[0]-res[:,3]*b[2]+e[1])-betta*res[:,4]
az=qm*(res[:,3]*b[1]-res[:,4]*b[0]+e[2])-betta*res[:,5]

fig=plt.figure()
p=fig.gca(projection="3d")
p.plot(res[:,0],res[:,1],res[:,2],'ro',color='black')
p.quiver(res[:,0],res[:,1],res[:,2],res[:,3],res[:,4],res[:,5],arrow_length_ratio=0.0001,length=5,color='red')
p.quiver(res[:,0],res[:,1],res[:,2],ax,ay,az,arrow_length_ratio=0.0001,length=40,color='green')

plt.grid()
plt.show()

print (res)