From 70406c0d76ea427eb939b9d3fd8b01c475938616 Mon Sep 17 00:00:00 2001
From: Andrej <andrej@baza>
Date: Sun, 6 Jun 2021 15:54:33 +0200
Subject: [PATCH] delujoc prototip delca v zemljinem magnetnem polju

---
 __pycache__/boris.cpython-37.pyc | Bin 0 -> 2886 bytes
 boris.py                         |  61 ++++++++++++++++++++++++-------
 zemlja.py                        |  32 ++++++++++++++++
 3 files changed, 80 insertions(+), 13 deletions(-)
 create mode 100644 __pycache__/boris.cpython-37.pyc
 create mode 100644 zemlja.py

diff --git a/__pycache__/boris.cpython-37.pyc b/__pycache__/boris.cpython-37.pyc
new file mode 100644
index 0000000000000000000000000000000000000000..434fcdb130acf84221517ccc00c9c87e7d2e435f
GIT binary patch
literal 2886
zcmai0O^h2!74E9;wtH;PjA!g*CO^z}et@+Yb|OVWO9(q?mVgwo0T#p}FOj#$Rd(8T
zcRSVX!FJlG$zDJ~xd!cFjJR+BX}NOXg0vv*0f}1)ArKPM+!i5$gpk7bx@RVtXr--I
zy?XD}tE%T$-+Seq#l@0=_V{PNdFR6g!}vQ1$E^tFE9mhbL4+aLz-Y2;<bm0=80q-H
z4hqczGh{(r7M8HPOcdVZKVg}sC>*iyh{LkfEQ)f|5tZhGSZtQW6|pR;4_UJ;R>bN<
zqgfF(QOCF_u8K9Wj+rI7EY9`0`cRw~7cf&57qf9ivV&FeDX}3oy5^iUYni+uUOJNN
z!fl?BWm%JTS(PjDjJUL4_yKR8ZLY~Bd3GP33HBYM@##bSD&tOrX|CGZ>O^6p?N}<m
zGhKc?3WlRZD))^@NEJ4?E_A%WA54^mHQGw#IQcGC7DllI-~IOc|6bESFW&5KT>kNs
zCy(E`(f#KKAN=yqlV5+)zdioT$|1w5wnOKk)zA@Y<8}1-A0Ua5+NpKH*k>~*VFx^Q
z1nclQn=PD@N~tj`r)*Y9t+Y@_=<7yWykN{0$L}R(>g=<x8h0PhmeR%4Ow0SMZp@a2
znJ&Z9!ai2IoG#D#nt>J7v;r%hTdf(hm2?Fb_A$-v(kfZcR@2o}ww1Ju^DA`&b7kyP
zN|%H!?00#G-!w3{@VtEg_FD!tr4C2H3!_!U!QF|JA4hwxH|T;(WZ+J{-Ja}vZrqRh
z9;A|s=7)*wDzEMLB3Jkwsbtvpy+J5l8G65av@1)Q+4!3KW%upz=45kv*S+EH$$k<k
z?9rA9Oo!3b#||<x-OEh(#+!SS&ApGazUL<@iW5JxzJ=X7qp*+V-lp64UE%el+w=Oc
zz`lc?=iYjG+{kQqKG}BX*j6Mr+_r}RgudI4l-G9;O<<slj`r8>LykcaR_-(k+C&&M
z4~N>GO8BO2<%M0TO%VaCZD`^~P4n?J&G)Wpc1yDz%|+6*5;0Pm-PbJ84p=YoqfoQ&
zYxX_O#+vOli}zbY<p&a2Cd%$<7Hc+Y@=@y-hPs5~|NH9Jy(o}dUMQ68ZFRhUyfqwn
zao}yme%Kwz&ksG7_-!yjMEh^;(#^gyoM=1a7O_d}f*7R=Ut&wlVivEOORU1yS(RPn
zRd7{S=jsY<+LR;IY;^khS$g~kj9@b+_>3pU0Ym4hA<UXFXGjn;Wo0Af!m1fa22<EI
zg141gurq}P%fdOD!3r8PP4Se*%;#vw%ntU%uhSY><~-Gzx&a6~)7x%<k}z<mG8pvy
zkYF2l-B5;r8~7X;-cUvT*bjW%mOFJbjAT|x=m$VEtfIi}p1<<@|2$YXw*4n>dWNDp
ztf(#=QUnx*BZEQ1QWfypitnpLoAD^nwx^UgY1rxlCYpQ)*Z<~JGu(JTySkiQrd7^P
zD|6L*1yY>VdXw&<&dxIRIq2GiH)=G}$ywaj!3dTjCKUA)@vk8kbCeuJ!MbyT*4%5l
zt_IKCRG-QDH4N%J5yH}Q4)M48uR;@N$1=4^vt<0WHl+bIa+02&%#kI;fn1myxNNh>
z1zdT=aK&7q>K#}D*PNI=8<%U(&6%C+NR!D0$kdCeJ!cr5v;g4nS@B5kyr6fE^a~hE
zG?%R^9<5qPtzJ1TrcSyb3ZnQ3ParfEKn35sT$mkxj)0R(&MC<vNh}<%I#?8?gQXPv
z&+#Co*1@tUpOUMh@=QhqR^|phxCsB^9nG%`ygntstS8%)dpTSva)*M*|DHa5O7GJ&
zQh2~J935Jim^xR-SKYTC%<h7loN&{eYpiN>I7sl=`Q4F{y6g$jibuP{L6qbmbOwGH
z4^dq;ANtyk=#kWB8<ojKQP5Tt${0xXr&3pF6;`VbnOh|35}wP}c;M~I0S<{&;2r5E
zxh|WZ=DZc(i+<QBsdI2!+cJz2iQhwM=V=oTYVM7di%E5n2&IB%6V0YCBF?6xuYU>i
zF<BbSLEW)f9kpj2-v(RdUjUDftJ}1T(`p5tAGTU2HHnh(BsZ!M2$BRvVNz{)J||Gh
zfBJ;lAGW*Kv-YR$<J|ZxRyJ7X>WfY~HwRv@E4<gf37dbR8|oV%r#v`G3d+tS58y*y
z!Rp4xctZOa#z)&3JU-cr?8t9ng?a<zD9Jzi?Gf`#v)z4A|JB$3Hhs8#Scf5g66G%T
zr$UeyEuJTJV_n<fC>Tz1_^jZ<={5%bF5Za@5`}zJ*RfQWp3N&N5Sjrw0bZ&mC@{m?
z*;|nnuDo;{Us&h`vei=b-%62Biejx+h^!JhLxiC8Vo}P<(`8I(=O#+~XdqvsljD~`
g2q+b{$|@WslRjR$<dht<RC53#2M<`4Tl}o`e;p3FK>z>%

literal 0
HcmV?d00001

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)
+