delujoc prototip delca v zemljinem magnetnem polju

boris.py

@@ -5,9 +5,12 @@ from matplotlib import rc
 rc('font',**{'family':'serif','serif':['Computer Modern']})
 rc('text', usetex=True)
 
-Q = 1.602176565e-19
+e = 1.602176565e-19 #osnovni naboj [C]
+m_pr = 1.672621777e-27 #masa protona [kg]
+m_el = 9.10938291e-31 #masa elektrona [kg]
+c = 299792458 #hitrost svetlobe [m/s]
 
-def boris(x0, v0, E, B, dt, q=1, m=1):
+def boris(x0, v0, E, B, dt, tdur, q=1, m=1):
     ''' Borisov algoritem zabjega skoka za
     integracijo diferencialne enačbe
 
@@ -17,7 +20,7 @@ def boris(x0, v0, E, B, dt, q=1, m=1):
     B : funkcija, ki daje jakost polja B(x)
     dt: casovni korak '''
 
-    duration = 1000
+    duration = int(tdur/dt)
 
     X = np.zeros((duration,3))
     V = np.zeros((duration,3))
@@ -40,17 +43,49 @@ def boris(x0, v0, E, B, dt, q=1, m=1):
 
     return X, V
 
-E = lambda x: np.array([0.0,0.0,0.])
-B = lambda x: np.array([0.,0.0,1.*x[0]])
+def B_zemlja(x):
+    ''' model zemljinega magnentnega polja aproksimiran z 
+    magnetnim dipolom '''
+    B0 = 3.07e-5 # [T]
+    Rz = 6378137 # radii zemlje [m]
 
-x0 = np.array([-1.,0.,0.])
-v0 = np.array([0.,1.,0.])
+    r = np.sqrt(np.sum(x**2))
+    k = -B0*Rz**3/r**5
+    return k*np.array([3*x[0]*x[2],3*x[1]*x[2],2*x[2]**2-x[0]**2-x[1]**2])
 
-dt = 1e-2
+def Wk(m,V):
+    return 0.5*m*np.sum(V**2,1)
 
-X, V = boris(x0,v0,E,B,dt)
+def plot3(X,enote='m',zemlja=False):
+    fig = plt.figure()
+    ax = fig.add_subplot(projection='3d')
+    if zemlja:
+        u = np.linspace(0, 2 * np.pi, 100)
+        v = np.linspace(0, np.pi, 100)
+        x = np.outer(np.cos(u), np.sin(v))
+        y = np.outer(np.sin(u), np.sin(v))
+        z = np.outer(np.ones(np.size(u)), np.cos(v))
+        ax.plot_surface(x,y,z,zorder=-4)
 
-ax = plt.figure().add_subplot(projection='3d')
-ax.plot(X[:,0],X[:,1],X[:,2])
-plt.xlabel(r'$x$ label')
-plt.show()
+    ax.plot(X[:,0],X[:,1],X[:,2],zorder=4)
+    ax.set_xlabel(r'$x$ [{}]'.format(enote))
+    ax.set_ylabel(r'$y$ [{}]'.format(enote))
+    ax.set_zlabel(r'$z$ [{}]'.format(enote))
+
+    plt.show()
+
+
+if __name__ == "__main__": # testni del kode
+    E = lambda x: np.array([0.0,0.2,0.3])
+    B = lambda x: np.array([0.,0.0,1.*x[0]])
+
+    x0 = np.array([-1.,0.,0.])
+    v0 = np.array([0.,1.,0.])
+
+    dt = 1e-2
+
+    X, V = boris(x0,v0,E,B,dt,20)
+
+    #ax.quiver(X[0:-1:20,0],X[0:-1:20,1],X[0:-1:20,2],V[0:-1:20,0],V[0:-1:20,1],V[0:-1:20,2], color='black',normalize=True,length=.5)
+    # quiver res naredi grdo sliko
+    plot3(X)

zemlja.py

@@ -0,0 +1,32 @@
+#!/usr/bin/python3
+
+### koda za generiranje rezultatov za seminar 
+
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib import rc
+rc('font',**{'family':'serif','serif':['Computer Modern']})
+rc('text', usetex=True)
+from boris import *
+
+Rz = 6378137 # radii zemlje [m]
+
+q = e 
+m = m_pr
+
+K = 1e7 # kineticna energija 10MeV
+K = K*e # v Jouleih
+
+get_v0 = lambda K: c*np.sqrt(1-1/(K/(m*c**2)+1)**2)
+
+x0 = np.array([4*Rz,0.,0.])
+vpad = 30
+v0 = get_v0(K)*np.array([0., np.sin(vpad*np.pi/180), np.cos(vpad*np.pi/180)])
+
+E = lambda x: np.array([0.,0.,0.]) # elektricno polje je nula
+dt = 1e-3
+tsim = 20
+
+X, V = boris(x0, v0, E, B_zemlja, dt, tsim, q=q,m=m)
+plot3(X/Rz,enote="Rz", zemlja=True)
+