diff --git a/__pycache__/boris.cpython-37.pyc b/__pycache__/boris.cpython-37.pyc new file mode 100644 index 0000000..434fcdb Binary files /dev/null and b/__pycache__/boris.cpython-37.pyc differ diff --git a/boris.py b/boris.py index c3619d3..028f6b9 100644 --- a/boris.py +++ b/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) diff --git a/zemlja.py b/zemlja.py new file mode 100644 index 0000000..f9fffb8 --- /dev/null +++ b/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) +